1. 前言

在這篇文章中，我會把我所接觸的數學建模的知識的代碼分享給大家，有的是我自己寫的，更多的是網上借鑒并修改為可執行可用的代碼

需要說明的是，我不太了解其中的數學實現的具體方法和原理，我只分享源碼和作用條件以及使用的經驗說明（更詳細見注釋），網上有不少文章對某些方法講得非常透徹，因此我沒必要贅述

各位只需要忘記代碼或使用條件時從我這里 copy 就好了

2. 灰色預測模型

這是用一種在數學上對數據進行累加和累減以取巧性地削弱預測數據與原始數據的關聯性

最好用于短期預測，只能用于遞增序列

2.1. GM（1，1）

我也不懂為什么叫 GM(1,1)，但是用得最多的就是它，又稱一階灰色方程

clc,clear;close all y = [0.6 1.1 1.6] yuceshu = 3; % 本程序主要用來計算根據灰色理論建立的模型的預測值 % 應用的數學模型是 GM(1,1) % 原始數據的處理方法是一次累加 % 渡辺曜改進了這個代碼,通常GM預測第一個數是不夠優秀的 % y 是應該輸入的數組，yuceshu是你想預測的后幾個數據 我們一般就預測兩個 % 示例 huiseyuce([1,2,3],2) n=length(y); yy=ones(n,1); yy(1)=y(1); for i=2:nyy(i)=yy(i-1)+y(i); end B=ones(n-1,2); for i=1:(n-1)B(i,1)=-(yy(i)+yy(i+1))/2;B(i,2)=1; end BT=B'; for j=1:n-1YN(j)=y(j+1); end YN=YN'; A=inv(BT*B)*BT*YN; a=A(1); u=A(2); t=u/a; i=1:n+yuceshu; yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1); for j=n+yuceshu:-1:2ys(j)=yys(j)-yys(j-1); end x=1:n; % 把下面這里的2改成1 就可以看到GM擬合的第一個數據 xs=2:n+yuceshu; yn=ys(2:n+yuceshu); plot(x,y,'^r',xs,yn,'*-b'); det=0;sum1=0; sumpe=0; for i=1:nsumpe=sumpe+y(i); end pe=sumpe/n; for i=1:nsum1=sum1+(y(i)-pe).^2; end s1=sqrt(sum1/n); sumce=0; for i=2:nsumce=sumce+(y(i)-yn(i)); end ce=sumce/(n-1); sum2=0; for i=2:nsum2=sum2+(y(i)-yn(i)-ce).^2; end s2=sqrt(sum2/(n-1)); c=(s2)/(s1); disp(['后驗差比值為：',num2str(c)]); if c<0.35disp('系統預測精度好') else if c<0.5disp('系統預測精度合格')else if c<0.65disp('系統預測精度勉強')elsedisp('系統預測精度不合格')endend enddisp(['下個擬合值為 ',num2str(ys(n+1))]); for i=2:yuceshudisp(['再下個擬合值為',num2str(ys(n+i))]); end

2.2. 分數階 GM(1,1)

分數階灰色方程的預測效果優于原始灰色模型，我只找到了 python 版本的可運行代碼

首先，將 GM(1,1)打包成一個類

GM.py

##渡辺筆記 ##傳統的灰色模型 import numpy as np import xlrd import xlwt import datetime class Grey_model(object):def __init__(self,input_value):##初始化過程，先建立好變量空間，即累加序列，背景值序列，B和Y矩陣，擬合序列，預測值序列等。self.input_value=input_valueself.accumulation_value=np.zeros(len(input_value))self.background_value=np.zeros(len(input_value)-1)self.y_matrix_value=np.mat(np.zeros((len(input_value)-1,1)))self.b_matrix_value=np.mat(np.ones((len(input_value)-1,2)))self.accumulation()##計算累加序列def accumulation(self):for i in range(len(self.input_value)):self.accumulation_value[i]=np.sum(self.input_value[0:i+1])##計算Z值def background_values(self):for i in range(len(self.accumulation_value)-1):self.background_value[i]=(self.accumulation_value[i]+self.accumulation_value[i+1])/2##計算B矩陣def b_matrix(self):self.background_values()for i in range (self.b_matrix_value.shape[0]):self.b_matrix_value[i,0]=-self.background_value[i]##計算Y矩陣def y_matrix(self):for i in range(self.y_matrix_value.shape[0]):self.y_matrix_value[i]=self.accumulation_value[i+1]-self.accumulation_value[i]##計算參數矩陣U：U=(B^T*Y)^-1*B^T*Ydef u_matrix(self):self.y_matrix()self.b_matrix()self.u_matrix_values=(self.b_matrix_value.T*self.b_matrix_value)**-1 * self.b_matrix_value.T * self.y_matrix_value##下面把矩陣格式轉化為數組格式，再轉化為列表格式self.u_matrix_array=np.array(self.u_matrix_values.T)[0]self.u_matrix_value=self.u_matrix_array.tolist()##計算預測值def predict(self,number_of_forecast):self.u_matrix()self.predicted_accumulation_value=[]##使用了float(%.2f% 來調整輸出值的小數的個數self.predicted_value=[float(%.2f%(self.input_value[0]))]for i in range(len(self.input_value)+int(number_of_forecast)):self.predicted_accumulation_value.append(((self.accumulation_value[0]-self.u_matrix_value[1]/self.u_matrix_value[0])*np.exp(-self.u_matrix_value[0]*(i)))+self.u_matrix_value[1]/self.u_matrix_value[0])##上式中e^-at,由于python中從0開始算，故t減去1for i in range(len(self.predicted_accumulation_value)-1):self.predicted_value.append(float(%.2f%(self.predicted_accumulation_value[i+1]-self.predicted_accumulation_value[i])))##計算誤差def test_error(self):MAPE_list=[]for i in range(len(self.input_value)):MAPE_list.append(abs((self.predicted_value[i]-self.input_value[i])/self.input_value[i]))self.MAPE=str(100*(np.mean(MAPE_list[1:])))[:5]+ '%'## 打開xls，只能用這種格式了 work_book = xlrd.open_workbook('first_shuju.xls') ## 按sheet表名稱獲取sheet對象，名稱分大小寫 sheet_1 = work_book.sheet_by_name('Sheet1') col_0_value = sheet_1.col_values(6,2,15)##輸入變量，order表示分數階的階數 input_value = col_0_value if __name__ == '__main__':input_value = col_0_valuenumber_of_forecast=3A=Grey_model(input_value)A.predict(number_of_forecast)A.test_error()print('預測值')print(A.predicted_value)print('\nMAPE值')print(A.MAPE)## xlwt_data = A.predicted_valuedef save_excel(target_list, output_file_name):將數據寫入xls文件if not output_file_name.endswith('.xls'):output_file_name += '.xls'workbook = xlwt.Workbook(encoding='utf-8')ws = workbook.add_sheet(sheet1)rows = len(target_list)lines = len(target_list[0])for i in range(rows):for j in range(lines):ws.write(i, j, target_list[i][j])workbook.save(output_file_name)## output_file_name = 'putong_GM.xls' ## xlwt_data = tuple(xlwt_data) ## data = [] ## data.append(xlwt_data) ## print(data) ## save_excel(data, output_file_name)

在這個類的基礎上，運用分數階進行優化

FGM.py

##渡辺筆記 ##分數階灰色模型 from GM import Grey_model import numpy as np from scipy.special import gamma import xlrd import xlwt##更新累加序列（式1） def updata_fractional_accumulation(input_value,order):accumulation_value=np.zeros(len(input_value))for i in range(len(input_value)):for j in range(i+1):tmp=gamma(i-j+order-1+1)/(gamma(i-j+1)*gamma(order-1+1))accumulation_value[i]+=tmp*input_value[j]return accumulation_value##更新還原方法（參考式2和3） def updata_predicted_results(predicted_accumulation_value,order):if order!=1:predicted_one_accumulation_value=updata_fractional_accumulation(predicted_accumulation_value,1-order)else:predicted_one_accumulation_value=predicted_accumulation_valuepredicted_value=[float(%.2f%(predicted_one_accumulation_value[0]))]for i in range(len(predicted_accumulation_value)-1):predicted_value.append(float(%.2f%(predicted_one_accumulation_value[i+1]-predicted_one_accumulation_value[i])))return predicted_value## 打開xls，只能用這種格式了 work_book = xlrd.open_workbook('first_shuju.xls') ## 按sheet表名稱獲取sheet對象，名稱分大小寫 sheet_1 = work_book.sheet_by_name('Sheet1') col_0_value = sheet_1.col_values(6,2,19)print(col_0_value)##輸入變量，order表示分數階的階數 input_value = col_0_value number_of_forecast=10 order=0.1 A=Grey_model(input_value)##引入GM模型的計算模塊 A.accumulation_value=updata_fractional_accumulation(input_value,order)##更新累加值 A.predict(number_of_forecast) A.predicted_value= updata_predicted_results(A.predicted_accumulation_value,order)##更新預測值的還原方式 A.test_error()##誤差用MAPE值來衡量 print('預測值') print(A.predicted_value) print('\nMAPE值') print(A.MAPE)def save_excel(target_list, output_file_name):將數據寫入xls文件if not output_file_name.endswith('.xls'):output_file_name += '.xls'workbook = xlwt.Workbook(encoding='utf-8')ws = workbook.add_sheet(sheet1)rows = len(target_list)lines = len(target_list[0])for i in range(rows):for j in range(lines):ws.write(i, j, target_list[i][j])workbook.save(output_file_name)xlwt_data = A.predicted_value output_file_name = 'fenjieshu_putong_GM.xls' xlwt_data = tuple(xlwt_data) data = [] data.append(xlwt_data) print(data) save_excel(data, output_file_name)

3. 最短路徑

3.1. Floyd

感覺大多都能用 floyd 解決

% 這是floyd算法，以下注釋為渡邊所加 % floyd不能解決多條最短路徑和負權回路的問題 % 解稠密圖效果最佳，邊權可正可負，是雙源解法 % 對于稠密圖，效率要高于執行|V|次Dijkstra算法，也要高于執行|V|次SPFA算法。 % 缺點：時間復雜度比較高，不適合計算大量數據，注意 % 輸入示例 % W=[0 50 inf 40 25 10; 50 0 15 20 inf 25; inf 15 0 10 20 inf; 40 20 10 0 10 25;25 inf 20 10 0 55; 10 25 inf 25 55 0]; % [D, P, dis, path] = floyd(W, 1, 4) function [D, P, dis, path] = floyd(W, start, stop) %start為指定起始結點，stop為指定終止結點 D = W; %最短距離矩陣賦初值 n = length(D); %n為結點個數 P = zeros(n,n); %路由矩陣賦初值for i = 1:nfor j = 1:nP(i,j) = j;end endfor k = 1:nfor i = 1:nfor j = 1:nif D(i,k) + D(k,j) < D(i,j)D(i,j) = D(i,k) + D(k,j); %核心代碼P(i,j) = P(i,k);endendend endif nargin ~= 3errordlg('參數個數輸入有誤！', 'Warning!') elsedis = D(start, stop); %指定兩結點間的最短距離m(1) = start;i = 1;while P(m(i),stop) ~= stopk = i + 1;m(k) = P(m(i),stop);i = i + 1;endm(i+1) = stop;path = m; %指定兩結點之間的最短路徑 end

3.2. Dijkstra

Dijkstra 運算速度快，但是用不了負權圖

% 這是 dijkstra 算法，以下注釋為渡邊所加 % 不能用于負權圖，但是由于時間復雜度低，適合在大數據中運算 % 是單源解法 % 輸入示例 % W=[0 50 inf 40 25 10; 50 0 15 20 inf 25; inf 15 0 10 20 inf; 40 20 10 0 10 25;25 inf 20 10 0 55; 10 25 inf 25 55 0]; % [dis, path] = dijkstra(W, 1, 4) function [dis,path] = dijkstra( W,st,e ) % W 權值矩陣 st 搜索的起點 e 搜索的終點 n=length(W);%節點數 D = W(st,:); visit= ones(1:n); visit(st)=0; parent = zeros(1,n);%記錄每個節點的上一個節點path =[];for i=1:n-1temp = [];%從起點出發，找最短距離的下一個點，每次不會重復原來的軌跡，設置visit判斷節點是否訪問for j=1:nif visit(j)temp =[temp D(j)];elsetemp =[temp inf];endend[value,index] = min(temp);visit(index) = 0;%更新 如果經過index節點，從起點到每個節點的路徑長度更小，則更新，記錄前趨節點，方便后面回溯循跡for k=1:nif D(k)>D(index)+W(index,k)D(k) = D(index)+W(index,k);parent(k) = index;endendenddistance = D(e);%最短距離 %回溯法 從尾部往前尋找搜索路徑 t = e; while t~=st && t>0path =[t,path];p=parent(t);t=p; end path =[st,path];%最短路徑 dis = distance; end

3.3. matlab 自帶最短路徑

clc,clear,close all; % 程序基于但不基于 Dijkstra，解決了負權問題 s = [1 1 1 2 2 3 3 4 5 5 6 7]; t = [2 4 8 3 7 4 6 5 6 8 7 8]; weights = [10 10 1 10 1 -10 1 1 12 12 12 12]; % names = {'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H'}; % G = digraph(s,t,weights,names); G = digraph(s,t,weights); p = plot(G,'Layout','force','EdgeLabel',G.Edges.Weight);[path1,d] = shortestpath(G,6,8); highlight(p,path1,'EdgeColor','r')

4. 層次分析法

感覺這個方法就像我們土木工程所謂的“拍腦袋”，先拍腦袋得出一個比較矩陣，待反復測試符合規范后，就可以拿去比較了

% 本代碼只需要在提示輸入處，輸入構成的比較矩陣 就會輸出各指標的權重值。 % 例子： 選拔干部考慮5個條件：品德x1，才能x2，資歷x3，年齡x4，群眾關系x5。某決策人用成對比較法，得到成對比較陣如下： % [1,2,7,5,5; % 1/2,1,4,3,3; % 1/7,1/4,1,1/2,1/3; % 1/5,1/3,2,1,1; % 1/5,1/3,3,1,1] % 其中的x1=1;x2=2;x3=7;x4=5;x5=5; 怎么來的呢？ % 其實這些x的值就是由決策人決定了，對應值也就是決策人認為的重要性了 clc; clear; % 判斷矩陣A A= [1 3 3 31/3 1 3 51/5 1/3 1 31/5 1/5 1/3 1]; [m,n]=size(A); %獲取指標個數 RI=[0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45 1.49 1.51]; R=rank(A); %求判斷矩陣的秩 [V,D]=eig(A); %求判斷矩陣的特征值和特征向量，V特征值，D特征向量； tz=max(D); B=max(tz); %最大特征值 [row, col]=find(D==B); %最大特征值所在位置 C=V(:,col); %對應特征向量 CI=(B-n)/(n-1); %計算一致性檢驗指標CI CR=CI/RI(1,n); if CR<0.10disp('CI=');disp(CI);disp('CR=');disp(CR);disp('對比矩陣A通過一致性檢驗，各向量權重向量Q為：');Q=zeros(n,1);for i=1:nQ(i,1)=C(i,1)/sum(C(:,1)); %特征向量標準化endQ' %輸出權重向量 elsedisp('對比矩陣A未通過一致性檢驗，需對對比矩陣A重新構造'); end sc = Q';

5. 聚類分析

這類方法對于我的專業毫無實際用處，我并不需要知道管道應該分成幾類，因為規范早就規定好了的

聚類分析方法豐富，要根據實際情況謹慎選取

5.1. k-mean

% k_mean 算法，渡邊筆記 % 對處理大型數據集，該算法保持可收縮性，高效性； % 當簇接近正態分布時，效果最好； % 若簇中含有異常點，將導致均值偏離嚴重；（即對噪聲、孤立點敏感）； % 不適用于發現非凸形狀的簇，或大小差別大的簇 % 中k是事先給定的，這個k值難以估計，很多時候并不知道數據集應該分成多少個類別最合適 clear; clc;data = [11 2 02 2 24 3 3]; % 輸入數據 k = 3; % 分類數X = mapminmax(data',0,1)'; % 按列最小最大規范化到[0,1] T = kmeans(X,k); % 直接調用kmeans函數 for i = 1:kkong = find(T == i);fprintf('第 %d 類 :',i);disp(data(kong)') end

5.2. 層次聚類

% 層次聚類，渡邊筆記 % 看圖說話，分類清晰明了 data =[17.0 0.31 0.42 0.74 1.16 1.81 3.79 18.0 0.1 0.32 0.31 0.69 0.76 2.68 19.0 0.18 0.35 0.47 0.77 1.57 3.93 20.0 0.23 0.15 0.5 1.04 0.99 4.11 21.0 0.13 0.15 0.2 1.42 1.53 2.51]';X=mapminmax(data,0,1)'; % 按行最小最大規范化到[0,1] Y = pdist(X); % 計算矩陣X中樣本兩兩之間的距離，但得到的Y是個行向量 D = squareform(Y); % 將行向量的Y轉換成方陣，方便觀察兩點距離(方陣的對角線元素都是0) Z = linkage(Y); % 產生層次聚類樹，默認采用最近距離作為類間距離的計算公式 dendrogram(Z); % 圖示層次聚類樹 title('層次聚類法');ylabel('標準距離');

5.3. 直接聚類

% 直接聚類，渡筆記 % 不管機理，不管距離，比較愚蠢但是有用且方便 clear; clc; data =[36.6 4.46 36.72 6.37 32.39 4.63 15.43 36.80 12.46 47.21 18.66 9.22 1.69 10.76 67.88 2.76 39.1 4.2 36.92 1.87 15.15 48.9 2.85 36.85 7.23 38.29 3.51 11.27 60.5 2.23 27.25 6.8 45 8.73 9.99 ]; X=mapminmax(data,0,1); % 按列最小最大規范化到[0,1]T1=clusterdata(X,0.2); % 如果0<cutoff<2，則當不一致系數大于cutoff時，分到不同類(簇)中 T2=clusterdata(X,3); % 如果cutoff是一個≥2的整數，則形成的不同類別數為cutoff k1 = max(T1); k2 = max(T2); for i = 1:k1kong = find(T1 == i);fprintf('第 %d 類 :',i);disp(data(kong)') end disp('--------------------------------'); % for i = 1:k2 % kong = find(T2 == i); % fprintf('第 %d 類 :',i); % disp(data(kong)') % end

6. 最小生成樹

好像是個示例，太久遠忘了干什么用的了

zhuixiao_shengchengshu.m

clc,clear,close all; % 使用 Kruskal 的算法來計算 無向圖的最小生成樹 s = [1 1 1 2 5 3 6 4 7 8 8 8]; t = [2 3 4 5 3 6 4 7 2 6 7 5];% 使用加權邊創建并繪制一個立方體圖 weights = [100 10 10 10 10 20 10 30 50 10 70 10]; G = graph(s,t,weights); % figure; p = plot(G,'EdgeLabel',G.Edges.Weight);% 計算并在圖上方繪制圖的最小生成樹。T 包含的節點與 G 相同，但包含的邊僅為后者的子集 % figure; [T,pred] = minspantree(G); % p = plot(G,'EdgeLabel',G.Edges.Weight); highlight(p,T) % 高亮顯示

7. 誤差檢驗

7.1. 一般檢驗

看著用吧，我覺得沒什么意義，反正計算機計算出的誤差都挺小的，最多能在論文中湊字數

% 渡辺筆記 % 各類檢驗，除決定系數是1最好，都是0最好 YReal = [1 2 3 4 5]; YPred = [1 2 3 4 5.1]; % 平均絕對百分比誤差（MAPE） mape = mean(abs((YReal - YPred)./YReal)); disp(mape); % 均方根誤差（RMSE） rmse = sqrt(mean((YPred-YReal).^2)); disp(rmse); % 殘差平方和（SSE） sse = sum((YReal - YPred).^2); disp(sse); % 均方誤差（MSE） mse = mean(sum((YReal - YPred).^2)); disp(mse); % 平均絕對誤差（MAE） mae = mean(abs(YReal - YPred)); disp(mae); % 決定系數（R2-R-Square） r2 = 1 - (sum((YPred - YReal).^2) / sum((YReal - mean(YReal)).^2)); disp(r2)

7.2. T 檢驗

% 渡邊筆記 t檢驗 % 待檢驗的數據應該近服從正態分布 clc,clear;close all; x = [1 2 2 3 3 3 3 4 4]; y = [1 2 2 3 3 3 4 4 4]; ALPHA = 0.05; [H,P,CI,STATS]=ttest2(x,y, ALPHA); % 其中，x，y均為行向量（維度必須相同），各表示一組數據 % ALPHA為可選參數，表示設置一個值作為t檢驗執行的顯著性水平 % 在不設置ALPHA的情況下默認ALPHA為0.05，即計算x和y在5%的顯著性水平下是否來自同一分布(假設是否被接受） % H = 0，P > 0.05，表明零假設被接受，即x，y在統計上可看做來自同一分布的數據 % H = 1，P < 0.05，表明零假設被拒絕，即x，y在統計上認為是來自不同分布的數據，即有區分度 disp(['H = ',num2str(H)]); disp(['P = ',num2str(P)]); if H == 0disp('x，y在統計上可看做來自同一分布的數據') elsedisp('x，y在統計上可看做來自不同分布的數據') end

7.3. 回歸檢驗

clc,clear;close all; % 回歸檢驗；渡邊筆記 x = [0 1 2 3 4];y = [1 1 4 10 16];x = x'; y = y'; % x,y 需要轉置才可以使用X = [ones(size(y)) x]; % 定義回歸方程[B,BINT,R,RINT,STATS] = regress(y,X); % 回歸 % [B,BINT,R,RINT,STATS] = REGRESS(Y,X,alpha) % B 為回歸系數向量 % BINT 回歸系數的估計區間 % R 殘差 % RINT 置信區間 % STATS 用于檢驗回歸模型的統計量。有4個數值：判定系數r2，F統計量觀測值，檢驗的p的值，誤差方差的估計 % r2越接近1，回歸方程越顯著；F>F1?α(k,n?k?1)時拒絕H0，F越大，回歸方程越顯著；p<α時拒絕H0 % ALPHA 顯著性水平（缺少時默認0.05） y2 = B(2).* x + B(1); %預測值 plot(x', y' ,'+', x, y2); %回歸效果圖 disp('回歸分析統計量為'); disp([ 'R^2 : ',num2str(STATS(1))]); % 模型由 這么多原始數據 解釋，最優為 1 disp([ 'F : ',num2str(STATS(2))]); disp([ '檢驗 p : ',num2str(STATS(3))]); % 日常小于 0.05 最好 disp([ '誤差方差 : ',num2str(STATS(4))]); % 越小越好

7.4. 交叉檢驗

% 交叉檢驗，渡邊筆記 % 提高數據的利用率 % 有效的避免過學習以及欠學習狀態的發生,最后得到的結果也比較具有說服性 % 常用的 K 值有 3,6,10 等 clear,clc;close all;% 導入數據 data = reshape(1:144,12,12); % 每行的是這個樣本屬性的數據，最后一個數據是標簽 [data_r, data_c] = size(data); % 將數據樣本隨機分割為 n 部分，做 n 折交叉檢驗 n = 5; indices = crossvalind('Kfold', data_r, n); for i = 1 : n% 獲取第i份測試數據的索引邏輯值test = (indices == i);% 取反，獲取第i份訓練數據的索引邏輯值train = ~test;% 測試用的數據test_data = data(test, 1 : data_c - 1);test_label = data(test, data_c);% 訓練用的數據train_data = data(train, 1 : data_c - 1);train_label = data(train, data_c);% 下面放入檢驗數據的代碼end

8. 差分

% 數據差分，渡邊筆記 % 差分就是后一個減去前一個，使得數據平穩化 % n 階就是做 n 次差分，具體直接看例子就好 clc,clear;close all; Y = [0 5 15 30 50 75 105]; disp(['未差分時為 : ',num2str(Y)]) Y_1= diff(Y, 1); disp(['一階差分結果為 : ',num2str(Y_1)]) Y_2 = diff(Y, 2); disp(['二階差分結果為 : ',num2str(Y_2)]) Y_3 = diff(Y, 3); disp(['三階差分結果為 : ',num2str(Y_3)])

9. 立法白噪音

% 立法數白噪聲檢驗，渡邊筆記 % 如果易知原始數據不是平穩的，則不能做隨機性檢驗 % 接下來要求差分,目的: 變成平穩的數據 % p 如果比 0.05 小就不是白噪聲序列，可以使用時間序列 % 某種現象的指標數值按照時間順序排列而成的數值序列 % 因為時間序列是某個指標數值長期變化的數值表現 % 所以時間序列數值變化背后必然蘊含著數值變換的規律性 % 這些規律性就是時間序列分析的切入點 % 如果原數據平穩，但是沒通過，可以直接差分 clc,clear;close all;x = []; % 時間，一般做題就是順序時間排列 yanchi=[6 12 18]; % 通常是做6 12 18 24步延遲，這個數據的選擇上限請根據報錯來調整 y = [970279 1259308 1127571 1163959 1169540 1076938 991350 953275 951508 904434 889381 864015 836236 ]';; [h1] = adftest(y); %檢驗是否平穩 if h1 == 1disp('數據是平穩的');y_1 = y; elsedisp('數據是不平穩的');i = 1;while 1y_1 = diff(y,i); % 在這里對數據進行 n 階差分[h1,p1,adf,ljz] = adftest(y_1); %檢驗是否平穩if h1 == 1disp(['差分后是平穩的，做了 ',num2str(i),' 階差分']);subplot(2,2,3)plot(y_1); % 一階差分數據時序圖title('一階差分數據時序圖')subplot(2,2,4)autocorr(y_1); % 一階差分數據的自相關系數圖title('一階差分數據自相關系數圖');breakendi = i + 1;end end % 隨時間的變化值 subplot(2,2,1) plot(y); % 原始數據時序圖 title('原始數據時序圖') subplot(2,2,2) autocorr(y); % 原始數據的自相關系數圖 title('原始數據自相關系圖像')[H,pValue,Qstat,CriticalValue]=lbqtest(y_1,'lags',yanchi); % H.結果，pValue.p值, Qstat.卡方統計量 fprintf('%15s%15s%15s','延遲階數','卡方統計量','p值'); fprintf('\n'); for i=1:length(yanchi) % i = 1,時候為6，i = 2時候為12fprintf('%18f%19f%19f',yanchi(i),Qstat(i),pValue(i));fprintf('\n'); endif sum(find(pValue > 0.05))disp('但是沒通過立法白噪聲檢驗');i = 1;while 1y_1 = diff(y,i); % 在這里對數據進行 n 階差分[h1,p1,adf,ljz] = adftest(y_1); %檢驗是否平穩if h1 == 1disp(['差分后是平穩的，做了 ',num2str(i),' 階差分']);subplot(2,2,3)plot(y_1); % 一階差分數據時序圖title('一階差分數據時序圖')subplot(2,2,4)autocorr(y_1); % 一階差分數據的自相關系數圖title('一階差分數據自相關系數圖');breakendi = i + 1;end end% 再來一次[H,pValue,Qstat,CriticalValue]=lbqtest(y_1,'lags',yanchi); % H.結果，pValue.p值, Qstat.卡方統計量 fprintf('%15s%15s%15s','延遲階數','卡方統計量','p值'); fprintf('\n'); for i=1:length(yanchi) % i = 1,時候為6，i = 2時候為12fprintf('%18f%19f%19f',yanchi(i),Qstat(i),pValue(i));fprintf('\n'); end if sum(find(pValue < 0.05))disp('數據通過立法白噪聲檢驗'); elsedisp('數據沒通過立法白噪聲檢驗'); end

10. 熵權法

% 渡邊筆記 熵權法 clc;clear;close all; % 熵權法的思想是：變量數值變化越大，變異程度越大，則其權重應該更大；反之權重則越小x = [ 2.41 52.59 0 9.78 1.17 1.42 53.21 0 6.31 1.63 4.71 35.16 1 9.17 3.02 14.69 15.16 2.13 10.35 7.97 0.94 72.99 0 7.39 0.61 1.43 72.62 0 8.16 0.51 2.21 67.5 0 9.84 0.85 3.79 51.21 0 12.95 1.43 1.23 85.09 3.97 4.08 0.13 1.71 82.07 2.88 4.97 0.33 3.63 66.9 3.18 8.57 0.71 5.72 49.77 3.44 10.52 1.83 1.49 79.51 6.53 2.58 0.27 1.66 81.44 5.18 2.74 0.36 2.41 76.32 5.88 4.13 0.54 4.42 59.65 7.64 8.38 1.02 3.27 88.42 3.36 2.85 0.14 11.27 70.05 5.77 6.07 0.19 13.18 62.45 5.66 7.85 0.74 15.83 56.28 2.92 9.97 1.14 11.59 80.23 1.04 3.64 0.2 26.67 55.7 2.02 8.13 0.38 28.51 51.07 2.12 9.66 1.46 3.69 87.26 0 3.12 0.18 3.27 84.43 0 5.43 0.31 3.98 79.99 0 6.62 0.57 1.59 86.5 0 6.14 0.14 4.31 82.26 0 4.71 0.2 4.6 72.79 0 8.27 0.52 4.99 81.93 0 7.52 0.16 4.66 75.09 0 10.24 0.33 5.08 61.02 1.57 15.7 0.53 12.49 83.06 0 1.2 1.06 4.67 92.77 0 0.33 0.58 5.8 90.32 0 0.91 0.8 97.76 0 0 0 2.14 94.75 0 0 1.42 2.83 93.76 0 0 1.18 3.24 3.48 81.43 7.45 1.33 0.14 4.2 80 5.3 2.21 0.18 8.83 71.28 5.34 2.9 0.43 5.39 79.6 6.87 2.64 0.31 7.67 74.74 5.91 3.4 0.66 19.65 55.4 4.87 6.14 1.2 2.63 90.74 3.18 1.42 0.14 2.8 89.7 2.85 1.96 0.14 4.07 85.12 3.43 3.52 0.25 5.7 83.4 0 4.48 0.1 4.03 81.35 0 6.18 0.19 4.11 73.45 0 9.71 0.45 2.78 89.53 0 4.23 0.2 3.92 83.2 0 7.59 0.32 5.21 71.37 3.09 10.29 0.72 18.98 76.81 0 1.05 0.31 19.79 73.56 0 0.88 0.42 19.86 70.07 0 1.72 0.74 16.61 67.57 3.77 3.15 1.16 6.91 82.18 4.19 0 0.1 2.93 83.06 1.93 5.14 0.32 8.47 78.11 4.04 4.02 0.31 12.29 70.48 3.89 4.32 0.69 3.98 84.81 4.76 1.97 0.18 7.67 78.13 4.22 4.57 0.35 14.04 66.89 4.41 6.27 0.47 14.62 59.29 5.28 8.35 0.77 1.97 85.16 4.87 3.27 0.23 2.16 86.83 3.82 2.25 0.15 4.81 74.9 5.05 5.97 0.5 7.44 57.98 6.75 10.73 1.04 2.04 86.01 4.79 2.95 0.13 3.49 79.79 5.67 4.28 0.15 6.47 68.02 6.71 5.74 0.2 7.94 59.12 7.14 5.93 1.42 ];[m,n]=size(x); lamda=ones(1,n); % 人為修權，1代表不修改計算后的指標權重 for i=1:nx(:,i)=(x(:,i)-min(x(:,i)))/(max(x(:,i))-min(x(:,i)))+1; % 對原始數據進行非負數化、歸一化處理，值介于1-2之間 end for i=1:mfor j=1:np(i,j)=x(i,j)/sum(x(:,j));end end k=1/log(m); for i=1:mfor j=1:nif p(i,j)~=0e(i,j)=p(i,j)*log(p(i,j));elsee(i,j)=0;endend end for j=1:nE(j)=-k*sum(e(:,j)); end d=1-E; for j=1:nw(j)=d(j)/sum(d);% 指標權重計算 end for j=1:nw(j)=w(j)*lamda(j)/sum(w.*lamda);% 修改指標權重 end for i=1:mscore(i,1)=sum(x(i,:).*w);% 計算綜合分數 end disp('各指標權重為：') disp(w) disp('各項綜合分數為：') disp(score) Out = mean (score,1)

11. 數據修復

% 判斷缺失值和異常值并修復，順便光滑噪音，渡邊筆記 clc,clear;close all; x = 0:0.06:10; y = sin(x)+0.2*rand(size(x)); y(22:34) = NaN; y(89:95) = 50; testdata = [x' y'];subplot(2,2,1); plot(testdata(:,1),testdata(:,2)); title('原始數據');%% 判斷數據中是否存在缺失值 if sum(isnan(testdata(:)))disp('數據存在缺失值'); elsedisp('數據不存在缺失值'); end%% 判斷數據中是否存在異常值 % 1.mean 三倍標準差法 2.median 離群值法 3.quartiles 非正態的離群值法 % 4.grubbs 正態的離群值法 5.gesd 多離群值相互掩蓋的離群值法 choice_1 = 5; yichangzhi_fa = char('mean', 'median', 'quartiles', 'grubbs','gesd'); yi_chang = isoutlier(y,strtrim(yichangzhi_fa(choice_1,:))); if sum(yi_chang)disp('數據存在異常值'); elsedisp('數據不存在異常值'); end%% 對異常值賦空值 F = find(yi_chang == 1); y(F) = NaN; % 令數據點缺失 testdata = [x' y'];subplot(2,2,2); plot(testdata(:,1),testdata(:,2)); title('去除差異值');%% 對數據進行補全 % 數據補全方法選擇 % 1.線性插值 linear 2.分段三次樣條插值 spline 3.保形分段三次樣條插值 pchip % 4.移動滑窗插補 movmean chazhi_fa = char('linear', 'spline', 'pchip', 'movmean'); choice_2 = 2; if choice_2 ~= 4testdata_1 = fillmissing(testdata,strtrim(chazhi_fa(choice_2,:))); % strtrim 是為了去除字符串組的空格 elsetestdata_1 = fillmissing(testdata,'movmean',10); % 窗口長度為 10 的移動均值 endsubplot(2,2,3); plot(testdata_1(:,1),testdata_1(:,2)); title('數據補全結果');%% 進行數據平滑處理 % 濾波器選擇 1.Savitzky-golay 2.rlowess 3.rloess choice_3 = 2; lvboqi = char('Savitzky-golay', 'rlowess', 'pchip', 'rloess'); % 通過求 n 元素移動窗口的中位數，來對數據進行平滑處理 windows = 8; testdata_2 = smoothdata(testdata_1(:,2),strtrim(lvboqi(choice_3,:)),windows) ;subplot(2,2,4); plot(x,testdata_2) title('數據平滑結果');

12. 相關性分析

clc,clear,close all % 相關性分析，渡邊筆記 % 列為指標，行為數據data = [1 22 93 4 ]; % Pearson相關系數 r1 = corr(data,'type','Pearson'); disp(r1); % Kendall相關系數 r2 = corr(data,'type','Kendall'); disp(r2); % Spearman相關系數 r3 = corr(data,'type','Spearman'); disp(r3);

13. 馬氏鏈

13.1. 一般馬氏鏈

% 馬氏鏈模型，渡邊筆記 % 系統在每個時期所處的狀態是隨機的 % 從一時期到下時期的狀態按一定概率轉移 % 下時期狀態只取決于本時期狀態和轉移概率。即已知現在，將來與過去無關（無后效性） % 過去不影響未來的判斷，是馬氏鏈的本質 % P_0 是初始分布，P_n 是轉移矩陣，則在 n 未來的概率為 % P_0 * ( P_n ^ n ) clc,clear;close all;% 第一次買鹽，概率為 P_0 % 已知 買鹽傾向 P_n % P_0 = [0.2 0.4 0.4]; % P_n = [0.8 0.1 0.1 % 0.5 0.1 0.4 % 0.5 0.3 0.2]; % P = P_0 * P_n ^ 3; % disp(P)% 一般馬氏鏈模型 a=[160 120 120180 90 30180 30 90]; % 輸入統計矩陣 P_0 = [0.4 0.3 0.4]; % 初始概率% 求狀態轉移矩陣 [m,~] = size(a); ni=sum(a,2); %計算矩陣f的行和 for i = 1:mfor j = 1:mZT(i,j) = a(i,j)./ni(i); %求狀態轉移的頻率end end disp('轉移矩陣為 : '); disp(ZT) fprintf('%c', 8); % 刪掉換行符 % 求 n 期概率 n = 2; P = P_0 * ( ZT ^n ); disp(['第 ',num2str(n),' 期概率為 ',num2str(P)]);% 概率平衡時，有 P = P * P_n 恒成立 X = ZT'; X(1,:) = []; for i = 2:mX(i-1,i) = X(i-1,i) - 1; end X(m,:) = ones(1,m); Y = zeros(m-1,1); Y(m,:) = ones(1); x = X\Y; x = x'; disp(['n 趨向于無窮的 X 的解為 : ',num2str(x)]);

13.2. 帶利潤的馬氏鏈

% 帶利潤的馬氏鏈模型，渡邊筆記 clc,clear;close all;a=[5 54 6 ]; % 輸入統計矩陣 R = [9 33 -7]; % 輸入利潤矩陣% 對統計矩陣，求狀態轉移矩陣 [m,~] = size(a); ni=sum(a,2); %計算矩陣f的行和 for i = 1:mfor j = 1:mZT(i,j) = a(i,j)./ni(i); %求狀態轉移的頻率end end disp('轉移矩陣為 : '); disp(ZT) fprintf('%c', 8); % 刪掉換行符% 求 n 期利潤 n = 3;% n = 1 時 for i = 1:mfor j = 1:mv(i,j) = R(i,j) * ZT(i,j);endSum = sum(v,2); end % n = n 時for i = 1:mfor j = 1:mV(i,j) = Sum(j).^(n-1) * ZT(i,j);endV_Sum = sum(V,2);V_Sum = V_Sum + Sum;if n > 1 % 邊界條件disp(['處于狀態 ',num2str(i),' 的 ',num2str(n),' 期利潤為 ',num2str(V_Sum(i,:))])elsedisp(['處于狀態 ',num2str(i),' 的 ',num2str(n),' 期利潤為 ',num2str(Sum(i,:))])end end% 帶利潤的馬氏鏈一般不存在平衡 % 因為錢是賺不完的

14. 蒙特卡羅

% 渡邊筆記 % Monte_Carlo 是一種方法的概述，不是什么特定的算法 % 主要思想是利用大量隨機數來實現真實值的擬合 %% 例1 計算 0-3 上 x^2 的積分，正解應該是 9 N=500; x=unifrnd(0,3,N,1); % 設定隨機數及其終止上限 M = mean(3*x.^2);disp(['M = ',num2str(M)]); % 將隨機數代入原式求平均值，特別的，蒙特卡洛需要式乘以積分區間,下同 %% 例2 求相交面積 %給定曲線y =2 – x2 和曲線y3 = x2，曲線的交點為：P1(–1，1)、P2(1，1),函數圖像如下 x_1 = -1.5:.1:1.5; y_1 = 2 - x_1.^2; x_2 = -2:.1:2; y_2 = x_2.^2.^(1/3); subplot(1,2,1); plot(x_1,y_1,'r',x_2,y_2,'b') %曲線圍成平面有限區域，用蒙特卡羅方法計算區域面積。 P=rand(1000,2); %隨機產生1000個點 %定義x y 的范圍 x=2*P(:,1)-1; % 這一步是巧妙設置在(-1,2) y=2*P(:,2); II=find(y<=2-x.^2&y.^3>=x.^2); %找出在函數范圍的數 %計算索引的長度 M=length(II); %計算面積 S = 4*M/1000;disp(['S = ',num2str(S)]); subplot(1,2,2); plot(x(II),y(II),'g.'); close all % 刪除這個就可以顯示圖像 %% 例3 求圓周率 n = 1000; %總的實驗次數 m = 0; %落在圓中點的次數 %循環實驗 for i = 1:nx = 2 * rand - 1;y = 2 * rand - 1;if (x^2 + y^2 <= 1)m = m + 1;end end PI=4*m/n; disp(['PI = ',num2str(PI)]);

15. 模糊評價

先打包成一個 function

fuzzy_zhpj.m

function[B]=fuzzy_zhpj(model,A,R) %模糊綜合評判 B=[]; [m,s1]=size(A); [s2,n]=size(R); if(s1~=s2)disp('A的列不等于R的行'); elseif(model==1) %主因素決定型for(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)x=0;if(A(i,k)<R(k,j))x=A(i,k);elsex=R(k,j);endif(B(i,j)<x)B(i,j)=x;endendendendelseif(model==2) %主因素突出型for(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)x=A(i,k)*R(k,j);if(B(i,j)<x)B(i,j)=x;endendendendelseif(model==3) %加權平均型for(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)B(i,j)=B(i,j)+A(i,k)*R(k,j);endendendelseif(model==4) %取小上界和型for(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)x=0;x=min(A(i,k),R(k,j));B(i,j)=B(i,j)+x;endB(i,j)=min(B(i,j),1);endendelseif(model==5) %均衡平均型C=[];C=sum(R);for(j=1:n)for(i=1:s2)R(i,j)=R(i,j)/C(j);endendfor(i=1:m)for(j=1:n)B(i,j)=0;for(k=1:s1)x=0;x=min(A(i,k),R(k,j));B(i,j)=B(i,j)+x;endendendelsedisp('模型賦值不當');end end end

使用方法

clc,clear; A1=[0.1 0.9]; % A1=[0.4 0.35 0.15 0.1]; R=[ 1 4 2 2];sco = fuzzy_zhpj(1,A1,R) % fuzzy_zhpj(2,A1,R) % fuzzy_zhpj(3,A1,R) % fuzzy_zhpj(1,A2,R)

16. 模擬退火

16.1. 運用數學方法

clc,clear % 模擬退火(求最小) temperature = 100; % 初始溫度 final_temperature = 1; % 最低溫度 time_temperature = 1; % 用于計算步數 time_tuihuo = 10; % 退火步長 cooling_rate = 0.95; % 降溫步數 % 初始化隨機數生成器，以使結果具備可重復性 % rng(0,'twister'); % 生成范圍a-b的隨機數 a = -10; b = 10; % rang_math = (b-a).*rand(1000,1) + a; rang_math = (b-a).*rand + a; f = @(x)(...x.^2 ... % 定義函數，求最小值); current_old = f(rang_math);while final_temperature < temperaturerang_math = (b-a).*rand + a; % 隨機數current_new = f(rang_math); % 生成新解diff = current_new - current_old;% Metropolis 準則if(diff<0)||(rand < exp(-diff/(temperature))) % 接受新解的條件route = rang_math;current_old = current_new;time_temperature = time_temperature + 1;endif time_temperature == time_tuihuotemperature = temperature * cooling_rate;time_temperature = 0;endendplot([a:b],f([a:b])); % 原函數圖像 hold on; plot(route,current_old,'ko', 'linewidth', 2); % 作解的圖像

16.2. 運用 matlab 自帶工具

######## 2.1. 一元示例

% 渡邊自己手打的 % 模擬退火解一元函數 clc,clear;close all;% 設置區間 x = -10:1:10;f = @(x)(... x.^4 ... % 定義函數，求最小值 );% f = @(x)f_xy(x(1),x(2)); x0 = rand(1,1); % 退火開始點 lb = []; ub = []; % 退火實施范圍，可以不設置 % 實時觀測 options = saoptimset('MaxIter',20,... % 迭代次數 'StallIterLim',300,... % 最高溫度 'TolFun',1e-100,... % 最低溫度 'AnnealingFcn',@annealingfast,'InitialTemperature',... 100,'TemperatureFcn',@temperatureexp,'ReannealInterval',500,'PlotFcns',... {@saplotbestx, @saplotbestf, @saplotx, @saplotf,@saplottemperature});[jie,fval] = simulannealbnd(f,x0,lb,ub,options); fprintf('最優解為 x = %.10f, y = %.10f\n', jie, fval ); figure; fun = x.^2 + 2*x; % 再定義一次 plot(x,fun); hold on; plot(jie,fval,'ko', 'linewidth', 2);

######## 2.2. 二元示例

% 渡邊自己手打的 % 模擬退火解二元函數 clc,clear;close all;% 設置區間 xx = -10:1:10; yy = -10:1:10; [x, y]=meshgrid(xx, yy);f_xy = @(x,y)(... x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20 ... % 定義函數，求最小值 ); fun = x.^2 + y.^2 - 10*cos(2*pi*x) - 10*cos(2*pi*y) + 20; % 再定義一次 f = @(x)f_xy(x(1),x(2)); x0 = rand(1,2); % 退火開始點 lb = []; ub = []; % 退火實施范圍，可以不設置 % 實時觀測 options = saoptimset('MaxIter',20,... % 迭代次數 'StallIterLim',300,... % 最高溫度 'TolFun',1e-100,... % 最低溫度 'AnnealingFcn',@annealingfast,'InitialTemperature',... 100,'TemperatureFcn',@temperatureexp,'ReannealInterval',500,'PlotFcns',... {@saplotbestx, @saplotbestf, @saplotx, @saplotf,@saplottemperature});[jie,fval] = simulannealbnd(f,x0,lb,ub,options); fprintf('最優解為 x = %.10f, y = %.10f\n', jie(1),jie(2) ); fprintf('最優值為 z = %.10f\n',fval); figure; surf(x,y,fun); hold on; plot3(jie(1),jie(2),fval,'ko', 'linewidth', 3);

17. 神經網絡

17.1. BP

% 渡辺筆記，BP 神經網絡 clc,clear,close all; %輸入數據矩陣,每行為對應的屬性 p = [1 1 12 2 23 3 3]; %目標（輸出）數據矩陣 t = [6 6 6];%對訓練集中的輸入數據矩陣和目標數據矩陣進行歸一化處理 [pn, inputStr] = mapminmax(p); [tn, outputStr] = mapminmax(t);%建立BP神經網絡 net = newff(pn, tn, [3 7 2], {'purelin', 'logsig', 'purelin'}); %每10輪回顯示一次結果 net.trainParam.show = 10; %最大訓練次數 net.trainParam.epochs = 5000; %網絡的學習速率 net.trainParam.lr = 0.05; %訓練網絡所要達到的目標誤差 net.trainParam.goal = 0.65 * 10^(-3); %網絡誤差如果連續6次迭代都沒變化，則matlab會默認終止訓練。為了讓程序繼續運行，用以下命令取消這條設置 net.divideFcn = ''; %開始訓練網絡 net = train(net, pn, tn); %使用訓練好的網絡，基于訓練集的數據對BP網絡進行仿真得到網絡輸出結果 %(因為輸入樣本（訓練集）容量較少，否則一般必須用新鮮數據進行仿真測試) answer = sim(net, pn); %反歸一化 answer1 = mapminmax('reverse', answer, outputStr);

17.2. 灰色

clc,clear,close all; % 渡邊筆記 灰色神經網絡filename_1 = 'first_shuju.xlsx'; sheet_1 = ''; range_3 = 'G3:G18'; for j = 1:16years(j) = j+3; end [num_3, ~,~] = xlsread(filename_1,sheet_1,range_3); [m,~] = size(num_3); IN=1:m; for i = 1:msr(i) = num_3(i); end OUT=sr; [X,minx,maxx,T,mint,maxt]=premnmx(IN,OUT); q=50; q1=0; q0=70; while(q1<q)q=q0;[M,N]=size(X);[L,N]=size(T);net=newff(minmax(X),[q,L],{'tansig','purelin'},'trainlm');net.trainParam.lr=0.05;net.trainParam.epochs=10000;net.trainParam.goal=1e-6;[net,tr]=train(net,X,T);Y=sim(net,X);Y=postmnmx(Y,mint,maxt);%灰色關聯分析，調整網絡隱層節點p=0.5;e=0.5;%此兩個系數的設定是根據一些論文，已經實驗的嘗試得出的an=repmat(net.b{1},1,N);op=tansig(net.iw{1,1}*X+an);op1=op';T0=T';T1=repmat(T0,1,q);DIF=abs(T1-op1);MIN=min(min(DIF));MAX=max(max(DIF));Si=(MIN+p*MAX)./(DIF+p*MAX);ri=sum(Si)/N;D=find(ri>=e);[q0,q1]=size(D);q0=q1; end q0; ri; D; q=q1;%進行測試 PRD=1:m; PRD=PRD'; P=tramnmx(PRD,minx,maxx); TNEW=sim(net,P'); TNEW=postmnmx(TNEW,mint,maxt);YY=OUT; YC=TNEW; figure plot(years,YY,'b+-',years,YC,'r') legend('灰色神經網絡','生態用水量'); xlabel('年份'); RES0=YC-YY; res0=RES0./YY; figure bar(years,res0); xlabel('訓練次數'); ylabel('訓練誤差');

17.3. 擬合

擬合神經網絡是什么我忘記了，反正能用就行了

shenjin.m

function [net,tr] = shenjin(choice,P,T) net=newff(minmax(P),[20,1],... % 隱單元的神經元數 20，輸出 1，可調 {'tansig','purelin'}); if(choice==1)% 采用 L-M 優化算法 TRAINLMnet.trainFcn='trainlm';% 設置訓練參數net.trainParam.epochs = 500;net.trainParam.goal = 1e-6;net=init(net);% 重新初始化 elseif(choice==2)% 采用貝葉斯正則化算法 TRAINBRnet.trainFcn='trainbr';% 設置訓練參數net.trainParam.epochs = 500;randn('seed',192736547);net = init(net);% 重新初始化 elseif(choice==3)% 采用動量梯度下降算法 TRAINGDM% 當前輸入層權值和閾值inputWeights=net.IW{1,1};inputbias=net.b{1};% 當前網絡層權值和閾值layerWeights=net.LW{2,1};layerbias=net.b{2};% 設置訓練參數net.trainParam.show = 50;net.trainParam.lr = 0.05;net.trainParam.mc = 0.9;net.trainParam.epochs = 1000;net.trainParam.goal = 1e-3; end % 調用相應算法訓練 BP 網絡 [net,tr]=train(net,P,T); end

調用方法

% L-M 優化算法（trainlm）和貝葉斯正則化算法（trainbr），用以訓練 BP 網絡 % 使其能夠擬合某一附加有白噪聲的正弦樣本數據，渡邊筆記 % 神經網絡的初始和常數設置，請在shenjin.m修改 clc,clear;close all;% 0.全選，如果數據不多時可以操作 % 1.L-M 優化算法 TRAINLM % 2.貝葉斯正則化算法 TRAINBR % 3.動量梯度下降算法 TRAINGDM % trainbfg---BFGS算法（擬牛頓反向傳播算法）訓練函數； % trainbr---貝葉斯歸一化法訓練函數； % traincgb---Powell-Beale共軛梯度反向傳播算法訓練函數； % traincgp---Polak-Ribiere變梯度反向傳播算法訓練函數； % traingd---梯度下降反向傳播算法訓練函數； % traingda---自適應調整學習率的梯度下降反向傳播算法訓練函數； % traingdm---附加動量因子的梯度下降反向傳播算法訓練函數； % traingdx---自適應調整學習率并附加動量因子的梯度下降反向傳播算法訓練函數； % trainlm---Levenberg-Marquardt反向傳播算法訓練函數； % trainoss---OSS（one step secant）反向傳播算法訓練函數； % trainrp---RPROP（彈性BP算法）反向傳播算法訓練函數； % trainscg---SCG（scaled conjugate gradient）反向傳播算法訓練函數； % trainb---以權值/閾值的學習規則采用批處理的方式進行訓練的函數； % trainc---以學習函數依次對輸入樣本進行訓練的函數； % trainr---以學習函數隨機對輸入樣本進行訓練的函數。 % 當調用train函數時，上述訓練函數被用于訓練網絡 L = 3; %共有三種方法,我覺得夠用了 choice = 0;% P 為輸入矢量 P = [1 2 3 4 5]; % T 為目標矢量 T = [1 2 3 4 5];% 繪制樣本數據點 % 3種全選，比較誰更合適 % choice = 0 if choice == 0for i = 1:L% 神經網絡[net,tr] = shenjin(i,P,T);% 對 BP 網絡進行仿真A = sim(net,P);title(['神經網絡 方法',num2str(i),' 計算結果']);% 計算仿真誤差E = T - A;MSE(i) = mse(E); % mse 函數用于計算均方誤差% 均方誤差（mean-square error, MSE）是反映估計量與被估計量之間差異程度的一種度量endfor i = 1:Ldisp(['神經網絡 方法',num2str(i),' 的_均方誤差_為 : ',num2str(MSE(i))])end[~,Min] = min(MSE);disp(['神經網絡 方法',num2str(Min),' 比較優']) elsesubplot(1,2,1)plot(P,T,'+');title('原始數據點');% 神經網絡[net,tr] = shenjin(i,P,T);% 對 BP 網絡進行仿真A = sim(net,P);title(['神經網絡 方法',num2str(choice),' 計算結果']);% 計算仿真誤差E = T - A;MSE = mse(E); % mse 函數用于計算均方誤差% 均方誤差（mean-square error, MSE）是反映估計量與被估計量之間差異程度的一種度量disp(['神經網絡 方法',num2str(choice),' 的_均方誤差_為 : ',num2str(MSE)]) end% 這個用于分析是否可以很好的進行線性回歸 % [m,b,r] = postreg(A, P); % 神經網絡的簡單使用方法 % P = [1 2 3]; % [Y] = sim(net,P); % 輸入即可輸出，維度應與元數據保存一致

18. 貪心算法

% 渡邊筆記 % 貪心算法是一種思想，即尋找當前最優解，至使結果最優解 % 貪心算法必須經過證明才可以使用 % 其證明圍繞著：整個問題的最優解一定由在貪心策略中存在的子問題的最優解得來的 % 一旦貪心算法得以證明，它就是最快最優的算法 % 實際上，貪心算法適用的情況很少 % 判斷一個問題是否適合用貪心法求解，目前還沒有一個通用的方法，在競賽中，需要憑個人的經驗來判斷 % 貪心算法畢竟是貪心算法，只能緩一時，而不是長久之計， % 問題的模型、參數對貪心算法求解結果具有決定性作用，這在某種程度上是不能接受的 % 于是聰明的人類就發明了各種智能算法（也叫啟發式算法） % 但在我看來所謂的智能算法本質上就是貪心算法和隨機化算法結合 % 例如傳統遺傳算法用的選擇策略就是典型的貪心選擇，正是這些貪心算法和隨機算法的結合，我們才看到今天各種各樣的智能算法 %% 例1 設有n個正整數，將它們連接成一排，組成一個最大的多位整數。 % n=3時，3個整數13，312，343，連成的最大整數為34331213。 % n=4時，4個整數7，13，4，246，連成的最大整數為7424613。 % 此題很容易想到使用貪心法，比如把整數按從大到小的順序連接起來，例子符合，但測試結果不全對 % 反例：12，121應該組成12121而非12112 % 正確的標準是：先把整數轉換成字符串，然后比較a+b和b+a，如果a+b>=b+a，就把a排在b的前面，反之則把a排在b的后面 %% 例2 最短路徑問題 clear,clc;close all; n=4; % 設置點數 x = zeros(1,n); % 產生一個行向量存儲坐標 y = zeros(1,n); luxian = 1:1:n; % 生成1到n的矩陣，存儲路線 paixu = 1:1:n;% 這一步是作出隨機點的圖 for i = 1 : nx(i) = randn * n ; %生成正態分布的隨機坐標點y(i) = randn * n ;hold ontext(x(i)+0.3,y(i)+0.3,num2str(i)) %用text做好標記， end% 這一步是計算點與點之間的距離 d = zeros(n) ; % 生成一個n*n的矩陣用來存儲點與點之間的距離，可刪去 for i = 1 : nfor j = 1 : nd(i,j) = sqrt( ( x(i) - x(j) ) ^ 2 + ( y(i) - y(j) ) ^ 2);end endluxian(1) = 1; %設置起點，路線的起點是1 num = 1; for a = 1:(n-2) %去除起點和終點需要n-2次判斷 paixu(:,1)=[]; %刪除上一個最優點 dis = zeros(1,(n-a)); %用來存剩下各個點的距離for b = 1:(n-a) %用來獲取剩下各個點的距離dis(b) = d (num, paixu(b));end num1 = find( dis == min(dis) ); %根據距離最小得到最優點位置 t = paixu(1); %通過中間變量互換位置paixu(1) = paixu(num1); %最優點位置換到第一個paixu(num1) = t;num = paixu(1); %獲取下次進行操作的數luxian(a+1) = paixu(1); %將最優點存入luxian數組（空出開頭） end luxian(n) = paixu(2); %補上最后一個點plot(x(luxian),y(luxian),'--o') ; %標出點用虛線相連得到路徑 grid on;title('貪心算法計算最優路徑');

19. 小波分析

%%初始化程序 clear,clc t1=clock;%% 載入噪聲信號數據，數據為.mat格式，并且和程序放置在同一個文件夾下 x = 0:0.06:10; YSJ = sin(x)+0.2*rand(size(x));plot(YSJ); title('原始數據');%% 數據預處理，數據可能是存儲在矩陣或者是EXCEL中的二維數據，銜接為一維的，如果數據是一維數據，此步驟也不會影響數據 [c,l]=size(YSJ); Y=[]; for i=1:cY=[Y,YSJ(i,:)]; end [c1,l1]=size(Y); X=[1:l1];%% 繪制噪聲信號圖像 figure(1); plot(X,Y); xlabel('橫坐標'); ylabel('縱坐標'); title('原始信號');%% 硬閾值處理 lev=3; xd=wden(Y,'heursure','h','one',lev,'db4');%硬閾值去噪處理后的信號序列 figure(2) plot(X,xd) xlabel('橫坐標'); ylabel('縱坐標'); title('硬閾值去噪處理') set(gcf,'Color',[1 1 1])%% 軟閾值處理 lev=3; xs=wden(Y,'heursure','s','one',lev,'db4');%軟閾值去噪處理后的信號序列 figure(3) plot(X,xs) xlabel('橫坐標'); ylabel('縱坐標'); title('軟閾值去噪處理') set(gcf,'Color',[1 1 1]) %% 固定閾值后的去噪處理 lev=3; xz=wden(Y,'sqtwolog','s','sln',lev,'db4');%固定閾值去噪處理后的信號序列 figure(4) plot(X,xz); xlabel('橫坐標'); ylabel('縱坐標'); title('固定閾值后的去噪處理') set(gcf,'Color',[1 1 1]) %% 計算信噪比SNR Psig=sum(Y*Y')/l1; Pnoi1=sum((Y-xd)*(Y-xd)')/l1; Pnoi2=sum((Y-xs)*(Y-xs)')/l1; Pnoi3=sum((Y-xz)*(Y-xz)')/l1; SNR1=10*log10(Psig/Pnoi1); SNR2=10*log10(Psig/Pnoi2); SNR3=10*log10(Psig/Pnoi3); %% 計算均方根誤差RMSE RMSE1=sqrt(Pnoi1); RMSE2=sqrt(Pnoi2); RMSE3=sqrt(Pnoi3); %% 輸出結果 disp('-------------三種閾值設定方式的降噪處理結果---------------'); disp(['硬閾值去噪處理的SNR=',num2str(SNR1),'，RMSE=',num2str(RMSE1)]); disp(['軟閾值去噪處理的SNR=',num2str(SNR2),'，RMSE=',num2str(RMSE2)]); disp(['固定閾值后的去噪處理SNR=',num2str(SNR3),'，RMSE=',num2str(RMSE3)]); t2=clock; tim=etime(t2,t1); disp(['------------------運行耗時',num2str(tim),'秒-------------------'])

20. 遺傳算法

我是非常不建議使用這個的，太復雜太麻煩了，每次換算遺傳因子都要想半天，而且代碼冗余，又長又臭

20.1. TSP 示例

cross.m

%交叉操作函數 cross.m function [A,B]=cross(A,B) L=length(A); if L<10W=L; elseif ((L/10)-floor(L/10))>=rand&&L>10W=ceil(L/10)+8; elseW=floor(L/10)+8; end %%W為需要交叉的位數 p=(L-W+1);%隨機產生一個交叉位置 %fprintf('p=%d ',p);%交叉位置 for i=1:Wx=find(A==B(1,p+i-1));y=find(B==A(1,p+i-1));[A(1,p+i-1),B(1,p+i-1)]=exchange(A(1,p+i-1),B(1,p+i-1));[A(1,x),B(1,y)]=exchange(A(1,x),B(1,y)); endend

exchange.m

%對調函數 exchange.mfunction [x,y]=exchange(x,y) temp=x; x=y; y=temp;end

fit.m

%適應度函數fit.m，每次迭代都要計算每個染色體在本種群內部的優先級別，類似歸一化參數。越大約好！ function fitness=fit(len,m,maxlen,minlen) fitness=len; for i=1:length(len)fitness(i,1)=(1-(len(i,1)-minlen)/(maxlen-minlen+0.0001)).^m; end

Mutation.m

%變異函數 Mutation.mfunction a=Mutation(A) index1=0;index2=0; nnper=randperm(size(A,2)); index1=nnper(1); index2=nnper(2); %fprintf('index1=%d ',index1); %fprintf('index2=%d ',index2); temp=0; temp=A(index1); A(index1)=A(index2); A(index2)=temp; a=A;end

myLength.m

%染色體的路程代價函數 mylength.m function len=myLength(D,p)%p是一個排列 [N,NN]=size(D); len=D(p(1,N),p(1,1)); for i=1:(N-1)len=len+D(p(1,i),p(1,i+1)); endend

plot_route.m

%連點畫圖函數 plot_route.mfunction plot_route(a,R) scatter(a(:,1),a(:,2),'rx'); hold on; plot([a(R(1),1),a(R(length(R)),1)],[a(R(1),2),a(R(length(R)),2)]); hold on; for i=2:length(R)x0=a(R(i-1),1);y0=a(R(i-1),2);x1=a(R(i),1);y1=a(R(i),2);xx=[x0,x1];yy=[y0,y1];plot(xx,yy);hold on; endend

調用方法

% 渡邊筆記 % 多種群遺傳算法解決 TSP 問題 % TSP問題即旅行商問題 % 經典的TSP可以描述為 % 一個商品推銷員要去若干個城市推銷商品，該推銷員從一個城市出發，需要經過所有城市后，回到出發地 % 應如何選擇行進路線，以使總的行程最短。從圖論的角度來看，該問題實質是在一個帶權完全無向圖中，找一個權值最小的哈密爾頓回路clear,clc;close all; % 輸入參數 M=100; %種群的個數 ITER=200; %迭代次數 m=2; %適應值歸一化淘汰加速指數 Pc=0.8; %交叉概率 Pmutation=0.05; %變異概率 % 城市坐標 pos=[0 01 1-1 -12 3] ;%生成城市之間距離矩陣 [N,~] = size(pos); %%城市的個數 D=zeros(N,N); for i=1:Nfor j=i+1:Ndis=(pos(i,1)-pos(j,1)).^2+(pos(i,2)-pos(j,2)).^2;D(i,j)=dis^(0.5);D(j,i)=D(i,j);end end%生成初始群體 popm=zeros(M,N); for i=1:Mpopm(i,:)=randperm(N);%隨機排列，比如[2 4 5 6 1 3] end %隨機選擇一個種群 R=popm(1,:);%初始化種群及其適應函數 fitness = zeros(M,1); len = zeros(M,1); for i=1:M % 計算每個染色體對應的總長度len(i,1)=myLength(D,popm(i,:)); end maxlen=max(len);%最大回路 minlen=min(len);%最小回路 fitness=fit(len,m,maxlen,minlen); rr=find(len==minlen);%找到最小值的下標，賦值為rr R=popm(rr(1,1),:);%提取該染色體，賦值為Rfitness=fitness/sum(fitness); distance_min=zeros(ITER+1,1); %%各次迭代的最小的種群的路徑總長 nn=M; iter=0; while iter<=ITER%選擇操作p=fitness./sum(fitness);q=cumsum(p);%累加for i=1:(M-1)len_1(i,1)=myLength(D,popm(i,:));r=rand;tmp=find(r<=q);popm_sel(i,:)=popm(tmp(1),:);end[fmax,indmax]=max(fitness);%求當代最佳個體popm_sel(M,:)=popm(indmax,:);%交叉操作nnper=randperm(M); % A=popm_sel(nnper(1),:); % B=popm_sel(nnper(2),:);for i=1:M*Pc*0.5A=popm_sel(nnper(i),:);B=popm_sel(nnper(i+1),:);[A,B]=cross(A,B); % popm_sel(nnper(1),:)=A; % popm_sel(nnper(2),:)=B;popm_sel(nnper(i),:)=A;popm_sel(nnper(i+1),:)=B;end %變異操作for i=1:Mpick=rand;while pick==0pick=rand;endif pick<=Pmutationpopm_sel(i,:)=Mutation(popm_sel(i,:));endend%%求適應度函數NN=size(popm_sel,1);len=zeros(NN,1);for i=1:NNlen(i,1)=myLength(D,popm_sel(i,:));endmaxlen=max(len);minlen=min(len);distance_min(iter+1,1)=minlen;fitness=fit(len,m,maxlen,minlen);rr=find(len==minlen);R=popm_sel(rr(1,1),:); % for i=1:N % fprintf('%d ',R(i)); % end % fprintf('\n');popm=[];popm=popm_sel;iter=iter+1;%pause(1); end%畫出所有城市坐標 subplot(1,3,1); scatter(pos(:,1),pos(:,2),'rx');fprintf('最短路程是 '); for i=1:Nfprintf('%d ',R(i));%把R順序打印出來 end fprintf('\n最短距離是 %d\n',minlen); title('點圖'); subplot(1,3,2); plot_route(pos,R); % 路程圖 title('路程圖'); subplot(1,3,3); plot(distance_min); % 最短距離隨迭代次數的變化 title('最短距離的變化');

21. 蟻群算法

21.1. TSP 示例

yiyu%% 蟻群算法求解TSP問題%% 清空環境變量 clear clc%% 導入數據 load citys_data.mat%% 計算城市間相互距離 n = size(citys,1); D = zeros(n,n); for i = 1:nfor j = 1:nif i ~= jD(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2));elseD(i,j) = 1e-4;endend end%% 初始化參數 m = 50; % 螞蟻數量 alpha = 1; % 信息素重要程度因子 beta = 5; % 啟發函數重要程度因子 rho = 0.1; % 信息素揮發因子 Q = 1; % 常系數 Eta = 1./D; % 啟發函數 Tau = ones(n,n); % 信息素矩陣 Table = zeros(m,n); % 路徑記錄表 iter = 1; % 迭代次數初值 iter_max = 200; % 最大迭代次數 Route_best = zeros(iter_max,n); % 各代最佳路徑 Length_best = zeros(iter_max,1); % 各代最佳路徑的長度 Length_ave = zeros(iter_max,1); % 各代路徑的平均長度%% 迭代尋找最佳路徑 while iter <= iter_max% 隨機產生各個螞蟻的起點城市start = zeros(m,1);for i = 1:mtemp = randperm(n);start(i) = temp(1);endTable(:,1) = start;% 構建解空間citys_index = 1:n;% 逐個螞蟻路徑選擇for i = 1:m% 逐個城市路徑選擇for j = 2:ntabu = Table(i,1:(j - 1)); % 已訪問的城市集合(禁忌表)allow_index = ~ismember(citys_index,tabu);allow = citys_index(allow_index); % 待訪問的城市集合P = allow;% 計算城市間轉移概率for k = 1:length(allow)P(k) = Tau(tabu(end),allow(k))^alpha * Eta(tabu(end),allow(k))^beta;endP = P/sum(P);% 輪盤賭法選擇下一個訪問城市Pc = cumsum(P);target_index = find(Pc >= rand);target = allow(target_index(1));Table(i,j) = target;endend% 計算各個螞蟻的路徑距離Length = zeros(m,1);for i = 1:mRoute = Table(i,:);for j = 1:(n - 1)Length(i) = Length(i) + D(Route(j),Route(j + 1));endLength(i) = Length(i) + D(Route(n),Route(1));end% 計算最短路徑距離及平均距離if iter == 1[min_Length,min_index] = min(Length);Length_best(iter) = min_Length;Length_ave(iter) = mean(Length);Route_best(iter,:) = Table(min_index,:);else[min_Length,min_index] = min(Length);Length_best(iter) = min(Length_best(iter - 1),min_Length);Length_ave(iter) = mean(Length);if Length_best(iter) == min_LengthRoute_best(iter,:) = Table(min_index,:);elseRoute_best(iter,:) = Route_best((iter-1),:);endend% 更新信息素Delta_Tau = zeros(n,n);% 逐個螞蟻計算for i = 1:m% 逐個城市計算for j = 1:(n - 1)Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i);endDelta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i);endTau = (1-rho) * Tau + Delta_Tau;% 迭代次數加1，清空路徑記錄表iter = iter + 1;Table = zeros(m,n); end%% 結果顯示 [Shortest_Length,index] = min(Length_best); Shortest_Route = Route_best(index,:); disp(['最短距離:' num2str(Shortest_Length)]); disp(['最短路徑:' num2str([Shortest_Route Shortest_Route(1)])]);%% 繪圖 figure(1) plot([citys(Shortest_Route,1);citys(Shortest_Route(1),1)],...[citys(Shortest_Route,2);citys(Shortest_Route(1),2)],'o-'); grid on for i = 1:size(citys,1)text(citys(i,1),citys(i,2),[' ' num2str(i)]); end text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),' 起點'); text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),' 終點'); xlabel('城市位置橫坐標') ylabel('城市位置縱坐標') title(['蟻群算法優化路徑(最短距離:' num2str(Shortest_Length) ')']) figure(2) plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r:') legend('最短距離','平均距離') xlabel('迭代次數') ylabel('距離') title('各代最短距離與平均距離對比')

citys_data.mat

1304 2312 3639 1315 4177 2244 3712 1399 3488 1535 3326 1556 3238 1229 4196 1004 4312 790 4386 570 3007 1970 2562 1756 2788 1491 2381 1676 1332 695 3715 1678 3918 2179 4061 2370 3780 2212 3676 2578 4029 2838 4263 2931 3429 1908 3507 2367 3394 2643 3439 3201 2935 3240 3140 3550 2545 2357 2778 2826 2370 2975

22. 元胞自動機

22.1. 森林火災示例

clc,clear; n=200; Pltg = 5e-6; Pgrw = 1e-2; NW = [n 1:n-1]; SE = [2:n 1]; veg = zeros(n); imh = image( cat(3,(veg == 1),(veg == 2),zeros(n)) );for i = 1:30num = (veg(NW,:)==1) + (veg(:,NW)==1) + (veg(:,SE)==1) + (veg(SE,:)==1);veg = 2*( (veg==2) | (veg==0 & rand(n)<Pgrw) ) - ...((veg==2) & (num>0|rand(n)<Pltg ) );set(imh,'cdata',cat(3,(veg==1),(veg==2),zeros(n)));drawnow end

22.2. 交通流示例

gaplength.m

function gaps = gaplength(x,L) ncar = length(x); gaps = inf*ones(1,ncar); if ncar>1gaps = x([2:end 1]) -x;gaps(gaps<0) = gaps(gaps<0) +L; end

ns.m

function [flux,vmean] = ns(rho,p,L,tmax)ncar = ceil(L*rho); vmax = 5; x = sort(randperm(L,ncar)); v = vmax*ones(1,ncar); flux = 0;vmean = 0;h = plotcirc(L,x,0.5);for t = 1:tmaxv = min(v+1,vmax);gaps = gaplength(x,L);v = min(v,gaps - 1);vdrops = (rand(1,ncar)<p);v = max(v - vdrops,0);x = x + v;passed = x>L;x(passed) = x(passed)-L;if t>tmax/2flux = flux +sum(v)/L;vmean = vmean+mean(v);endplotcirc(L,x,0.5,h); endflux = flux/(tmax/2); vmean = vmean/(tmax/2);

plotcirc.m

function h=plotcirc(L,x,dt,h) W = 0.05;R=1; ncar = length(x);theta = 0:2*pi/L:2*pi; xc = cos(theta);yc = sin(theta); xin = (R-W/2)*xc;yin = (R-W/2)*yc; xot = (R+W/2)*xc;yot = (R+W/2)*yc;xi = [xin(x);xin(x+1);xot(x+1);xot(x)]; yi = [yin(x);yin(x+1);yot(x+1);yot(x)];if nargin == 3color = randperm(ncar);h = fill(xi,yi,color);hold onplot([xin;xot],[yin;yot],'k','linewidth',1.5);plot([xin;xot]',[yin;yot]','k','linewidth',1.5); elsefor i=1:ncarset(h(i),'xdata',xi(:,i),'ydata',yi(:,i));end end pause(dt)

調用方法

clc; clear; [flux,vmean] = ns(0.2,1e-9,25,30);

23. 主成分分析

% 渡邊筆記 主成分分析 % 當研究的問題涉及多個變量，并且變量間相關性明顯，即包含的信息有所重疊時 % 可以考慮用主成分分析的方法，這樣更容易抓住事物的主要矛盾，使問題簡化 clc,clear;close all; % 與主成分分析相關的 MATLAB 函數有pcacov、princomp函數% pcacov函數用來根據 協方差矩陣 或 相關系數矩陣 進行主成分分析 % [COEFF,latent,explained]=pcacov（V） % V 是總體或樣本的協方差矩陣或相關系數矩陣，對于p維總體，V是pxp的矩陣 % 輸出參數 COEFF 是p個主成分的系數矩陣，它是pxp的矩陣，它的第i列是第i個主成分的系數向量 % 輸出參數 latent 是p個主成分的方差構成的向量，即V的p個特征值的大小(從大到小)構成的向量 % 輸出參數 explained 是p個主成分的貢獻率向量，已經轉化為百分比 V = [1 0.79 0.36 0.76 0.25 0.510.79 1 0.31 0.55 0.17 0.350.36 0.31 1 0.35 0.64 0.580.76 0.55 0.35 1 0.16 0.380.25 0.17 0.64 0.16 1 0.630.51 0.35 0.58 0.38 0.63 1 ]; [a,~] = size(V); [COEFF_1,latent,explained] = pcacov(V); result_1(1,:) = {'特征值', '差值', '貢獻率', '累積貢獻率'}; result_1(2: a+1 ,1) = num2cell(latent); % diff函數式用于求導數和差分的. result_1(2: a ,2) = num2cell(-diff(latent)); % cumsum函數通常用于計算一個數組各行的累加值。 result_1(2: a+1 ,3:4) = num2cell([explained, cumsum(explained)]); % disp(result_1);% princomp函數用來根據樣本觀測值矩陣進行主成分分析，其調用格式如下： % [COEFF,SCORE,latent,tsquare] = princomp(X,'econ') % 根據樣本觀測值矩陣X進行主成分分析。輸入參數X是n行p列的矩陣，每一行對應一個觀測(樣品)，每一列對應一個變量 % 輸出參數COEFF是p個主成分析的系數矩陣，是pxp的矩陣，它的第i列對應第i個主成分的系數向量 % 輸出參數SCORE是n個樣品的p個主成分得分矩陣，它是n行p列的矩陣 % 每一行對應一個觀測，每一列對應一個主成分，第i行第j列元素表示第i個樣品的第j個主成分得分 % SCORE與X是一一對應的關系，是X在新坐標系中的數據，可以通過X*系數矩陣得到 % 返回樣本協方差矩陣的特征值向量latent，它是由p個特征值構成的列向量，其中特征值按降序排列。 % 返回一個包含n個元素的列向量tsquare，它的第i個元素的第i個觀測對應的霍特林T^2統計量 % 表述了第i個觀測與數據集(樣本觀測矩陣)的中心之間的距離，可用來尋找遠離中心的極端數據 % 通過設置參數‘econ’參數，使得當 n <= p 時 % 只返回latent中的前n-1個元素(去掉不必要的0元素)及COEFF和SCORE矩陣中相應的列 X = [2000 1000 3000 900 951900 990 3000 700 662000 900 3400 800 56]; % 行為樣本，列為特征量 XZ = zscore(X); % 數據標準化 [COEFF_2,SCORE,latent,tsquare] = pca(XZ); [m, n] = size(X); % 求X的行數和列數 result_2 = cell(m, 4); % 定義一個n+1行，4列的元胞數組 result_2(1,:) = {'特征值', '差值', '貢獻率', '累積貢獻率'}; explained = 100*latent/sum(latent);% 計算貢獻率 %result_2中第1列的第2行到最后一行存放的數據(latent)特征值 result_2(2:end,1) = num2cell(latent); %result1中第2列的第2行到倒數第2行存放的數據(latent的方差，特征值的方差) result_2(2:end-1,2) = num2cell(-diff(latent)); %result1中第3列和第4列的第2行到最后一行分別存放主成分的貢獻率和累積貢獻率 result_2(2:end,3:4) = num2cell([explained, cumsum(explained)]); disp(result_2); [b_1,b_2] = size(X); [b_3,b_4] = size(SCORE); disp(SCORE);

24. 主成分回歸

這是我論文的代碼，找不到原生的主成分回歸了

clc,clear,close all; % 渡辺論文代碼 m_1 % 文件的路徑 filename_1 = 'first_shuju.xlsx'; sheet_1 = ''; % Excle 的 工作表 名稱 range_1 = 'B43:Q58'; range_2 = 'B1:T1'; range_3 = 'G3:G18'; Year_0 = 2004; years(1) = Year_0; for j = 2:13years(j) = years(j-1)+1; end[num_1, ~,~] = xlsread(filename_1,sheet_1,range_1); [~, num_2,~] = xlsread(filename_1,sheet_1,range_2); [num_3, ~,~] = xlsread(filename_1,sheet_1,range_3); data = num_1(1:13,:);X = num_1(1:13,:); XZ = zscore(X); % 數據標準化 [COEFF_2,SCORE,latent,tsquare] = pca(XZ); [m, n] = size(X); % 求X的行數和列數 result_2 = cell(m, 4); % 定義一個n+1行，4列的元胞數組 result_2(1,:) = {'特征值', '差值', '貢獻率', '累積貢獻率'}; explained = 100*latent/sum(latent);% 計算貢獻率 %result_2中第1列的第2行到最后一行存放的數據(latent)特征值 result_2(2:end,1) = num2cell(latent); %result1中第2列的第2行到倒數第2行存放的數據(latent的方差，特征值的方差) result_2(2:end-1,2) = num2cell(-diff(latent)); %result1中第3列和第4列的第2行到最后一行分別存放主成分的貢獻率和累積貢獻率 result_2(2:end,3:4) = num2cell([explained, cumsum(explained)]); disp(result_2); disp(COEFF_2); xlswrite('result.xlsx',result_2); xlswrite('SCORE.xlsx',SCORE);y = num_3(1:13); % x,y 需要轉置才可以使用 num_1_2 = num_1(1:13,:); ppp = 0; for i=1:13for ii = 1:6pp(i,ii) = num_1_2(i,:) * COEFF_2(:,ii);end end SCORE = pp; X = [ones(size(y)) SCORE(:,1) SCORE(:,2) SCORE(:,3) SCORE(:,4) ...SCORE(:,5) SCORE(:,6)]; % 定義回歸方程[B,BINT,R,RINT,STATS] = regress(y,X); % 回歸 % [B,BINT,R,RINT,STATS] = REGRESS(Y,X,alpha) % b 為回歸系數向量 % BINT 回歸系數的估計區間 % R 殘差 % RINT 置信區間 % STATS 用于檢驗回歸模型的統計量。有4個數值：判定系數r2，F統計量觀測值，檢驗的p的值，誤差方差的估計 % r2越接近1，回歸方程越顯著；F>F1?α(k,n?k?1)時拒絕H0，F越大，回歸方程越顯著；p<α時拒絕H0 % ALPHA 顯著性水平（缺少時默認0.05） y2 = B(2).* SCORE(:,1)' + B(3).*SCORE(:,2)' + B(4).*SCORE(:,3)' + B(5).*SCORE(:,4)' ...+ B(6).*SCORE(:,5)' + B(7).*SCORE(:,6)' + B(1)'; %預測值 plot([2004:2019],num_3','ro-',years,y2,'b*'); legend('生態用水量', '主成分回歸'); xlabel('時間/年'); ylabel('城市生態用水量/億立方米'); disp('回歸分析統計量為'); disp([ 'R^2 : ',num2str(STATS(1))]); % 模型由 這么多原始數據 解釋，最優為 1 disp([ 'F : ',num2str(STATS(2))]); disp([ '檢驗 p : ',num2str(STATS(3))]); % 日常小于 0.05 最好 disp([ '誤差方差 : ',num2str(STATS(4))]); % 越小越好 xlswrite('zhuchengf_huigui.xlsx',B); hold on;num_1_2 = num_1(14:16,:); ppp = 0; for n = 1:3for ii = 1:6pp_1(n,ii) = num_1_2(n,:) * COEFF_2(:,ii);end end SCORE = pp_1; y3 = B(2).* SCORE(:,1)' + B(3).*SCORE(:,2)' + B(4).*SCORE(:,3)' + B(5).*SCORE(:,4)' ...+ B(6).*SCORE(:,5)' + B(7).*SCORE(:,6)' + B(1)'; %預測值 plot([2017:2019],y3,'m*'); line([2016.5,2016.5],[0,17],'linestyle','--'); axis([2004 2019 0 17]); legend('城市生態用水量-原始數據', '主成分回歸在數據集的結果','主成分回歸在測試集的結果'); xlswrite('save_4.xlsx',y3);

25. 分治算法

25.1. 快速排序

QuickSort.m

function a=QuickSort(a,leftIndex,rightIndex) % a是待排序序列 %leftIndex和rightIndex是開始的左右兩個邊界 %%----------------------------------------------------------%% % if leftIndex>rightIndex % break; % end if leftIndex<rightIndexi=leftIndex;j=rightIndex;temp=a(i);%選定的這個數挖掉，相當于挖坑while i<jwhile (i<j)&&(a(j)>=temp)%從右往左，找到第一個小于設定的數，j=j-1;enda(i)=a(j);%將找到的第一個小于設定的數填坑到最開始挖的坑里面去while (i<j)&&(a(i)<=temp)%從做到由，找到第一個大于選定的數i=i+1;enda(j)=a(i);%將找到的第一個大于選定的數填入上一步挖的坑里面去 % if i==j % a(j)=temp; % endenda(j)=temp;%最后，i=j，將選定的數再填到上一步挖的坑里面去a=QuickSort(a,leftIndex,j-1);%對左邊序列進行遞歸a=QuickSort(a,i+1,rightIndex);%對右邊序列進行遞歸 end end

調用方法

% 分治算法，渡邊筆記 % 分治算法是一種思路，把大問題歸納為小問題來解決 % 使用分治法，問題應該滿足 % 該問題的規模縮小到一定的程度就可以容易地解決 % 該問題可以分解為若干個規模較小的相同問題，即該問題具有最優子結構性質 % 利用該問題分解出的子問題的解可以合并為該問題的解 % 該問題所分解出的各個子問題是相互獨立的，即子問題之間不包含公共的子子問題 % 第一條特征是絕大多數問題都可以滿足的，因為問題的計算復雜性一般是隨著問題規模的增加而增加 % 第二條特征是應用分治法的前提它也是大多數問題可以滿足的，此特征反映了遞歸思想的應用 % 第三條特征是關鍵，能否利用分治法完全取決于問題是否具有第三條特征，如果具備了第一條和第二條特征，而不具備第三條特征，則可以考慮用貪心法或動態規劃法。 % 第四條特征涉及到分治法的效率，如果各子問題是不獨立的則分治法要做許多不必要的工作， % 重復地解公共的子問題，此時雖然可用分治法，但一般用動態規劃法較好 clc,clear;close all; a=[72 57 88 60 42 83 73 48 85]; % 輸入排序數leftIndex=1; [~, rightIndex]=size(a); disp(['未排序的序列為：',num2str(a)]) a=QuickSort(a,leftIndex,rightIndex); disp(['未排序的序列為：',num2str(a)])

25.2. 最近點對

cloest.m

%分治求最近點對 function [d,x1,x2] = cloest(S,low,high) P=[]; if(high - low == 1)[d,x1,x2] = deal(pdist2(S(low,:),S(high,:)),S(low,:),S(high,:)); elseif(high - low == 2)d1 = pdist2(S(low,:),S(low + 1,:));d2 = pdist2(S(low + 1,:),S(high,:));d3 = pdist2(S(low,:),S(high,:));if((d1 < d2) && (d1 < d3))[d,x1,x2] = deal(d1,S(low,:),S(low + 1,:));elseif(d2 < d3)[d,x1,x2] = deal(d2,S(low+1,:),S(high,:));else[d,x1,x2] = deal(d3,S(low,:),S(high,:));end elsemid = floor((low+high)/2);[d1,x11,x12] = cloest(S,low,mid);[d2,x21,x22] = cloest(S,mid+1,high);if(d1 <= d2)[d,x1,x2] = deal(d1,x11,x12);else[d,x1,x2] = deal(d2,x21,x22);endindex = 0;for i = mid:-1:lowif(S(mid,1) - S(i,1) >= d)break;endindex = index + 1;P(index,:) = S(i,:);endfor i = mid+1:highif(S(i,1) - S(mid,1) >= d)break;endindex = index + 1;P(index,:) = S(i,:);endsortrows(P,2);for i = 1:indexfor j = i+1:indexif(P(j,2) - P(i,2) >= d)break;elsed3 = pdist2(P(i,:),P(j,:));if(d3 < d)[d,x1,x2] = deal(d3,P(i,:),P(j,:));endendendend end

調用方法

% 分治法求最近點，渡邊筆記 clear;clc n = 20; % 隨機生成20個點 A=rand(n,2)*10; % 將20個點按橫坐標升序排列 A = sortrows(A,1); % 分治法求隨機點的最近點對 [mindist1,y1,y2] = cloest(A,1,n); % 使用紅色的點標記隨機點 x=A(:,1); y=A(:,2); for i = 1:nplot(x,y,'r.'); hold ontext(A(i,1),A(i,2),num2str(i));hold on end %使用綠色十字標記分治法的最近點對 plot(y1(1),y1(2),'go', 'linewidth', 2);hold on plot(y2(1),y2(2),'go', 'linewidth', 2);hold on

26. 雜項

26.1. 有向圖

clc,clear,close all; s = [1 1 1 2 2 3 3 4 5 5 6 7]; t = [2 4 8 3 7 4 6 5 6 8 7 8]; weights = [10 10 1 10 1 10 1 1 12 12 12 12]; % names = {'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H'}; % G = digraph(s,t,weights,names); G = digraph(s,t,weights); plot(G,'Layout','force','EdgeLabel',G.Edges.Weight)

26.2. 數據擬合

懶得修改了，反正方法也就那樣

clc,clear % 產生數據 lamuda_0 = xlsread('shuju.xlsx','','B3:U23'); lamuda_test = xlsread('shuju.xlsx','','A3:A23'); lamuda_test = lamuda_test'; x = xlsread('shuju.xlsx','','V3:V23'); % 擬合 for j = 1:21lamuda = 355 - lamuda_test(j);p = fittype('x * exp( S*lamuda )');f = fit(x,lamuda_0(j,:)',p); % plot(f,x,lamuda_0(j,:)');S(j) = f.S; end

26.3. 標準化

% 標準化處理，渡邊筆記 % 歸一化的依據非常簡單，不同變量往往量綱不同 % 歸一化可以消除量綱對最終結果的影響，使不同變量具有可比性 % 如兩個人體重差10KG，身高差0.02M % 在衡量兩個人的差別時體重的差距會把身高的差距完全掩蓋，歸一化之后就不會有這樣的問題 % 數據歸一化應該針對對于的 y，而不是針對每條數據，針對每條數據是完全沒有意義的，因為只是等比例縮放，對之后的分類沒有任何作用 clc,clear;close all; x = [1 2 31 2 3]; % Min-max 極大極小標準化 % 當有新數據加入時，可能導致xmax和xmin的變化，需要重新計算 X_min_max = mapminmax(x, 0, 1); % 默認標準化區間為 0-1 disp('極大極小標準化為 :') disp(X_min_max);% z-score 標準化 % 其結果沒有任何意義，只能用于比較 % X_z_score = zscore(x); % disp('z-score 標準化 :') % disp(X_z_score);

26.4. excel 數據轉置

clc,clear % 省會城市 順序表 % 北京 天津 石家莊 太原 呼和浩特 沈陽 長春 哈爾濱 上海 南京 % 杭州 合肥 福州 南昌 濟南 鄭州 武漢 長沙 廣州 南寧 % ?？?重慶 成都 貴陽 昆明 拉薩 西安 蘭州 西寧 銀川 % 烏魯木齊 filename_lujin_1 = 'C:\Users\99791\Desktop\生態用水\數據\國家統計局\分省\'; % 文件的前面的路徑名稱 filename_name_1 = '分省_城市綠地面積.xls'; % 文件的名字 filename_1 = [filename_lujin_1,filename_name_1]; sheet_1 = ''; % Excle 的 工作表 名稱 range_1 = 'B5:Q35'; % 讀取的范圍 [num_1, tex_1, raw_1] = xlsread(filename_1,sheet_1,range_1); range_1_1 = 'A2'; [num_2, tex_2, raw_2] = xlsread(filename_1,sheet_1,range_1_1); num_1 = num_1'; dubian_bianliang_1 = num_1; % 同上 filename_lujin_2 = 'C:\Users\99791\Desktop\生態用水\數據\'; filename_name_2 = 'dubian_shuju.xlsx'; filename_2 = [filename_lujin_2,filename_name_2]; sheet_2 = ''; t_1 = 21; t_2 = 'W'; % 改變寫入 列 range_2 = [t_2,'3:',t_2,'18']; for i = 1:31dubian_bianliang_2 = dubian_bianliang_1(:,i);xlswrite(filename_2,dubian_bianliang_2,sheet_2,range_2)range_2 = [t_2,num2str(3+i*t_1),':',t_2,num2str(18+i*t_1)]; end range_2_1 = [t_2,'1']; xlswrite(filename_2,tex_2,sheet_2,range_2_1)

26.5. 線性規劃

lingo 還是比 matlab 強，下面的方法好久不用了

% function 函數解多元方程 最小值 % function 基本語法 % [x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) % x的返回值是決策向量x的取值，fval的返回值是目標函數f(x)的取值 % fun是用 M 文件定義的函數f(x),代表了(非)線性目標函數 % x0是x的初始值 % A,b,Aeq,beq 定義了線性約束 ,如果沒有線性約束，則A=[],b=[],Aeq=[],beq=[] % lb和ub是變量x的下界和上界，如果下界和上界沒有約束，則lb=[],ub=[],也可以寫成lb的各分量都為 -inf,ub的各分量都為inf % nonlcon是用 M 文件定義的非線性向量函數約束 % options 定義了優化參數，不填寫表示使用 Matla b默認的參數設置 clc,clear;close all; options = optimset('Algorithm','active-set'); % 可以不設置 [X,Y,EXITFLAG] = fmincon(@(x)( x(1)^3 + x(2)^3 + x(3)^3 +x(4)^3 ),... [0 0 0 0],[],[],[],[],... [-2 -2 -2 -2],[-1 -1 -1 -1]... % 設置范圍 ,[],options); disp(['函數的滿足最小值的 X 解為 : ',num2str(X)]); disp(['函數的最小值的 Y 解為 : ',num2str(Y)]);

總結

以上是生活随笔為你收集整理的数学建模方法总结(matlab)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。