Machine learning algorithm implemented in Python

• All the code

Def computerCost(X,y,theta): m = len(y) J = 0 J = (np. Transpose (X*theta-y))*(X*theta-y)/(2*mCopy the code

• Notice that X is the real number with a 1 in front of it, because it has theta of 0.

• The goal is to scale the data to a range that is easy to use the gradient descent algorithm
• 
• Where is the average of all this Feture data
• It can be the maximum-minimum value, or the standard deviation of the data corresponding to this feature
• Implementation code:

Feature def featureNormaliza(X): X_norm = np.array(X) # Zeros ((1, x.shape [1])) sigma = Np.zeros ((1,X.shape[1])) mu = np.mean(X_norm,0) # Sigma = np.std(X_norm,0) # For I in range(x.shape [1]): # traversal column X_norm [, I] = (X_norm [, I] - mu [I])/sigma [I] # normalized return X_norm, mu, the sigmaCopy the code

• Note that you also need mean normalization data for forecasting

• Import packages

From sklearn import linear_model from sklearn. Preprocessing import StandardScaler #Copy the code

• The normalized

Scaler = StandardScaler() scaler.fit(X) x_train = scaler.transform(X) x_test = Scaler. Transform (np) array ([1650]))Copy the code

• Linear model fitting

LinearRegression() model.linear regression () model.fit(x_train, y)Copy the code

• To predict

Result = model. Predict (x_test)Copy the code

• All the code

You can see that whenTend to be1.y=1, the price paid at this time is consistent with the predicted valuecostTend to be0If,Tend to be0.y=1, the cost at this timecostIt’s a very large value, and our ultimate goal is to minimize the generation value

• In the same wayIs shown below (y=0) :
  
  

• The goal is to prevent overfitting
• Add a term to the cost function
• Note that j starts with the increment of 1, because theta(0) is a constant term and the first column of X will add 1 column 1, so the product is still Theta (0). Feature does not matter, and regularization is not necessary
• Cost after regularization:

Def costFunction(initial_theta,X,y,inital_lambda): M = len(y) J = 0 h = sigmoid(np.dot(X,initial_theta)) # calculate h(z) theta1 = initial_thet.copy () # The pre-theta (0) value is 0 theTA1 [0] = 0 Temp = NP.dot (NP.transpose (theta1),theta1) J = (-np.dot(np.transpose(y),np.log(h))-np.dot(Np.transpose (1-y),np.log(1-h))+temp*inital_lambda/2)/m # regularized cost equation return JCopy the code

• The regularized gradient of the cost

Def gradient(initial_theta,X,y,inital_lambda): M = len(y) grad = Np.zeros ((initial_thet.shape [0])) h = sigmoid(np.dot(X, initial_Theta))# calculate h(z) theta1 = Initial_thet.copy () theta1[0] = 0 grad = np.dot(NP.transpose (X), H-Y)/m+inital_lambda/m*theta1 # regularized gradient return gradCopy the code

• Implementation code:

Def sigmoid(z): h = np.zeros((len(z),1)) h = 1.0/(1.0+np.exp(-z)) return hCopy the code

• Gradient descent servicescipyIn theoptimizeIn thefmin_bfgsfunction
• The optimization algorithm fMIN_BFGS (quasi-Newton method Broyden-Fletcher-Goldfarb-Shanno) in SCIPY was called
• CostFunction is a self-implemented costFunction,
• Initial_theta stands for the initialized value,
• Fprime specifies the gradient of the costFunction
• Args is the remaining test parameters, passed in as a tuple, and will eventually return theta that minimizes costFunction

    result = optimize.fmin_bfgs(costFunction, initial_theta, fprime=gradient, args=(X,y,initial_lambda))    
Copy the code

• Import packages

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.cross_validation import train_test_split
import numpy as np
Copy the code

• Divide training set and test set

X_train,x_test,y_train,y_test = train_test_split(X,y,test_size=0.2)Copy the code

• The normalized

Scaler = StandardScaler() fit(x_train) x_train = scaler. Fit_transform (x_train) x_test = scaler.fit_transform(x_test)Copy the code

• Logistic regression

Model = LogisticRegression() model. Fit (x_train,y_train)Copy the code

• To predict

# predict = model. Predict (x_test) right = sum(predict == y_test) predict = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 % % % % f '(right * 100.0 / predict shape [0])) # the accuracy of the calculation on the test setCopy the code

• All the code

• How to use logistic regression to solve the problem of multiple classification, OneVsAll is to treat a certain category as a class, all other categories as a class, so there is a dichotomous problem
• As shown in the figure below, the data on the way are divided into three categories. First, the red ones are regarded as one category, and the others as another category for logistic regression. Then, the blue ones are regarded as one category, and the others are regarded as one category, and so on…
• It can be seen that when the number of classes is greater than 2, logistic regression classification should be carried out as many times as there are classes

All_theta def oneVsAll(X,y,num_labels,Lambda); X = np.hstack((np.ones((m,1)),X)) # X = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 Initial_theta = np.zeros((n+1,1)) # Initializes a category of theta # mapping y for I in range(num_labels): Class_y [, I] = np. Int32 (y = = I). Reshape # (1, 1) attention to reshape (1, 1) to assignment # np. Savetxt (" class_y. CSV, "class_y [00 0-6. :], Delimiter =',') "" traverses each category and calculates the corresponding theta value" "for I in range(num_labels): result = optimize.fmin_bfgs(costFunction, initial_theta, fprime=gradient, 0 Args =(X,class_y[:, I],Lambda) # calling gradient descent optimization method ALL_theta [:, I] = Result. 0 (1,-1) # Putting all_theta = np.transpose(all_theta) return all_thetaCopy the code

• 10 times classification, accuracy on training set:

• 1. Import the package

from scipy import io as spio
import numpy as np
from sklearn import svm
from sklearn.linear_model import LogisticRegression
Copy the code

• 2. Load data

Data = loadmat_data("data_digits. Mat ") X = data['X'] # 20x20px y = data['y'] # shape=(5000, 1) y = np.ravel(y)Copy the code

• 3. Fitting model

Model = LogisticRegression() model.fit(X, y) #Copy the code

• 4, forecasting

Predict = model. Predict (X) # predict = model. Predict (X) # predict = model. Predict (X) # predict = model.Copy the code

• 5. Output results (accuracy on training set)

Def nnCostFunction(nn_params,input_layer_size,hidden_layer_size,num_labels,X,y,Lambda): Shape [0] # restore theta1 and theta2 theta1 = nn_params[0:hidden_layer_size*(input_layer_size+1)].reshape(hidden_layer_size,input_layer_size+1) Theta2 = nn_params[hidden_layer_size*(input_layer_size+1):length].reshape(num_labels,hidden_layer_size+1) # Np.savetxt (" theta1.csv ",Theta1,delimiter=',') m = x.shape [0] class_y = Np.zeros ((m,num_labels)) # The y of the data corresponds to 0-9, # mapping y for I in range(num_labels): 0 0 class_Y [:, I] = np.int32(y== I).0 Because "" Theta1_colCount = theta1. shape[1] Theta1_x = Theta1[:,1:Theta1_colCount] Theta2_colCount = theta2.shape [1] Theta2_x = Theta2[:,1:Theta2_colCount] # regularize to theta^2 term = Np. Dot (np) transpose (np) vstack ((Theta1_x. Reshape (1, 1), Theta2_x. Reshape (1, 1)))), np. Vstack ((Theta1_x. Reshape (1, 1), Theta2 _x. Reshape (1, 1)))) ' ' 'forward travel, every time you need to fill in a column 1 offset bias,' 'a1 = np. Hstack ((np) ones ((m, 1)), X)) z2 = np. Dot (a1, np. Transpose (Theta1)) Hstack ((np.ones((m,1)),a2) z3 = Np.dot (A2, nP.transpose (Theta2)) h = sigmoid(z3) J = - (np) dot (np. Transpose (class_y. Reshape (1, 1)), np. The log (practice eshape (1, 1))) + np. Dot (np) transpose (1 - class_y. Reshape (1, 1)), np, lo G (1 - practice eshape (1, 1))) - Lambda * term / 2)/m return np ravel (J)Copy the code

Def nnGradient(nn_params,input_layer_size,hidden_layer_size,num_labels,X,y,Lambda): length = nn_params.shape[0] Theta1 = Nn_params [0:hidden_layer_size*(input_layer_size+1)].0 (hidden_layer_size,input_layer_size+1).copy() Otherwise let's change the value of Theta, Nn_params also modifies Theta2 = nn_params[hidden_layer_size*(input_layer_size+1):length].reshape(num_labels,hidden_layer_size+1).copy() m = X.shape[0] Class_y = np.zeros((m,num_labels)) # y for I in range(num_labels): 0 0 class_Y [:, I] = np.int32(y== I).0 Because "" Theta1_colCount = theta1. shape[1] Theta1_x = Theta1[:,1:Theta1_colCount] Theta2_colCount = theta2.shape [1] Theta2_x = Theta2[:,1:Theta2_colCount] Theta1_grad = Np.zeros ((theta1.shape)) # Theta2_x = Theta2[:,1:Theta2_colCount] Theta1_grad = Np.zeros ((theta1.shape)) # Np. zeros((theta2.shape)) # Hstack ((Np.ones ((m,1)),X)) z2 = Np.dot (A1, NP.transpose (Theta1)) A2 = sigmoid(z2) A2 = Hstack ((np.ones((m,1)),a2)) z3 = np.dot(A2, nP.transpose (Theta2)) h = sigmoid(z3) ''' delta3 = np.zeros((m,num_labels)) delta2 = np.zeros((m,hidden_layer_size)) for i in range(m): # delta3 [I:] = (h [I:] - class_y [I:]) * sigmoidGradient (z3 [I:]) # error of mean square error (mse) delta3 [I:] [I:] - class_y = h [I:] # cross entropy error rate Theta2_grad = Theta2_grad+np.dot(np.transpose(delta3[i,:].reshape(1,-1)),a2[i,:].reshape(1,-1)) delta2[i,:] = np.dot(delta3[i,:].reshape(1,-1),Theta2_x)*sigmoidGradient(z2[i,:]) Theta1_grad = Theta1_grad+np.dot(np.transpose(delta2[i,:].reshape(1,-1)),a1[i,:].reshape(1,-1)) Theta1[:,0] = 0 Theta2[:,0] = 0 "Gradient" "grad" = (np) vstack ((Theta1_grad. Reshape (1, 1), Theta2_grad. Reshape (1, 1))) + Lambda * np vstack ((Theta1. Reshape (1, 1), Theta2. Reshape ( - 1, 1))))/m return np ravel (grad)Copy the code

• Neural networks cannot be initialized like logistic regressionthetafor0, because if the weight of each edge is 0 and each neuron has the same output, it will get the same gradient in the back propagation, and only one result will be predicted.
• So it should be initialized to something close to zero
• The implementation code

Theta def randInitializeWeights(L_in,L_out): W = np. Zeros ((1 + L_in, L_out)) # corresponding to the weight of theta epsilon_init = (6.0 / (L_out + L_in)) * * W = 0.5 Rand (L_out,1+L_in)*2*epsilon_init-epsilon_init # np.random.rand(L_out,1+L_in) generates a random matrix of L_out*(1+L_in) size return WCopy the code

• Gradient check:
• Random display of 100 handwritten digits
• Show theta1 weight
• Training set prediction accuracy
• Prediction accuracy of training set after normalization

• In logistic regression, our cost is:, where:,
• As shown, ify=1.costThe cost function is shown in the figure
  
  
  
  We want to make, i.e.,z>>0In this casecostThe cost function tends to be minimal (that’s what we want), so the utilityredThe function ofInstead of cost in logistic regression
• wheny=0As well asInstead of
• The resulting cost function is: finally we want to
• The cost function in our logistic regression was:
  
  
  
  You can think of it hereIt’s just a matter of form hereCThe greater the value of the SVM decision boundarymarginThe larger it is, as will be explained below

Loadmat ('data1.mat') X = data1['X'] y = data1['y'] y = np.ravel(y) plot_data(X,y) model = Svm.svc (C=1.0,kernel='linear').fit(X,y) #Copy the code

• Nonlinear separable, the default kernel function isrbf

Data2 = IO. Loadmat ('data2.mat') X = data2['X'] y = data2['y'] y = np.ravel(y) PLT = plot_data(X,y) Plt.show () model = svm.svc (gamma=100).fit(X,y) # gammaCopy the code

• All the code

Def findClosestCentroids(X,initial_centroids): M = x.shape [0] # number of data items K = initial_centroids.shape[0] # number of classes dis = Np.zeros ((m,K)) # Storage calculation of the distance from each point to K classes IDx = For I in range(m): for j in range(K): for j in range(K): for j in range(K): Dis [I, j] = np. Dot ((X [I:] - initial_centroids [j:]), reshape (1, 1), (X [I:] - initial_centroids [j:]), reshape (1, 1)) Returns the column number corresponding to the minimum value of each row of dis, 0 0 -np.min (dis, Axis =1) returns the minimum value of each line -np.min where(dis == Np.min (dis, Axis =1).0 There may be multiple coordinates corresponding to the minimum value, which can be found by WHERE. Therefore, when returning, the first m coordinates need to be returned (because for multiple minimum values, Dummy,idx = np.min where(dis == np.min(dis, axis=1).0 (-1,1)) return idx[0:dis.shape[0]] #Copy the code

• Computing class center implementation code:

Def computerCentroids(X,idx,K): def computerCentroids(X,idx,K): n = x.shape [1] centroids = np.zeros((K,n)) for I in range(K): Centroids [I:] = np. The mean (X [np. Ravel (independence idx = = I), :], axis = 0). Reshape index (1, 1) # if a one-dimensional, axis = 0 for each column, independence idx = = once I find out what kind of, Then calculate the mean return centroidsCopy the code

• The label of Y is not known for clustering, so the real number of clustering is not known
• The Elbow Method
• Cost functionJandKIf there is an inflection point, as shown in the figure below,KTake the value at the inflection point, right hereK=3
• If very smooth is not clear, artificial choice.
• The second is human observational selection

• Import packages

    from sklearn.cluster import KMeans
Copy the code

• Use models to fit data

Model = KMeans(n_clusters=3).fit(X) # n_clusters specifies 3 classes and fits the dataCopy the code

• The clustering center

Centroids = model.cluster_centers_ # Cluster centersCopy the code

• All the code

• As shown in the figure below, all data points can be projected onto a straight line, minimizing the sum of squares (projection error) of projected distances
• Pay attention to data needsThe normalizedTo deal with
• Idea is to find1aThe vector u, all data are projected onto it to minimize the projection distance
• thennD-->kDIs to look forkA vector, all data are projected onto it to minimize the projection error
• Eg :3D–>2D,2 vectors represent a plane, the projection error of all points to this plane is minimum

• Data preprocessing (mean normalization)
• Formula:
• That is, subtract the mean value of the corresponding feature and divide by the standard deviation of the corresponding feature (also can be the maximum – minimum value).
• Implementation code:

Def featureNormalize(X, X, X, X, X, X, X, X, X, X, X, X, X, X, X); sigma = np.zeros((1,n)) mu = np.mean(X,axis=0) sigma = np.std(X,axis=0) for i in range(n): X[:,i] = (X[:,i]-mu[i])/sigma[i] return X,mu,sigmaCopy the code

• To calculateSigma covariance matrix(Covariance Matrix) :
• Notice hereΣIt’s not the summation notation
• Covariance matrixSymmetric positive definite(If you don’t understand positive definite, look at line substitution.)
• Size of thenxn.nforfeatureThe dimensions of
• Implementation code:

Sigma = Np.dot (Np.transpose (X_norm),X_norm)/m # find SigmaCopy the code

• To calculateΣEigenvalues and eigenvectors of
• Can be usedsvdSingular value decomposition function:U, S, V = SVD (Σ)
• And is returnedΣDiagonal matrix of the same sizeS(byΣOf) [Pay attention to:matlabThe function returns a diagonal matrix atpythonIs a vector, saving space.]
• There are also two unitary matrices U and V, and
• 
• Pay attention to:svdDelta functionSIs in descending order of eigenvalue, if not usedsvd, need to pressThe eigenvalueSize rearrangementU
• Dimension reduction
• selectUIn the firstKColumn (suppose to drop toKD)
• 
• ZThe corresponding data after dimension reduction
• Implementation code:

Def projectData(X_norm,U,K): Z = np.zeros((X_norm. Shape [0],K)) U_reduce = U[:,0:K] # take first K Z = np.dot(X_norm,U_reduce) return ZCopy the code

• Process summary:
• Sigma = X'*X/m
• U,S,V = svd(Sigma)
• Ureduce = U[:,0:k]
• Z = Ureduce'*x

Def recoverData(Z,U,K); X_rec = Np.zeros ((Z.shape[0],U.s shape[0])) U_recude = U[:,0:K] X_rec = Np.dot (Z, NP.transpose (U_recude)) # Restore data (approximate) return  X_recCopy the code

• The 2-d data has been reduced to 1-d
• The direction of the projection
• 2D reduced to 1D and the corresponding relationship
• Face data dimensionality reduction
• The original data
• Visual sectionUThe matrix information
  
  
• Restore data

• Import required packages:

#-*- coding: Utf-8 -*- # Author: Bob # Date:2016.12.22 Import numpy as NP from matplotlib import pyplot as PLT from scipy import IO as utF-8 -*- # Author: Bob # Date:2016.12.22 Import numpy as NP from matplotlib import pyplot as PLT from scipy import IO as  spio from sklearn.decomposition import pca from sklearn.preprocessing import StandardScalerCopy the code

• Normalized data

Scaler = StandardScaler() scaler.fit(X) x_train = scaler.transform(X)Copy the code

• PCA model was used to fit data and reduce dimension
• n_componentsCorresponds to the dimension to be transferred

Model = pca.pca (n_components=K).fit(x_train) # N_components defines the dimension to be reduced Z = model.transform(x_train) # transform will perform the dimension reduction operationCopy the code

• Data recovery
• model.components_It will be used for dimensionality reductionUmatrix

Data recovery and plotting Ureduce = model.components_ # Get Ureduce x_rec = np.dot(Z,Ureduce) # Data recoveryCopy the code

• Mismeasurement of skew data

• Since the data can be very skewed (i.e. the number of y=1 is very small (y=1 indicates an exception)), F1Score can be calculated using Precision/Recall (on the CV cross validation set).

• For example, in predicting cancer, if the model can get 99% correct prediction and 1% error rate, but the actual probability of cancer is very small, only 0.5%, then we always predict that there is no cancer, but y=0 can get a smaller error rate. Using error rate for evaluation is not scientific.

• As shown below:
  
  

• , that is:Predict positive samples correctly/all predict positive samples correctly

• , that is:Correctly predicted positive sample/true value is positive sample

• Always make y=1(fewer classes) and calculate Precision and Recall

• 

• Taking cancer prediction as an example, it is assumed that all predictions are no-cancer, TN=199, FN=1, TP=0, FP=0, so Precision=0/0, Recall=0/1=0, although accuracy=199/200=99.5%, it is not credible.

• The selection of epsilon.

• Try multiple ε values to make F1Score high

• The implementation code

Def selectThreshold(yval,pval): Initialize the required variable bestEpsilon = 0.bestf1 = 0.f1 = 0.step = (np.max(pval)-np.min(pval))/1000 "" calculate" for epsilon in np.arange(np.min(pval),np.max(pval),step): CvPrecision = pval<epsilon tp = np.sum((cvPrecision == 1) & (yval == 1)). Astype (float) # sum is int, Float fp = np.sum((cvPrecision == 1) & (yval == 0)). Astype (float) fn = np.sum((cvPrecision == 1) & (yval == 0)). Astype (float) precision = TP /(TP +fp) # recision = TP /(TP +fn) # recall rate F1 = (2*precision*recision)/(precision+recision) # F1Score if F1 > bestF1: BestF1 = F1 bestEpsilon = epsilon return bestEpsilon,bestF1Copy the code

• Problems existing in element Gaussian distribution
• The red dots are outliers and the rest are normal (such as CPU and memory changes).
• The gaussian distribution of X1 is as follows:
• The gaussian distribution of X2 is as follows:
• It can be seen that the corresponding values of P (x1) and P (x2) do not change much, so they are not considered abnormal
• Since we believe that features are independent of each other, the above figure is extended in a circular way
• Multivariate Gaussian distribution
• It’s not buildingp(x1),p(x2)... p(xn)It’s about unityp(x)
• Parameters:.ΣforCovariance matrix
• 
• In the same way,| Σ |The smaller,p(x)The more sharp
• For example,, represents the positive correlation between x1 and x2, that is, the larger x1 is, the larger X2 will be, as shown in the figure below, the red anomaly points can be detected. If:, represents the negative correlation between x1 and x2
• Implementation code:

Def multivariateGaussian(X,mu,Sigma2): k = len(mu) if (sigma2.shape [0]>1): Sigma2 = np. Diag (Sigma2) ' ' ' 'multivariate gaussian distribution function' X = X - mu argu = (np. 2 * PI * * (k / 2) * np linalg. Det (Sigma2) * * p = (0.5) Argu *np.exp(-0.5*np.sum(np.dot(X,np.linalg.inv(Sigma2))*X,axis=1)) # axis indicates that each line returns pCopy the code

• Display the data
• contour
• Outlier labeling