机器学习算法Python实现

目录

一、线性回归

1、代价函数

# 计算代价函数
def computerCost(X,y,theta):
 m = len(y)
 J = 0
 
 J = (np.transpose(X*theta-y))*(X*theta-y)/(2*m) #计算代价J
 return J

2、梯度下降算法

# 梯度下降算法
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_history 

3、均值归一化

# 归一化feature
def 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、最终运行结果

5、使用scikit-learn库中的线性模型实现

from sklearn import linear_model
from 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、代价函数

2、梯度

3、正则化

# 代价函数
def costFunction(initial_theta,X,y,inital_lambda):
 m = len(y)
 J = 0
 
 h = sigmoid(np.dot(X,initial_theta)) # 计算h(z)
 theta1 = initial_theta.copy() # 因为正则化j=1从1开始,不包含0,所以复制一份,前theta(0)值为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 # 正则化的代价方程
 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] = 0
 grad = np.dot(np.transpose(X),h-y)/m+inital_lambda/m*theta1 #正则化的梯度
 return grad 

4、S型函数(即{{h_\theta }(x)})

# S型函数 
def sigmoid(z):
 h = np.zeros((len(z),1)) # 初始化,与z的长度一置
 
 h = 1.0/(1.0+np.exp(-z))
 return h

5、映射为多项式

# 映射为多项式 
def mapFeature(X1,X2):
 degree = 3; # 映射的最高次方
 out = np.ones((X1.shape[0],1)) # 映射后的结果数组(取代X)
 '''
 这里以degree=2为例,映射为1,x1,x2,x1^2,x1,x2,x2^2
 '''
 for i in np.arange(1,degree+1): 
 for j in range(i+1):
 temp = X1**(i-j)*(X2**j) #矩阵直接乘相当于matlab中的点乘.*
 out = np.hstack((out, temp.reshape(-1,1)))
 return out

6、使用scipy的优化方法

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

7、运行结果

8、使用scikit-learn库中的逻辑回归模型实现

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

 # 划分为训练集和测试集
 x_train,x_test,y_train,y_test = train_test_split(X,y,test_size=0.2)

 # 归一化
 scaler = StandardScaler()
 scaler.fit(x_train)
 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 predict
 print ('测试集准确率:%f%%'%(right*100.0/predict.shape[0])) #计算在测试集上的准确度

逻辑回归_手写数字识别_OneVsAll

1、随机显示100个数字

# 显示100个数字
def display_data(imgData):
 sum = 0
 '''
 显示100个数(若是一个一个绘制将会非常慢,可以将要画的数字整理好,放到一个矩阵中,显示这个矩阵即可)
 - 初始化一个二维数组
 - 将每行的数据调整成图像的矩阵,放进二维数组
 - 显示即可
 '''
 pad = 1
 display_array = -np.ones((pad+10*(20+pad),pad+10*(20+pad)))
 for i in range(10):
 for j in range(10):
 display_array[pad+i*(20+pad):pad+i*(20+pad)+20,pad+j*(20+pad):pad+j*(20+pad)+20] = (imgData[sum,:].reshape(20,20,order="F")) # order=F指定以列优先,在matlab中是这样的,python中需要指定,默认以行
 sum += 1
 
 plt.imshow(display_array,cmap='gray') #显示灰度图像
 plt.axis('off')
 plt.show()

2、OneVsAll

3、手写数字识别

# 求每个分类的theta,最后返回所有的all_theta 
def oneVsAll(X,y,num_labels,Lambda):
 # 初始化变量
 m,n = X.shape
 all_theta = np.zeros((n+1,num_labels)) # 每一列对应相应分类的theta,共10列
 X = np.hstack((np.ones((m,1)),X)) # X前补上一列1的偏置bias
 class_y = np.zeros((m,num_labels)) # 数据的y对应0-9,需要映射为0/1的关系
 initial_theta = np.zeros((n+1,1)) # 初始化一个分类的theta
 
 # 映射y
 for i in range(num_labels):
 class_y[:,i] = np.int32(y==i).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, fprime=gradient, args=(X,class_y[:,i],Lambda)) # 调用梯度下降的优化方法
 all_theta[:,i] = result.reshape(1,-1) # 放入all_theta中
 
 all_theta = np.transpose(all_theta) 
 return all_theta

4、预测

# 预测
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最前面加一列1
 
 h = sigmoid(np.dot(X,np.transpose(all_theta))) #预测
 '''
 返回h中每一行最大值所在的列号
 - np.max(h, axis=1)返回h中每一行的最大值(是某个数字的最大概率)
 - 最后where找到的最大概率所在的列号(列号即是对应的数字)
 '''
 p = np.array(np.where(h[0,:] == np.max(h, axis=1)[0])) 
 for i in np.arange(1, m):
 t = np.array(np.where(h[i,:] == np.max(h, axis=1)[i]))
 p = np.vstack((p,t))
 return p

5、运行结果

6、使用scikit-learn库中的逻辑回归模型实现

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

 data = loadmat_data("data_digits.mat") 
 X = data['X'] # 获取X数据,每一行对应一个数字20x20px
 y = data['y'] # 这里读取mat文件y的shape=(5000, 1)
 y = np.ravel(y) # 调用sklearn需要转化成一维的(5000,)

 model = LogisticRegression()
 model.fit(X, y) # 拟合

 predict = model.predict(X) #预测
 
 print u"预测准确度为:%f%%"%np.mean(np.float64(predict == y)*100)

三、BP神经网络

1、神经网络model

2、代价函数

3、正则化

# 代价函数
def nnCostFunction(nn_params,input_layer_size,hidden_layer_size,num_labels,X,y,Lambda):
 length = nn_params.shape[0] # theta的中长度
 # 还原theta1和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)) # 数据的y对应0-9,需要映射为0/1的关系
 # 映射y
 for i in range(num_labels):
 class_y[:,i] = np.int32(y==i).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^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))))
 
 '''正向传播,每次需要补上一列1的偏置bias'''
 a1 = 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)/m 
 
 return np.ravel(J)

4、反向传播BP

# 梯度
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)].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)
 m = X.shape[0]
 class_y = np.zeros((m,num_labels)) # 数据的y对应0-9,需要映射为0/1的关系 
 # 映射y
 for i in range(num_labels):
 class_y[:,i] = np.int32(y==i).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)) #第二层到第三层的权重
 
 Theta1[:,0] = 0;
 Theta2[:,0] = 0;
 '''正向传播,每次需要补上一列1的偏置bias'''
 a1 = 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,:]
 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))
 
 '''梯度'''
 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)

5、BP可以求梯度的原因

6、梯度检查

# 检验梯度是否计算正确
# 检验梯度是否计算正确
def checkGradient(Lambda = 0):
 '''构造一个小型的神经网络验证,因为数值法计算梯度很浪费时间,而且验证正确后之后就不再需要验证了'''
 input_layer_size = 3
 hidden_layer_size = 5
 num_labels = 3
 m = 5
 initial_Theta1 = debugInitializeWeights(input_layer_size,hidden_layer_size); 
 initial_Theta2 = debugInitializeWeights(hidden_layer_size,num_labels)
 X = debugInitializeWeights(input_layer_size-1,m)
 y = 1+np.transpose(np.mod(np.arange(1,m+1), num_labels))# 初始化y
 
 y = y.reshape(-1,1)
 nn_params = np.vstack((initial_Theta1.reshape(-1,1),initial_Theta2.reshape(-1,1))) #展开theta 
 '''BP求出梯度'''
 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-4
 for i in range(nn_params.shape[0]):
 step[i] = e
 loss1 = nnCostFunction(nn_params-step.reshape(-1,1), input_layer_size, hidden_layer_size, 
 num_labels, X, y, 
 Lambda)
 loss2 = nnCostFunction(nn_params+step.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 res

7、权重的随机初始化

# 随机初始化权重theta
def randInitializeWeights(L_in,L_out):
 W = np.zeros((L_out,1+L_in)) # 对应theta的权重
 epsilon_init = (6.0/(L_out+L_in))**0.5
 W = np.random.rand(L_out,1+L_in)*2*epsilon_init-epsilon_init # np.random.rand(L_out,1+L_in)产生L_out*(1+L_in)大小的随机矩阵
 return W

8、预测

# 预测
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, axis=1)返回h中每一行的最大值(是某个数字的最大概率)
 - 最后where找到的最大概率所在的列号(列号即是对应的数字)
 '''
 #np.savetxt("h2.csv",h2,delimiter=',')
 p = np.array(np.where(h2[0,:] == np.max(h2, axis=1)[0])) 
 for i in np.arange(1, m):
 t = np.array(np.where(h2[i,:] == np.max(h2, axis=1)[i]))
 p = np.vstack((p,t))
 return p 

9、输出结果

四、SVM支持向量机

1、代价函数

2、Large Margin

3、SVM Kernel(核函数)

4、使用scikit-learn中的SVM模型代码

 '''data1——线性分类'''
 data1 = spio.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) # 指定核函数为线性核函数

 '''data2——非线性分类'''
 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(gamma=100).fit(X,y) # gamma为核函数的系数,值越大拟合的越好

5、运行结果

五、K-Means聚类算法

1、聚类过程

# 找到每条数据距离哪个类中心最近 
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, axis=1)返回每一行的最小值
 - np.where(dis == np.min(dis, axis=1).reshape(-1,1)) 返回对应最小值的坐标
 - 注意:可能最小值对应的坐标有多个,where都会找出来,所以返回时返回前m个需要的即可(因为对于多个最小值,属于哪个类别都可以)
 ''' 
 dummy,idx = np.where(dis == np.min(dis, axis=1).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(idx==i),:], axis=0).reshape(1,-1) # 索引要是一维的,axis=0为每一列,idx==i一次找出属于哪一类的,然后计算均值
 return centroids

2、目标函数

3、聚类中心的选择

# 初始化类中心--随机取K个点作为聚类中心
def kMeansInitCentroids(X,K):
 m = X.shape[0]
 m_arr = np.arange(0,m) # 生成0-m-1
 centroids = np.zeros((K,X.shape[1]))
 np.random.shuffle(m_arr) # 打乱m_arr顺序 
 rand_indices = m_arr[:K] # 取前K个
 centroids = X[rand_indices,:]
 return centroids

4、聚类个数K的选择

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'%(i+1)
 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 # 返回聚类中心和数据属于哪个类

6、使用scikit-learn库中的线性模型实现聚类

 from sklearn.cluster import KMeans

 model = KMeans(n_clusters=3).fit(X) # n_clusters指定3类,拟合数据

 centroids = model.cluster_centers_ # 聚类中心

7、运行结果

六、PCA主成分分析(降维)

1、用处

2、2D-->1D,nD-->kD

3、主成分分析PCA与线性回归的区别

4、PCA降维过程

 # 归一化数据
 def featureNormalize(X):
 '''(每一个数据-当前列的均值)/当前列的标准差'''
 n = X.shape[1]
 mu = np.zeros((1,n));
 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,sigma

Sigma = np.dot(np.transpose(X_norm),X_norm)/m # 求Sigma

 # 映射数据
 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

5、数据恢复

 # 恢复数据 
 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_rec

6、主成分个数的选择(即要降的维度)

7、使用建议

8、运行结果

9、使用scikit-learn库中的PCA实现降维

#-*- coding: 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 StandardScaler

 '''归一化数据并作图'''
 scaler = StandardScaler()
 scaler.fit(X)
 x_train = scaler.transform(X)

 '''拟合数据'''
 K=1 # 要降的维度
 model = pca.PCA(n_components=K).fit(x_train) # 拟合数据,n_components定义要降的维度
 Z = model.transform(x_train) # transform就会执行降维操作

 '''数据恢复并作图'''
 Ureduce = model.components_ # 得到降维用的Ureduce
 x_rec = np.dot(Z,Ureduce) # 数据恢复

七、异常检测 Anomaly Detection

1、高斯分布(正态分布)Gaussian distribution

2、异常检测算法

# 参数估计函数(就是求均值和方差)
def estimateGaussian(X):
 m,n = X.shape
 mu = np.zeros((n,1))
 sigma2 = np.zeros((n,1))
 
 mu = np.mean(X, axis=0) # axis=0表示列,每列的均值
 sigma2 = np.var(X,axis=0) # 求每列的方差
 return mu,sigma2

3、评价p(x)的好坏,以及ε的选取

# 选择最优的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 = pval<epsilon
 tp = np.sum((cvPrecision == 1) & (yval == 1)).astype(float) # sum求和是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) # 召回率
 F1 = (2*precision*recision)/(precision+recision) # F1Score计算公式
 if F1 > bestF1: # 修改最优的F1 Score
 bestF1 = F1
 bestEpsilon = epsilon
 return bestEpsilon,bestF1

4、选择使用什么样的feature(单元高斯分布)

5、多元高斯分布

# 多元高斯分布函数 
def multivariateGaussian(X,mu,Sigma2):
 k = len(mu)
 if (Sigma2.shape[0]>1):
 Sigma2 = np.diag(Sigma2)
 '''多元高斯分布函数''' 
 X = X-mu
 argu = (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,axis=1)) # axis表示每行
 return p

6、单元和多元高斯分布特点

7、程序运行结果