數(shù)學(xué)建模模型
綜合評價，權(quán)重決策分析
步驟：
構(gòu)建層次結(jié)構(gòu)模型
構(gòu)造判斷矩陣/成對比較矩陣
層次單排序及其一致性檢驗(yàn)
層次總排序及其一致性檢驗(yàn)
缺點(diǎn)：主觀性比較強(qiáng)
# 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? #導(dǎo)入所需包并將其命名為np
def ConsisTest(X):? #函數(shù)接收一個如上述A似的矩陣
#計(jì)算權(quán)重
? #方法一：算術(shù)平均法
? ? ## 第一步：將判斷矩陣按照列歸一化（每個元素除以其所在列的和）
? ? X = np.array(X)? #將X轉(zhuǎn)換為np.array對象
? ? sum_X = X.sum(axis=0)? #計(jì)算X每列的和
? ? (n,n) = X.shape? #X為方陣，行和列相同，所以用一個n來接收
? ? sum_X = np.tile(sum_X,(n,1))? #將和向量重復(fù)n行組成新的矩陣
? ? stand_X = X/sum_X? #標(biāo)準(zhǔn)化X（X中每個元素除以其所在列的和）
? ??
? ? ## 第二步：將歸一化矩陣每一行求和
? ? sum_row = stand_X.sum(axis=1)
? ? ## 第三步：將相加后得到的向量中每個元素除以n即可得到權(quán)重向量
? ? print("算數(shù)平均法求權(quán)重的結(jié)果為：")
? ? print(sum_row/n)
? ??
? #方法二：特征值法
? ? ## 第一步：找出矩陣X的最大特征值以及其對應(yīng)的特征向量
? ? V,E = np.linalg.eig(X)? #V是特征值，E是特征值對應(yīng)的特征向量
? ? max_value = np.max(V)? #最大特征值
? ? #print("最大特征值是：",max_value)
? ? max_v_index = np.argmax(V)? #返回最大特征值所在位置
? ? max_eiv = E[:,max_v_index]? #最大特征值對應(yīng)的特征向量
? ??
? ? ## 第二步：對求出的特征向量進(jìn)行歸一化處理即可得到權(quán)重
? ? stand_eiv = max_eiv/max_eiv.sum()
? ? print("特征值法求權(quán)重的結(jié)果為：")
? ? print(stand_eiv)
? ? print("———————————————————————————————")
#一致性檢驗(yàn)
? ? ## 第一步：計(jì)算一致性指標(biāo)CI
? ? CI = (max_value-n)/(n-1)
? ? ## 第二步：查找對應(yīng)的平均隨機(jī)一致性指標(biāo)RI
? ? RI = np.array([0, 0, 0.58, 0.9, 1.12, 1.24, 1.32, 1.41, 1.45, 1.49])
? ? ## 第三步：計(jì)算一致性比例CR
? ? CR = CI/RI[n]
? ? if CR < 0.1:
? ? ? ? print("CR=",CR,"，小于0.1，通過一致性檢驗(yàn)")
? ? else:
? ? ? ? print("CR=",CR,"，大于等于0.1，沒有通過一致性檢驗(yàn)，請修改判斷矩陣")
? ? return None
#ConsisTest(A)
https://zhuanlan.zhihu.com/p/613113565
https://blog.csdn.net/ddjhpxs/article/details/105407103
https://zhuanlan.zhihu.com/p/266405027
一般是輸入時間序列,適合數(shù)據(jù)量少的時候
示例代碼
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
plt.rcParams['font.sans-serif'] = ['SimHei']? # 用來正常顯示中文標(biāo)簽
plt.rcParams['axes.unicode_minus'] = False? # 用來正常顯示負(fù)號
?
data = pd.read_excel("E:/桌面/灰色預(yù)測.xlsx",sheet_name=1)
list1 = data["序列"]
tlist1 = [2014,2015,2016,2017,2018]
fy = []
def GM11(x,n):
? ? '''
? ? 灰色預(yù)測
? ? x：序列，numpy對象
? ? n:需要往后預(yù)測的個數(shù)
? ? '''
? ? 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為發(fā)展系數(shù) b為灰色作用量
? ? [[a],[b]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)#計(jì)算待估參數(shù)??
? ? result = (x[0]-b/a)*np.exp(-a*(n-1))-(x[0]-b/a)*np.exp(-a*(n-2))? #預(yù)測方程
? ? 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)? ? #預(yù)測值
? ? ? ? fy.append(predict)
? ? print("后驗(yàn)差檢驗(yàn)")
? ? S2_2 = np.array(e).var()#殘差方差
? ? C = S2_2/S1_2#后驗(yàn)差比
? ? if C<=0.35:
? ? ? ? assess = '后驗(yàn)差比<=0.35，模型精度等級為好'
? ? elif C<=0.5:
? ? ? ? assess = '后驗(yàn)差比<=0.5，模型精度等級為合格'
? ? elif C<=0.65:
? ? ? ? assess = '后驗(yàn)差比<=0.65，模型精度等級為勉強(qiáng)'
? ? else:
? ? ? ? assess = '后驗(yàn)差比>0.65，模型精度等級為不合格'
? ? #預(yù)測數(shù)據(jù)
? ? 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':'發(fā)展灰數(shù)'},
? ? ? ? ? ? 'b':{'value':b,'desc':'控制灰數(shù)'},
? ? ? ? ? ? 'predict':{'value':result,'desc':'第%d個預(yù)測值'%n},
? ? ? ? ? ? 'C':{'value':C,'desc':assess},
? ? ? ? ? ? 'predict':{'value':predict,'desc':'往后預(yù)測%d個的序列'%(n)},
? ? ? ? ? ? }
? ??
?
? ??
if __name__ == "__main__":
? ? data = np.array(list1)
? ? x = data[0:5]#輸入數(shù)據(jù)
? ? #y = data[0:]#需要預(yù)測的數(shù)據(jù)
? ? result = GM11(x,1)
? ? predict = result['predict']['value']
? ? predict = np.round(predict,1)
? ? #print('真實(shí)值:',x)
? ? print('預(yù)測值:',predict)
? ? print(result)
?
print("殘差檢驗(yàn)")
a = []
for i in range(5):
? ? a.append(abs(fy[i]-list1[i]))
? ? print('%.5f%%' % ((abs(fy[i]-list1[i]))/list1[i]))
print("關(guān)聯(lián)度檢驗(yàn)")
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 關(guān)聯(lián)度為：",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='真實(shí)值曲線')? ?#真實(shí)值
plt.plot(x2,y2,'b+-',label='預(yù)測值曲線')? ?#預(yù)測值
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)? # 鄰接矩陣維度，即節(jié)點(diǎn)個數(shù)
? ? fmax = float('inf')
? ? w = [[0 for _ in range(n)] for j in range(n)]? # 鄰接矩陣轉(zhuǎn)化成維度矩陣，即0→max
? ? book = [0 for _ in range(n)]? # 是否已經(jīng)是最小的標(biāo)記列表
? ? dis = [fmax for i in range(n)]? # s到其他節(jié)點(diǎn)的最小距離
? ? book[s - 1] = 1? # 節(jié)點(diǎn)編號從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:? # 直連的節(jié)點(diǎn)最小距離就是network[i][j]
? ? ? ? ? dis[j] = network[i][j]
? ? for i in range(n - 1):? # n-1次遍歷，除了s節(jié)點(diǎn)
? ? ? 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直連的節(jié)點(diǎn)遍歷一遍
? ? ? ? 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)? # 因?yàn)榇鎯Φ氖巧弦惶?，所以先加入目的?jié)點(diǎn)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)
#運(yùn)行結(jié)果
#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
# 創(chuàng)建帶權(quán)圖
G = nx.Graph()
G.add_weighted_edges_from([(1, 2, 1), (1, 3, 2), (2, 3, 1), (2, 4, 3), (3, 4, 1)])
# 計(jì)算最短路徑
path = nx.dijkstra_path(G, 1, 4)
distance = nx.dijkstra_path_length(G, 1, 4)
# 輸出結(jié)果
print('Shortest path:', path)
print('Shortest distance:', distance)
#結(jié)果如下:
#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]? # 更新父結(jié)點(diǎn)
# 打印路徑
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
# 構(gòu)圖
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)]? # 關(guān)鍵地方，i-->j 的父結(jié)點(diǎn)初始化都為i
for u, v, c in datas:
? ? graph[u][v] = c # 因?yàn)槭怯邢驁D，邊權(quán)只賦給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
# 定義弗洛伊德算法函數(shù)
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
# 定義帶權(quán)圖
graph = np.array([
? ? [0, 1, 2, np.inf],
? ? [1, 0, 1, 3],
? ? [2, 1, 0, 1],
? ? [np.inf, 3, 1, 0]
])
# 計(jì)算最短路徑
dist = floyd(graph)
# 輸出結(jié)果
print('Shortest distance matrix:')
print(dist)
#得到的是所有節(jié)點(diǎn)之間的最短路徑矩陣
Shortest distance matrix:
[[0. 1. 2. 3.]
?[1. 0. 1. 2.]
?[2. 1. 0. 1.]
?[3. 2. 1. 0.]]
https://pythonjishu.com/znmehrtcqshdjef/
弗洛伊德和迪杰斯特拉算法的區(qū)別:
https://blog.csdn.net/weixin_43872728/article/details/100662957
import math
from random import random
import matplotlib.pyplot as plt
def func(x, y):? ? ? ? ? ? ? ? ? #函數(shù)優(yōu)化問題
? ? 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? ? ? ? ?#內(nèi)循環(huán)迭代次數(shù),即為L =100
? ? ? ? self.alpha = alpha? ? ? ?#降溫系數(shù)，alpha=0.99
? ? ? ? self.T0 = T0? ? ? ? ? ? ?#初始溫度T0為100
? ? ? ? self.Tf = Tf? ? ? ? ? ? ?#溫度終值Tf為0.01
? ? ? ? self.T = T0? ? ? ? ? ? ? #當(dāng)前溫度
? ? ? ? self.x = [random() * 11 -5? for i in range(iter)] #隨機(jī)生成100個x的值
? ? ? ? self.y = [random() * 11 -5? for i in range(iter)] #隨機(jī)生成100個y的值
? ? ? ? self.most_best =[]
? ? ? ? """
? ? ? ? random()這個函數(shù)取0到1之間的小數(shù)
? ? ? ? 如果你要取0-10之間的整數(shù)（包括0和10）就寫成 (int)random()*11就可以了，11乘以零點(diǎn)多的數(shù)最大是10點(diǎn)多，最小是0點(diǎn)多
? ? ? ? 該實(shí)例中x1和x2的絕對值不超過5（包含整數(shù)5和-5），（random() * 11 -5）的結(jié)果是-6到6之間的任意值（不包括-6和6）
? ? ? ? （random() * 10 -5）的結(jié)果是-5到5之間的任意值（不包括-5和5），所有先乘以11，取-6到6之間的值，產(chǎn)生新解過程中，用一個if條件語句把-5到5之間（包括整數(shù)5和-5）的篩選出來。
? ? ? ? """
? ? ? ? self.history = {'f': [], 'T': []}
? ? def generate_new(self, x, y):? ?#擾動產(chǎn)生新解的過程
? ? ? ? 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? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? #重復(fù)得到新解，直到產(chǎn)生的新解滿足約束條件
? ? ? ? return x_new, y_new?
? ? def Metrospolis(self, f, f_new):? ?#Metropolis準(zhǔn)則
? ? ? ? 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):? ? #獲取最優(yōu)目標(biāo)函數(shù)值
? ? ? ? f_list = []? ? #f_list數(shù)組保存每次迭代之后的值
? ? ? ? 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次之后目標(biāo)函數(shù)的最優(yōu)解和最優(yōu)解的下標(biāo)
? ? def run(self):
? ? ? ? count = 0
? ? ? ? #外循環(huán)迭代，當(dāng)前溫度小于終止溫度的閾值
? ? ? ? while self.T > self.Tf:? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ? #內(nèi)循環(huán)迭代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]) #產(chǎn)生新解
? ? ? ? ? ? ? ? f_new = self.func(x_new, y_new)? ? ? ? ? ? ? ? ? ? ? ? #產(chǎn)生新值
? ? ? ? ? ? ? ? if self.Metrospolis(f, f_new):? ? ? ? ? ? ? ? ? ? ? ? ?#判斷是否接受新值
? ? ? ? ? ? ? ? ? ? self.x[i] = x_new? ? ? ? ? ? ?#如果接受新值，則把新值的x,y存入x數(shù)組和y數(shù)組
? ? ? ? ? ? ? ? ? ? self.y[i] = y_new
? ? ? ? ? ? # 迭代L次記錄在該溫度下最優(yōu)解
? ? ? ? ? ? ft, _ = self.best()
? ? ? ? ? ? self.history['f'].append(ft)
? ? ? ? ? ? self.history['T'].append(self.T)
? ? ? ? ? ? #溫度按照一定的比例下降（冷卻）
? ? ? ? ? ? self.T = self.T * self.alpha? ? ? ??
? ? ? ? ? ? count += 1
? ? ? ? ? ??
? ? ? ? ? ? # 得到最優(yōu)解
? ? ? ? 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? ? ? ? ? ? ? ? ? ? ? ? ?# 導(dǎo)入模塊 math
import random? ? ? ? ? ? ? ? ? ? ? ?# 導(dǎo)入模塊 random
import pandas as pd? ? ? ? ? ? ? ? ?# 導(dǎo)入模塊 pandas 并簡寫成 pd
import numpy as np? ? ? ? ? ? ? ? ? # 導(dǎo)入模塊 numpy 并簡寫成 np YouCans
import matplotlib.pyplot as plt? ? ?# 導(dǎo)入模塊 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
# 子程序：初始化模擬退火算法的控制參數(shù)
def initParameter():
? ? # custom function initParameter():
? ? # Initial parameter for simulated annealing algorithm
? ? tInitial = 100.0? ? ? ? ? ? ? ? # 設(shè)定初始退火溫度(initial temperature)
? ? tFinal? = 1? ? ? ? ? ? ? ? ? ? ?# 設(shè)定終止退火溫度(stop temperature)
? ? nMarkov = 1000? ? ? ? ? ? ? ? # Markov鏈長度，也即內(nèi)循環(huán)運(yùn)行次數(shù)
? ? alfa? ? = 0.98? ? ? ? ? ? ? ? ?# 設(shè)定降溫參數(shù)，T(k)=alfa*T(k-1)
? ? return tInitial,tFinal,alfa,nMarkov
# 子程序：讀取TSPLib數(shù)據(jù)
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')? ? ? # 數(shù)據(jù)格式：i Xi Yi
? ? coordinates = loadData[:,1::]
? ? return coordinates
# 子程序：計(jì)算各城市間的距離，得到距離矩陣
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 求向量的范數(shù)（默認(rèn)求 二范數(shù)），得到 i、j 間的距離
? ? ? ? ? ? distMat[i][j] = distMat[j][i] = round(np.linalg.norm(coordinates[i]-coordinates[j]))
? ? return distMat? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 城市間距離取整（四舍五入）
# 子程序：計(jì)算 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')? ? ? # 繪制城市坐標(biāo)點(diǎn) 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')? # 繪制城市坐標(biāo)點(diǎn) 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))? # 設(shè)置圖形標(biāo)題
? ? 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) 產(chǎn)生 2個不相等的隨機(jī)整數(shù) i,j
? ? i = np.random.randint(nCities)? ? ? ? ? # 產(chǎn)生第一個 [0,n) 區(qū)間內(nèi)的隨機(jī)整數(shù)
? ? while True:
? ? ? ? j = np.random.randint(nCities)? ? ? # 產(chǎn)生一個 [0,n) 區(qū)間內(nèi)的隨機(jī)整數(shù)
? ? ? ? if i!=j: break? ? ? ? ? ? ? ? ? ? ? # 保證 i, j 不相等
? ? tourSwap = tourGiven.copy()? ? ? ? ? ? ?# 將給定路徑復(fù)制給新路徑 tourSwap
? ? tourSwap[i],tourSwap[j] = tourGiven[j],tourGiven[i] # 交換 城市 i 和 j 的位置————簡潔的實(shí)現(xiàn)方法
? ? return tourSwap
def main():
? ? # 主程序
? ? # # 讀取旅行城市位置的坐標(biāo)
? ? 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"? ? ? ? ? ? ? ? ? ? ? # 數(shù)據(jù)文件的地址和文件名
? ? # coordinates = read_TSPLib(fileName)? ? ? ? ? ? ? ? ?# 調(diào)用子程序，讀取城市坐標(biāo)數(shù)據(jù)文件
? ? # 模擬退火算法參數(shù)設(shè)置
? ? tInitial,tFinal,alfa,nMarkov = initParameter()? ? ? # 調(diào)用子程序，獲得設(shè)置參數(shù)
? ? nCities = coordinates.shape[0]? ? ? ? ? ? ? # 根據(jù)輸入的城市坐標(biāo) 獲得城市數(shù)量 nCities
? ? distMat = getDistMat(nCities, coordinates)? # 調(diào)用子程序，計(jì)算城市間距離矩陣
? ? nMarkov = nCities? ? ? ? ? ? ? ? ? ? ? ? ? ?# Markov鏈 的初值設(shè)置
? ? tNow? ? = tInitial? ? ? ? ? ? ? ? ? ? ? ? ? # 初始化 當(dāng)前溫度(current temperature)
? ? # 初始化準(zhǔn)備
? ? tourNow? ?= np.arange(nCities)? ?# 產(chǎn)生初始路徑，返回一個初值為0、步長為1、長度為n 的排列
? ? valueNow? = calTourMileage(tourNow,nCities,distMat) # 計(jì)算當(dāng)前路徑的總長度 valueNow
? ? tourBest? = tourNow.copy()? ? ? ? ? ? ? ? ? ? ? ? ? # 初始化最優(yōu)路徑，復(fù)制 tourNow
? ? valueBest = valueNow? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 初始化最優(yōu)路徑的總長度，復(fù)制 valueNow
? ? recordBest = []? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?# 初始化 最優(yōu)路徑記錄表
? ? recordNow? = []? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?# 初始化 最優(yōu)路徑記錄表
? ? # 開始模擬退火優(yōu)化過程
? ? iter = 0? ? ? ? ? ? ? ? ? ? ? ? # 外循環(huán)迭代次數(shù)計(jì)數(shù)器
? ? while tNow >= tFinal:? ? ? ? ? ?# 外循環(huán)，直到當(dāng)前溫度達(dá)到終止溫度時結(jié)束
? ? ? ? # 在當(dāng)前溫度下，進(jìn)行充分次數(shù)(nMarkov)的狀態(tài)轉(zhuǎn)移以達(dá)到熱平衡
? ? ? ? for k in range(nMarkov):? ? # 內(nèi)循環(huán)，循環(huán)次數(shù)為Markov鏈長度
? ? ? ? ? ? # 產(chǎn)生新解：
? ? ? ? ? ? tourNew = mutateSwap(tourNow, nCities)? ? ? # 通過 交換操作 產(chǎn)生新路徑
? ? ? ? ? ? # tourNew,deltaE = mutateSwapE(tourNow,nCities,distMat)? ?# 通過 交換操作 產(chǎn)生新路徑（計(jì)算 deltaE）
? ? ? ? ? ? valueNew = calTourMileage(tourNew,nCities,distMat) # 計(jì)算當(dāng)前路徑的總長度
? ? ? ? ? ? deltaE = valueNew - valueNow
? ? ? ? ? ? # 接受判別：按照 Metropolis 準(zhǔn)則決定是否接受新解
? ? ? ? ? ? if deltaE < 0:? ? ? ? ? ? ? ? ? ? ? ? ? # 更優(yōu)解：如果新解的目標(biāo)函數(shù)好于當(dāng)前解，則接受新解
? ? ? ? ? ? ? ? accept = True
? ? ? ? ? ? ? ? if valueNew < valueBest:? ? ? ? ? ? # 如果新解的目標(biāo)函數(shù)好于最優(yōu)解，則將新解保存為最優(yōu)解
? ? ? ? ? ? ? ? ? ? tourBest[:] = tourNew[:]
? ? ? ? ? ? ? ? ? ? valueBest = valueNew
? ? ? ? ? ? else:? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?# 容忍解：如果新解的目標(biāo)函數(shù)比當(dāng)前解差，則以一定概率接受新解
? ? ? ? ? ? ? ? pAccept = math.exp(-deltaE/tNow)? ? # 計(jì)算容忍解的狀態(tài)遷移概率
? ? ? ? ? ? ? ? if pAccept > random.random():
? ? ? ? ? ? ? ? ? ? accept = True
? ? ? ? ? ? ? ? else:
? ? ? ? ? ? ? ? ? ? accept = False
? ? ? ? ? ? # 保存新解
? ? ? ? ? ? if accept == True:? ? ? ? ? ? ? ? ? ? ? # 如果接受新解，則將新解保存為當(dāng)前解
? ? ? ? ? ? ? ? tourNow[:] = tourNew[:]
? ? ? ? ? ? ? ? valueNow = valueNew
? ? ? ? # 平移當(dāng)前路徑，以解決變換操作避開 0,（n-1）所帶來的問題，并未實(shí)質(zhì)性改變當(dāng)前路徑。
? ? ? ? tourNow = np.roll(tourNow,2)? ? ? ? ? ? ? ? # 循環(huán)移位函數(shù)，沿指定軸滾動數(shù)組元素
? ? ? ? # 完成當(dāng)前溫度的搜索，保存數(shù)據(jù)和輸出
? ? ? ? recordBest.append(valueBest)? ? ? ? ? ? ? ? # 將本次溫度下的最優(yōu)路徑長度追加到 最優(yōu)路徑記錄表
? ? ? ? recordNow.append(valueNow)? ? ? ? ? ? ? ? ? # 將當(dāng)前路徑長度追加到 當(dāng)前路徑記錄表
? ? ? ? print('i:{}, t(i):{:.2f}, valueNow:{:.1f}, valueBest:{:.1f}'.format(iter,tNow,valueNow,valueBest))
? ? ? ? # 緩慢降溫至新的溫度，
? ? ? ? iter = iter + 1
? ? ? ? tNow = tNow * alfa? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 指數(shù)降溫曲線：T(k)=alfa*T(k-1)
? ? # 結(jié)束模擬退火過程
? ? # 圖形化顯示優(yōu)化結(jié)果
? ? figure1 = plt.figure()? ? ?# 創(chuàng)建圖形窗口 1
? ? plot_tour(tourBest, valueBest, coordinates)
? ? figure2 = plt.figure()? ? ?# 創(chuàng)建圖形窗口 2
? ? plt.title("Optimization result of TSP{:d}".format(nCities)) # 設(shè)置圖形標(biāo)題
? ? plt.plot(np.array(recordBest),'b-', label='Best')? ? ? ? ? ?# 繪制 recordBest曲線
? ? plt.plot(np.array(recordNow),'g-', label='Now')? ? ? ? ? ? ?# 繪制 recordNow曲線
? ? plt.xlabel("iter")? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 設(shè)置 x軸標(biāo)注
? ? plt.ylabel("mileage of tour")? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?# 設(shè)置 y軸標(biāo)注
? ? plt.legend()? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? # 顯示圖例
? ? plt.show()
? ? print("Tour verification successful!")
? ? print("Best tour: \n", tourBest)
? ? print("Best value: {:.1f}".format(valueBest))
? ? exit()
if __name__ == '__main__':
? ? main()
運(yùn)行結(jié)果:
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"旅行商問題"應(yīng)用領(lǐng)域包括:如何規(guī)劃最合理高效的道路交通,以減少擁堵;如何更好地規(guī)劃物流,以減少運(yùn)輸成本;在互聯(lián)網(wǎng)環(huán)境中如何更好地設(shè)置節(jié)點(diǎn),以更好地讓信息流動
# 種群競爭模型
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? # 正常顯示正負(fù)號
# 兩個物種的參數(shù)
type_x1 = [500, 1.2, 15000, 0.0004]? # [初始數(shù)量, 自然增長率, 環(huán)境容量, 資源消耗]
type_x2 = [2560, 1.3, 3000, 0.0001]? # 上面資源消耗的意思是：對于可以供養(yǎng)x2的資源，單位數(shù)量的x1的消耗為單位數(shù)量x2消耗的倍數(shù)
# 阻滯作用
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()
# 運(yùn)行
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)
# 這個是繪制兩個物種之間的關(guān)系
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
應(yīng)用:不同企業(yè)推出的類似產(chǎn)品的競爭問題
很多個指標(biāo)最終變成幾個指標(biāo)
示例
# 數(shù)據(jù)處理
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球狀檢驗(yàn)
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檢驗(yàn)
# 檢查變量間的相關(guān)性和偏相關(guān)性，取值在0-1之間；KOM統(tǒng)計(jì)量越接近1，變量間的相關(guān)性越強(qiáng)，偏相關(guān)性越弱，因子分析的效果越好。
# 通常取值從0.6開始進(jìn)行因子分析
from factor_analyzer.factor_analyzer import calculate_kmo
?
kmo_all, kmo_model = calculate_kmo(df)
print(kmo_all)
?
# #標(biāo)準(zhǔn)化
?
# #所需庫
# from sklearn import preprocessing
# #進(jìn)行標(biāo)準(zhǔn)化
# df = preprocessing.scale(df)
# print(df)
?
# #求解系數(shù)相關(guān)矩陣
# covX = np.around(np.corrcoef(df.T),decimals=3)
# print(covX)
?
# #求解特征值和特征向量
# featValue, featVec=? np.linalg.eig(covX.T)? #求解系數(shù)相關(guān)矩陣的特征值和特征向量
# print(featValue, featVec)
?
?
#不標(biāo)準(zhǔn)化
#均值
def meanX(dataX):
? ? return np.mean(dataX,axis=0)#axis=0表示依照列來求均值。假設(shè)輸入list,則axis=1
average = meanX(df)
print(average)
?
#查看列數(shù)和行數(shù)
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)
?
#協(xié)方差陣
covX = np.cov(data_adjust.T)? ?#計(jì)算協(xié)方差矩陣
print(covX)
?
#計(jì)算協(xié)方差陣的特征值和特征向量
featValue, featVec=? np.linalg.eig(covX)? #求解協(xié)方差矩陣的特征值和特征向量
print(featValue, featVec)
?
####下面沒有區(qū)分標(biāo)準(zhǔn)化還是沒有標(biāo)準(zhǔn)化#######
?
#對特征值進(jìn)行排序并輸出 降序
featValue = sorted(featValue)[::-1]
print(featValue)
?
#繪制散點(diǎn)圖和折線圖
# 同樣的數(shù)據(jù)繪制散點(diǎn)圖和折線圖
plt.scatter(range(1, df.shape[1] + 1), featValue)
plt.plot(range(1, df.shape[1] + 1), featValue)
?
# 顯示圖的標(biāo)題和xy軸的名字
# 最好使用英文，中文可能亂碼
plt.title("Scree Plot")
plt.xlabel("Factors")
plt.ylabel("Eigenvalue")
?
plt.grid()? # 顯示網(wǎng)格
plt.show()? # 顯示圖形
?
#求特征值的貢獻(xiàn)度
gx = featValue/np.sum(featValue)
print(gx)
?
#求特征值的累計(jì)貢獻(xiàn)度
lg = np.cumsum(gx)
print(lg)
?
#選出主成分
k=[i for i in range(len(lg)) if lg[i]<0.85]
k = list(k)
print(k)
?
#選出主成分對應(yīng)的特征向量矩陣
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")
?
# 設(shè)置y軸字體大小
ax.yaxis.set_tick_params(labelsize=15)
plt.title("Factor Analysis", fontsize="xx-large")
?
# 設(shè)置y軸標(biāo)簽
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():
? ? # 讀取數(shù)據(jù)文件
? ? readPath = "../data/education2015.xlsx"? # 數(shù)據(jù)文件的地址和文件名
? ? dfFile = pd.read_excel(readPath, header=0)? # 首行為標(biāo)題行
? ? dfFile = dfFile.dropna()? # 刪除含有缺失值的數(shù)據(jù)
? ? # print(dfFile.dtypes)? # 查看 df 各列的數(shù)據(jù)類型
? ? # print(dfFile.shape)? # 查看 df 的行數(shù)和列數(shù)
? ? print(dfFile.head())
? ? # 數(shù)據(jù)準(zhǔn)備
? ? z_scaler = lambda x:(x-np.mean(x))/np.std(x)? # 定義數(shù)據(jù)標(biāo)準(zhǔn)化函數(shù)
? ? dfScaler = dfFile[['x1','x2','x3','x4','x5','x6','x7','x8','x9','x10']].apply(z_scaler)? # 數(shù)據(jù)歸一化
? ? dfData = pd.concat([dfFile[['地區(qū)']], 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)? # 準(zhǔn)備 sklearn.cluster.KMeans 模型數(shù)據(jù)
? ? print("Shape of cluster data:", X.shape)
? ? # KMeans 聚類分析(sklearn.cluster.KMeans)
? ? nCluster = 4
? ? kmCluster = KMeans(n_clusters=nCluster).fit(X)? # 建立模型并進(jìn)行聚類，設(shè)定 K=2
? ? print("Cluster centers:\n", kmCluster.cluster_centers_)? # 返回每個聚類中心的坐標(biāo)
? ? print("Cluster results:\n", kmCluster.labels_)? # 返回樣本集的分類結(jié)果
? ? # 整理聚類結(jié)果
? ? listName = dfData['地區(qū)'].tolist()? # 將 dfData 的首列 '地區(qū)' 轉(zhuǎn)換為 listName
? ? dictCluster = dict(zip(listName,kmCluster.labels_))? # 將 listName 與聚類結(jié)果關(guān)聯(lián)，組成字典
? ? 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聚類分析結(jié)果(分為{}類):".format(nCluster))? # 返回樣本集的分類結(jié)果
? ? for k in range(nCluster):
? ? ? ? print("第 {} 類：{}".format(k, listCluster[k]))? # 顯示第 k 類的結(jié)果
? ? return
if __name__ == '__main__':
? ? main()
https://blog.csdn.net/youcans/article/details/116596017
K-Means算法 示例
from sklearn.cluster import KMeans? # 導(dǎo)入 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)? # 建立模型并進(jìn)行聚類，設(shè)定 K=2
print(kmCluster.cluster_centers_)? # 返回每個聚類中心的坐標(biāo)
#[[10., 2.], [ 1., 2.]]? # print 顯示聚類中心坐標(biāo)
print(kmCluster.labels_)? # 返回樣本集的分類結(jié)果
#[1, 1, 1, 0, 0, 0]? # print 顯示分類結(jié)果
print(kmCluster.predict([[0, 0], [12, 3]]))? # 根據(jù)模型聚類結(jié)果進(jìn)行預(yù)測判斷
#[1, 0]? # print顯示判斷結(jié)果：樣本屬于哪個類別
https://blog.csdn.net/youcans/article/details/116596017
自變量很多，變量之間不完全相互獨(dú)立。逐步回歸，找出對因變量影響大的因素
示例
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 #得到結(jié)果中所有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（值得一看）
整體參考：小石老師數(shù)學(xué)建模14講（全）_數(shù)學(xué)建模比賽快速入門_13個常用數(shù)學(xué)模型
