数学建模模型
综合评价，权重决策分析
步骤：
构建层次结构模型
构造判断矩阵/成对比较矩阵
层次单排序及其一致性检验
层次总排序及其一致性检验
缺点：主观性比较强
# A = [[1,1,4,1/3,3],
# [1,1,4,1/3,3],
# [1/4,1/4,1,1/3,1/2],
# [3,3,3,1,3],
# [1/3,1/3,2,1/3,1]]
import numpy as np #导入所需包并将其命名为np
def ConsisTest(X): #函数接收一个如上述A似的矩阵
#计算权重
#方法一：算术平均法
## 第一步：将判断矩阵按照列归一化（每个元素除以其所在列的和）
X = np.array(X) #将X转换为np.array对象
sum_X = X.sum(axis=0) #计算X每列的和
(n,n) = X.shape #X为方阵，行和列相同，所以用一个n来接收
sum_X = np.tile(sum_X,(n,1)) #将和向量重复n行组成新的矩阵
stand_X = X/sum_X #标准化X（X中每个元素除以其所在列的和）
## 第二步：将归一化矩阵每一行求和
sum_row = stand_X.sum(axis=1)
## 第三步：将相加后得到的向量中每个元素除以n即可得到权重向量
print("算数平均法求权重的结果为：")
print(sum_row/n)
#方法二：特征值法
## 第一步：找出矩阵X的最大特征值以及其对应的特征向量
V,E = np.linalg.eig(X) #V是特征值，E是特征值对应的特征向量
max_value = np.max(V) #最大特征值
#print("最大特征值是：",max_value)
max_v_index = np.argmax(V) #返回最大特征值所在位置
max_eiv = E[:,max_v_index] #最大特征值对应的特征向量
## 第二步：对求出的特征向量进行归一化处理即可得到权重
stand_eiv = max_eiv/max_eiv.sum()
print("特征值法求权重的结果为：")
print(stand_eiv)
print("———————————————————————————————")
#一致性检验
## 第一步：计算一致性指标CI
CI = (max_value-n)/(n-1)
## 第二步：查找对应的平均随机一致性指标RI
RI = np.array([0, 0, 0.58, 0.9, 1.12, 1.24, 1.32, 1.41, 1.45, 1.49])
## 第三步：计算一致性比例CR
CR = CI/RI[n]
if CR < 0.1:
print("CR=",CR,"，小于0.1，通过一致性检验")
else:
print("CR=",CR,"，大于等于0.1，没有通过一致性检验，请修改判断矩阵")
return None
#ConsisTest(A)
https://zhuanlan.zhihu.com/p/613113565
https://blog.csdn.net/ddjhpxs/article/details/105407103
https://zhuanlan.zhihu.com/p/266405027
一般是输入时间序列,适合数据量少的时候
示例代码
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
data = pd.read_excel("E:/桌面/灰色预测.xlsx",sheet_name=1)
list1 = data["序列"]
tlist1 = [2014,2015,2016,2017,2018]
fy = []
def GM11(x,n):
'''
灰色预测
x：序列，numpy对象
n:需要往后预测的个数
'''
x1 = x.cumsum()#一次累加
z1 = (x1[:len(x1) - 1] + x1[1:])/2.0#紧邻均值
z1 = z1.reshape((len(z1),1))
B = np.append(-z1,np.ones_like(z1),axis=1)
Y = x[1:].reshape((len(x) - 1,1))
#a为发展系数 b为灰色作用量
[[a],[b]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)#计算待估参数
result = (x[0]-b/a)*np.exp(-a*(n-1))-(x[0]-b/a)*np.exp(-a*(n-2)) #预测方程
S1_2 = x.var()#原序列方差
e = list()#残差序列
for index in range(1,x.shape[0]+1):
predict = (x[0]-b/a)*np.exp(-a*(index-1))-(x[0]-b/a)*np.exp(-a*(index-2))
e.append(x[index-1]-predict)
print(predict) #预测值
fy.append(predict)
print("后验差检验")
S2_2 = np.array(e).var()#残差方差
C = S2_2/S1_2#后验差比
if C<=0.35:
assess = '后验差比<=0.35，模型精度等级为好'
elif C<=0.5:
assess = '后验差比<=0.5，模型精度等级为合格'
elif C<=0.65:
assess = '后验差比<=0.65，模型精度等级为勉强'
else:
assess = '后验差比>0.65，模型精度等级为不合格'
#预测数据
predict = list()
for index in range(x.shape[0]+1,x.shape[0]+n+1):
predict.append((x[0]-b/a)*np.exp(-a*(index-1))-(x[0]-b/a)*np.exp(-a*(index-2)))
predict = np.array(predict)
return {
'a':{'value':a,'desc':'发展灰数'},
'b':{'value':b,'desc':'控制灰数'},
'predict':{'value':result,'desc':'第%d个预测值'%n},
'C':{'value':C,'desc':assess},
'predict':{'value':predict,'desc':'往后预测%d个的序列'%(n)},
}
if __name__ == "__main__":
data = np.array(list1)
x = data[0:5]#输入数据
#y = data[0:]#需要预测的数据
result = GM11(x,1)
predict = result['predict']['value']
predict = np.round(predict,1)
#print('真实值:',x)
print('预测值:',predict)
print(result)
print("残差检验")
a = []
for i in range(5):
a.append(abs(fy[i]-list1[i]))
print('%.5f%%' % ((abs(fy[i]-list1[i]))/list1[i]))
print("关联度检验")
c= []
for i in range(5):
b = (min(a)+0.5*max(a))/(abs(a[i])+0.5*max(a))
c.append(b)
print("ρ = 0.5 关联度为：",np.mean(c))
#print(y)
print(predict)
fy.append(3.8)
#作图
import numpy as np
import matplotlib.pyplot as plt
x1 = np.array([2014,2015,2016,2017,2018])
y1 = np.array(x)
x2 = np.array([2014,2015,2016,2017,2018,2019])
y2 = np.array(fy)
plt.plot(x1,y1,'r*-',label='真实值曲线') #真实值
plt.plot(x2,y2,'b+-',label='预测值曲线') #预测值
plt.xlabel('年份')
plt.ylabel('值')
plt.legend()
plt.plot()
plt.show()
https://blog.csdn.net/weixin_53597801/article/details/121879893
示例1
def Dijkstra(network, s, d): # 迪杰斯特拉算法算s-d的最短路径，并返回该路径和值
print("Start Dijstra Path……")
path = [] # 用来存储s-d的最短路径
n = len(network) # 邻接矩阵维度，即节点个数
fmax = float('inf')
w = [[0 for _ in range(n)] for j in range(n)] # 邻接矩阵转化成维度矩阵，即0→max
book = [0 for _ in range(n)] # 是否已经是最小的标记列表
dis = [fmax for i in range(n)] # s到其他节点的最小距离
book[s - 1] = 1 # 节点编号从1开始，列表序号从0开始
midpath = [-1 for i in range(n)] # 上一跳列表
for i in range(n):
for j in range(n):
if network[i][j] != 0:
w[i][j] = network[i][j] # 0→max
else:
w[i][j] = fmax
if i == s - 1 and network[i][j] != 0: # 直连的节点最小距离就是network[i][j]
dis[j] = network[i][j]
for i in range(n - 1): # n-1次遍历，除了s节点
min = fmax
for j in range(n):
if book[j] == 0 and dis[j] < min: # 如果未遍历且距离最小
min = dis[j]
u = j
book[u] = 1
for v in range(n): # u直连的节点遍历一遍
if dis[v] > dis[u] + w[u][v]:
dis[v] = dis[u] + w[u][v]
midpath[v] = u + 1 # 上一跳更新
j = d - 1 # j是序号
path.append(d) # 因为存储的是上一跳，所以先加入目的节点d，最后倒置
while (midpath[j] != -1):
path.append(midpath[j])
j = midpath[j] - 1
path.append(s)
path.reverse() # 倒置列表
print("path:",path)
# print(midpath)
print("dis:",dis)
# return path
network = [[0, 1, 0, 2, 0, 0],
[1, 0, 2, 4, 3, 0],
[0, 2, 0, 0, 1, 4],
[2, 4, 0, 0, 6, 0],
[0, 3, 1, 6, 0, 2],
[0, 0, 4, 0, 2, 0]]
Dijkstra(network, 1, 6)
#运行结果
#Start Dijstra Path……
#path: [1, 2, 5, 6](这个是最短路径)
#dis: [2, 1, 3, 2, 4, 6]
https://blog.csdn.net/time_money/article/details/109899692
示例2
import networkx as nx
# 创建带权图
G = nx.Graph()
G.add_weighted_edges_from([(1, 2, 1), (1, 3, 2), (2, 3, 1), (2, 4, 3), (3, 4, 1)])
# 计算最短路径
path = nx.dijkstra_path(G, 1, 4)
distance = nx.dijkstra_path_length(G, 1, 4)
# 输出结果
print('Shortest path:', path)
print('Shortest distance:', distance)
#结果如下:
#Shortest path: [1, 2, 3, 4]
#Shortest distance: 3
示例1
import sys
sys.setrecursionlimit(100000000)
# 弗洛伊德算法
def floyd():
n = len(graph)
for k in range(n):
for i in range(n):
for j in range(n):
if graph[i][k] + graph[k][j] < graph[i][j]:
graph[i][j] = graph[i][k] + graph[k][j]
parents[i][j] = parents[k][j] # 更新父结点
# 打印路径
def print_path(i, j):
if i != j:
print_path(i, parents[i][j])
print(j, end='-->')
# Data [u, v, cost]
datas = [
[0, 1, 2],
[0, 2, 6],
[0, 3, 4],
[1, 2, 3],
[2, 0, 7],
[2, 3, 1],
[3, 0, 5],
[3, 2, 12],
]
n = 4
# 无穷大
inf = 9999999999
# 构图
graph = [[(lambda x: 0 if x[0] == x[1] else inf)([i, j]) for j in range(n)] for i in range(n)]
parents = [[i] * n for i in range(4)] # 关键地方，i-->j 的父结点初始化都为i
for u, v, c in datas:
graph[u][v] = c # 因为是有向图，边权只赋给graph[u][v]
#graph[v][u] = c # 如果是无向图，要加上这条。
floyd()
print('Costs:')
for row in graph:
for e in row:
print('∞' if e == inf else e, end='\t')
print()
print('\nPath:')
for i in range(n):
for j in range(n):
print('Path({}-->{}): '.format(i, j), end='')
print_path(i, j)
print(' cost:', graph[i][j])
# 最终的graph
# graph[i][j]表示i-->j的最短路径
# 0 2 5 4
# 9 0 3 4
# 6 8 0 1
# 5 7 10 0
#
# Path:
# Path(0-->0): 0--> cost: 0
# Path(0-->1): 0-->1--> cost: 2
# Path(0-->2): 0-->1-->2--> cost: 5
# Path(0-->3): 0-->3--> cost: 4
# Path(1-->0): 1-->2-->3-->0--> cost: 9
# Path(1-->1): 1--> cost: 0
# Path(1-->2): 1-->2--> cost: 3
# Path(1-->3): 1-->2-->3--> cost: 4
# Path(2-->0): 2-->3-->0--> cost: 6
# Path(2-->1): 2-->3-->0-->1--> cost: 8
# Path(2-->2): 2--> cost: 0
# Path(2-->3): 2-->3--> cost: 1
# Path(3-->0): 3-->0--> cost: 5
# Path(3-->1): 3-->0-->1--> cost: 7
# Path(3-->2): 3-->0-->1-->2--> cost: 10
# Path(3-->3): 3--> cost: 0
https://juejin.cn/post/6987592888318165028
示例2
import numpy as np
# 定义弗洛伊德算法函数
def floyd(graph):
n = len(graph)
dist = np.copy(graph)
for k in range(n):
for i in range(n):
for j in range(n):
if dist[i][j] > dist[i][k] + dist[k][j]:
dist[i][j] = dist[i][k] + dist[k][j]
return dist
# 定义带权图
graph = np.array([
[0, 1, 2, np.inf],
[1, 0, 1, 3],
[2, 1, 0, 1],
[np.inf, 3, 1, 0]
])
# 计算最短路径
dist = floyd(graph)
# 输出结果
print('Shortest distance matrix:')
print(dist)
#得到的是所有节点之间的最短路径矩阵
Shortest distance matrix:
[[0. 1. 2. 3.]
[1. 0. 1. 2.]
[2. 1. 0. 1.]
[3. 2. 1. 0.]]
https://pythonjishu.com/znmehrtcqshdjef/
弗洛伊德和迪杰斯特拉算法的区别:
https://blog.csdn.net/weixin_43872728/article/details/100662957
import math
from random import random
import matplotlib.pyplot as plt
def func(x, y): #函数优化问题
res= 4*x**2-2.1*x**4+x**6/3+x*y-4*y**2+4*y**4
return res
#x为公式里的x1,y为公式里面的x2
class SA:
def __init__(self, func, iter=100, T0=100, Tf=0.01, alpha=0.99):
self.func = func
self.iter = iter #内循环迭代次数,即为L =100
self.alpha = alpha #降温系数，alpha=0.99
self.T0 = T0 #初始温度T0为100
self.Tf = Tf #温度终值Tf为0.01
self.T = T0 #当前温度
self.x = [random() * 11 -5 for i in range(iter)] #随机生成100个x的值
self.y = [random() * 11 -5 for i in range(iter)] #随机生成100个y的值
self.most_best =[]
"""
random()这个函数取0到1之间的小数
如果你要取0-10之间的整数（包括0和10）就写成 (int)random()*11就可以了，11乘以零点多的数最大是10点多，最小是0点多
该实例中x1和x2的绝对值不超过5（包含整数5和-5），（random() * 11 -5）的结果是-6到6之间的任意值（不包括-6和6）
（random() * 10 -5）的结果是-5到5之间的任意值（不包括-5和5），所有先乘以11，取-6到6之间的值，产生新解过程中，用一个if条件语句把-5到5之间（包括整数5和-5）的筛选出来。
"""
self.history = {'f': [], 'T': []}
def generate_new(self, x, y): #扰动产生新解的过程
while True:
x_new = x + self.T * (random() - random())
y_new = y + self.T * (random() - random())
if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):
break #重复得到新解，直到产生的新解满足约束条件
return x_new, y_new
def Metrospolis(self, f, f_new): #Metropolis准则
if f_new <= f:
return 1
else:
p = math.exp((f - f_new) / self.T)
if random() < p:
return 1
else:
return 0
def best(self): #获取最优目标函数值
f_list = [] #f_list数组保存每次迭代之后的值
for i in range(self.iter):
f = self.func(self.x[i], self.y[i])
f_list.append(f)
f_best = min(f_list)
idx = f_list.index(f_best)
return f_best, idx #f_best,idx分别为在该温度下，迭代L次之后目标函数的最优解和最优解的下标
def run(self):
count = 0
#外循环迭代，当前温度小于终止温度的阈值
while self.T > self.Tf:
#内循环迭代100次
for i in range(self.iter):
f = self.func(self.x[i], self.y[i]) #f为迭代一次后的值
x_new, y_new = self.generate_new(self.x[i], self.y[i]) #产生新解
f_new = self.func(x_new, y_new) #产生新值
if self.Metrospolis(f, f_new): #判断是否接受新值
self.x[i] = x_new #如果接受新值，则把新值的x,y存入x数组和y数组
self.y[i] = y_new
# 迭代L次记录在该温度下最优解
ft, _ = self.best()
self.history['f'].append(ft)
self.history['T'].append(self.T)
#温度按照一定的比例下降（冷却）
self.T = self.T * self.alpha
count += 1
# 得到最优解
f_best, idx = self.best()
print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")
sa = SA(func)
sa.run()
plt.plot(sa.history['T'], sa.history['f'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()
https://blog.csdn.net/weixin_48241292/article/details/109468947
# -*- coding: utf-8 -*-
import math # 导入模块 math
import random # 导入模块 random
import pandas as pd # 导入模块 pandas 并简写成 pd
import numpy as np # 导入模块 numpy 并简写成 np YouCans
import matplotlib.pyplot as plt # 导入模块 matplotlib.pyplot 并简写成 plt
np.set_printoptions(precision=4)
pd.set_option('display.max_rows', 20)
pd.set_option('expand_frame_repr', False)
pd.options.display.float_format = '{:,.2f}'.format
# 子程序：初始化模拟退火算法的控制参数
def initParameter():
# custom function initParameter():
# Initial parameter for simulated annealing algorithm
tInitial = 100.0 # 设定初始退火温度(initial temperature)
tFinal = 1 # 设定终止退火温度(stop temperature)
nMarkov = 1000 # Markov链长度，也即内循环运行次数
alfa = 0.98 # 设定降温参数，T(k)=alfa*T(k-1)
return tInitial,tFinal,alfa,nMarkov
# 子程序：读取TSPLib数据
def read_TSPLib(fileName):
# custom function read_TSPLib(fileName)
# Read datafile *.dat from TSPlib
# return coordinates of each city by YouCans, XUPT
res = []
with open(fileName, 'r') as fid:
for item in fid:
if len(item.strip())!=0:
res.append(item.split())
loadData = np.array(res).astype('int') # 数据格式：i Xi Yi
coordinates = loadData[:,1::]
return coordinates
# 子程序：计算各城市间的距离，得到距离矩阵
def getDistMat(nCities, coordinates):
# custom function getDistMat(nCities, coordinates):
# computer distance between each 2 Cities
distMat = np.zeros((nCities,nCities)) # 初始化距离矩阵
for i in range(nCities):
for j in range(i,nCities):
# np.linalg.norm 求向量的范数（默认求 二范数），得到 i、j 间的距离
distMat[i][j] = distMat[j][i] = round(np.linalg.norm(coordinates[i]-coordinates[j]))
return distMat # 城市间距离取整（四舍五入）
# 子程序：计算 TSP 路径长度
def calTourMileage(tourGiven, nCities, distMat):
# custom function caltourMileage(nCities, tour, distMat):
# to compute mileage of the given tour
mileageTour = distMat[tourGiven[nCities-1], tourGiven[0]] # dist((n-1),0)
for i in range(nCities-1): # dist(0,1),...dist((n-2)(n-1))
mileageTour += distMat[tourGiven[i], tourGiven[i+1]]
return round(mileageTour) # 路径总长度取整（四舍五入）
# 子程序：绘制 TSP 路径图
def plot_tour(tour, value, coordinates):
# custom function plot_tour(tour, nCities, coordinates)
num = len(tour)
x0, y0 = coordinates[tour[num - 1]]
x1, y1 = coordinates[tour[0]]
plt.scatter(int(x0), int(y0), s=15, c='r') # 绘制城市坐标点 C(n-1)
plt.plot([x1, x0], [y1, y0], c='b') # 绘制旅行路径 C(n-1)~C(0)
for i in range(num - 1):
x0, y0 = coordinates[tour[i]]
x1, y1 = coordinates[tour[i + 1]]
plt.scatter(int(x0), int(y0), s=15, c='r') # 绘制城市坐标点 C(i)
plt.plot([x1, x0], [y1, y0], c='b') # 绘制旅行路径 C(i)~C(i+1)
plt.xlabel("Total mileage of the tour:{:.1f}".format(value))
plt.title("Optimization tour of TSP{:d}".format(num)) # 设置图形标题
plt.show()
# 子程序：交换操作算子
def mutateSwap(tourGiven, nCities):
# custom function mutateSwap(nCities, tourNow)
# produce a mutation tour with 2-Swap operator
# swap the position of two Cities in the given tour
# 在 [0,n) 产生 2个不相等的随机整数 i,j
i = np.random.randint(nCities) # 产生第一个 [0,n) 区间内的随机整数
while True:
j = np.random.randint(nCities) # 产生一个 [0,n) 区间内的随机整数
if i!=j: break # 保证 i, j 不相等
tourSwap = tourGiven.copy() # 将给定路径复制给新路径 tourSwap
tourSwap[i],tourSwap[j] = tourGiven[j],tourGiven[i] # 交换 城市 i 和 j 的位置————简洁的实现方法
return tourSwap
def main():
# 主程序
# # 读取旅行城市位置的坐标
coordinates = np.array([[565.0, 575.0], [25.0, 185.0], [345.0, 750.0], [945.0, 685.0], [845.0, 655.0],
[880.0, 660.0], [25.0, 230.0], [525.0, 1000.0], [580.0, 1175.0], [650.0, 1130.0],
[1605.0, 620.0], [1220.0, 580.0], [1465.0, 200.0], [1530.0, 5.0], [845.0, 680.0],
[725.0, 370.0], [145.0, 665.0], [415.0, 635.0], [510.0, 875.0], [560.0, 365.0],
[300.0, 465.0], [520.0, 585.0], [480.0, 415.0], [835.0, 625.0], [975.0, 580.0],
[1215.0, 245.0], [1320.0, 315.0], [1250.0, 400.0], [660.0, 180.0], [410.0, 250.0],
[420.0, 555.0], [575.0, 665.0], [1150.0, 1160.0], [700.0, 580.0], [685.0, 595.0],
[685.0, 610.0], [770.0, 610.0], [795.0, 645.0], [720.0, 635.0], [760.0, 650.0],
[475.0, 960.0], [95.0, 260.0], [875.0, 920.0], [700.0, 500.0], [555.0, 815.0],
[830.0, 485.0], [1170.0, 65.0], [830.0, 610.0], [605.0, 625.0], [595.0, 360.0],
[1340.0, 725.0], [1740.0, 245.0]])
# fileName = "../data/eil76.dat" # 数据文件的地址和文件名
# coordinates = read_TSPLib(fileName) # 调用子程序，读取城市坐标数据文件
# 模拟退火算法参数设置
tInitial,tFinal,alfa,nMarkov = initParameter() # 调用子程序，获得设置参数
nCities = coordinates.shape[0] # 根据输入的城市坐标 获得城市数量 nCities
distMat = getDistMat(nCities, coordinates) # 调用子程序，计算城市间距离矩阵
nMarkov = nCities # Markov链 的初值设置
tNow = tInitial # 初始化 当前温度(current temperature)
# 初始化准备
tourNow = np.arange(nCities) # 产生初始路径，返回一个初值为0、步长为1、长度为n 的排列
valueNow = calTourMileage(tourNow,nCities,distMat) # 计算当前路径的总长度 valueNow
tourBest = tourNow.copy() # 初始化最优路径，复制 tourNow
valueBest = valueNow # 初始化最优路径的总长度，复制 valueNow
recordBest = [] # 初始化 最优路径记录表
recordNow = [] # 初始化 最优路径记录表
# 开始模拟退火优化过程
iter = 0 # 外循环迭代次数计数器
while tNow >= tFinal: # 外循环，直到当前温度达到终止温度时结束
# 在当前温度下，进行充分次数(nMarkov)的状态转移以达到热平衡
for k in range(nMarkov): # 内循环，循环次数为Markov链长度
# 产生新解：
tourNew = mutateSwap(tourNow, nCities) # 通过 交换操作 产生新路径
# tourNew,deltaE = mutateSwapE(tourNow,nCities,distMat) # 通过 交换操作 产生新路径（计算 deltaE）
valueNew = calTourMileage(tourNew,nCities,distMat) # 计算当前路径的总长度
deltaE = valueNew - valueNow
# 接受判别：按照 Metropolis 准则决定是否接受新解
if deltaE < 0: # 更优解：如果新解的目标函数好于当前解，则接受新解
accept = True
if valueNew < valueBest: # 如果新解的目标函数好于最优解，则将新解保存为最优解
tourBest[:] = tourNew[:]
valueBest = valueNew
else: # 容忍解：如果新解的目标函数比当前解差，则以一定概率接受新解
pAccept = math.exp(-deltaE/tNow) # 计算容忍解的状态迁移概率
if pAccept > random.random():
accept = True
else:
accept = False
# 保存新解
if accept == True: # 如果接受新解，则将新解保存为当前解
tourNow[:] = tourNew[:]
valueNow = valueNew
# 平移当前路径，以解决变换操作避开 0,（n-1）所带来的问题，并未实质性改变当前路径。
tourNow = np.roll(tourNow,2) # 循环移位函数，沿指定轴滚动数组元素
# 完成当前温度的搜索，保存数据和输出
recordBest.append(valueBest) # 将本次温度下的最优路径长度追加到 最优路径记录表
recordNow.append(valueNow) # 将当前路径长度追加到 当前路径记录表
print('i:{}, t(i):{:.2f}, valueNow:{:.1f}, valueBest:{:.1f}'.format(iter,tNow,valueNow,valueBest))
# 缓慢降温至新的温度，
iter = iter + 1
tNow = tNow * alfa # 指数降温曲线：T(k)=alfa*T(k-1)
# 结束模拟退火过程
# 图形化显示优化结果
figure1 = plt.figure() # 创建图形窗口 1
plot_tour(tourBest, valueBest, coordinates)
figure2 = plt.figure() # 创建图形窗口 2
plt.title("Optimization result of TSP{:d}".format(nCities)) # 设置图形标题
plt.plot(np.array(recordBest),'b-', label='Best') # 绘制 recordBest曲线
plt.plot(np.array(recordNow),'g-', label='Now') # 绘制 recordNow曲线
plt.xlabel("iter") # 设置 x轴标注
plt.ylabel("mileage of tour") # 设置 y轴标注
plt.legend() # 显示图例
plt.show()
print("Tour verification successful!")
print("Best tour: \n", tourBest)
print("Best value: {:.1f}".format(valueBest))
exit()
if __name__ == '__main__':
main()
运行结果:
Tour verification successful!
Best tour:
[18 7 40 2 16 17 31 38 39 36 24 27 26 11 50 3 5 4 23 47 37 14 42 9
8 32 10 51 13 12 25 46 28 29 1 6 41 20 30 21 0 48 35 34 33 43 45 15
49 19 22 44]
Best value: 9544.0
https://www.cnblogs.com/youcans/p/14728931.html
TSP"旅行商问题"应用领域包括:如何规划最合理高效的道路交通,以减少拥堵;如何更好地规划物流,以减少运输成本;在互联网环境中如何更好地设置节点,以更好地让信息流动
# 种群竞争模型
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint # 用以求常微分
plt.rcParams['figure.dpi'] = 100 # 绘图的dpi
plt.rcParams['font.sans-serif'] = ['SimHei'] # 正常显示中文
plt.rcParams['axes.unicode_minus'] = False # 正常显示正负号
# 两个物种的参数
type_x1 = [500, 1.2, 15000, 0.0004] # [初始数量, 自然增长率, 环境容量, 资源消耗]
type_x2 = [2560, 1.3, 3000, 0.0001] # 上面资源消耗的意思是：对于可以供养x2的资源，单位数量的x1的消耗为单位数量x2消耗的倍数
# 阻滞作用
def propagate(init, time, x1, x2):
ix1, ix2 = init
rx1 = x1[1]*ix1*(1-ix1/x1[2])-x1[3]*ix1*ix2
rx2 = x2[1]*ix2*(1-ix2/x2[2])-x2[3]*ix1*ix2
rx = np.array([rx1, rx2])
return rx
# 画图
def ploter(time, numer):
plt.xlabel('时间')
plt.ylabel('物种量')
plt.plot(time, numer[:,0], "b-", label="物种$x_1$")
plt.plot(time, numer[:,1], "r-", label="物种$x_2$")
plt.legend()
plt.show()
# 运行
time = np.linspace(0, 200, 1000) # 时间为200个单位，均分为1000份
init = np.array([type_x1[0], type_x2[0]])
numer = odeint(propagate, y0=init, t=time, args=(type_x1, type_x2))
ploter(time, numer)
# 这个是绘制两个物种之间的关系
plt.subplot(1, 2, 2)
plt.plot(number[:, 0], number[:, 1], 'k', linewidth=2)
plt.xlabel('x1')
plt.ylabel('x2')
plt.grid()
https://blog.csdn.net/qq_39798423/article/details/89400206
https://blog.csdn.net/qq_25601345/article/details/107757565
应用:不同企业推出的类似产品的竞争问题
很多个指标最终变成几个指标
示例
# 数据处理
import pandas as pd
import numpy as np
# 绘图
import seaborn as sns
import matplotlib.pyplot as plt
df = pd.read_csv(r"D:\桌面\aa.csv", encoding='gbk', index_col=0).reset_index(drop=True)
print(df)
# Bartlett's球状检验
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity
chi_square_value, p_value = calculate_bartlett_sphericity(df)
print(chi_square_value, p_value)
# KMO检验
# 检查变量间的相关性和偏相关性，取值在0-1之间；KOM统计量越接近1，变量间的相关性越强，偏相关性越弱，因子分析的效果越好。
# 通常取值从0.6开始进行因子分析
from factor_analyzer.factor_analyzer import calculate_kmo
kmo_all, kmo_model = calculate_kmo(df)
print(kmo_all)
# #标准化
# #所需库
# from sklearn import preprocessing
# #进行标准化
# df = preprocessing.scale(df)
# print(df)
# #求解系数相关矩阵
# covX = np.around(np.corrcoef(df.T),decimals=3)
# print(covX)
# #求解特征值和特征向量
# featValue, featVec= np.linalg.eig(covX.T) #求解系数相关矩阵的特征值和特征向量
# print(featValue, featVec)
#不标准化
#均值
def meanX(dataX):
return np.mean(dataX,axis=0)#axis=0表示依照列来求均值。假设输入list,则axis=1
average = meanX(df)
print(average)
#查看列数和行数
m, n = np.shape(df)
print(m,n)
#均值矩阵
data_adjust = []
avgs = np.tile(average, (m, 1))
print(avgs)
#去中心化
data_adjust = df - avgs
print(data_adjust)
#协方差阵
covX = np.cov(data_adjust.T) #计算协方差矩阵
print(covX)
#计算协方差阵的特征值和特征向量
featValue, featVec= np.linalg.eig(covX) #求解协方差矩阵的特征值和特征向量
print(featValue, featVec)
####下面没有区分标准化还是没有标准化#######
#对特征值进行排序并输出 降序
featValue = sorted(featValue)[::-1]
print(featValue)
#绘制散点图和折线图
# 同样的数据绘制散点图和折线图
plt.scatter(range(1, df.shape[1] + 1), featValue)
plt.plot(range(1, df.shape[1] + 1), featValue)
# 显示图的标题和xy轴的名字
# 最好使用英文，中文可能乱码
plt.title("Scree Plot")
plt.xlabel("Factors")
plt.ylabel("Eigenvalue")
plt.grid() # 显示网格
plt.show() # 显示图形
#求特征值的贡献度
gx = featValue/np.sum(featValue)
print(gx)
#求特征值的累计贡献度
lg = np.cumsum(gx)
print(lg)
#选出主成分
k=[i for i in range(len(lg)) if lg[i]<0.85]
k = list(k)
print(k)
#选出主成分对应的特征向量矩阵
selectVec = np.matrix(featVec.T[k]).T
selectVe=selectVec*(-1)
print(selectVec)
#主成分得分
finalData = np.dot(data_adjust,selectVec)
print(finalData)
#绘制热力图
plt.figure(figsize = (14,14))
ax = sns.heatmap(selectVec, annot=True, cmap="BuPu")
# 设置y轴字体大小
ax.yaxis.set_tick_params(labelsize=15)
plt.title("Factor Analysis", fontsize="xx-large")
# 设置y轴标签
plt.ylabel("Sepal Width", fontsize="xx-large")
# 显示图片
plt.show()
# 保存图片
# plt.savefig("factorAnalysis", dpi=500)
https://blog.csdn.net/qq_25990967/article/details/121366143
https://blog.csdn.net/qq_45722196/article/details/127584340
K-均值 示例
# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans, MiniBatchKMeans
def main():
# 读取数据文件
readPath = "../data/education2015.xlsx" # 数据文件的地址和文件名
dfFile = pd.read_excel(readPath, header=0) # 首行为标题行
dfFile = dfFile.dropna() # 删除含有缺失值的数据
# print(dfFile.dtypes) # 查看 df 各列的数据类型
# print(dfFile.shape) # 查看 df 的行数和列数
print(dfFile.head())
# 数据准备
z_scaler = lambda x:(x-np.mean(x))/np.std(x) # 定义数据标准化函数
dfScaler = dfFile[['x1','x2','x3','x4','x5','x6','x7','x8','x9','x10']].apply(z_scaler) # 数据归一化
dfData = pd.concat([dfFile[['地区']], dfScaler], axis=1) # 列级别合并
df = dfData.loc[:,['x1','x2','x3','x4','x5','x6','x7','x8','x9','x10']] # 基于全部 10个特征聚类分析
# df = dfData.loc[:,['x1','x2','x7','x8','x9','x10']] # 降维后选取 6个特征聚类分析
X = np.array(df) # 准备 sklearn.cluster.KMeans 模型数据
print("Shape of cluster data:", X.shape)
# KMeans 聚类分析(sklearn.cluster.KMeans)
nCluster = 4
kmCluster = KMeans(n_clusters=nCluster).fit(X) # 建立模型并进行聚类，设定 K=2
print("Cluster centers:\n", kmCluster.cluster_centers_) # 返回每个聚类中心的坐标
print("Cluster results:\n", kmCluster.labels_) # 返回样本集的分类结果
# 整理聚类结果
listName = dfData['地区'].tolist() # 将 dfData 的首列 '地区' 转换为 listName
dictCluster = dict(zip(listName,kmCluster.labels_)) # 将 listName 与聚类结果关联，组成字典
listCluster = [[] for k in range(nCluster)]
for v in range(0, len(dictCluster)):
k = list(dictCluster.values())[v] # 第v个城市的分类是 k
listCluster[k].append(list(dictCluster.keys())[v]) # 将第v个城市添加到 第k类
print("\n聚类分析结果(分为{}类):".format(nCluster)) # 返回样本集的分类结果
for k in range(nCluster):
print("第 {} 类：{}".format(k, listCluster[k])) # 显示第 k 类的结果
return
if __name__ == '__main__':
main()
https://blog.csdn.net/youcans/article/details/116596017
K-Means算法 示例
from sklearn.cluster import KMeans # 导入 sklearn.cluster.KMeans 类
import numpy as np
X = np.array([[1,2], [1,4], [1,0], [10,2], [10,4], [10,0]])
kmCluster = KMeans(n_clusters=2).fit(X) # 建立模型并进行聚类，设定 K=2
print(kmCluster.cluster_centers_) # 返回每个聚类中心的坐标
#[[10., 2.], [ 1., 2.]] # print 显示聚类中心坐标
print(kmCluster.labels_) # 返回样本集的分类结果
#[1, 1, 1, 0, 0, 0] # print 显示分类结果
print(kmCluster.predict([[0, 0], [12, 3]])) # 根据模型聚类结果进行预测判断
#[1, 0] # print显示判断结果：样本属于哪个类别
https://blog.csdn.net/youcans/article/details/116596017
自变量很多，变量之间不完全相互独立。逐步回归，找出对因变量影响大的因素
示例
import numpy as np
import pandas as pd
import statsmodels.api as sm
file = r'C:\Users\data.xlsx'
data = pd.read_excel(file)
data.columns = ['y', 'x1', 'x2', 'x3', 'x4', 'x5', 'x6', 'x7', 'x8', 'x9']
def looper(limit):
cols = ['x1', 'x2', 'x3', 'x5', 'x6', 'x7', 'x8', 'x9']
for i in range(len(cols)):
data1 = data[cols]
x = sm.add_constant(data1) #生成自变量
y = data['y'] #生成因变量
model = sm.OLS(y, x) #生成模型
result = model.fit() #模型拟合
pvalues = result.pvalues #得到结果中所有P值
pvalues.drop('const',inplace=True) #把const取得
pmax = max(pvalues) #选出最大的P值
if pmax>limit:
ind = pvalues.idxmax() #找出最大P值的index
cols.remove(ind) #把这个index从cols中删除
else:
return result
result = looper(0.05)
result.summary()
https://blog.csdn.net/bf02jgtrs00xktcx/article/details/108231365（值得一看）
