用v9做的网站上传服务器,免费收录软文网站,正规漫画网站开发流程,flash企业网站源码机器学习算法Python实现目录1、代价函数其中#xff1a;下面就是要求出theta#xff0c;使代价最小#xff0c;即代表我们拟合出来的方程距离真实值最近共有m条数据#xff0c;其中代表我们要拟合出来的方程到真实值距离的平方#xff0c;平方的原因是因为可能有负值…机器学习算法Python实现目录1、代价函数其中下面就是要求出theta使代价最小即代表我们拟合出来的方程距离真实值最近共有m条数据其中代表我们要拟合出来的方程到真实值距离的平方平方的原因是因为可能有负值正负可能会抵消前面有系数2的原因是下面求梯度是对每个变量求偏导2可以消去实现代码# 计算代价函数def computerCost(X,y,theta):m len(y)J 0J (np.transpose(X*theta-y))*(X*theta-y)/(2*m) #计算代价Jreturn J注意这里的X是真实数据前加了一列1因为有theta(0)2、梯度下降算法代价函数对求偏导得到所以对theta的更新可以写为其中为学习速率控制梯度下降的速度一般取0.01,0.03,0.1,0.3.....为什么梯度下降可以逐步减小代价函数假设函数f(x)泰勒展开f(x△x)f(x)f(x)*△xo(△x)令△x-α*f(x) ,即负梯度方向乘以一个很小的步长α将△x代入泰勒展开式中f(x△x)f(x)-α*[f(x)]²o(△x)可以看出α是取得很小的正数[f(x)]²也是正数所以可以得出f(x△x)f(x)所以沿着负梯度放下函数在减小多维情况一样。实现代码# 梯度下降算法def gradientDescent(X,y,theta,alpha,num_iters):m len(y)n len(theta)temp np.matrix(np.zeros((n,num_iters))) # 暂存每次迭代计算的theta转化为矩阵形式J_history np.zeros((num_iters,1)) #记录每次迭代计算的代价值for i in range(num_iters): # 遍历迭代次数h np.dot(X,theta) # 计算内积matrix可以直接乘temp[:,i] theta - ((alpha/m)*(np.dot(np.transpose(X),h-y))) #梯度的计算theta temp[:,i]J_history[i] computerCost(X,y,theta) #调用计算代价函数print .,return theta,J_history3、均值归一化目的是使数据都缩放到一个范围内便于使用梯度下降算法其中 为所有此feture数据的平均值可以是最大值-最小值也可以是这个feature对应的数据的标准差实现代码# 归一化featuredef featureNormaliza(X):X_norm np.array(X) #将X转化为numpy数组对象才可以进行矩阵的运算#定义所需变量mu np.zeros((1,X.shape[1]))sigma np.zeros((1,X.shape[1]))mu np.mean(X_norm,0) # 求每一列的平均值(0指定为列1代表行)sigma np.std(X_norm,0) # 求每一列的标准差for i in range(X.shape[1]): # 遍历列X_norm[:,i] (X_norm[:,i]-mu[i])/sigma[i] # 归一化return X_norm,mu,sigma注意预测的时候也需要均值归一化数据4、最终运行结果代价随迭代次数的变化导入包from sklearn import linear_modelfrom sklearn.preprocessing import StandardScaler #引入缩放的包归一化# 归一化操作scaler StandardScaler()scaler.fit(X)x_train scaler.transform(X)x_test scaler.transform(np.array([1650,3]))线性模型拟合# 线性模型拟合model linear_model.LinearRegression()model.fit(x_train, y)预测#预测结果result model.predict(x_test)1、代价函数可以综合起来为其中为什么不用线性回归的代价函数表示因为线性回归的代价函数可能是非凸的对于分类问题使用梯度下降很难得到最小值上面的代价函数是凸函数的图像如下即y1时可以看出当趋于1y1,与预测值一致此时付出的代价cost趋于0若趋于0y1,此时的代价cost值非常大我们最终的目的是最小化代价值同理的图像如下(y0)2、梯度同样对代价函数求偏导可以看出与线性回归的偏导数一致推到过程3、正则化目的是为了防止过拟合在代价函数中加上一项注意j是重1开始的因为theta(0)为一个常数项X中最前面一列会加上1列1所以乘积还是theta(0),feature没有关系没有必要正则化正则化后的代价# 代价函数def costFunction(initial_theta,X,y,inital_lambda):m len(y)J 0h sigmoid(np.dot(X,initial_theta)) # 计算h(z)theta1 initial_theta.copy() # 因为正则化j1从1开始不包含0所以复制一份前theta(0)值为0theta1[0] 0temp 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 # 正则化的代价方程return J正则化后的代价的梯度# 计算梯度def gradient(initial_theta,X,y,inital_lambda):m len(y)grad np.zeros((initial_theta.shape[0]))h sigmoid(np.dot(X,initial_theta))# 计算h(z)theta1 initial_theta.copy()theta1[0] 0grad np.dot(np.transpose(X),h-y)/minital_lambda/m*theta1 #正则化的梯度return grad4、S型函数(即)实现代码# S型函数def sigmoid(z):h np.zeros((len(z),1)) # 初始化与z的长度一置h 1.0/(1.0np.exp(-z))return h5、映射为多项式因为数据的feture可能很少导致偏差大所以创造出一些feture结合eg:映射为2次方的形式:实现代码# 映射为多项式def mapFeature(X1,X2):degree 3; # 映射的最高次方out np.ones((X1.shape[0],1)) # 映射后的结果数组(取代X)这里以degree2为例映射为1,x1,x2,x1^2,x1,x2,x2^2for i in np.arange(1,degree1):for j in range(i1):temp X1**(i-j)*(X2**j) #矩阵直接乘相当于matlab中的点乘.*out np.hstack((out, temp.reshape(-1,1)))return out6、使用scipy的优化方法梯度下降使用scipy中optimize中的fmin_bfgs函数调用scipy中的优化算法fmin_bfgs(拟牛顿法Broyden-Fletcher-Goldfarb-ShannocostFunction是自己实现的一个求代价的函数initial_theta表示初始化的值,fprime指定costFunction的梯度args是其余测参数以元组的形式传入最后会将最小化costFunction的theta返回result optimize.fmin_bfgs(costFunction, initial_theta, fprimegradient, args(X,y,initial_lambda))7、运行结果data1决策边界和准确度data2决策边界和准确度导入包from sklearn.linear_model import LogisticRegressionfrom sklearn.preprocessing import StandardScalerfrom sklearn.cross_validation import train_test_splitimport numpy as np划分训练集和测试集# 划分为训练集和测试集x_train,x_test,y_train,y_test train_test_split(X,y,test_size0.2)归一化# 归一化scaler StandardScaler()x_train scaler.fit_transform(x_train)x_test scaler.fit_transform(x_test)逻辑回归#逻辑回归model LogisticRegression()model.fit(x_train,y_train)预测# 预测predict model.predict(x_test)right sum(predict y_test)predict np.hstack((predict.reshape(-1,1),y_test.reshape(-1,1))) # 将预测值和真实值放在一块好观察print predictprint (测试集准确率%f%%%(right*100.0/predict.shape[0])) #计算在测试集上的准确度1、随机显示100个数字我没有使用scikit-learn中的数据集像素是20*20px彩色图如下灰度图实现代码# 显示100个数字def display_data(imgData):sum 0显示100个数(若是一个一个绘制将会非常慢可以将要画的数字整理好放到一个矩阵中显示这个矩阵即可)- 初始化一个二维数组- 将每行的数据调整成图像的矩阵放进二维数组- 显示即可pad 1display_array -np.ones((pad10*(20pad),pad10*(20pad)))for i in range(10):for j in range(10):display_array[padi*(20pad):padi*(20pad)20,padj*(20pad):padj*(20pad)20] (imgData[sum,:].reshape(20,20,orderF)) # orderF指定以列优先在matlab中是这样的python中需要指定默认以行sum 1plt.imshow(display_array,cmapgray) #显示灰度图像plt.axis(off)plt.show()2、OneVsAll如何利用逻辑回归解决多分类的问题OneVsAll就是把当前某一类看成一类其他所有类别看作一类这样有成了二分类的问题了如下图把途中的数据分成三类先把红色的看成一类把其他的看作另外一类进行逻辑回归然后把蓝色的看成一类其他的再看成一类以此类推...可以看出大于2类的情况下有多少类就要进行多少次的逻辑回归分类3、手写数字识别共有0-910个数字需要10次分类由于数据集y给出的是0,1,2...9的数字而进行逻辑回归需要0/1的label标记所以需要对y处理说一下数据集前500个是0,500-1000是1,...,所以如下图处理后的y前500行的第一列是1其余都是0,500-1000行第二列是1其余都是0....然后调用梯度下降算法求解theta实现代码# 求每个分类的theta最后返回所有的all_thetadef oneVsAll(X,y,num_labels,Lambda):# 初始化变量m,n X.shapeall_theta np.zeros((n1,num_labels)) # 每一列对应相应分类的theta,共10列X np.hstack((np.ones((m,1)),X)) # X前补上一列1的偏置biasclass_y np.zeros((m,num_labels)) # 数据的y对应0-9需要映射为0/1的关系initial_theta np.zeros((n1,1)) # 初始化一个分类的theta# 映射yfor i in range(num_labels):class_y[:,i] np.int32(yi).reshape(1,-1) # 注意reshape(1,-1)才可以赋值#np.savetxt(class_y.csv, class_y[0:600,:], delimiter,)遍历每个分类计算对应的theta值for i in range(num_labels):result optimize.fmin_bfgs(costFunction, initial_theta, fprimegradient, args(X,class_y[:,i],Lambda)) # 调用梯度下降的优化方法all_theta[:,i] result.reshape(1,-1) # 放入all_theta中all_theta np.transpose(all_theta)return all_theta4、预测之前说过预测的结果是一个概率值利用学习出来的theta代入预测的S型函数中每行的最大值就是是某个数字的最大概率所在的列号就是预测的数字的真实值,因为在分类时所有为0的将y映射在第一列为1的映射在第二列依次类推实现代码# 预测def predict_oneVsAll(all_theta,X):m X.shape[0]num_labels all_theta.shape[0]p np.zeros((m,1))X np.hstack((np.ones((m,1)),X)) #在X最前面加一列1h sigmoid(np.dot(X,np.transpose(all_theta))) #预测返回h中每一行最大值所在的列号- np.max(h, axis1)返回h中每一行的最大值(是某个数字的最大概率)- 最后where找到的最大概率所在的列号(列号即是对应的数字)p np.array(np.where(h[0,:] np.max(h, axis1)[0]))for i in np.arange(1, m):t np.array(np.where(h[i,:] np.max(h, axis1)[i]))p np.vstack((p,t))return p5、运行结果10次分类在训练集上的准确度1、导入包from scipy import io as spioimport numpy as npfrom sklearn import svmfrom sklearn.linear_model import LogisticRegression2、加载数据data loadmat_data(data_digits.mat)X data[X] # 获取X数据每一行对应一个数字20x20pxy data[y] # 这里读取mat文件y的shape(5000, 1)y np.ravel(y) # 调用sklearn需要转化成一维的(5000,)3、拟合模型model LogisticRegression()model.fit(X, y) # 拟合4、预测predict model.predict(X) #预测print u预测准确度为%f%%%np.mean(np.float64(predict y)*100)5、输出结果(在训练集上的准确度)三、BP神经网络1、神经网络model先介绍个三层的神经网络如下图所示输入层(input layer)有三个units(为补上的bias通常设为1)表示第j层的第i个激励也称为为单元unit为第j层到第j1层映射的权重矩阵就是每条边的权重所以可以得到隐含层输出层 其中S型函数也成为激励函数可以看出 为3x4的矩阵为1x4的矩阵 》j1的单元数x(j层的单元数1)2、代价函数假设最后输出的即代表输出层有K个单元 其中代表第i个单元输出与逻辑回归的代价函数差不多就是累加上每个输出(共有K个输出)3、正则化L--所有层的个数--第l层unit的个数正则化后的代价函数为共有L-1层然后是累加对应每一层的theta矩阵注意不包含加上偏置项对应的theta(0)正则化后的代价函数实现代码# 代价函数def nnCostFunction(nn_params,input_layer_size,hidden_layer_size,num_labels,X,y,Lambda):length nn_params.shape[0] # theta的中长度# 还原theta1和theta2Theta1 nn_params[0:hidden_layer_size*(input_layer_size1)].reshape(hidden_layer_size,input_layer_size1)Theta2 nn_params[hidden_layer_size*(input_layer_size1):length].reshape(num_labels,hidden_layer_size1)# np.savetxt(Theta1.csv,Theta1,delimiter,)m X.shape[0]class_y np.zeros((m,num_labels)) # 数据的y对应0-9需要映射为0/1的关系# 映射yfor i in range(num_labels):class_y[:,i] np.int32(yi).reshape(1,-1) # 注意reshape(1,-1)才可以赋值去掉theta1和theta2的第一列因为正则化时从1开始Theta1_colCount Theta1.shape[1]Theta1_x Theta1[:,1:Theta1_colCount]Theta2_colCount Theta2.shape[1]Theta2_x Theta2[:,1:Theta2_colCount]# 正则化向theta^2term 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))))正向传播,每次需要补上一列1的偏置biasa1 np.hstack((np.ones((m,1)),X))z2 np.dot(a1,np.transpose(Theta1))a2 sigmoid(z2)a2 np.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.log(h.reshape(-1,1)))np.dot(np.transpose(1-class_y.reshape(-1,1)),np.log(1-h.reshape(-1,1)))-Lambda*term/2)/mreturn np.ravel(J)4、反向传播BP上面正向传播可以计算得到J(θ),使用梯度下降法还需要求它的梯度BP反向传播的目的就是求代价函数的梯度假设4层的神经网络,记为--l层第j个单元的误差《》(向量化)没有因为对于输入没有误差因为S型函数的导数为所以上面的和可以在前向传播中计算出来反向传播计算梯度的过程为(是大写的)for i1-m:--正向传播计算(l2,3,4...L)-反向计算、...--最后即得到代价函数的梯度实现代码# 梯度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_size1)].reshape(hidden_layer_size,input_layer_size1).copy() # 这里使用copy函数否则下面修改Theta的值nn_params也会一起修改Theta2 nn_params[hidden_layer_size*(input_layer_size1):length].reshape(num_labels,hidden_layer_size1).copy()m X.shape[0]class_y np.zeros((m,num_labels)) # 数据的y对应0-9需要映射为0/1的关系# 映射yfor i in range(num_labels):class_y[:,i] np.int32(yi).reshape(1,-1) # 注意reshape(1,-1)才可以赋值去掉theta1和theta2的第一列因为正则化时从1开始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_grad np.zeros((Theta2.shape)) #第二层到第三层的权重正向传播每次需要补上一列1的偏置biasa1 np.hstack((np.ones((m,1)),X))z2 np.dot(a1,np.transpose(Theta1))a2 sigmoid(z2)a2 np.hstack((np.ones((m,1)),a2))z3 np.dot(a2,np.transpose(Theta2))h sigmoid(z3)反向传播delta为误差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,:]) # 均方误差的误差率delta3[i,:] h[i,:]-class_y[i,:] # 交叉熵误差率Theta2_grad Theta2_gradnp.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_gradnp.dot(np.transpose(delta2[i,:].reshape(1,-1)),a1[i,:].reshape(1,-1))Theta1[:,0] 0Theta2[:,0] 0梯度grad (np.vstack((Theta1_grad.reshape(-1,1),Theta2_grad.reshape(-1,1)))Lambda*np.vstack((Theta1.reshape(-1,1),Theta2.reshape(-1,1))))/mreturn np.ravel(grad)5、BP可以求梯度的原因实际是利用了链式求导法则因为下一层的单元利用上一层的单元作为输入进行计算大体的推导过程如下最终我们是想预测函数与已知的y非常接近求均方差的梯度沿着此梯度方向可使代价函数最小化。可对照上面求梯度的过程。求误差更详细的推导过程6、梯度检查检查利用BP求的梯度是否正确利用导数的定义验证求出来的数值梯度应该与BP求出的梯度非常接近验证BP正确后就不需要再执行验证梯度的算法了实现代码# 检验梯度是否计算正确# 检验梯度是否计算正确def checkGradient(Lambda 0):构造一个小型的神经网络验证因为数值法计算梯度很浪费时间而且验证正确后之后就不再需要验证了input_layer_size 3hidden_layer_size 5num_labels 3m 5initial_Theta1 debugInitializeWeights(input_layer_size,hidden_layer_size);initial_Theta2 debugInitializeWeights(hidden_layer_size,num_labels)X debugInitializeWeights(input_layer_size-1,m)y 1np.transpose(np.mod(np.arange(1,m1), num_labels))# 初始化yy y.reshape(-1,1)nn_params np.vstack((initial_Theta1.reshape(-1,1),initial_Theta2.reshape(-1,1))) #展开thetaBP求出梯度grad nnGradient(nn_params, input_layer_size, hidden_layer_size,num_labels, X, y, Lambda)使用数值法计算梯度num_grad np.zeros((nn_params.shape[0]))step np.zeros((nn_params.shape[0]))e 1e-4for i in range(nn_params.shape[0]):step[i] eloss1 nnCostFunction(nn_params-step.reshape(-1,1), input_layer_size, hidden_layer_size,num_labels, X, y,Lambda)loss2 nnCostFunction(nn_paramsstep.reshape(-1,1), input_layer_size, hidden_layer_size,num_labels, X, y,Lambda)num_grad[i] (loss2-loss1)/(2*e)step[i]0# 显示两列比较res np.hstack((num_grad.reshape(-1,1),grad.reshape(-1,1)))print res7、权重的随机初始化神经网络不能像逻辑回归那样初始化theta为0,因为若是每条边的权重都为0每个神经元都是相同的输出在反向传播中也会得到同样的梯度最终只会预测一种结果。所以应该初始化为接近0的数实现代码# 随机初始化权重thetadef randInitializeWeights(L_in,L_out):W np.zeros((L_out,1L_in)) # 对应theta的权重epsilon_init (6.0/(L_outL_in))**0.5W np.random.rand(L_out,1L_in)*2*epsilon_init-epsilon_init # np.random.rand(L_out,1L_in)产生L_out*(1L_in)大小的随机矩阵return W8、预测正向传播预测结果实现代码# 预测def predict(Theta1,Theta2,X):m X.shape[0]num_labels Theta2.shape[0]#p np.zeros((m,1))正向传播预测结果X np.hstack((np.ones((m,1)),X))h1 sigmoid(np.dot(X,np.transpose(Theta1)))h1 np.hstack((np.ones((m,1)),h1))h2 sigmoid(np.dot(h1,np.transpose(Theta2)))返回h中每一行最大值所在的列号- np.max(h, axis1)返回h中每一行的最大值(是某个数字的最大概率)- 最后where找到的最大概率所在的列号(列号即是对应的数字)#np.savetxt(h2.csv,h2,delimiter,)p np.array(np.where(h2[0,:] np.max(h2, axis1)[0]))for i in np.arange(1, m):t np.array(np.where(h2[i,:] np.max(h2, axis1)[i]))p np.vstack((p,t))return p9、输出结果梯度检查随机显示100个手写数字显示theta1权重训练集预测准确度归一化后训练集预测准确度四、SVM支持向量机1、代价函数在逻辑回归中我们的代价为其中如图所示如果y1cost代价函数如图所示我们想让即z0这样的话cost代价函数才会趋于最小(这是我们想要的)所以用途中红色的函数代替逻辑回归中的cost当y0时同样用代替最终得到的代价函数为最后我们想要之前我们逻辑回归中的代价函数为可以认为这里的只是表达形式问题这里C的值越大SVM的决策边界的margin也越大下面会说明2、Large Margin如下图所示,SVM分类会使用最大的margin将其分开先说一下向量内积表示u的欧几里得范数(欧式范数)向量V在向量u上的投影的长度记为p则向量内积根据向量夹角公式推导一下即可前面说过当C越大时margin也就越大我们的目的是最小化代价函数J(θ),当margin最大时C的乘积项要很小所以近似为我们最后的目的就是求使代价最小的θ由可以得到p即为x在θ上的投影如下图所示假设决策边界如图找其中的一个点到θ上的投影为p,则或者若是p很小则需要很大这与我们要求的θ使最小相违背所以最后求的是large margin3、SVM Kernel(核函数)对于线性可分的问题使用线性核函数即可对于线性不可分的问题在逻辑回归中我们是将feature映射为使用多项式的形式SVM中也有多项式核函数但是更常用的是高斯核函数也称为RBF核高斯核函数为假设如图几个点令...可以看出若是x与距离较近》(即相似度较大)若是x与距离较远》(即相似度较低)高斯核函数的σ越小f下降的越快如何选择初始的训练集选择对于给出的x计算f,令所以最小化J求出θ如果》预测y14、使用scikit-learn中的SVM模型代码线性可分的,指定核函数为lineardata1——线性分类data1 spio.loadmat(data1.mat)X data1[X]y data1[y]y np.ravel(y)plot_data(X,y)model svm.SVC(C1.0,kernellinear).fit(X,y) # 指定核函数为线性核函数非线性可分的默认核函数为rbfdata2——非线性分类data2 spio.loadmat(data2.mat)X data2[X]y data2[y]y np.ravel(y)plt plot_data(X,y)plt.show()model svm.SVC(gamma100).fit(X,y) # gamma为核函数的系数值越大拟合的越好5、运行结果线性可分的决策边界线性不可分的决策边界五、K-Means聚类算法1、聚类过程聚类属于无监督学习不知道y的标记分为K类K-Means算法分为两个步骤第一步簇分配随机选K个点作为中心计算到这K个点的距离分为K个簇第二步移动聚类中心重新计算每个簇的中心移动中心重复以上步骤。如下图所示随机分配的聚类中心重新计算聚类中心移动一次最后10步之后的聚类中心计算每条数据到哪个中心最近实现代码# 找到每条数据距离哪个类中心最近def findClosestCentroids(X,initial_centroids):m X.shape[0] # 数据条数K initial_centroids.shape[0] # 类的总数dis np.zeros((m,K)) # 存储计算每个点分别到K个类的距离idx np.zeros((m,1)) # 要返回的每条数据属于哪个类计算每个点到每个类中心的距离for i in range(m):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))返回dis每一行的最小值对应的列号即为对应的类别- np.min(dis, axis1)返回每一行的最小值- np.where(dis np.min(dis, axis1).reshape(-1,1)) 返回对应最小值的坐标- 注意可能最小值对应的坐标有多个where都会找出来所以返回时返回前m个需要的即可(因为对于多个最小值属于哪个类别都可以)dummy,idx np.where(dis np.min(dis, axis1).reshape(-1,1))return idx[0:dis.shape[0]] # 注意截取一下计算类中心实现代码# 计算类中心def computerCentroids(X,idx,K):n X.shape[1]centroids np.zeros((K,n))for i in range(K):centroids[i,:] np.mean(X[np.ravel(idxi),:], axis0).reshape(1,-1) # 索引要是一维的,axis0为每一列idxi一次找出属于哪一类的然后计算均值return centroids2、目标函数也叫做失真代价函数最后我们想得到其中表示第i条数据距离哪个类中心最近其中即为聚类的中心3、聚类中心的选择随机初始化从给定的数据中随机抽取K个作为聚类中心随机一次的结果可能不好可以随机多次最后取使代价函数最小的作为中心实现代码(这里随机一次)# 初始化类中心--随机取K个点作为聚类中心def kMeansInitCentroids(X,K):m X.shape[0]m_arr np.arange(0,m) # 生成0-m-1centroids np.zeros((K,X.shape[1]))np.random.shuffle(m_arr) # 打乱m_arr顺序rand_indices m_arr[:K] # 取前K个centroids X[rand_indices,:]return centroids4、聚类个数K的选择聚类是不知道y的label的所以不知道真正的聚类个数肘部法则(Elbow method)作代价函数J和K的图若是出现一个拐点如下图所示K就取拐点处的值下图此时K3若是很平滑就不明确人为选择。第二种就是人为观察选择5、应用——图片压缩将图片的像素分为若干类然后用这个类代替原来的像素值执行聚类的算法代码# 聚类算法def runKMeans(X,initial_centroids,max_iters,plot_process):m,n X.shape # 数据条数和维度K initial_centroids.shape[0] # 类数centroids initial_centroids # 记录当前类中心previous_centroids centroids # 记录上一次类中心idx np.zeros((m,1)) # 每条数据属于哪个类for i in range(max_iters): # 迭代次数print u迭代计算次数%d%(i1)idx findClosestCentroids(X, centroids)if plot_process: # 如果绘制图像plt plotProcessKMeans(X,centroids,previous_centroids) # 画聚类中心的移动过程previous_centroids centroids # 重置centroids computerCentroids(X, idx, K) # 重新计算类中心if plot_process: # 显示最终的绘制结果plt.show()return centroids,idx # 返回聚类中心和数据属于哪个类导入包from sklearn.cluster import KMeans使用模型拟合数据model KMeans(n_clusters3).fit(X) # n_clusters指定3类拟合数据聚类中心centroids model.cluster_centers_ # 聚类中心7、运行结果二维数据类中心的移动图片压缩六、PCA主成分分析(降维)1、用处数据压缩(Data Compression),使程序运行更快可视化数据例如3D--2D等......2、2D--1DnD--kD如下图所示所有数据点可以投影到一条直线是投影距离的平方和(投影误差)最小注意数据需要归一化处理思路是找1个向量u,所有数据投影到上面使投影距离最小那么nD--kD就是找k个向量所有数据投影到上面使投影误差最小eg:3D--2D,2个向量就代表一个平面了所有点投影到这个平面的投影误差最小即可3、主成分分析PCA与线性回归的区别线性回归是找x与y的关系然后用于预测yPCA是找一个投影面最小化data到这个投影面的投影误差4、PCA降维过程数据预处理(均值归一化)公式就是减去对应feature的均值然后除以对应特征的标准差(也可以是最大值-最小值)实现代码# 归一化数据def featureNormalize(X):(每一个数据-当前列的均值)/当前列的标准差n X.shape[1]mu np.zeros((1,n));sigma np.zeros((1,n))mu np.mean(X,axis0)sigma np.std(X,axis0)for i in range(n):X[:,i] (X[:,i]-mu[i])/sigma[i]return X,mu,sigma计算协方差矩阵Σ(Covariance Matrix)注意这里的Σ和求和符号不同协方差矩阵对称正定(不理解正定的看看线代)大小为nxn,n为feature的维度实现代码Sigma np.dot(np.transpose(X_norm),X_norm)/m # 求Sigma计算Σ的特征值和特征向量可以是用svd奇异值分解函数U,S,V svd(Σ)返回的是与Σ同样大小的对角阵S(由Σ的特征值组成)[注意matlab中函数返回的是对角阵在python中返回的是一个向量节省空间]还有两个酉矩阵U和V且注意svd函数求出的S是按特征值降序排列的若不是使用svd,需要按特征值大小重新排列U降维选取U中的前K列(假设要降为K维)Z就是对应降维之后的数据实现代码# 映射数据def projectData(X_norm,U,K):Z np.zeros((X_norm.shape[0],K))U_reduce U[:,0:K] # 取前K个Z np.dot(X_norm,U_reduce)return Z过程总结Sigma X*X/mU,S,V svd(Sigma)Ureduce U[:,0:k]Z Ureduce*x5、数据恢复因为所以 (注意这里是X的近似值)又因为Ureduce为正定矩阵【正定矩阵满足所以】所以这里实现代码# 恢复数据def recoverData(Z,U,K):X_rec np.zeros((Z.shape[0],U.shape[0]))U_recude U[:,0:K]X_rec np.dot(Z,np.transpose(U_recude)) # 还原数据(近似)return X_rec6、主成分个数的选择(即要降的维度)如何选择投影误差(project error)总变差(total variation):若误差率(error ratio)则称99%保留差异性误差率一般取1%5%10%等如何实现若是一个个试的话代价太大之前U,S,V svd(Sigma),我们得到了S这里误差率error ratio:可以一点点增加K尝试。7、使用建议不要使用PCA去解决过拟合问题Overfitting还是使用正则化的方法(如果保留了很高的差异性还是可以的)只有在原数据上有好的结果但是运行很慢才考虑使用PCA8、运行结果2维数据降为1维要投影的方向2D降为1D及对应关系人脸数据降维原始数据可视化部分U矩阵信息恢复数据导入需要的包#-*- coding: utf-8 -*-# Author:bob# Date:2016.12.22import numpy as npfrom matplotlib import pyplot as pltfrom scipy import io as spiofrom sklearn.decomposition import pcafrom sklearn.preprocessing import StandardScaler归一化数据归一化数据并作图scaler StandardScaler()scaler.fit(X)x_train scaler.transform(X)使用PCA模型拟合数据并降维n_components对应要将的维度拟合数据K1 # 要降的维度model pca.PCA(n_componentsK).fit(x_train) # 拟合数据n_components定义要降的维度Z model.transform(x_train) # transform就会执行降维操作数据恢复model.components_会得到降维使用的U矩阵数据恢复并作图Ureduce model.components_ # 得到降维用的Ureducex_rec np.dot(Z,Ureduce) # 数据恢复七、异常检测 Anomaly Detection1、高斯分布(正态分布)Gaussian distribution分布函数其中u为数据的均值σ为数据的标准差σ越小对应的图像越尖参数估计(parameter estimation)2、异常检测算法例子训练集,其中假设相互独立建立model模型过程选择具有代表异常的feature:xi参数估计计算p(x),若是P(x)这里只是单元高斯分布假设了feature之间是独立的下面会讲到多元高斯分布会自动捕捉到feature之间的关系参数估计实现代码# 参数估计函数(就是求均值和方差)def estimateGaussian(X):m,n X.shapemu np.zeros((n,1))sigma2 np.zeros((n,1))mu np.mean(X, axis0) # axis0表示列每列的均值sigma2 np.var(X,axis0) # 求每列的方差return mu,sigma23、评价p(x)的好坏以及ε的选取对偏斜数据的错误度量因为数据可能是非常偏斜的(就是y1的个数非常少(y1表示异常))所以可以使用Precision/Recall计算F1Score(在CV交叉验证集上)例如预测癌症假设模型可以得到99%能够预测正确1%的错误率但是实际癌症的概率很小只有0.5%那么我们始终预测没有癌症y0反而可以得到更小的错误率。使用error rate来评估就不科学了。如下图记录 即正确预测正样本/所有预测正样本 即正确预测正样本/真实值为正样本总是让y1(较少的类)计算Precision和Recall还是以癌症预测为例假设预测都是no-cancerTN199FN1TP0FP0所以Precision0/0Recall0/10尽管accuracy199/20099.5%但是不可信。ε的选取尝试多个ε值使F1Score的值高实现代码# 选择最优的epsilon即使F1Score最大def selectThreshold(yval,pval):初始化所需变量bestEpsilon 0.bestF1 0.F1 0.step (np.max(pval)-np.min(pval))/1000计算for epsilon in np.arange(np.min(pval),np.max(pval),step):cvPrecision pvaltp np.sum((cvPrecision 1) (yval 1).ravel()).astype(float) # sum求和是int型的需要转为floatfp np.sum((cvPrecision 1) (yval 0).ravel()).astype(float)fn np.sum((cvPrecision 0) (yval 1).ravel()).astype(float)precision tp/(tpfp) # 精准度recision tp/(tpfn) # 召回率F1 (2*precision*recision)/(precisionrecision) # F1Score计算公式if F1 bestF1: # 修改最优的F1 ScorebestF1 F1bestEpsilon epsilonreturn bestEpsilon,bestF14、选择使用什么样的feature(单元高斯分布)如果一些数据不是满足高斯分布的可以变化一下数据例如log(xC),x^(1/2)等如果p(x)的值无论异常与否都很大可以尝试组合多个feature,(因为feature之间可能是有关系的)5、多元高斯分布单元高斯分布存在的问题如下图红色的点为异常点其他的都是正常点(比如CPU和memory的变化)x1对应的高斯分布如下x2对应的高斯分布如下可以看出对应的p(x1)和p(x2)的值变化并不大就不会认为异常因为我们认为feature之间是相互独立的所以如上图是以正圆的方式扩展多元高斯分布并不是建立p(x1),p(x2)...p(xn)而是统一建立p(x)其中参数,Σ为协方差矩阵同样|Σ|越小p(x)越尖例如表示x1,x2正相关即x1越大x2也就越大如下图也就可以将红色的异常点检查出了若表示x1,x2负相关实现代码# 多元高斯分布函数def multivariateGaussian(X,mu,Sigma2):k len(mu)if (Sigma2.shape[0]1):Sigma2 np.diag(Sigma2)多元高斯分布函数X X-muargu (2*np.pi)**(-k/2)*np.linalg.det(Sigma2)**(-0.5)p argu*np.exp(-0.5*np.sum(np.dot(X,np.linalg.inv(Sigma2))*X,axis1)) # axis表示每行return p6、单元和多元高斯分布特点单元高斯分布人为可以捕捉到feature之间的关系时可以使用计算量小多元高斯分布自动捕捉到相关的feature计算量大因为mn或Σ可逆时可以使用。(若不可逆可能有冗余的x因为线性相关不可逆或者就是m7、程序运行结果显示数据等高线异常点标注
