网站导航
> 技术文档 > 数学建模常用30个算法——Python代码（一）_python学习资料及常用模型算法代码

数学建模常用30个算法——Python代码（一）_python学习资料及常用模型算法代码

网友: 技术文档 2025-08-08 02:53:44

数学建模常用算法

  • 数学建模常用30个算法——Python代码
  • 1 时间序列分析
    • ARIMA时间序列预测模型
    • 灰色预测模型
  • 2. 机器学习/深度学习
    • 监督学习
      • BP神经网络模型
      • 卷积神经网络（CNN）模型
      • 决策树分类模型
      • 逻辑回归模型
      • 随机森林分类模型
      • 支持向量机（SVM）模型
    • 无监督学习
      • K-means聚类模型
      • 主成分分析（PCA）算法

数学建模常用30个算法——Python代码

1 时间序列分析

ARIMA时间序列预测模型

用于分析和预测时间序列数据，结合自回归（AR）、差分（I）和移动平均（MA）来建模趋势和季节性。

 import pandas # 读取数据，指定日期为索引列 data = pandas.read_csv( \'D:\\\\DATA\\\\pycase\\\\number2\\\\9.3\\\\Data.csv\' , index_col=\'日期\') # 绘图过程中 import matplotlib.pyplot as plt # 用来正常显示中文标签 plt.rcParams[\'font.sans-serif\']=[\'SimHei\'] # 用来正常显示负号 plt.rcParams[\'axes.unicode_minus\'] = False # 查看趋势图data.plot() #有增长趋势，不平稳 # 附加：查看自相关系数合片自相关系数（查分之后），可以用于平稳性的检测，也可用于定阶系数预估 #自相关图（） from statsmodels.graphics.tsaplots import plot_acf plot_acf(data).show() #自相关图既不是拖尾也不是截尾。以上的图的自相关是一个三角对称的形式，这种趋势是单调趋势的典型图形，说明这个序列不是平稳序列 # 1 平稳性检测 from statsmodels.tsa.stattools import adfuller as ADF def tagADF(t): result = pandas.DataFrame(index=[ \"Test Statistic Value\", \"p-value\", \"Lags Used\", \"Number of Observations Used\", \"Critical Value(1%)\", \"Critical Value(5%)\", \"Critical Value(10%)\" ], columns=[\'销量\'] ); result[\'销量\'][\'Test Statistic Value\'] = t[0] result[\'销量\'][\'p-value\'] = t[1] result[\'销量\'][\'Lags Used\'] = t[2] result[\'销量\'][\'Number of Observations Used\'] = t[3] result[\'销量\'][\'Critical Value(1%)\'] = t[4][\'1%\'] result[\'销量\'][\'Critical Value(5%)\'] = t[4][\'5%\'] result[\'销量\'][\'Critical Value(10%)\'] = t[4][\'10%\'] return result; print(\'原始序列的ADF检验结果为:\',tagADF(ADF(data[u\'销量\']))) # 添加标签后展现 # 平稳判断：得到统计量大于三个置信度(1%,5%,10%)临界统计值，p值显著大于0.05，该序列为非平稳序列。# 备注：得到的统计量显著小于3个置信度（1%，5%，10%）的临界统计值时，为平稳 此时p值接近于0 此处不为0，尝试增加数据量，原数据太少 # 2 进行数据差分，一般一阶差分就可以 D_data = data.diff(1).dropna()D_data.columns = [u\'销量差分\'] #差分图趋势查看 D_data.plot() plt.show() # 附加：查看自相关系数合片自相关系数（查分之后），可以用于平稳性的检测，也可用于定阶系数预估 #自相关图 plot_acf(D_data).show() plt.show() #偏自相关图 from statsmodels.graphics.tsaplots import plot_pacf plot_pacf(D_data).show() # 3 平稳性检测 print(u\'差分序列的ADF检验结果为：\', tagADF(ADF(D_data[u\'销量差分\']))) # 解释：Test Statistic Value值小于两个水平值，p值显著小于0.05，一阶差分后序列为平稳序列。 # 4 白噪声检验from statsmodels.stats.diagnostic import acorr_ljungbox #返回统计量和p值 print(u\'差分序列的白噪声检验结果为：\', acorr_ljungbox(D_data, lags=1)) # 分别为stat值（统计量）和P值 # P值小于0.05，所以一阶差分后的序列为平稳非白噪声序列。 # 5 p，q定阶 from statsmodels.tsa.arima_model import ARIMA #一般阶数不超过length/10 pmax = int(len(D_data)/10) #一般阶数不超过length/10 qmax = int(len(D_data)/10) #bic矩阵 bic_matrix = [] for p in range(pmax+1): tmp = [] for q in range(qmax+1):#存在部分报错，所以用try来跳过报错。 try: tmp.append(ARIMA(data, (p,1,q)).fit().bic) except: tmp.append(None) bic_matrix.append(tmp) #从中可以找出最小值 bic_matrix = pandas.DataFrame(bic_matrix) #先用stack展平，然后用idxmin找出最小值位置。 p,q = bic_matrix.stack().idxmin() print(u\'BIC最小的p值和q值为：%s、%s\' %(p,q))# 取BIC信息量达到最小的模型阶数，结果p为0，q为1，定阶完成。 # 6 建立模型和预测 model = ARIMA(data, (p,1,q)).fit() #给出一份模型报告 model.summary2() #作为期5天的预测，返回预测结果、标准误差、置信区间。 model.forecast(5)

灰色预测模型

适用于小样本数据预测，通过灰色系统理论对不完全信息进行建模，常用于短期趋势分析。

# -*- coding: utf-8 -*-\"\"\"Spyder EditorThis is a temporary script file.\"\"\"import numpy as npimport mathhistory_data = [724.57,746.62,778.27,800.8,827.75,871.1,912.37,954.28,995.01,1037.2]n = len(history_data)X0 = np.array(history_data)#累加生成history_data_agg = [sum(history_data[0:i+1]) for i in range(n)]X1 = np.array(history_data_agg)#计算数据矩阵B和数据向量YB = np.zeros([n-1,2])Y = np.zeros([n-1,1])for i in range(0,n-1): B[i][0] = -0.5*(X1[i] + X1[i+1]) B[i][1] = 1 Y[i][0] = X0[i+1]#计算GM(1,1)微分方程的参数a和u#A = np.zeros([2,1])A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)a = A[0][0]u = A[1][0]#建立灰色预测模型XX0 = np.zeros(n)XX0[0] = X0[0]for i in range(1,n): XX0[i] = (X0[0] - u/a)*(1-math.exp(a))*math.exp(-a*(i));#模型精度的后验差检验e = 0 #求残差平均值for i in range(0,n): e += (X0[i] - XX0[i])e /= n#求历史数据平均值aver = 0; for i in range(0,n): aver += X0[i]aver /= n#求历史数据方差s12 = 0; for i in range(0,n): s12 += (X0[i]-aver)**2;s12 /= n#求残差方差s22 = 0; for i in range(0,n): s22 += ((X0[i] - XX0[i]) - e)**2;s22 /= n#求后验差比值C = s22 / s12 #求小误差概率cout = 0for i in range(0,n): if abs((X0[i] - XX0[i]) - e) < 0.6754*math.sqrt(s12): cout = cout+1 else: cout = coutP = cout / nif (C < 0.35 and P > 0.95): #预测精度为一级 m = 10 #请输入需要预测的年数 #print(\'往后m各年负荷为：\') f = np.zeros(m) for i in range(0,m): f[i] = (X0[0] - u/a)*(1-math.exp(a))*math.exp(-a*(i+n)) else: print(\'灰色预测法不适用\')

2. 机器学习/深度学习

监督学习

BP神经网络模型

一种多层前馈神经网络，通过反向传播算法优化权重，适用于分类和回归任务。

# -*- coding: utf-8 -*-\"\"\"Created on Mon Oct 1 22:15:54 2018@author: Heisenberg\"\"\"import numpy as npimport mathimport randomimport stringimport matplotlib as mplimport matplotlib.pyplot as plt#random.seed(0) #当我们设置相同的seed，每次生成的随机数相同。如果不设置seed，则每次会生成不同的随机数  #参考https://blog.csdn.net/jiangjiang_jian/article/details/79031788#生成区间[a,b]内的随机数def random_number(a,b): return (b-a)*random.random()+a#生成一个矩阵，大小为m*n,并且设置默认零矩阵def makematrix(m, n, fill=0.0): a = [] for i in range(m): a.append([fill]*n) return a#函数sigmoid(),这里采用tanh，因为看起来要比标准的sigmoid函数好看def sigmoid(x): return math.tanh(x)#函数sigmoid的派生函数def derived_sigmoid(x): return 1.0 - x**2#构造三层BP网络架构class BPNN: def __init__(self, num_in, num_hidden, num_out): #输入层，隐藏层，输出层的节点数 self.num_in = num_in + 1 #增加一个偏置结点 self.num_hidden = num_hidden + 1 #增加一个偏置结点 self.num_out = num_out  #激活神经网络的所有节点（向量） self.active_in = [1.0]*self.num_in self.active_hidden = [1.0]*self.num_hidden self.active_out = [1.0]*self.num_out  #创建权重矩阵 self.wight_in = makematrix(self.num_in, self.num_hidden) self.wight_out = makematrix(self.num_hidden, self.num_out)  #对权值矩阵赋初值 for i in range(self.num_in):  for j in range(self.num_hidden):  self.wight_in[i][j] = random_number(-0.2, 0.2) for i in range(self.num_hidden):  for j in range(self.num_out):  self.wight_out[i][j] = random_number(-0.2, 0.2) #最后建立动量因子（矩阵） self.ci = makematrix(self.num_in, self.num_hidden) self.co = makematrix(self.num_hidden, self.num_out)  #信号正向传播 def update(self, inputs): if len(inputs) != self.num_in-1:  raise ValueError(\'与输入层节点数不符\')  #数据输入输入层 for i in range(self.num_in - 1):  #self.active_in[i] = sigmoid(inputs[i]) #或者先在输入层进行数据处理  self.active_in[i] = inputs[i] #active_in[]是输入数据的矩阵  #数据在隐藏层的处理 for i in range(self.num_hidden - 1):  sum = 0.0  for j in range(self.num_in):  sum = sum + self.active_in[i] * self.wight_in[j][i]  self.active_hidden[i] = sigmoid(sum) #active_hidden[]是处理完输入数据之后存储，作为输出层的输入数据  #数据在输出层的处理 for i in range(self.num_out):  sum = 0.0  for j in range(self.num_hidden):  sum = sum + self.active_hidden[j]*self.wight_out[j][i]  self.active_out[i] = sigmoid(sum) #与上同理  return self.active_out[:] #误差反向传播 def errorbackpropagate(self, targets, lr, m): #lr是学习率， m是动量因子 if len(targets) != self.num_out:  raise ValueError(\'与输出层节点数不符！\')  #首先计算输出层的误差 out_deltas = [0.0]*self.num_out for i in range(self.num_out):  error = targets[i] - self.active_out[i]  out_deltas[i] = derived_sigmoid(self.active_out[i])*error  #然后计算隐藏层误差 hidden_deltas = [0.0]*self.num_hidden for i in range(self.num_hidden):  error = 0.0  for j in range(self.num_out):  error = error + out_deltas[j]* self.wight_out[i][j]  hidden_deltas[i] = derived_sigmoid(self.active_hidden[i])*error  #首先更新输出层权值 for i in range(self.num_hidden):  for j in range(self.num_out):  change = out_deltas[j]*self.active_hidden[i]  self.wight_out[i][j] = self.wight_out[i][j] + lr*change + m*self.co[i][j]  self.co[i][j] = change#然后更新输入层权值 for i in range(self.num_in):  for i in range(self.num_hidden):  change = hidden_deltas[j]*self.active_in[i]  self.wight_in[i][j] = self.wight_in[i][j] + lr*change + m* self.ci[i][j]  self.ci[i][j] = change#计算总误差 error = 0.0 for i in range(len(targets)):  error = error + 0.5*(targets[i] - self.active_out[i])**2 return error #测试 def test(self, patterns): for i in patterns:  print(i[0], \'->\', self.update(i[0])) #权重 def weights(self): print(\"输入层权重\") for i in range(self.num_in):  print(self.wight_in[i]) print(\"输出层权重\") for i in range(self.num_hidden):  print(self.wight_out[i])  def train(self, pattern, itera=100000, lr = 0.1, m=0.1): for i in range(itera):  error = 0.0  for j in pattern:  inputs = j[0]  targets = j[1]  self.update(inputs)  error = error + self.errorbackpropagate(targets, lr, m)  if i % 100 == 0:  print(\'误差 %-.5f\' % error) #实例def demo(): patt = [  [[1,2,5],[0]],  [[1,3,4],[1]],  [[1,6,2],[1]],  [[1,5,1],[0]],  [[1,8,4],[1]]  ] #创建神经网络，3个输入节点，3个隐藏层节点，1个输出层节点 n = BPNN(3, 3, 1) #训练神经网络 n.train(patt) #测试神经网络 n.test(patt) #查阅权重值 n.weights() if __name__ == \'__main__\': demo()

卷积神经网络（CNN）模型

主要用于图像识别和计算机视觉，通过卷积层提取空间特征，具有局部连接和权值共享的特点。

#卷积神经网络手写字体识别代码from tensorflow.examples.tutorials.mnist import input_dataimport tensorflow as tf #初始化权重函数def weight_variable(shape): initial = tf.truncated_normal(shape,stddev=0.1);#生成维度是shape标准差是0.1的正态分布数 return tf.Variable(initial) #初始化偏置项def bias_variable(shape): initial = tf.constant(0.1,shape=shape)#生成维度为shape的全为0.1的向量 return tf.Variable(initial) #定义卷积函数def conv2d(x,w): return tf.nn.conv2d(x,w,strides=[1,1,1,1],padding=\'SAME\') #strides： 卷积时在图像每一维的步长，这是一个一维的向量， #[ 1, strides, strides, 1]，第一位和最后一位固定必须是1 #padding参数是string类型，值为“SAME” 和 “VALID”，表示的是卷积的形式。 #设置为\"SAME\"时，会在原图像边缘外增加几圈0来使卷积后的矩阵和原图像矩阵的维度相同 #设置为\"VALID\"则不考虑这一点，卷积后的矩阵维度会相应减少，例如原图像如果是5*5，卷积核是3*3 #那么卷积过后的输出矩阵回是3*3的 #定义一个2*2的最大池化层def max_pool_2_2(x): return tf.nn.max_pool(x,ksize=[1,2,2,1],strides=[1,2,2,1],padding=\'SAME\') #第一个参数value：需要池化的输入，一般池化层接在卷积层后面， #依然是[batch, height, width, channels]这样的shape #第二个参数ksize：池化窗口的大小，取一个四维向量，一般是[1, height, width, 1]， #因为我们不想在batch和channels上做池化，所以这两个维度设为了1 #第三个参数strides：和卷积类似，窗口在每一个维度上滑动的步长， #一般也是[1, stride,stride, 1] #第四个参数padding：和卷积类似，可以取\'VALID\' 或者\'SAME\' if __name__ == \"__main__\": #定义输入变量 x = tf.placeholder(\"float\",shape=[None,784])#占位 #浮点型变量，行数不定，列数为784(每个图像是一个长度为784（28*28）的向量) #定义输出变量 y_ = tf.placeholder(\"float\",shape=[None,10])#占位 #浮点型变量，行数不定，列数为10,输出一个长度为10的向量来表示每个数字的可能性 #初始化权重,第一层卷积，32的意思代表的是输出32个通道 # 其实，也就是设置32个卷积，每一个卷积都会对图像进行卷积操作 w_conv1 = weight_variable([5,5,1,32])###生成了32个5*5的矩阵 #初始化偏置项 b_conv1 = bias_variable([32]) x_image = tf.reshape(x,[-1,28,28,1]) #将输入的x转成一个4D向量，第2、3维对应图片的宽高，最后一维代表图片的颜色通道数 # 输入的图像为灰度图，所以通道数为1，如果是RGB图，通道数为3 # tf.reshape(x,[-1,28,28,1])的意思是将x自动转换成28*28*1的数组 # -1的意思是代表不知道x的shape，它会按照后面的设置进行转换 # 卷积并激活 h_conv1 = tf.nn.relu(conv2d(x_image,w_conv1) + b_conv1) #池化 h_pool1 = max_pool_2_2(h_conv1) #第二层卷积 #初始权重 w_conv2 = weight_variable([5,5,32,64]) #在32个第一层卷积层上每个再用一个5*5的卷积核在做特征提取，并输出到第二层卷积层， #第二层设置了64个卷积层 #初始化偏置项 b_conv2 = bias_variable([64]) #将第一层卷积池化后的结果作为第二层卷积的输入加权求和后激活 h_conv2 = tf.nn.relu(conv2d(h_pool1,w_conv2) + b_conv2) #池化 h_pool2 = max_pool_2_2(h_conv2) # 设置全连接层的权重 w_fc1 = weight_variable([7*7*64,1024]) #28*28的原图像经过两次池化后变为7*7，设置了1024个输出单元 # 设置全连接层的偏置 b_fc1 = bias_variable([1024]) # 将第二层卷积池化后的结果，转成一个7*7*64的数组 h_pool2_flat = tf.reshape(h_pool2,[-1,7*7*64]) # 通过全连接之后并激活 h_fc1 = tf.nn.relu(tf.matmul(h_pool2_flat,w_fc1) + b_fc1) # 防止过拟合 keep_prob = tf.placeholder(\"float\")#占位 h_fc1_drop = tf.nn.dropout(h_fc1,keep_prob) #设置每个单元保留的概率来随机放弃一些单元来防止过拟合 #输出层 w_fc2 = weight_variable([1024,10]) b_fc2 = bias_variable([10]) #加权求和并激活 y_conv = tf.nn.softmax(tf.matmul(h_fc1_drop,w_fc2) + b_fc2) #日志输出，每迭代100次输出一次日志 #定义交叉熵为损失函数 cross_entropy = -tf.reduce_sum(y_ * tf.log(y_conv)) #最小化交叉熵 train_step = tf.train.AdamOptimizer(1e-4).minimize(cross_entropy) #计算准确率 correct_prediction = tf.equal(tf.argmax(y_conv,1),tf.argmax(y_,1)) accuracy = tf.reduce_mean(tf.cast(correct_prediction,\"float\")) sess = tf.Session() sess.run(tf.initialize_all_variables()) #上面的两行是在为tf的输出变量做准备 # 下载minist的手写数字的数据集 mnist = input_data.read_data_sets(\"MNIST_data/\", one_hot=True) for i in range(20000):#迭代20000次 batch = mnist.train.next_batch(50)#设置batch，即每次用来训练模型的数据个数 if i % 100 == 0:#每100次迭代输出一次精度 train_accuracy = accuracy.eval(session=sess,  feed_dict={x:batch[0],y_:batch[1],keep_prob:1.0}) #喂给之前占位的x和y_本次训练的batch个数据中第一个数据的图像矩阵和标签，不考虑过拟合 #计算当前的精度 print(\"step %d,training accuracy %g\"%(i,train_accuracy)) train_step.run(session = sess,feed_dict={x:batch[0],y_:batch[1],keep_prob:0.5}) #当i不能被100整除时的训练过程，考虑过拟合，单元保留的概率为0.5 print(\"test accuracy %g\" % accuracy.eval(session=sess,feed_dict={ x: mnist.test.images, y_: mnist.test.labels, keep_prob: 1.0})) #输出测试集的精度

决策树分类模型

基于树状结构进行决策，通过特征分裂实现分类，直观易解释，如ID3、C4.5和CART算法。

逻辑回归模型

用于二分类问题，通过Sigmoid函数将线性回归结果映射到概率，适用于线性可分数据。

import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport ospath=\'data\'+os.sep+\'Logireg_data.txt\'pdData=pd.read_csv(path,header=None,names=[\'Exam1\',\'Exam2\',\'Admitted\'])pdData.head()print(pdData.head())print(pdData.shape)positive=pdData[pdData[\'Admitted\']==1]#定义正nagative=pdData[pdData[\'Admitted\']==0]#定义负fig,ax=plt.subplots(figsize=(10,5))ax.scatter(positive[\'Exam1\'],positive[\'Exam2\'],s=30,c=\'b\',marker=\'o\',label=\'Admitted\')ax.scatter(nagative[\'Exam1\'],nagative[\'Exam2\'],s=30,c=\'r\',marker=\'x\',label=\'not Admitted\')ax.legend()ax.set_xlabel(\'Exam 1 score\')ax.set_ylabel(\'Exam 2 score\')plt.show()#画图##实现算法 the logistics regression 目标建立一个分类器 设置阈值来判断录取结果##sigmoid 函数def sigmoid(z): return 1/(1+np.exp(-z))#画图nums=np.arange(-10,10,step=1)fig,ax=plt.subplots(figsize=(12,4))ax.plot(nums,sigmoid(nums),\'r\')#画图定义plt.show()#按照理论实现预测函数def model(X,theta): return sigmoid(np.dot(X,theta.T)) pdData.insert(0,\'ones\',1)#插入一列orig_data=pdData.as_matrix()cols=orig_data.shape[1]X=orig_data[:,0:cols-1]y=orig_data[:,cols-1:cols]theta=np.zeros([1,3])print(X[:5])print(X.shape,y.shape,theta.shape)##损失函数def cost(X,y,theta): left=np.multiply(-y,np.log(model(X,theta))) right=np.multiply(1-y,np.log(1-model(X,theta))) return np.sum(left-right)/(len(X))print(cost(X,y,theta)) #计算梯度def gradient(X, y, theta): grad = np.zeros(theta.shape) error = (model(X, theta) - y).ravel() for j in range(len(theta.ravel())): # for each parmeter term = np.multiply(error, X[:, j]) grad[0, j] = np.sum(term) / len(X) return grad##比较3种不同梯度下降方法STOP_ITER=0STOP_COST=1STOP_GRAD=2 def stopCriterion(type,value,threshold): if type==STOP_ITER: return value>threshold elif type==STOP_COST: return abs(value[-1]-value[-2])<threshold elif type==STOP_GRAD: return np.linalg.norm(value)<threshold import numpy.random#打乱数据洗牌def shuffledata(data): np.random.shuffle(data) cols=data.shape[1] X=data[:,0:cols-1] y=data[:,cols-1:] return X,y import time def descent(data, theta, batchSize, stopType, thresh, alpha): # 梯度下降求解 init_time = time.time() i = 0 # 迭代次数 k = 0 # batch X, y = shuffledata(data) grad = np.zeros(theta.shape) # 计算的梯度 costs = [cost(X, y, theta)] # 损失值 while True: grad = gradient(X[k:k + batchSize], y[k:k + batchSize], theta) k += batchSize # 取batch数量个数据 if k >= n: k = 0 X, y = shuffledata(data) # 重新洗牌 theta = theta - alpha * grad # 参数更新 costs.append(cost(X, y, theta)) # 计算新的损失 i += 1 if stopType == STOP_ITER: value = i elif stopType == STOP_COST: value = costs elif stopType == STOP_GRAD: value = grad if stopCriterion(stopType, value, thresh): break return theta, i - 1, costs, grad, time.time() - init_time#选择梯度下降def runExpe(data, theta, batchSize, stopType, thresh, alpha): #import pdb; pdb.set_trace(); theta, iter, costs, grad, dur = descent(data, theta, batchSize, stopType, thresh, alpha) name = \"Original\" if (data[:,1]>2).sum() > 1 else \"Scaled\" name += \" data - learning rate: {} - \".format(alpha) if batchSize==n: strDescType = \"Gradient\" elif batchSize==1: strDescType = \"Stochastic\" else: strDescType = \"Mini-batch ({})\".format(batchSize) name += strDescType + \" descent - Stop: \" if stopType == STOP_ITER: strStop = \"{} iterations\".format(thresh) elif stopType == STOP_COST: strStop = \"costs change < {}\".format(thresh) else: strStop = \"gradient norm < {}\".format(thresh) name += strStop print (\"***{}\\nTheta: {} - Iter: {} - Last cost: {:03.2f} - Duration: {:03.2f}s\".format( name, theta, iter, costs[-1], dur)) fig, ax = plt.subplots(figsize=(12,4)) ax.plot(np.arange(len(costs)), costs, \'r\') ax.set_xlabel(\'Iterations\') ax.set_ylabel(\'Cost\') ax.set_title(name.upper() + \' - Error vs. Iteration\') return thetan= 100runExpe(orig_data,theta,n,STOP_ITER,thresh=5000,alpha=0.000001)plt.show()runExpe(orig_data,theta,n,STOP_GRAD,thresh=0.05,alpha=0.001)plt.show()runExpe(orig_data,theta,n,STOP_COST,thresh=0.000001,alpha=0.001)plt.show()#对比runExpe(orig_data, theta, 1, STOP_ITER, thresh=5000, alpha=0.001)plt.show()runExpe(orig_data, theta, 1, STOP_ITER, thresh=15000, alpha=0.000002)plt.show()runExpe(orig_data, theta, 16, STOP_ITER, thresh=15000, alpha=0.001)plt.show()##对数据进行标准化 将数据按其属性(按列进行)减去其均值，然后除以其方差。#最后得到的结果是，对每个属性/每列来说所有数据都聚集在0附近，方差值为1 from sklearn import preprocessing as pp scaled_data = orig_data.copy()scaled_data[:, 1:3] = pp.scale(orig_data[:, 1:3]) runExpe(scaled_data, theta, n, STOP_ITER, thresh=5000, alpha=0.001)#设定阈值def predict(X, theta): return [1 if x >= 0.5 else 0 for x in model(X, theta)] # if __name__==\'__main__\': scaled_X = scaled_data[:, :3]y = scaled_data[:, 3]predictions = predict(scaled_X, theta)correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]accuracy = (sum(map(int, correct)) % len(correct))print (\'accuracy = {0}%\'.format(accuracy))

随机森林分类模型

集成学习方法，通过多棵决策树投票提高泛化能力，抗过拟合且适用于高维数据。

随机森林需要调整的参数有：（1） 决策树的个数（2） 特征属性的个数（3） 递归次数（即决策树的深度）\'\'\'from numpy import inffrom numpy import zerosimport numpy as npfrom sklearn.model_selection import train_test_split #生成数据集。数据集包括标签，全包含在返回值的dataset上def get_Datasets(): from sklearn.datasets import make_classification dataSet,classLabels=make_classification(n_samples=200,n_features=100,n_classes=2) #print(dataSet.shape,classLabels.shape) return np.concatenate((dataSet,classLabels.reshape((-1,1))),axis=1) #切分数据集，实现交叉验证。可以利用它来选择决策树个数。但本例没有实现其代码。#原理如下：#第一步，将训练集划分为大小相同的K份；#第二步，我们选择其中的K-1分训练模型，将用余下的那一份计算模型的预测值，#这一份通常被称为交叉验证集；第三步，我们对所有考虑使用的参数建立模型#并做出预测，然后使用不同的K值重复这一过程。#然后是关键，我们利用在不同的K下平均准确率最高所对应的决策树个数#作为算法决策树个数 def splitDataSet(dataSet,n_folds): #将训练集划分为大小相同的n_folds份； fold_size=len(dataSet)/n_folds data_split=[] begin=0 end=fold_size for i in range(n_folds): data_split.append(dataSet[begin:end,:]) begin=end end+=fold_size return data_split#构建n个子集def get_subsamples(dataSet,n): subDataSet=[] for i in range(n): index=[] #每次都重新选择k个 索引 for k in range(len(dataSet)): #长度是k index.append(np.random.randint(len(dataSet))) #(0,len(dataSet)) 内的一个整数 subDataSet.append(dataSet[index,:]) return subDataSet # subDataSet=get_subsamples(dataSet,10)############################################################################# #根据某个特征及值对数据进行分类def binSplitDataSet(dataSet,feature,value): mat0=dataSet[np.nonzero(dataSet[:,feature]>value)[0],:] mat1=dataSet[np.nonzero(dataSet[:,feature]<value)[0],:] return mat0,mat1 \'\'\' feature=2 value=1 dataSet=get_Datasets() mat0,mat1= binSplitDataSet(dataSet,2,1)\'\'\' #计算方差，回归时使用def regErr(dataSet): return np.var(dataSet[:,-1])*np.shape(dataSet)[0] #计算平均值，回归时使用def regLeaf(dataSet): return np.mean(dataSet[:,-1]) def MostNumber(dataSet): #返回多类 #number=set(dataSet[:,-1]) len0=len(np.nonzero(dataSet[:,-1]==0)[0]) len1=len(np.nonzero(dataSet[:,-1]==1)[0]) if len0>len1: return 0 else: return 1 #计算基尼指数 一个随机选中的样本在子集中被分错的可能性 是被选中的概率乘以被分错的概率 def gini(dataSet): corr=0.0 for i in set(dataSet[:,-1]):  #i 是这个特征下的 某个特征值 corr+=(len(np.nonzero(dataSet[:,-1]==i)[0])/len(dataSet))**2 return 1-corr def select_best_feature(dataSet,m,alpha=\"huigui\"): f=dataSet.shape[1]#拿过这个数据集，看这个数据集有多少个特征，即f个 index=[] bestS=inf; bestfeature=0;bestValue=0; if alpha==\"huigui\": S=regErr(dataSet) else: S=gini(dataSet) for i in range(m): index.append(np.random.randint(f)) #在f个特征里随机，注意是随机！选择m个特征，然后在这m个特征里选择一个合适的分类特征。 for feature in index: for splitVal in set(dataSet[:,feature]):  #set() 函数创建一个无序不重复元素集，用于遍历这个特征下所有的值 mat0,mat1=binSplitDataSet(dataSet,feature,splitVal)  if alpha==\"huigui\": newS=regErr(mat0)+regErr(mat1) #计算每个分支的回归方差 else: newS=gini(mat0)+gini(mat1) #计算被分错率 if bestS>newS: bestfeature=feature bestValue=splitVal bestS=newS if (S-bestS)<0.001 and alpha==\"huigui\":# 对于回归来说，方差足够了，那就取这个分支的均值 return None,regLeaf(dataSet) elif (S-bestS)<0.001: #print(S,bestS) return None,MostNumber(dataSet) #对于分类来说，被分错率足够下了，那这个分支的分类就是大多数所在的类。 #mat0,mat1=binSplitDataSet(dataSet,feature,splitVal) return bestfeature,bestValue def createTree(dataSet,alpha=\"huigui\",m=20,max_level=10): #实现决策树，使用20个特征，深度为10， bestfeature,bestValue=select_best_feature(dataSet,m,alpha=alpha) if bestfeature==None: return bestValue retTree={} max_level-=1 if max_level<0: #控制深度 return regLeaf(dataSet) retTree[\'bestFeature\']=bestfeature retTree[\'bestVal\']=bestValue lSet,rSet=binSplitDataSet(dataSet,bestfeature,bestValue) #lSet是根据特征bestfeature分到左边的向量，rSet是根据特征bestfeature分到右边的向量 retTree[\'right\']=createTree(rSet,alpha,m,max_level) retTree[\'left\']=createTree(lSet,alpha,m,max_level) #每棵树都是二叉树，往下分类都是一分为二。 #print(\'retTree:\',retTree) return retTree def RondomForest(dataSet,n,alpha=\"huigui\"): #树的个数 #dataSet=get_Datasets() Trees=[] # 设置一个空树集合 for i in range(n): X_train, X_test, y_train, y_test = train_test_split(dataSet[:,:-1], dataSet[:,-1], test_size=0.33, random_state=42) X_train=np.concatenate((X_train,y_train.reshape((-1,1))),axis=1) Trees.append(createTree(X_train,alpha=alpha)) return Trees # 生成好多树################################################################### #预测单个数据样本，重头！！如何利用已经训练好的随机森林对单个样本进行 回归或分类！def treeForecast(trees,data,alpha=\"huigui\"): if alpha==\"huigui\": if not isinstance(trees,dict): #isinstance() 函数来判断一个对象是否是一个已知的类型 return float(trees) if data[trees[\'bestFeature\']]>trees[\'bestVal\']: # 如果数据的这个特征大于阈值，那就调用左支 if type(trees[\'left\'])==\'float\':  #如果左支已经是节点了，就返回数值。如果左支还是字典结构，那就继续调用， 用此支的特征和特征值进行选支。  return trees[\'left\'] else: return treeForecast(trees[\'left\'],data,alpha) else: if type(trees[\'right\'])==\'float\': return trees[\'right\'] else: return treeForecast(trees[\'right\'],data,alpha) else: if not isinstance(trees,dict):#分类和回归是同一道理 return int(trees) if data[trees[\'bestFeature\']]>trees[\'bestVal\']: if type(trees[\'left\'])==\'int\': return trees[\'left\'] else: return treeForecast(trees[\'left\'],data,alpha) else: if type(trees[\'right\'])==\'int\': return trees[\'right\'] else: return treeForecast(trees[\'right\'],data,alpha) #随机森林 对 数据集打上标签 0、1 或者是 回归值def createForeCast(trees,test_dataSet,alpha=\"huigui\"): cm=len(test_dataSet) yhat=np.mat(zeros((cm,1))) for i in range(cm):  # yhat[i,0]=treeForecast(trees,test_dataSet[i,:],alpha) # return yhat #随机森林预测def predictTree(Trees,test_dataSet,alpha=\"huigui\"): #trees 是已经训练好的随机森林 调用它！ cm=len(test_dataSet) yhat=np.mat(zeros((cm,1))) for trees in Trees: yhat+=createForeCast(trees,test_dataSet,alpha) #把每次的预测结果相加 if alpha==\"huigui\": yhat/=len(Trees) #如果是回归的话，每棵树的结果应该是回归值，相加后取平均 else: for i in range(len(yhat)):  #如果是分类的话，每棵树的结果是一个投票向量，相加后， #看每类的投票是否超过半数，超过半数就确定为1 if yhat[i,0]>len(Trees)/2: yhat[i,0]=1 else: yhat[i,0]=0 return yhat if __name__ == \'__main__\' : dataSet=get_Datasets() print(dataSet[:,-1].T)  #打印标签，与后面预测值对比 .T其实就是对一个矩阵的转置 RomdomTrees=RondomForest(dataSet,4,alpha=\"fenlei\") #这里我训练好了 很多树的集合，就组成了随机森林。一会一棵一棵的调用。 print(\"---------------------RomdomTrees------------------------\") #print(RomdomTrees[0]) test_dataSet=dataSet #得到数据集和标签 yhat=predictTree(RomdomTrees,test_dataSet,alpha=\"fenlei\") # 调用训练好的那些树。综合结果，得到预测值。 print(yhat.T)#get_Datasets() print(dataSet[:,-1].T-yhat.T)

支持向量机（SVM）模型

通过寻找最优超平面进行分类，适用于小样本、非线性及高维数据（借助核函数）。

from numpy import *import randomimport matplotlib.pyplot as pltimport numpy def kernelTrans(X,A,kTup):  # 核函数(此例未使用) m,n=shape(X) K = mat(zeros((m,1))) if kTup[0] ==\'lin\': K=X*A.T elif kTup[0]==\'rbf\': for j in range(m): deltaRow = X[j,:]-A K[j]=deltaRow*deltaRow.T  # ||w||^2 = w^T * w K =exp(K/(-1*kTup[1]**2))  # K = e^(||x-y||^2 / (-2*sigma^2)) else: raise NameError(\"Houston we Have a problem --\") return K class optStruct: def __init__(self,dataMain,classLabel,C,toler,kTup): self.X = dataMain  # 样本矩阵 self.labelMat = classLabel self.C = C # 惩罚因子 self.tol = toler# 容错率 self.m = shape(dataMain)[0]  # 样本点个数 self.alphas = mat(zeros((self.m,1))) # 产生m个拉格郎日乘子，组成一个m×1的矩阵 self.b =0 # 决策面的截距 self.eCache = mat(zeros((self.m,2))) # 产生m个误差 E=f(x)-y ,设置成m×2的矩阵，矩阵第一列是标志位，标志为1就是E计算好了，第二列是误差E # self.K = mat(zeros((self.m,self.m))) # for i in range(self.m):  # K[,]保存的是任意样本之间的相似度(用高斯核函数表示的相似度) # self.K[:,i]=kernelTrans(self.X,self.X[i,:],kTup) def loadDataSet(filename):  # 加载数据 dataMat = [] labelMat = [] fr = open(filename) for line in fr.readlines(): lineArr = line.split() dataMat.append([float(lineArr[0]),float(lineArr[1])]) labelMat.append(float(lineArr[2])) # 一维列表 return dataMat, labelMat def selectJrand(i, m): # 随机选择一个不等于i的下标 j =i while(j==i): j = int(random.uniform(0,m)) return j def clipAlpha(aj, H,L): if aj>H:# 如果a^new 大于上限值，那么就把上限赋给它 aj = H if L>aj:# 如果a^new 小于下限值，那么就把下限赋给它 aj = L return aj def calcEk(oS, k):  # 计算误差E, k代表第k个样本点，它是下标，oS是optStruct类的实例 # fXk = float(multiply(oS.alphas,oS.labelMat).T * oS.K[:,k] + oS.b) # 公式f(x)=sum(ai*yi*xi^T*x)+b fXk = float(multiply(oS.alphas,oS.labelMat).T * (oS.X*oS.X[k,:].T)) +oS.b Ek = fXk - float(oS.labelMat[k]) # 计算误差 E=f(x)-y return Ek def selectJ(i, oS, Ei): # 选择两个拉格郎日乘子,在所有样本点的误差计算完毕之后，寻找误差变化最大的那个样本点及其误差 maxK = -1 # 最大步长的因子的下标 maxDeltaE = 0 # 最大步长 Ej = 0  # 最大步长的因子的误差 oS.eCache[i] = [1,Ei] valiEcacheList = nonzero(oS.eCache[:,0].A)[0] # nonzero结果是两个array数组，第一个数组是不为0的元素的x坐标，第二个数组是该位置的y坐标 # 此处寻找误差矩阵第一列不为0的数的下标 print(\"valiEcacheList is {}\".format(valiEcacheList)) if (len(valiEcacheList))>1: for k in valiEcacheList: # 遍历所有计算好的Ei的下标，valiEcacheLIst保存了所有样本点的E，计算好的有效位置是1,没计算好的是0 if k == i: continue Ek = calcEk(oS,k) deltaE = abs(Ei-Ek) # 距离第一个拉格朗日乘子a1绝对值最远的作为第二个朗格朗日乘子a2 if deltaE>maxDeltaE: maxK = k  # 记录选中的这个乘子a2的下标 maxDeltaE = deltaE # 记录他俩的绝对值 Ej = Ek  # 记录a2此时的误差 return maxK, Ej else: # 如果是第一次循环，随机选择一个alphas j = selectJrand(i, oS.m) # j = 72 Ej = calcEk(oS, j) return j,Ej def updateEk(oS, k): Ek = calcEk(oS, k) oS.eCache[k] = [1,Ek] # 把第k个样本点的误差计算出来，并存入误差矩阵，有效位置设为1 def innerL(i, oS): Ei = calcEk(oS, i)  # KKT条件， 若yi*(w^T * x +b)-1<0 则 ai=C 若yi*(w^T * x +b)-1>0 则 ai=0 print(\"i is {0},Ei is {1}\".format(i,Ei)) if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)): j,Ej = selectJ(i,oS,Ei) print(\"第二个因子的坐标{}\".format(j)) alphaIold = oS.alphas[i].copy() # 用了浅拷贝， alphaIold 就是old a1,对应公式 alphaJold = oS.alphas[j].copy() if oS.labelMat[i] != oS.labelMat[j]: # 也是根据公式来的，y1 不等于 y2时 L = max(0,oS.alphas[j] - oS.alphas[i]) H = min(oS.C, oS.C+oS.alphas[j]-oS.alphas[i]) else: L = max(0,oS.alphas[j]+oS.alphas[i]-oS.C) H = min(oS.C,oS.alphas[j]+oS.alphas[i]) if L==H: # 如果这个j让L=H,i和j这两个样本是同一类别，且ai=aj=0或ai=aj=C，或者不同类别，aj=C且ai=0 # 当同类别时 ai+aj = 常数 ai是不满足KKT的，假设ai=0,需增大它，那么就得减少aj，aj已经是0了，不能最小了，所以此情况不允许发生 # 当不同类别时 ai-aj=常数，ai是不满足KKT的，ai=0,aj=C,ai需增大，它则aj也会变大，但是aj已经是C的不能再大了，故此情况不允许 print(\"L=H\") return 0 # eta = 2.0*oS.K[i,j]-oS.K[i,i]-oS.K[j,j] # eta=K11+K22-2*K12 eta = 2.0*oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T if eta >= 0:  # 这里跟公式正好差了一个负号，所以对应公式里的 K11+K22-2*K12 <=0，即开口向下，或为0成一条直线的情况不考虑 print(\"eta>=0\") return 0 oS.alphas[j]-=oS.labelMat[j]*(Ei-Ej)/eta # a2^new = a2^old+y2(E1-E2)/eta print(\"a2 归约之前是{}\".format(oS.alphas[j])) oS.alphas[j]=clipAlpha(oS.alphas[j],H,L) # 根据公式，看看得到的a2^new是否在上下限之内 print(\"a2 归约之后is {}\".format(oS.alphas[j])) # updateEk(oS,j)  # 把更新后的a2^new的E更新一下 if abs(oS.alphas[j]-alphaJold)<0.00001: print(\"j not moving enough\") return 0 oS.alphas[i] +=oS.labelMat[j]*oS.labelMat[i]*(alphaJold-oS.alphas[j]) # 根据公式a1^new = a1^old+y1*y2*(a2^old-a2^new) print(\"a1更新之后是{}\".format(oS.alphas[i])) # updateEk(oS,i) # b1^new = b1^old+(a1^old-a1^new)y1*K11+(a2^old-a2^new)y2*K12-E1 # b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i]-oS.labelMat[j]*\\ # (oS.alphas[j]-alphaJold)*oS.K[i,j] b1 = oS.b-Ei-oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T-oS.labelMat[j]* \\ (oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T # b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]-oS.labelMat[j]* \\ # (oS.alphas[j]-alphaJold)*oS.K[j,j] b2 = oS.b-Ej-oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T-oS.labelMat[j]* \\ (oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T updateEk(oS,j) # 个人认为更新误差应在更新b之后，因为公式算出的b的公式使用的是以前的Ei updateEk(oS,i) # b2^new=b2^old+(a1^old-a1^new)y1*K12+(a2^old-a2^new)y2*K22-E2 if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1.A[0][0] elif (0<oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2.A[0][0] else: oS.b = (b1+b2)/2.0 print(\"b is {}\".format(oS.b)) return 1 else: return 0 def smoP(dataMatIn, classLabels, C,toler,maxIter,kTup=(\'lin\',)): oS = optStruct(mat(dataMatIn), mat(classLabels).transpose(),C,toler,kTup) iter = 0 entireSet = True  # 两种遍历方式交替 alphaPairsChanged = 0 while (iter<maxIter) and ((alphaPairsChanged>0) or (entireSet)): alphaPairsChanged = 0 if entireSet: for i in range(oS.m): alphaPairsChanged += innerL(i,oS) print(\"fullSet, iter:%d i: %d pairs changed %d\"%(iter,i ,alphaPairsChanged)) iter+=1 print(\"第一种遍历alphaRairChanged is {}\".format(alphaPairsChanged)) print(\"-----------eCache is {}\".format(oS.eCache)) print(\"***********alphas is {}\".format(oS.alphas)) print(\"---------------------------------------\") else: nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0] # 这时数组相乘，里面其实是True 和False的数组，得出来的是 # 大于0并且小于C的alpha的下标 for i in nonBoundIs: alphaPairsChanged += innerL(i,oS) print(\"non-bound, iter: %d i:%d, pairs changed %d\"%(iter,i,alphaPairsChanged)) print(\"第二种遍历alphaPairChanged is {}\".format(alphaPairsChanged)) iter+=1 if entireSet: entireSet = False # 当第二种遍历方式alpha不再变化，那么继续第一种方式扫描，第一种方式不再变化，此时alphachanged为0且entireSet为false，退出循环 elif (alphaPairsChanged==0): entireSet=True print(\"iteration number: %d\"%iter) return oS.b,oS.alphas def calcWs(alphas,dataArr,classLabels): # 通过alpha来计算w X = mat(dataArr) labelMat = mat(classLabels).transpose() m,n = shape(X) w = zeros((n,1)) for i in range(m): w += multiply(alphas[i]*labelMat[i], X[i,:].T) # w = sum(ai*yi*xi) return w def draw_points(dataArr,classlabel, w,b,alphas): myfont = FontProperties(fname=\'/usr/share/fonts/simhei.ttf\') # 显示中文 plt.rcParams[\'axes.unicode_minus\'] = False # 防止坐标轴的‘-’变为方块 m = len(classlabel) red_points_x=[] red_points_y =[] blue_points_x=[] blue_points_y =[] svc_points_x =[] svc_points_y =[] # print(type(alphas)) svc_point_index = nonzero((alphas.A>0) * (alphas.A <0.8))[0] svc_points = array(dataArr)[svc_point_index] svc_points_x = [x[0] for x in list(svc_points)] svc_points_y = [x[1] for x in list(svc_points)] print(\"svc_points_x\",svc_points_x) print(\"svc_points_y\",svc_points_y) for i in range(m): if classlabel[i] ==1: red_points_x.append(dataArr[i][0]) red_points_y.append(dataArr[i][1]) else: blue_points_x.append(dataArr[i][0]) blue_points_y.append(dataArr[i][1]) fig = plt.figure()  # 创建画布 ax = fig.add_subplot(111) ax.set_title(\"SVM-Classify\")  # 设置图片标题 ax.set_xlabel(\"x\")  # 设置坐标名称 ax.set_ylabel(\"y\") ax1=ax.scatter(red_points_x, red_points_y, s=30,c=\'red\', marker=\'s\') #s是shape大小，c是颜色，marker是形状，\'s\'代表是正方形,默认\'o\'是圆圈 ax2=ax.scatter(blue_points_x, blue_points_y, s=40,c=\'green\') # ax.set_ylim([-6,5]) print(\"b\",b) print(\"w\",w) x = arange(-4.0, 4.0, 0.1)  # 分界线x范围,步长为0.1 # x = arange(-2.0,10.0) if isinstance(b,numpy.matrixlib.defmatrix.matrix): b = b.A[0][0] y = (-b-w[0][0]*x)/w[1][0] # 直线方程 Ax + By + C = 0 ax3,=plt.plot(x,y, \'k\') ax4=plt.scatter(svc_points_x,svc_points_y,s=50,c=\'orange\',marker=\'p\') plt.legend([ax1, ax2,ax3,ax4], [\"red points\",\"blue points\", \"decision boundary\",\"support vector\"], loc=\'lower right\') # 标注 plt.show() dataArr,labelArr = loadDataSet(\'/home/zhangqingfeng/test/svm_test_data\')b,alphas = smoP(dataArr,labelArr,0.8,0.001,40)w=calcWs(alphas,dataArr,labelArr)draw_points(dataArr,labelArr,w,b,alphas)

无监督学习

K-means聚类模型

K均值聚类算法是先随机选取K个对象作为初始的聚类中心。然后计算每个对象与各个种子聚类中心之间的距离，把每个对象分配给距离它最近的聚类中心。聚类中心以及分配给它们的对象就代表一个聚类。每分配一个样本，聚类的聚类中心会根据聚类中现有的对象被重新计算。这个过程将不断重复直到满足某个终止条件。终止条件可以是没有（或最小数目）对象被重新分配给不同的聚类，没有（或最小数目）聚类中心再发生变化，误差平方和局部最小。

  1. 便于理解，首先创建一个明显分为2类20*2的例子（每一列为一个变量共2个变量，每一行为一个样本共20个样本）：
import numpy as npc1x=np.random.uniform(0.5,1.5,(1,10))c1y=np.random.uniform(0.5,1.5,(1,10))c2x=np.random.uniform(3.5,4.5,(1,10))c2y=np.random.uniform(3.5,4.5,(1,10))x=np.hstack((c1x,c2x))y=np.hstack((c2y,c2y))X=np.vstack((x,y)).Tprint(X) 结果：[[1.4889993 4.18741329] [0.73017615 4.07842216] [1.15522846 4.05744838] [1.40768457 3.76674812] [1.376212 3.95063903] [1.20821055 4.34138767] [0.73898392 3.55026013] [0.97116627 3.65432314] [0.98267302 4.16731561] [1.06346541 4.44383585] [4.10945954 4.18741329] [3.75288064 4.07842216] [4.29638229 4.05744838] [3.95221785 3.76674812] [4.09826192 3.95063903] [4.04840874 4.34138767] [4.29594009 3.55026013] [3.56931245 3.65432314] [3.57962941 4.16731561] [3.65208848 4.44383585]]
  1. 引用Python库将样本分为两类（k=2），并绘制散点图：
#只需将X修改即可进行其他聚类分析import matplotlib.pyplot as pltfrom sklearn.cluster import KMeans kemans=KMeans(n_clusters=2)result=kemans.fit_predict(X) #训练及预测print(result) #分类结果 plt.rcParams[\'font.family\'] = [\'sans-serif\']plt.rcParams[\'font.sans-serif\'] = [\'SimHei\'] #散点图标签可以显示中文 x=[i[0] for i in X]y=[i[1] for i in X]plt.scatter(x,y,c=result,marker=\'o\')plt.xlabel(\'x\')plt.ylabel(\'y\')plt.show() 结果：[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1]
  1. 如果K值未知，可采用肘部法选择K值（假设最大分类数为9类，分别计算分类结果为1-9类的平均离差，离差的提升变化下降最抖时的值为最优聚类数K）：
import matplotlib.pyplot as pltfrom sklearn.cluster import KMeansfrom scipy.spatial.distance import cdistK=range(1,10)meanDispersions=[]for k in K: kemans=KMeans(n_clusters=k) kemans.fit(X) #计算平均离差 m_Disp=sum(np.min(cdist(X,kemans.cluster_centers_,\'euclidean\'),axis=1))/X.shape[0] meanDispersions.append(m_Disp)plt.rcParams[\'font.family\'] = [\'sans-serif\']plt.rcParams[\'font.sans-serif\'] = [\'SimHei\'] #使折线图显示中文plt.plot(K,meanDispersions,\'bx-\')plt.xlabel(\'k\')plt.ylabel(\'平均离差\')plt.title(\'用肘部方法选择K值\')plt.show()

主成分分析（PCA）算法

通过线性变换降维，保留最大方差的主成分，用于数据压缩和特征提取。

import matplotlib.pyplot as plt  #加载matplotlib用于数据的可视化from sklearn.decomposition import PCA  #加载PCA算法包from sklearn.datasets import load_iris data=load_iris()y=data.targetx=data.datapca=PCA(n_components=2) #加载PCA算法，设置降维后主成分数目为2reduced_x=pca.fit_transform(x)#对样本进行降维red_x,red_y=[],[]blue_x,blue_y=[],[]green_x,green_y=[],[] for i in range(len(reduced_x)): if y[i] ==0: red_x.append(reduced_x[i][0]) red_y.append(reduced_x[i][1]) elif y[i]==1: blue_x.append(reduced_x[i][0]) blue_y.append(reduced_x[i][1]) else: green_x.append(reduced_x[i][0]) green_y.append(reduced_x[i][1])#可视化plt.scatter(red_x,red_y,c=\'r\',marker=\'x\')plt.scatter(blue_x,blue_y,c=\'b\',marker=\'D\')plt.scatter(green_x,green_y,c=\'g\',marker=\'.\')plt.show()
网络标签：

相关问题