<jupyter_start><jupyter_text>"In this exercise, you will implement one-vs-all logistic regression and neural networks to recognize hand-written digits.""For this exercise, you will use logistic regression and neural networks to recognize handwritten digits (from 0 to 9). Automated handwritten digit recognition is widely used today - from recognizing zip codes (postal codes) on mail envelopes to recognizing amounts written on bank checks. This exercise will show you how the methods you’ve learned can be used for this classification task."<jupyter_code>#import the relevant libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
np.random.seed(79728)
%matplotlib inline
#%whos<jupyter_output><empty_output><jupyter_text>I added the line magic whos here. It shows the variables and functions you have defined as a list when you want to get information about what you definedWe need to import a matlab file. scipy has a library for this.<jupyter_code>import scipy.io as sio
data = sio.loadmat('../data/ex3data1.mat')
data<jupyter_output><empty_output><jupyter_text>This is a dictionary file and X, y are what we are looking for. Let's define and extract them.<jupyter_code>X = np.float64(data["X"])
y = np.float64(data["y"])
X.shape, y.shape<jupyter_output><empty_output><jupyter_text>It is a good idea to check the type of your data and convert it to float apparently. The reason for this is, when you preserve y as an integer vector, taking the negative of it will do interesting tricks on your data. Python will convert y to a 8bit integer vector, which has no representation for negative numbers. Then when you take negative of say 2, it will convert it to 256 - 2 = 254! What we have here is 5000 samples. Each row is an array of numbers that represent brightness level of pixels. Each square consist of $20x20=400$ pixel. We will insert a column of one as the constant value.<jupyter_code>X = np.insert(X, 0, 1, axis=1)
X.shape
#We'll have 401 columns as a result<jupyter_output><empty_output><jupyter_text>Let's check how many there are of each element:<jupyter_code>pd.Series(y.flatten()).value_counts()<jupyter_output><empty_output><jupyter_text>Now let's look at how the data is structured very briefly. As expected the edges of boxes are empty. The zigzags in the middle mostly represent the top and the bottom pixels which are empty. They happen around the multiples of 20.<jupyter_code>plt.figure(figsize=(12, 8))

ser_mean = np.asarray([np.mean([X[:, i]]) for i in range(1, 401)])
ser_sd = np.asarray([np.std([X[:, i]]) for i in range(1, 401)])
plt.plot(np.arange(1, 401), ser_mean, label = "Mean")
plt.plot(np.arange(1, 401), ser_sd, label = "Standart Deviation")
plt.legend()
plt.title("Mean and Std Dev of X columns");<jupyter_output><empty_output><jupyter_text>Finally, let's illustrate some of the data we have. As we know from the description of the dataset, each entry in X consist of 400 numbers, each number illustrating the degree of darkness of each pixel. We can actually map this to a graph, or pic if you'd like through imshow function in matplotlib library. We'll first draw 36 random entry to show. Then transform the one dimensional arrays of info into 20x20 square matrix. Finally we'll plot the each data row with a for loop.<jupyter_code>num_show = np.random.randint(0, 5001, 36)

plt.figure(figsize=(8, 8))

for i in range(1, 37):
    some_digit = X[num_show[i-1]][1:401]
    some_digit_image = some_digit.reshape(20, 20).T
    #We appearantly need to transpose the matrix
    
    plt.subplot(6, 6, i)
    plt.imshow(some_digit_image, cmap = plt.cm.binary,
           interpolation="nearest")
    plt.axis("off")
    plt.subplots_adjust(left=0, right = 1, bottom=0, top=1, wspace=0, hspace=0)<jupyter_output><empty_output><jupyter_text>**Cost Function**Let's try to write down the cost function that we will try to minimize for logistic regression. We will use the cost function from "Logistic Regression with regularization". We already vectorized our logistic cost function there.First define the Hypothesis function<jupyter_code>def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

def h(theta, X):
    """
    Hypothesis Function where
    X is an n x k dimensional array of explanatory variables
    theta is an one dimensional k elements vector
    Result will be one dimensional vector of n variables
    """
    return sigmoid(np.dot(X, theta))<jupyter_output><empty_output><jupyter_text>Now the cost function. One thing to note here is that we keep our theta variable as a flat, one dimensional array. That's because we need to use scipy optimization libraries. It needs flattened input.<jupyter_code>def cost_function(theta, Lambda, X, y):
    """
    This is a cost function that returns the cost of theta given X and y
    X is an n x k dimensional array of explanatory variables
    y is an array with n elements of explained variables. I flatten it no matter what its shape is
    theta is a vector with k elements of multipliers for linear function
    """
    #Let's write this function without the regularization
    item1 =  - y.flatten() * np.log(h(theta, X))
    item2 = (1 - y.flatten()) * np.log(1 - h(theta, X))
    item3 = (Lambda/(2*X.shape[0]) ) * np.sum(np.power(theta, 2)[1:]) 
    #Starting from the second element as the first theta will not be punished 
    return np.sum(item1 - item2) / (X.shape[0]) + item3

Lambda = 0.0
#theta = np.random.rand(X.shape[1]) - 0.5
theta = np.zeros(X.shape[1])<jupyter_output><empty_output><jupyter_text>How about gradient function? We already took care of it too. It has already been vectorized.<jupyter_code>def gradient(theta, Lambda, X, y):
    """
    This function will take in theta value and calculate the gradient 
    descent values.
    X is an n x k matrix of explanatory variables
    y is a n sized array of explained variables
    theta is a vector with k elements of multipliers for linear function
    """
    #Let's try this function without regularization
    errors = h(theta, X) - y.flatten()
    #errors will be calculated more than once, so let's do it once and store it
    correction2 = (Lambda/X.shape[0]) * theta
    correction2[0] = 0.0 #First element will not be part of regularization
    correction = np.sum(errors.reshape((X.shape[0], 1)) * X, axis=0) * (1.0 / X.shape[0])
    return correction + correction2<jupyter_output><empty_output><jupyter_text>**One vs all Classification**We will apply one vs all classification. What it does is that it estimates the probability function separately for each unique value of y. We will actually estimate the probability through hypothesis function $h_\theta(X)$ and assign to highest value. This means transforming y vector for every value of y. For example for value 5, we will call y = 1 if y = 5 originally, y = 0 otherwise. We will do this for each value between 0-9. Let's try it for y = 5.<jupyter_code>y_flattened = y.flatten() #easier to work with one dimensional arrays
#Find where y != 5 and mark them equal to zero: 
y_flattened[np.where(y_flattened != 5.0)] = 0.0
#Find where y == 5 and mark them equal to zero: 
y_flattened[np.where(y_flattened == 5.0)] = 1.0<jupyter_output><empty_output><jupyter_text>Now let's see if we had a successful transformation:<jupyter_code>pd.Series(y_flattened).value_counts()<jupyter_output><empty_output><jupyter_text>Let's try to optimize the function and find the $\theta$s. I have tried many different methods and the one that works well is SLSQP. You can google it to see how it works.<jupyter_code>Lambda = 0.0
theta = np.zeros(X.shape[1])
result = opt.minimize(fun=cost_function, method='SLSQP', 
                      jac = gradient, x0=theta, args=(Lambda, X, y_flattened.reshape((X.shape[0], 1))), 
                      options={'maxiter':1000})<jupyter_output><empty_output><jupyter_text>Let's write the function that would optimize our cost function and find the relevant $\theta$s for all values of y.<jupyter_code>def onevsall(X, y, Lambda = 0.0):
    results = []
    thetas = []
    successes = []
    
    for i in range(1, 11):
        theta = np.zeros(X.shape[1])
        #define the initial theta
        y_flattened = y.flatten()#easier to work with one dimensional arrays
        #Find where y != i and mark them equal to zero: 
        y_flattened[np.where(y_flattened != float(i))] = 0.0
        #Find where y == i and mark them equal to zero: 
        y_flattened[np.where(y_flattened == float(i))] = 1.0
        #optimization routine
        result = opt.minimize(fun=cost_function, method='SLSQP', 
                      jac = gradient, x0=theta, args=(Lambda, X, y_flattened.reshape((X.shape[0], 1))), 
                      options={'maxiter':1000})
        results.append(result)
        thetas.append(result['x'])
        successes.append(result['success'])
    
    return results, thetas, successes
res, thet, succes = onevsall(X, y)<jupyter_output>/Users/alperen/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: divide by zero encountered in log
  # Remove the CWD from sys.path while we load stuff.
/Users/alperen/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:10: RuntimeWarning: invalid value encountered in multiply
  # Remove the CWD from sys.path while we load stuff.<jupyter_text>'thet' contains all the multipliers for logistic regressionLet's see if all the loops successfully converged:<jupyter_code>succes<jupyter_output><empty_output><jupyter_text>Let's calculate the probabilities for each class.<jupyter_code>#convert theta list into numpy array
thetas = np.array(thet).T
thetas.shape
Probabilities = h(thetas, X)
Probabilities.shape<jupyter_output><empty_output><jupyter_text>Probabilites is an 2d array where you have 10 entries for each 5000 entries. You may think that this is like a matrix, each row correspond to an entry and each column correspond to the probability of that entry is that particular given number. But my advice is not to think arrays like matrix, but stacks of arrays, which are one dimensional vectors. Now let's locate each of the argument maximizing location of probability within its probability array.<jupyter_code>Prob_argmax = np.argmax(Probabilities, axis=1)<jupyter_output><empty_output><jupyter_text>Let's check the distribution of values as well:<jupyter_code>pd.Series(Prob_argmax.flatten()).value_counts()<jupyter_output><empty_output><jupyter_text>It is looking pretty good but it does not mean that our guesses correspond to the correct values. We need to make sure that our guesses are correct. Because of the way I have arranged the sequence of calculating $\theta$, here I have the following situation: 0 in images, which correspond to 10 in the dataset here correspond to 9! 1 in the dataset correspond to 0; 2 to 1; 3 to 2 and so on. Good news is that, I can add one to each number and it will solve my problem<jupyter_code>Prob_argmax += 1<jupyter_output><empty_output><jupyter_text>Finally let's check the prediction error:<jupyter_code>total_corrects = np.sum( (y.flatten() == np.float64(Prob_argmax)) )
total_dpoints = X.shape[0]

accuracy_rate = total_corrects/total_dpoints
accuracy_rate<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08860533595074507, 'fraction_numerical': 0.018738288569643973, 'mean_word_length': 4.40144578313253, 'pattern_counts': {'":': 0, '<': 80, '<?xml version=': 0, '>': 80, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '595813', 'license': 'MIT', 'n_tokens_mistral': 3597, 'n_tokens_neox': 3357, 'n_words': 1350, 'path': 'Ex 3 - Neural Networks/Multiclass-Logistic Regression.ipynb', 'repo_name': 'alperenevrin/MachineLearniNg', 'size': '11207'}