数模代码整理

绘图类

python绘图

调用 python matplotlib.pyplot库进行绘制，除如下列举的例子外，其余可参考：Matplotlib绘图详解

折线图

import matplotlib.pyplot as plt
import numpy as np

xpoints = np.array([1, 8])
ypoints = np.array([3, 10])

plt.plot(xpoints, ypoints)
plt.show()

散点图

import matplotlib.pyplot as plt
import numpy as np

xpoints = np.array([1, 2, 6, 8])
ypoints = np.array([3, 8, 1, 10])

plt.plot(xpoints, ypoints, 'o')
plt.show()

饼图

import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif']=['Microsoft YaHei']  #显示中文标签,处理中文乱码问题
plt.rcParams['axes.unicode_minus']=False  #坐标轴负号的处理
plt.axes(aspect='equal')  #将横、纵坐标轴标准化处理，确保饼图是一个正圆，否则为椭圆

#构造数据
edu = [0.2515, 0.3724, 0.3336, 0.0368, 0.0057]
labels = ['中专','大专','本科','硕士','其他']
explode = [0, 0.1, 0, 0, 0]  #生成数据，用于凸显大专学历人群
colors = ['#9999ff', '#ff9999', '#7777aa', '#2442aa', '#dd5555']  #自定义颜色

plt.pie(x=edu,  #绘图数据
        explode=explode, #指定饼图某些部分的突出显示，即呈现爆炸式
        labels=labels,  #添加教育水平标签
        colors=colors,
        autopct='%.2f%%',  #设置百分比的格式，这里保留两位小数
        pctdistance=0.8,  #设置百分比标签与圆心的距离
        labeldistance=1.1,  #设置教育水平标签与圆心的距离
        startangle=180,  #设置饼图的初始角度
        radius=1.2,  #设置饼图的半径
        counterclock=False,  #是否逆时针，这里设置为顺时针方向
        wedgeprops={'linewidth':1.5, 'edgecolor':'green'},  #设置饼图内外边界的属性值
        textprops={'fontsize':10, 'color':'black'},  #设置文本标签的属性值
        )

#添加图标题
plt.title('失信用户的受教育水平分布')
#显示图形
plt.show()

数据处理和算法内容

在《按赛题类型划分算法及代码（84种）》文件夹内搜索即可。下面列举的是课件内容。

3-1 熵权法

import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
import matplotlib.pyplot as plt
plt.rcParams['font.family']='SimHei'
x=pd.read_excel('待选材料参数.xlsx',engine='openpyxl').values
x[:,1]=-x[:,1]
x[:,3]=-x[:,3]
x=np.delete(x,-1,axis=0)
x=MinMaxScaler().fit(x).transform(x)
p=x/np.sum(x,axis=0)
e=np.zeros(x.shape[1])
for i in range(x.shape[1]):
    sum0=0
    for j in range(x.shape[0]):
        if p[j,i]==0:
            sum0+=0
        else:
            sum0+=p[j,i]*np.log(p[j,i])
    e[i]=-1/np.log(3)*sum0
w=(1-e)/np.sum(1-e)
defen=x@w.reshape(-1,1)
defen=defen.flatten()
a=np.sort(defen)

3-3 改进模拟退火算法

import numpy as np
from numpy import loadtxt,radians,sin,cos,inf,exp
from numpy import array,r_,c_,arange,savetxt
from numpy.lib.scimath import arccos
from numpy.random import shuffle,randint,rand
from matplotlib.pyplot import plot, show, rc
d=np.loadtxt('各个城市的最短路径.txt')

N=d.shape[0]

path=arange(N); L=inf
for j in range(1000):
    path0=arange(1,N-1); shuffle(path0)
    path0=r_[0,path0,N-1]; L0=d[0,path0[1]]  #初始化
    for i in range(1,N-1):L0+=d[path0[i],path0[i+1]]
    if L0<L: path=path0; L=L0
print(path,'\n',L)        
e=0.1**30; M=20000; at=0.999; T=1
gl=np.zeros(N-2)
for i in range(N-2):
    gl[i]=d[path[i],path[i+1]]
while T>e:
    pl=gl/gl.sum()
    c=np.random.choice(a=np.arange(1,N-1),size=2,replace=False,p=pl)
    c.sort()
    c1=c[0]; c2=c[1]
    df=d[path[c1-1],path[c2]]+d[path[c1],path[c2+1]]-\
    d[path[c1-1],path[c1]]-d[path[c2],path[c2+1]]  #续行
    if df<0:
        path=r_[path[0],path[1:c1],path[c2:c1-1:-1],path[c2+1:N+1]]; L=L+df
        for i in range(N-2):
            gl[i]=d[path[i],path[i+1]]
    else:
        if exp(-df/T)>=rand(1):
            path=r_[path[0],path[1:c1],path[c2:c1-1:-1],path[c2+1:N+1]]
            L=L+df
            for i in range(N-2):
                gl[i]=d[path[i],path[i+1]]
    T=T*at
print(path,'\n',L)  #输出巡航路径及路径长度

3-3 旅行商问题

from gurobipy import *
import numpy as np
mat=np.loadtxt('各个城市的最短路径.txt')
n=mat.shape[0]
m=Model()
x=m.addVars(n,n,vtype=GRB.BINARY,name='x')
u=m.addVars(n)
m.setObjective(sum(mat[i,j]*x[i,j] for i in range(n) for j in range(n) if i!=j),GRB.MINIMIZE)
m.addConstrs(sum(x[i,j] for i in range(n) if i!=j)==1 for j in range(n))
m.addConstrs(sum(x[i,j] for j in range(n) if i!=j)==1 for i in range(n))
m.addConstrs(u[i]-u[j]+n*x[i,j]<=n-1 for i in range(1,n) for j in range(1,n) if i!=j)
m.optimize()

3-3 粒子群

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
plt.rc('font',family='SimHei')
plt.rc('axes',unicode_minus=False)
def fitness_func(X):
    A=10;pi=np.pi
    x=X[:,0];y=X[:,1]
    return 2*A+x**2-A*np.cos(2*pi*x)+y**2-A*np.cos(2*pi*y)

def velocity_update(V,X,pbest,gbest,c1,c2,w,max_val):
    size=X.shape[0]
    r1=np.random.random((size,1))
    r2=np.random.random((size,1))
    V=w*V+c1*r1*(pbest-X)+c2*r2*(gbest-X)
    V[V<-max_val]=-max_val
    V[V>max_val]=max_val
    return V
def position_update(X,V):
    return X+V
w=1;c1=2;c2=2;r1=None;r2=None;dim=2;size=20;iter_num=1000
max_val=0.5;best_fitness=float(9e10)
fitness_value_list=[]
X=np.random.uniform(-5,5,size=(size,dim))
V=np.random.uniform(-0.5,0.5,size=(size,dim))
p_fitness=fitness_func(X)
g=p_fitness.min()
fitness_value_list.append(g)
pbest=X
gbest=X[p_fitness.argmin()]
for i in range(1,iter_num):
    V=velocity_update(V,X,pbest,gbest,c1,c2,w,max_val)
    X=position_update(X,V)
    p_fitness2=fitness_func(X)
    g_fitness2=p_fitness2.min()
    for j in range(size):
        if p_fitness[j]>p_fitness2[j]:
            pbest[j]=X[j]
            p_fitness[j]=p_fitness2[j]
    if g>g_fitness2:
        gbest=X[p_fitness2.argmin()]
        g=g_fitness2
    fitness_value_list.append(g)
print(fitness_value_list[-1])
print(gbest)
plt.plot(fitness_value_list,color='r')
plt.title('迭代过程')
plt.show()

3-3 遗传算法

import numpy as np
from numpy.random import randint, rand, shuffle
from matplotlib.pyplot import plot, show, rc
a=np.loadtxt("各个城市的最短路径.txt")
d=a;N=len(d)
w=50; g=10  #w为种群的个数，g为进化的代数
J=[]; 
for i in np.arange(w):
    c=np.arange(1,N-1); shuffle(c)
    c1=np.r_[0,c,101]; flag=1
    while flag>0:
        flag=0
        for m in np.arange(1,N-3):
            for n in np.arange(m+1,N-2):
                if d[c1[m],c1[n]]+d[c1[m+1],c1[n+1]]<\
                    d[c1[m],c1[m+1]]+d[c1[n],c1[n+1]]:
                    c1[m+1:n+1]=c1[n:m:-1]; flag=1
    c1[c1]=np.arange(N); J.append(c1)
J=np.array(J)/(N-1)
for k in np.arange(g):
    A=J.copy()
    c1=np.arange(w); shuffle(c1) 
    c2=randint(2,100,w)  #交叉点的数据
    for i in np.arange(0,w,2):
        temp=A[c1[i],c2[i]:N-1]  #保存中间变量
        A[c1[i],c2[i]:N-1]=A[c1[i+1],c2[i]:N-1]
        A[c1[i+1],c2[i]:N-1]=temp
    B=A.copy()
    by=[]  #初始化变异染色体的序号
    while len(by)<1: by=np.where(rand(w)<0.1)
    by=by[0]; B=B[by,:]
    G=np.r_[J,A,B]
    ind=np.argsort(G,axis=1)  #把染色体翻译成0,1，…，101
    NN=G.shape[0]; L=np.zeros(NN)
    for j in np.arange(NN):
        for i in np.arange(101):
            L[j]=L[j]+d[ind[j,i],ind[j,i+1]]
    ind2=np.argsort(L)
    J=G[ind2,:]
path=ind[ind2[0],:]; zL=L[ind2[0]]
print("所求的巡航路径长度为：",zL)

3-4 图论与路径

最大流

import numpy as np
import networkx as nx
import pylab as plt
L=[(1,2,5),(1,3,3),(2,4,2),(3,2,1),(3,5,4),(4,3,1),(4,5,3),(4,6,2),(5,6,5)]
G=nx.DiGraph()
for k in range(len(L)):
    G.add_edge(L[k][0]-1,L[k][1]-1, capacity=L[k][2])
value, flow_dict= nx.maximum_flow(G, 0, 5)
print("最大流的流量为：",value)
print("最大流为：", flow_dict)
n = len(flow_dict)
adj_mat = np.zeros((n, n), dtype=int)
for i, adj in flow_dict.items():
    for j, weight in adj.items():
        adj_mat[i,j] = weight
print("最大流的邻接矩阵为：\n",adj_mat)
ni,nj=np.nonzero(adj_mat)  #非零弧的两端点编号
key=range(n)
s=['v'+str(i+1) for i in range(n)]
s=dict(zip(key,s)) #构造用于顶点标注的字符字典
plt.rc('font',size=16)
pos=nx.shell_layout(G)  #设置布局
w=nx.get_edge_attributes(G,'capacity')
nx.draw(G,pos,font_weight='bold',labels=s,node_color='r')
nx.draw_networkx_edge_labels(G,pos,edge_labels=w)
path_edges=list(zip(ni,nj))
nx.draw_networkx_edges(G,pos,edgelist=path_edges,
            edge_color='r',width=3)
plt.show()

最小流

import numpy as np
import networkx as nx
L=[(1,2,5,3),(1,3,3,6),(2,4,2,8),(3,2,1,2),(3,5,4,2),
   (4,3,1,1),(4,5,3,4),(4,6,2,10),(5,6,5,2)]
G=nx.DiGraph()
for k in range(len(L)):
    G.add_edge(L[k][0]-1,L[k][1]-1, capacity=L[k][2], weight=L[k][3])
mincostFlow=nx.max_flow_min_cost(G,0,5)
print("所求流为：",mincostFlow)
mincost=nx.cost_of_flow(G, mincostFlow)
print("最小费用为：", mincost)
flow_mat=np.zeros((6,6),dtype=int)
for i,adj in mincostFlow.items():
    for j,f in adj.items():
        flow_mat[i,j]=f
print("最小费用最大流的邻接矩阵为：\n",flow_mat)

最小生成树

import numpy as np
import networkx as nx
import pylab as plt
L=[(1,2,8),(1,3,4),(1,5,2),(2,3,4),(3,4,2),(3,5,1),(4,5,5)]
b=nx.Graph()
b.add_nodes_from(range(1,6))
b.add_weighted_edges_from(L)
T=nx.minimum_spanning_tree(b)  #返回可迭代对象
w=nx.get_edge_attributes(T,'weight') #提取字典数据
TL=sum(w.values())  #计算最小生成树的长度
print("最小生成树为:",w)
print("最小生成树的长度为：",TL)
pos=nx.shell_layout(b)
nx.draw(T,pos,node_size=280,with_labels=True,node_color='r')
nx.draw_networkx_edge_labels(T,pos,edge_labels=w)
plt.show()

最短路径

import numpy as np
from numpy import inf
import pandas as pd
import time
def folyd(mat):
    n=len(mat)
    for k in range(n):
        for i in range(n):
            for j in range(n):
                if mat[i][k]+mat[k][j]<mat[i][j]:
                    mat[i][j]=mat[i][k]+mat[k][j]
    return mat
def dijkstra(mat,start):
    n=len(mat)
    dis=[];temp=[];
    dis.extend(mat[start]);temp.extend(mat[start])
    temp[start]=inf
    parents=[start]*n;visited=[start];
    for i in range(n-1):
        i=temp.index(min(temp))
        temp[i]=inf
        visited.append(i)
        print(visited)
        for j in range(n):
            if j not in visited:
                if mat[i][j]+dis[i]<dis[j]:
                    temp[j]=dis[j]=mat[i][j]+dis[i]
                    parents[j]=i
        path=[i]
        k=i
        while k !=start:
            path.append(parents[k])
            k=parents[k]
        path.reverse()
        print(str(path))
        print(str(i),':','->'.join(str(path).split()))
    print(dis)
a=np.loadtxt('ali535.txt')
time0=time.time()
mat=folyd(a)
time1=time.time()
print('已完成:',time1-time0)
time0=time.time()
for i in range(mat.shape[0]):
    dijkstra(a,i)
time1=time.time()
print('已完成:',time1-time0)
pd.DataFrame(mat).to_excel('最短路径.xlsx')

3-5 回归模型

一元线性回归

% 清除环境变量
clear; clc; close all;

x=[143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164]';
X=[ones(16,1) x] ;
y=[88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102]';
n=16;

% 绘制散点图
figure;
scatter(x, y, 'filled');
title('散点图');
xlabel('X值');
ylabel('Y值');
grid on;

% 线性回归拟合
[b, bint, r, rint, stats] = regress(y, [X]);


% 计算R平方和调整R平方
y_pred = b(1) + b(2) * x;
SS_res = sum((y - y_pred).^2);
SS_tot = sum((y - mean(y)).^2);
R2 = 1 - SS_res / SS_tot;

b
stats



% 打印结果
fprintf('截距(beta_0): %f\n', b(1));
fprintf('斜率(beta_1): %f\n', b(2));
fprintf('R平方: %f\n', R2);
% fprintf('调整R平方: %f\n', R2_adj);

% 绘制拟合曲线
hold on;
x_fit = linspace(min(x), max(x), 100);
y_fit = b(1) + b(2) * x_fit;
plot(x_fit, y_fit, 'r-', 'LineWidth', 2);

%%
hold on;
rcoplot(r,rint)
hold off;

%%

% 残差分析
residuals = y - y_pred;
figure;
subplot(2, 1, 1);
plot(residuals);
title('残差图');
xlabel('观测序号');
ylabel('残差');

% 残差的直方图
subplot(2, 1, 2);
histogram(residuals, 20);
title('残差直方图');
xlabel('残差');
ylabel('频率');

多元线性回归

clc;
clear;
x1=[3.5 5.3 5.1 5.8 4.2 6.0 6.8 5.5 3.1 7.2 4.5 4.9 8.0 6.5 6.5 3.7 6.2 7.0 4.0 4.5 5.9 5.6 4.8 3.9];
x2=[9 20 18 33 31 13 25 30 5 47 25 11 23 35 39 21 7 40 35 23 33 27 34 15];
x3=[6.1 6.4 7.4 6.7 7.5 5.9 6.0 4.0 5.8 8.3 5.0 6.4 7.6 7.0 5.0 4.0 5.5 7.0 6.0 3.5 4.9 4.3 8.0 5.0];
Y=[33.2 40.3 38.7 46.8 41.4 37.5 39.0 40.7 30.1 52.9 38.2 31.8 43.3 44.1 42.5 33.6 34.2 48.0 38.0 35.9 40.4 36.8 45.2 35.1];
n=24; m=3;
X=[ones(n,1),x1',x2',x3'];
[b,bint,r,rint,stats]=regress(Y',X,0.05);
% b,bint,r,rint,stats
b, stats

figure;
subplot(1,3,1),
plot(x1,Y,'r*'),
title('x1与y的散点图情况'); % 添加标题
grid on
subplot(1,3,2),
plot(x2,Y,'g+'),
title('x2与y的散点图情况'); % 添加标题
grid on
subplot(1,3,3),
plot(x3,Y,'bp'),
title('x3与y的散点图情况'); % 添加标题
grid on

figure;
rcoplot(r,rint)

非线性回归

clc;
clear;
x1 = [1000 600 1200 500 300 400 1300 1100 1300 300] ;
x2 = [5 7 6 6 8 7 5 4 3 9 ]; 
y = [100 75 80 70 50 65 90 100 110 60 ]';
x=[x1',x2'];
rstool(x , y, 'purequadratic')


%%
clc;
clear;
x1 = [1000 600 1200 500 300 400 1300 1100 1300 300] ;
x2 = [5 7 6 6 8 7 5 4 3 9 ]; 
y = [100 75 80 70 50 65 90 100 110 60 ]';
X = [ones(10,1) x1' x2' (x1.^2)' (x2.^2)'];
[b,bint,r,rint,stats]=regress(y,X,0.05);

3-6 有监督分类模型

KNN算法

from numpy import *
import operator
 
def createDataSet():
    group = array([[1.0, 1.1], [1.0, 1.0], [0, 0], [0, 0.1]])
    labels = ['A', 'A', 'B', 'B']
    return group, labels
 
group, labels = createDataSet()
print(group)
print(labels)
import matplotlib.pyplot as plt
 
def createDataSet():
    group = np.array([[1.0, 1.1], [1.0, 1.0], [0, 0], [0, 0.1]])
    labels = ['A', 'A', 'B', 'B']
    return group, labels
 
def plotDataSet(group, labels):
    label_dict = {'A': 'red', 'B': 'blue'}
    colors = [label_dict[label] for label in labels]
    fig, ax = plt.subplots()
    ax.scatter(group[:, 0], group[:, 1], c=colors)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title('Data Set')
    plt.show()
 
group, labels = createDataSet()
plotDataSet(group, labels)

def classify0(inX, dataSet, labels, k):
    #numpy函数shape[0]返回dataSet的行数
    dataSetSize = dataSet.shape[0]
    #在列向量方向上重复inX共1次(横向),行向量方向上重复inX共dataSetSize次(纵向)
    diffMat = np.tile(inX, (dataSetSize, 1)) - dataSet
    #二维特征相减后平方
    sqDiffMat = diffMat**2
    #sum()所有元素相加,sum(0)列相加,sum(1)行相加
    sqDistances = sqDiffMat.sum(axis=1)
    #开方,计算出距离
    distances = sqDistances**0.5
    #返回distances中元素从小到大排序后的索引值
    sortedDistIndices = distances.argsort()
    #定一个记录类别次数的字典
    classCount = {}
    for i in range(k):
        #取出前k个元素的类别
        voteIlabel = labels[sortedDistIndices[i]]
        #dict.get(key,default=None),字典的get()方法,返回指定键的值,如果值不在字典中返回默认值。
        #计算类别次数
        classCount[voteIlabel] = classCount.get(voteIlabel,0) + 1
    #python3中用items()替换python2中的iteritems()
    #key=operator.itemgetter(1)根据字典的值进行排序
    #key=operator.itemgetter(0)根据字典的键进行排序
    #reverse降序排序字典
    sortedClassCount = sorted(classCount.items(),key=operator.itemgetter(1),reverse=True)
    # import pdb
    # pdb.set_trace()
    #返回次数最多的类别,即所要分类的类别
    return sortedClassCount[0][0]

SVM算法、朴素贝叶斯算法

参考代码文件夹内容

3-7 无监督分类算法

DBSCAN

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from sklearn import metrics
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler

# 设置matplotlib支持中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
plt.rcParams['axes.unicode_minus'] = False  # 正常显示负号

# 生成样本数据
centers = [[1, 1], [-1, -1], [1, -1]]
X, labels_true = make_blobs(n_samples=750, centers=centers, cluster_std=0.4,
                            random_state=0)

X = StandardScaler().fit_transform(X)

# 计算k-距离图
min_samples = 10
k_distances = []
for i in range(len(X)):
    distances = np.linalg.norm(X - X[i], axis=1)
    distances = np.sort(distances)
    k_distance = distances[min_samples - 1]  # 第min_samples个最近邻的距离
    k_distances.append(k_distance)

# 排序k-距离
sorted_k_distances = sorted(k_distances)

# 使用选定的eps值计算DBSCAN
eps = 0.3  # 从k-距离图中选择的eps值
db = DBSCAN(eps=eps, min_samples=min_samples).fit(X)
core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
core_samples_mask[db.core_sample_indices_] = True
labels = db.labels_

# 计算标签中的聚类数量，忽略噪声点
n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
n_noise_ = list(labels).count(-1)

print('估计的聚类数量: %d' % n_clusters_)
print('估计的噪声点数量: %d' % n_noise_)
print("同质性: %0.3f" % metrics.homogeneity_score(labels_true, labels))
print("完整性: %0.3f" % metrics.completeness_score(labels_true, labels))
print("V-测量: %0.3f" % metrics.v_measure_score(labels_true, labels))
print("调整后的Rand指数: %0.3f"
      % metrics.adjusted_rand_score(labels_true, labels))
print("调整后的互信息: %0.3f"
      % metrics.adjusted_mutual_info_score(labels_true, labels))
print("轮廓系数: %0.3f"
      % metrics.silhouette_score(X, labels))

# 创建画布
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

# 绘制k-距离图
ax1.plot(range(1, len(sorted_k_distances) + 1), sorted_k_distances, marker='o')
ax1.set_xlabel('数据点索引')
ax1.set_ylabel('k-距离')
ax1.set_title('K-距离图')
ax1.grid(True)

# 绘制聚类结果
# 黑色用于表示噪声点
unique_labels = set(labels)
colors = ['#FF0000', '#00FF00', '#0000FF', '#FFFF00', '#00FFFF', '#FF00FF', '#800080']  # 自定义颜色列表
colors = colors[:len(unique_labels)]  # 确保颜色列表长度与聚类数量一致
for k, col in zip(unique_labels, colors):
    if k == -1:
        # 黑色用于噪声点
        col = 'black'

    class_member_mask = (labels == k)

    # 核心样本
    xy = X[class_member_mask & core_samples_mask]
    ax2.scatter(xy[:, 0], xy[:, 1], c=col, edgecolors='k', s=100, marker='s')  # 使用正方形标记

    # 非核心样本
    xy = X[class_member_mask & ~core_samples_mask]
    ax2.scatter(xy[:, 0], xy[:, 1], c=col, edgecolors='k', s=60, marker='o')  # 使用圆形标记

ax2.set_title('估计的聚类数量: %d' % n_clusters_)
plt.tight_layout()  # 调整子图间距
plt.show()

K-means

# #!/usr/bin/env python
# # -*- coding: utf-8 -*-
# # @Time    : 2024/8/9 15:21
# # @Author  : iszhou
# # @Email   :
# # @File    : demo.py
# # @purpose :
# # @software: PyCharm
# # @install packages: pip install -i https://pypi.tuna.tsinghua.edu.cn/simple bag_name
#


import numpy as np
from sklearn.datasets import make_blobs
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
import matplotlib.pyplot as plt

# 生成模拟数据
X, y_true = make_blobs(n_samples=300, centers=5, cluster_std=0.50, random_state=12345)

# 数据标准化
scaler = StandardScaler()
X = scaler.fit_transform(X)

# 肘部法则确定K值
distortions = []
silhouette_scores = []  # 添加轮廓系数列表
K = range(2, 10)  # K 的范围从 2 开始，因为 K=1 没有意义

for k in K:
    kmeanModel = KMeans(n_clusters=k).fit(X)
    distortions.append(kmeanModel.inertia_)
    silhouette_scores.append(silhouette_score(X, kmeanModel.labels_))  # 计算轮廓系数

# 创建一个大的图形窗口
fig, axs = plt.subplots(1, 3, figsize=(18, 6))

# 绘制肘部图
axs[0].plot(K, distortions, 'bx-')
axs[0].set_xlabel('k')
axs[0].set_ylabel('Distortion')
axs[0].set_title('The Elbow Method showing the optimal k')

# 绘制轮廓系数图
axs[1].plot(K, silhouette_scores, 'gx-')
axs[1].set_xlabel('k')
axs[1].set_ylabel('Silhouette Score')
axs[1].set_title('Silhouette Scores for Different k Values')

# 自动选择最佳K值
optimal_k = silhouette_scores.index(max(silhouette_scores)) + 2  # 找到最大轮廓系数对应的 K 值
print(f"Optimal number of clusters: {optimal_k}")

# 使用确定的K值进行聚类
kmeans = KMeans(n_clusters=optimal_k, random_state=0)
kmeans.fit(X)
labels = kmeans.labels_
centroids = kmeans.cluster_centers_

# 可视化聚类结果
axs[2].scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis')
axs[2].scatter(centroids[:, 0], centroids[:, 1], c='red', s=300, alpha=0.5)
axs[2].set_title('K-means Clustering with Optimal K')
axs[2].set_xlabel('Feature 1')
axs[2].set_ylabel('Feature 2')

# 显示整个图形
plt.tight_layout()
plt.show()

层次聚类

# 综合分类数据集
import numpy as np
from sklearn import metrics
from sklearn.cluster import AgglomerativeClustering
from sklearn.datasets import make_blobs
from matplotlib import pyplot as plt

# 创建数据集
# X为样本特征，Y为样本簇类别， 共1000个样本，每个样本2个特征，共4个簇，
# 簇中心在[-1,-1], [0,0],[1,1], [2,2]， 簇方差分别为[0.4, 0.2, 0.2, 0.2]
X, y = make_blobs(n_samples=1000, n_features=2,
                  centers=[[-1, -1], [0, 0], [1, 1], [2, 2]],
                  cluster_std=[0.4, 0.3, 0.2, 0.1],
                  random_state=9)
# 数据可视化
plt.scatter(X[:, 0], X[:, 1], marker='o')
plt.show()

# 定义颜色列表
colors = ['#FF0000', '#00FF00', '#0000FF', '#FFFF00', '#00FFFF', '#FF00FF', '#800080']

for index, metric in enumerate(["cosine", "euclidean", "cityblock"]):
    model = AgglomerativeClustering(
        n_clusters=4, linkage="average", affinity=metric
    )
    model.fit(X)
    print("%s Silhouette Coefficient: %0.3f"
          % (metric, metrics.silhouette_score(X, model.labels_, metric='sqeuclidean')))

    plt.figure()
    plt.axes([0, 0, 1, 1])
    for l, c in zip(np.arange(model.n_clusters), colors[:model.n_clusters]):
        row_ix = np.where(l == model.labels_)
        plt.scatter(X[row_ix, 0], X[row_ix, 1], c=c)
    plt.axis("tight")
    plt.axis("off")
    plt.suptitle("AgglomerativeClustering(affinity=%s)" % metric, size=20)

plt.show()


import scipy.cluster.hierarchy as shc

plt.figure(figsize =(8, 8))
plt.title('Visualising the data')
Dendrogram = shc.dendrogram((shc.linkage(X, method ='ward')))
plt.show()

时间序列模型

ARIMA

clc;clear;
%% 1. 读取数据 - 请将'B.xlsx'替换为您的数据文件名，并将'data(:,2)'根据要预测的列确定
data = readmatrix('sj2.xlsx');
time_series_data = data(:,2);
%data（x,y） x是行数 y是列数
%data（:,y）就是指y列对应的所有行的值组成的一个向量
%data(:,[y1:y2]) 就是指y1到y2列的对应的所有行的值组成的一个矩阵

%% 2. 划分训练集和测试集 - 这里使用80%的数据作为训练集，可以根据需要调整比例
train_size = round(length(time_series_data) * 0.8); 
%round(x)将x的每个元素四舍五入为最近的整数,length(x)就是求向量的长度,size(x)输出行列数
train_data = time_series_data(1:train_size);
test_data = time_series_data(train_size+1:end);
 
%% 3. 初始化最小AIC和BIC以及最优参数 - 选择模型参数的范围（p、d、q的最大值）
max_p = 5;
max_d = 2;%在ARIMA模型中，超参数d最常见的取值是0、1、2这些很小的数字。
max_q = 5;%在ARIMA模型当中，p和q的值往往取值不高，一般是[1,5]以内的正整数。
min_aic = Inf;
min_bic = Inf;
best_p = 0;
best_d = 0;
best_q = 0;
%实践中更常用的方法是从最小值p=1、q=1的方式开始进行尝试，不断改变p和q的取值，直到模型通过检验或达到我们需要的精度要求。
%% 4. 循环遍历不同的p, d, q值，尝试拟合ARIMA模型，并计算AIC和BIC
for p = 0:max_p
    for d = 0:max_d
        for q = 0:max_q
            % 创建ARIMA模型
            Mdl = arima(p, d, q);%创建不同p,d,q值下的arima模型，目的是为了找到最佳p,d,q值
            % 拟合模型，并计算AIC和BIC
            try   %里面语句报错时停止执行，直接跳到catch
                [EstMdl,~,logL] = estimate(Mdl, train_data, 'Display', 'off');%用时间序列对象的训练集拟合arima模型
                %[EstMdl，EstParamCov，logL，info] = estimate(model,data)model为指定模型，data为数据样本。返回估计参数，EstParamCov—关联的估计方差-协方差矩阵、logL—优化的对数似然目标函数。
                [aic, bic] = aicbic(logL, p + q + 1, length(train_data));
                %[aic, bic] = aicbic(logL,k,n)logL是模型的对数似然函数值,k为模型的参数数量,n为观测数据的数量
                %AIC和BIC是模型选择的常用准则，用于评估不同模型的拟合能力和复杂度。
                %AIC和BIC的值越小，表示模型的拟合效果越好且更简洁。
            catch
                continue;%会终止当前循环，进行下一次循环
            end
 
            % 更新最优参数
            if bic < min_bic
                %AIC是一个惩罚项，用于防止过度拟合，它的值越小表示模型越好。
                %BIC也类似，但是在惩罚项的计算中对参数数量的贡献更大，因此BIC更倾向于选择更简单的模型。
                min_aic = aic;
                min_bic = bic;
                best_p = p;
                best_d = d;
                best_q = q;
            end
        end
    end
end
 
%% 5. 使用最优参数创建ARIMA模型
best_mdl = arima(best_p, best_d, best_q);

%% 6. 拟合模型
EstMdl = estimate(best_mdl, train_data);
%返回指定的ARIMA模型EstMdl  
%% 7. 对测试集数据后的值进行预测 - 设定预测步长
num_steps =  10; % 预测测试集之后的10天数据
[forecast,forecast_RMSE] = forecast(EstMdl, num_steps, 'Y0', train_data);
%预测未来n个时间步长的值，其中Y0为数据初始化模型 
%forecast预测的未来值向量,forecast_RMSE是包含预测区间的置信区间的矩阵

%% 7. 输出预测结果
disp(['预测结果（', num2str(num_steps), '个步长）:']);
disp(forecast);
%% 8. 可视化预测结果
figure;
hold on;
plot(data(:,1), data(:,2),'k', 'LineWidth', 1);hold on
%plot(time_series_data, 'k', 'LineWidth', 1);hold on
%plot(train_size+1:train_size+length(test_data), test_data, 'b', 'LineWidth', 1); hold on% 绘制测试集数据
plot(data(train_size+1,1):data(train_size+length(test_data),1), test_data, 'b', 'LineWidth', 1); hold on% 绘制测试集数据
plot(data(train_size+1,1):data(train_size+num_steps,1), forecast, 'r', 'LineWidth', 1);hold on % 绘制预测数据
 
xlim([data(1,1), data(length(time_series_data),1)]);
title('ARIMA 时间序列预测');
xlabel('时间');
ylabel('值');
legend('实际数据', '测试集数据', '预测', 'Location', 'best');
 
% 10. 输出模型参数
disp(['最优模型参数: p = ', num2str(best_p), ', d = ', num2str(best_d), ', q = ', num2str(best_q)]);
disp(['最小 AIC: ', num2str(min_aic)]);
disp(['最小 BIC: ', num2str(min_bic)]);

灰色预测（GM）

clc;clear;
%% 1.将原始数据进行累加，形成有规律的数据列
syms a u;%创建符号标量变量a和u
c=[a,u]';%构成矩阵
A=[15 16.1 17.3 18.4 18.7 19.6 19.9 21.3 22.5];%输入数据，可以修改
Ago=cumsum(A);%原始数据一次累加,得到1-AGO序列xi(1)。
%如果A是一个向量，cumsum(A)返回一个向量，该向量中第m行的元素是A中第1行到第m行的所有元素累加和；
%如果A是一个矩阵，cumsum(A)返回一个和A同行同列的矩阵，矩阵中第m行第n列元素是A中第1行到第m行的所有第n列元素的累加和；

%% 2.对累加数列进行平滑，生成均值序列，取α为0.5
n=length(A);%原始数据个数
for k=1:(n-1)
    Z(k)=(Ago(k)+Ago(k+1))/2; %Z(i)为xi(1)的紧邻均值生成序列
end

%% 3.利用最小二乘法求解参数a和u
Yn = A;%Yn为常数项向量
Yn(1)=[]; %从第二个数开始，即x(2),x(3)...
Yn=Yn';
E=[-Z;ones(1,n-1)]';%累加生成数据做均值
c=(E'*E)\(E'*Yn);%利用公式求出a，u
%左除：a\b表示矩阵a的逆乘以b。右除：a/b表示矩阵a乘以矩阵b的逆
c= c';
a=c(1)%得到a的值
u=c(2)%得到u的值

%% 4.求解累加序列
F=[];
F(1)=A(1);
for k=2:n
    F(k)=(A(1)-u/a)/exp(a*(k-1))+u/a;%求出GM(1,1)模型公式
end

%% 5.还原成原序列即得到预测函数
G=[];
G(1)=A(1);
for k=2:n
    G(k)=F(k)-F(k-1);%两者做差还原原序列，得到预测数据
end

%% 6.可视化结果
t1=1:n;
t2=1:n;
plot(t1,A,'bo--');
hold on;
plot(t2,G,'r*-'); 
title('预测结果');
legend('真实值','预测值');

%% 7.模型检验
e=A-G;%误差
q=e/A;%相对误差
s1=var(A);%方差
s2=var(e);%残差的方差
c=s2/s1%后验差比值
len=length(e);
p=0;  %小误差概率
for i=1:len
    if(abs(e(i))<0.6745*s1)
        p=p+1;
    end
end
p=p/len
% 好 P>0.95 C<0.35
% 合格 P>0.80 C<0.45
% 勉强合格 P>0.70 C<0.50
% 不合格 P<=0.70 C>=0.65

3-9 机器学习模型算法

Logistic回归

clear
clc
%% 1.数据准备
%二分类 随机生成数据。  200个数据  每个数据2个特征
data=1*rand(300,2);%生成300个0~1之间（开环，不包含0和1两个数）均匀分布的伪随机数
label=zeros(300,1);
%label(sqrt(data(:,1).^2+data(:,2).^2)<8)=1;
label((data(:,2)+data(:,1)>1))=1;
%在data上加常数特征项；
data=[data,ones(size(data,1),1)];%size(A,1)该语句返回的时矩阵A的行数
 
%打乱顺序
randIndex = randperm(size(data,1));%返回一行包含从1到n的整数
data_new=data(randIndex,:);
label_new=label(randIndex,:);
 
%80%训练  20%测试
k=0.8*size(data,1);
X1=data_new(1:k,:);
Y1=label_new(1:k,:);
X2=data_new(k+1:end,:);
Y2=label_new(k+1:end,:);
 
[m1,n1] = size(X1);
[m2,n2] = size(X2);
Features=size(data,2);%特征个数
%% 2.开始训练
%设定学习率为1
delta=1;  
theta1=rand(1,Features);
%theta1=[.5,.5];
%%训练模型
 
%梯度下降算法求解theta（每次都是对全部的数据进行训练）
num = 300; %最大迭代次数
L=[];
while(num)
    dt=zeros(1,Features);
    loss=0;
    for i=1:m1    %训练集数据量
        xx=X1(i,1:Features); 
        yy=Y1(i,1);
        h=1/(1+exp(-(theta1 * xx')));  %得到预测值
        dt=dt+(h-yy) * xx;   %计算梯度
        loss=loss+ yy*log(h)+(1-yy)*log(1-h);
    end
    loss=-loss/m1;%计算风险函数的值
    L=[L,loss];
     
    theta2=theta1 - delta*dt/m1;
    theta1=theta2;
    num = num - 1;
     
    if loss<0.01   %误差小到一定程度跳出循环
        break;
    end
end
theta1
%% 3.可视化结果
figure
subplot(1,2,1)
plot(L)
title('loss')
 
subplot(1,2,2)
x=0:0.1:10;
y=(-theta1(1)*x-theta1(3))/theta1(2);
plot(x,y,'linewidth',2)
hold on
plot(data(label==1,1),data(label==1,2),'ro')
hold on
plot(data(label==0,1),data(label==0,2),'go')
axis([0 1 0 1])
 
%% 4.测试数据
acc=0;
for i=1:m2
    xx=X2(i,1:Features)';
    yy=Y2(i);
    finil=1/(1+exp(-theta2 * xx));
    if finil>0.5 && yy==1
        acc=acc+1;
    end
    if finil<=0.5 && yy==0
        acc=acc+1;
    end
end
acc/m2

BP神经网络

close all; clearvars; clear; %清空工作环境
data = readmatrix('sj1.xlsx');
P = data(:,1)';
T = data(:,2)';
%由于feedforwardnet函数自动对样本进行归一化和划分训练、验证、测试集,所以就不用手动将数据进行归一化处理，但不知道有没有打乱顺序
% net = feedforwardnet(5, 'traingd'); %是'5'是指隐含层有5个神经元，这里只有一个隐含层，多个隐含层神经元的个数设置为[5,3,...]
net = feedforwardnet([10,5], 'traingd');%两层隐含层，相应神经元个数分别为10和5
%traingd---梯度下降反向传播算法训练函数; traingda---自适应调整学习率的梯度下降反向传播算法训练函数；trainlm--- 中型网络，,内存需求最大,收敛速度最快
net.trainParam.lr = 0.01; %学习速率
net.trainParam.epochs = 1000; %最大训练次数
net.trainParam.goal = 1e-6; %最小误差，达到该精度，停止训练
net.trainParam.show = 50; %每50次展示训练结果
net = train(net, P, T); %net：需要训练的神经网络，P：网络输入，T：网络期望输出，并输出训练好的神经网络
Y = net(P); %计算预测值
perf = perform(net, Y, T);%计算误差性能
plot(P, T, 'r-');hold on
plot(P, Y ,'b-');hold on

LSTM预测

close all; clearvars; clear; %清空工作环境
load XTrain;
load XTest;
load YTrain;
load YTest;
%% 1.数据标准化
XTrain_mu = mean([XTrain{:}],2); %求均值
XTrain_sig = std([XTrain{:}],0,2); %求标准差
XTest_mu = mean([XTest{:}],2);
XTest_sig = std([XTest{:}],0,2);
YTrain_mu = mean([YTrain{:}],2);
YTrain_sig = std([YTrain{:}],0,2);
YTest_mu = mean([YTest{:}],2);
YTest_sig = std([YTest{:}],0,2);

for i = 1:numel(XTrain)  %返回数组A中元素个数
    XTrain{i} = (XTrain{i} - XTrain_mu) ./ XTrain_sig ;
    YTrain{i}=(YTrain{i} - YTrain_mu) ./ YTrain_sig;
end

for i = 1:numel(XTest)
    XTest{i}=(XTest{i} - XTest_mu) ./ XTest_sig;
    YTest{i}=(YTest{i} - YTest_mu) ./ YTest_sig;
end
%% 2.定义网络结构
%创建一个 LSTM 网络，该网络包含一个具有 200 个隐藏单元的 LSTM 层，
%然后是一个大小为 50 的全连接层和一个丢弃概率为 0.5 的丢弃层。
numResponses = size(YTrain{1},1);
numHiddenUnits = 200;
numFeatures=15;
layers = [ ...
    sequenceInputLayer(numFeatures)
    lstmLayer(numHiddenUnits,'OutputMode','sequence') %序列到序列回归
    fullyConnectedLayer(50)
    dropoutLayer(0.5)
    fullyConnectedLayer(numResponses)
    regressionLayer];
% 层级设置=[
% 序列输入层（特征数量）
% lstm层（隐藏单元个数,'输出模式','多对一'）
% 全连接层（目标个数）%目标个数即分类的类别数
% 创建丢弃层，丢弃层以给定的概率将输入元素随机设置为零
% 回归层计算回归任务的半均方误差损失]
%% 3.设定参数
%使用求解器 "adam"进行 90 轮训练。
%指定学习率为 0.01。要防止梯度爆炸，将梯度阈值设置为1。
%要使序列保持按长度排序，将 Shuffle 选项设置为 "never"。
%在图中显示训练进度并监控均方根误差 (RMSE) 度量
maxEpochs = 90;
options = trainingOptions('adam', ...
    'MaxEpochs',maxEpochs, ...
    'InitialLearnRate',0.01, ...
    'GradientThreshold',1, ...
    'Shuffle','never', ...
    'Plots','training-progress',...
    'Verbose',0);
% 选项 = 训练选项设置('累加器', ...
%     '最大训练轮次',maxEpochs, ...
%     '最小步距',miniBatchSize, ...
%     '序列长度','整个序列', ...
%     '乱序','否', ...
%     '画图','训练过程', ...
%     '在命令窗口中显示训练进度信息',否);

%% 4.模型训练
net = trainNetwork(XTrain,YTrain,layers,options);

% net = trainNetwork(imds,layers,options) 为图像分类问题训练网络。图像数据存储区 imds存储输入的图像数据， layers定义网络体系结构，并 options定义训练选项。
% net = trainNetwork(ds,layers,options) 使用数据存储训练网络ds。对于具有多个输入的网络，请将此语法与组合或转换后的数据存储区结合使用。
% net = trainNetwork(X,Y,layers,options) 为图像分类和回归问题训练网络。数字数组X包含预测变量，并 Y包含分类标签或数字响应。
% net = trainNetwork(sequences,Y,layers,options) 训练网络以解决序列分类和回归问题（例如LSTM或BiLSTM网络），其中sequences 包含序列或时间序列预测变量并Y包含响应。对于分类问题，Y是分类向量或分类序列的单元格数组。对于回归问题，Y是目标矩阵或数字序列的单元格数组。
% net = trainNetwork(tbl,layers,options) 为分类和回归问题训练网络。该表 tbl包含数字数据或数据的文件路径。预测变量必须位于的第一列中tbl。有关目标或响应变量的信息，请参见tbl。
% net = trainNetwork(tbl,responseName,layers,options) 为分类和回归问题训练网络。预测变量必须位于的第一列中tbl。该 responseName参数指定在响应变量tbl。
% [net,info] = trainNetwork(___) 还可以使用先前语法中的任何输入参数返回有关训练的信息。

%% 5.回归预测
YPred = predict(net,XTest);%将训练好的模型net和即将要预测的数据XTest输入，得到预测结果
%% 6.输出可视化
idx = randperm(numel(YPred),4); %随机输出4组预测数据
figure
for i = 1:numel(idx)
    subplot(2,2,i)
    plot(YTest{idx(i)},'--')
    hold on
    plot(YPred{idx(i)},'.-')
    hold off
    title("Test Observation " + idx(i))
    xlabel("Time Step")
    ylabel("青霉素浓度")
    rmse = sqrt(mean((YPred{i} - YTest{i}).^2)) %求解均方根误差评估模型性能
end
legend(["True" "Predicted"],'Location','southeast')

使用时需导入Xtrain.mat、Ytrain.mat

3-10

随机森林回归

close all  %关闭开启图窗
clear  %清空变量
clc    %清空命令行
% rng(2222)  %随机数种子固定结果
res =readmatrix('sj3.xlsx');

%% 数据归一化
X = res(:,1:end-1);
Y = res(:,end);
[x,psin] = mapminmax(X',0,1);%mapminmax对每一行进行规划，所以需要对X进行转置
[y,psout] = mapminmax(Y',0,1);

%% 划分训练集和测试集
num = size(res,1);%总样本数
state =randperm(num);%打乱样本，提高泛化性
ratio = 0.7;
train_num = floor(num*ratio);

x_train = x(:,state(1:train_num))';%工具箱需要每个特征竖着分布，因此需要转置
y_train = y(state(1:train_num))';
x_test = x(:,state(train_num+1:end))';
y_test = y(state(train_num+1:end))';

%% 训练模型
trees = 100;%决策树数目
leaf = 3;%定义叶子节点个数，最小叶子数，过小容易拟合
wuc = 'on';%打开误差图
Importance = 'on';%计算特征重要性
%regression——回归
net = TreeBagger(trees,x_train,y_train,'OOBPredictorImportance',Importance,...
    'Method','regression','OOBPrediction',wuc,'minleaf',leaf);
import = net.OOBPermutedPredictorDeltaError;

%% 预测
re1 = predict(net , x_train);
re2 = predict(net , x_test);

%% 数据反归一化
%实际值
Y_train = Y(state(1:train_num));
Y_test = Y(state(train_num+1:end));

%预测值
pre1 = mapminmax('reverse',re1,psout);
pre2 = mapminmax('reverse',re2,psout);
%%
%均方根误差
error1 = sqrt(mean(pre1 - Y_train).^2);
error2 = sqrt(mean(pre2 - Y_test).^2);

%相关指标计算
%R2
R1 = 1 - norm(Y_train - pre1)^2 / norm(Y_train -  mean(Y_train))^2;
R2 = 1 - norm(Y_test - pre2)^2 / norm(Y_test -  mean(Y_test))^2;

%MAE
mae1 = mean(abs(Y_train - pre1));
mae2 = mean(abs(pre2 - Y_test));

disp('训练集预测精度指标如下：')
disp(['训练集的R2：',num2str(R1)])
disp(['训练集的MAE：',num2str(mae1)])
disp(['训练集的RMSE：',num2str(error1)])
disp('测试集预测精度指标如下：')
disp(['测试集的R2：',num2str(R2)])
disp(['测试集的MAE：',num2str(mae2)])
disp(['测试集的RMSE：',num2str(error2)])
%%
figure
plot(1:train_num,Y_train,'r-O',1:train_num,pre1,'b-+','LineWidth',1)
legend('真实值','预测值')
xlabel('样本点')
ylabel('预测值')
title('训练集预测结果对比')
figure
plot(1:num-train_num,Y_test,'r-O',1:num-train_num,pre2,'b-+','LineWidth',1)
legend('真实值','预测值')
xlabel('样本点')
ylabel('预测值')
title('测试集预测结果对比')

%% 绘制误差曲线
figure
plot(1:trees,oobError(net),'r--O','LineWidth',1)
legend('误差迭代曲线')
xlabel('决策树（迭代次数）')
ylabel('误差')
grid
%% 绘制重要性
figure
bar(import,'green')
yticks([])
xlabel('特征')
ylabel('重要性')

%% 新数据预测
% newdata = xlsread("新的多输入");
% newy = newpre(newdata,psin,psout,net);
% figure
% plot(newy)
% xlabel('样本点')
% ylabel('预测值')
% xlswrite("新的输出",newy)
% 
% function y = newpre(newdata,psin,psout,net)
%    x = mapminmax('apply',newdata',psin);
%    pre = predict(net,x');
%    y = mapminmax('reverse',pre,psout);
% end