数学建模方法总结(matlab)
生活随笔
收集整理的這篇文章主要介紹了
数学建模方法总结(matlab)
小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.
1. 前言
在這篇文章中，我會把我所接觸的數(shù)學(xué)建模的知識的代碼分享給大家，有的是我自己寫的，更多的是網(wǎng)上借鑒并修改為可執(zhí)行可用的代碼
需要說明的是，我不太了解其中的數(shù)學(xué)實現(xiàn)的具體方法和原理，我只分享源碼和作用條件以及使用的經(jīng)驗說明（更詳細(xì)見注釋），網(wǎng)上有不少文章對某些方法講得非常透徹，因此我沒必要贅述
各位只需要忘記代碼或使用條件時從我這里 copy 就好了
2. 灰色預(yù)測模型
這是用一種在數(shù)學(xué)上對數(shù)據(jù)進行累加和累減以取巧性地削弱預(yù)測數(shù)據(jù)與原始數(shù)據(jù)的關(guān)聯(lián)性
最好用于短期預(yù)測，只能用于遞增序列
2.1. GM（1，1）
我也不懂為什么叫 GM(1,1)，但是用得最多的就是它，又稱一階灰色方程
clc,clear;close all y = [0.6 1.1 1.6] yuceshu = 3; % 本程序主要用來計算根據(jù)灰色理論建立的模型的預(yù)測值 % 應(yīng)用的數(shù)學(xué)模型是 GM(1,1) % 原始數(shù)據(jù)的處理方法是一次累加 % 渡辺曜改進了這個代碼,通常GM預(yù)測第一個數(shù)是不夠優(yōu)秀的 % y 是應(yīng)該輸入的數(shù)組，yuceshu是你想預(yù)測的后幾個數(shù)據(jù) 我們一般就預(yù)測兩個 % 示例 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擬合的第一個數(shù)據(jù) 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('系統(tǒng)預(yù)測精度好') else if c<0.5disp('系統(tǒng)預(yù)測精度合格')else if c<0.65disp('系統(tǒng)預(yù)測精度勉強')elsedisp('系統(tǒng)預(yù)測精度不合格')endend enddisp(['下個擬合值為 ',num2str(ys(n+1))]); for i=2:yuceshudisp(['再下個擬合值為',num2str(ys(n+i))]); end2.2. 分?jǐn)?shù)階 GM(1,1)
分?jǐn)?shù)階灰色方程的預(yù)測效果優(yōu)于原始灰色模型，我只找到了 python 版本的可運行代碼
首先，將 GM(1,1)打包成一個類
GM.py
##渡辺筆記 ##傳統(tǒng)的灰色模型 import numpy as np import xlrd import xlwt import datetime class Grey_model(object):def __init__(self,input_value):##初始化過程，先建立好變量空間，即累加序列，背景值序列，B和Y矩陣，擬合序列，預(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]##計算參數(shù)矩陣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##下面把矩陣格式轉(zhuǎn)化為數(shù)組格式，再轉(zhuǎn)化為列表格式self.u_matrix_array=np.array(self.u_matrix_values.T)[0]self.u_matrix_value=self.u_matrix_array.tolist()##計算預(yù)測值def predict(self,number_of_forecast):self.u_matrix()self.predicted_accumulation_value=[]##使用了float(%.2f% 來調(diào)整輸出值的小數(shù)的個數(shù)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表示分?jǐn)?shù)階的階數(shù) 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('預(yù)測值')print(A.predicted_value)print('\nMAPE值')print(A.MAPE)## xlwt_data = A.predicted_valuedef save_excel(target_list, output_file_name):將數(shù)據(jù)寫入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)在這個類的基礎(chǔ)上，運用分?jǐn)?shù)階進行優(yōu)化
FGM.py
##渡辺筆記 ##分?jǐn)?shù)階灰色模型 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表示分?jǐn)?shù)階的階數(shù) 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)##更新預(yù)測值的還原方式 A.test_error()##誤差用MAPE值來衡量 print('預(yù)測值') print(A.predicted_value) print('\nMAPE值') print(A.MAPE)def save_excel(target_list, output_file_name):將數(shù)據(jù)寫入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不能解決多條最短路徑和負(fù)權(quán)回路的問題 % 解稠密圖效果最佳，邊權(quán)可正可負(fù)，是雙源解法 % 對于稠密圖，效率要高于執(zhí)行|V|次Dijkstra算法，也要高于執(zhí)行|V|次SPFA算法。 % 缺點：時間復(fù)雜度比較高，不適合計算大量數(shù)據(jù)，注意 % 輸入示例 % 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為指定起始結(jié)點，stop為指定終止結(jié)點 D = W; %最短距離矩陣賦初值 n = length(D); %n為結(jié)點個數(shù) 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('參數(shù)個數(shù)輸入有誤！', 'Warning!') elsedis = D(start, stop); %指定兩結(jié)點間的最短距離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; %指定兩結(jié)點之間的最短路徑 end3.2. Dijkstra
Dijkstra 運算速度快，但是用不了負(fù)權(quán)圖
% 這是 dijkstra 算法，以下注釋為渡邊所加 % 不能用于負(fù)權(quán)圖，但是由于時間復(fù)雜度低，適合在大數(shù)據(jù)中運算 % 是單源解法 % 輸入示例 % 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 權(quán)值矩陣 st 搜索的起點 e 搜索的終點 n=length(W);%節(jié)點數(shù) D = W(st,:); visit= ones(1:n); visit(st)=0; parent = zeros(1,n);%記錄每個節(jié)點的上一個節(jié)點path =[];for i=1:n-1temp = [];%從起點出發(fā)，找最短距離的下一個點，每次不會重復(fù)原來的軌跡，設(shè)置visit判斷節(jié)點是否訪問for j=1:nif visit(j)temp =[temp D(j)];elsetemp =[temp inf];endend[value,index] = min(temp);visit(index) = 0;%更新 如果經(jīng)過index節(jié)點，從起點到每個節(jié)點的路徑長度更小，則更新，記錄前趨節(jié)點，方便后面回溯循跡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; end3.3. matlab 自帶最短路徑
clc,clear,close all; % 程序基于但不基于 Dijkstra，解決了負(fù)權(quán)問題 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. 層次分析法
感覺這個方法就像我們土木工程所謂的“拍腦袋”，先拍腦袋得出一個比較矩陣，待反復(fù)測試符合規(guī)范后，就可以拿去比較了
% 本代碼只需要在提示輸入處，輸入構(gòu)成的比較矩陣 就會輸出各指標(biāo)的權(quán)重值。 % 例子： 選拔干部考慮5個條件：品德x1，才能x2，資歷x3，年齡x4，群眾關(guān)系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的值就是由決策人決定了，對應(yīng)值也就是決策人認(rèn)為的重要性了 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); %獲取指標(biāo)個數(shù) 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); %對應(yīng)特征向量 CI=(B-n)/(n-1); %計算一致性檢驗指標(biāo)CI CR=CI/RI(1,n); if CR<0.10disp('CI=');disp(CI);disp('CR=');disp(CR);disp('對比矩陣A通過一致性檢驗，各向量權(quán)重向量Q為：');Q=zeros(n,1);for i=1:nQ(i,1)=C(i,1)/sum(C(:,1)); %特征向量標(biāo)準(zhǔn)化endQ' %輸出權(quán)重向量 elsedisp('對比矩陣A未通過一致性檢驗，需對對比矩陣A重新構(gòu)造'); end sc = Q';5. 聚類分析
這類方法對于我的專業(yè)毫無實際用處，我并不需要知道管道應(yīng)該分成幾類，因為規(guī)范早就規(guī)定好了的
聚類分析方法豐富，要根據(jù)實際情況謹(jǐn)慎選取
5.1. k-mean
% k_mean 算法，渡邊筆記 % 對處理大型數(shù)據(jù)集，該算法保持可收縮性，高效性； % 當(dāng)簇接近正態(tài)分布時，效果最好； % 若簇中含有異常點，將導(dǎo)致均值偏離嚴(yán)重；（即對噪聲、孤立點敏感）； % 不適用于發(fā)現(xiàn)非凸形狀的簇，或大小差別大的簇 % 中k是事先給定的，這個k值難以估計，很多時候并不知道數(shù)據(jù)集應(yīng)該分成多少個類別最合適 clear; clc;data = [11 2 02 2 24 3 3]; % 輸入數(shù)據(jù) k = 3; % 分類數(shù)X = mapminmax(data',0,1)'; % 按列最小最大規(guī)范化到[0,1] T = kmeans(X,k); % 直接調(diào)用kmeans函數(shù) for i = 1:kkong = find(T == i);fprintf('第 %d 類 :',i);disp(data(kong)') end5.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)'; % 按行最小最大規(guī)范化到[0,1] Y = pdist(X); % 計算矩陣X中樣本兩兩之間的距離，但得到的Y是個行向量 D = squareform(Y); % 將行向量的Y轉(zhuǎn)換成方陣，方便觀察兩點距離(方陣的對角線元素都是0) Z = linkage(Y); % 產(chǎn)生層次聚類樹，默認(rèn)采用最近距離作為類間距離的計算公式 dendrogram(Z); % 圖示層次聚類樹 title('層次聚類法');ylabel('標(biāo)準(zhǔn)距離');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); % 按列最小最大規(guī)范化到[0,1]T1=clusterdata(X,0.2); % 如果0<cutoff<2，則當(dāng)不一致系數(shù)大于cutoff時，分到不同類(簇)中 T2=clusterdata(X,3); % 如果cutoff是一個≥2的整數(shù)，則形成的不同類別數(shù)為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)') % end6. 最小生成樹
好像是個示例，太久遠(yuǎn)忘了干什么用的了
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];% 使用加權(quán)邊創(chuàng)建并繪制一個立方體圖 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 包含的節(jié)點與 G 相同，但包含的邊僅為后者的子集 % figure; [T,pred] = minspantree(G); % p = plot(G,'EdgeLabel',G.Edges.Weight); highlight(p,T) % 高亮顯示7. 誤差檢驗
7.1. 一般檢驗
看著用吧，我覺得沒什么意義，反正計算機計算出的誤差都挺小的，最多能在論文中湊字?jǐn)?shù)
% 渡辺筆記 % 各類檢驗，除決定系數(shù)是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); % 決定系數(shù)（R2-R-Square） r2 = 1 - (sum((YPred - YReal).^2) / sum((YReal - mean(YReal)).^2)); disp(r2)7.2. T 檢驗
% 渡邊筆記 t檢驗 % 待檢驗的數(shù)據(jù)應(yīng)該近服從正態(tài)分布 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均為行向量（維度必須相同），各表示一組數(shù)據(jù) % ALPHA為可選參數(shù)，表示設(shè)置一個值作為t檢驗執(zhí)行的顯著性水平 % 在不設(shè)置ALPHA的情況下默認(rèn)ALPHA為0.05，即計算x和y在5%的顯著性水平下是否來自同一分布(假設(shè)是否被接受） % H = 0，P > 0.05，表明零假設(shè)被接受，即x，y在統(tǒng)計上可看做來自同一分布的數(shù)據(jù) % H = 1，P < 0.05，表明零假設(shè)被拒絕，即x，y在統(tǒng)計上認(rèn)為是來自不同分布的數(shù)據(jù)，即有區(qū)分度 disp(['H = ',num2str(H)]); disp(['P = ',num2str(P)]); if H == 0disp('x，y在統(tǒng)計上可看做來自同一分布的數(shù)據(jù)') elsedisp('x，y在統(tǒng)計上可看做來自不同分布的數(shù)據(jù)') end7.3. 回歸檢驗
clc,clear;close all; % 回歸檢驗；渡邊筆記 x = [0 1 2 3 4];y = [1 1 4 10 16];x = x'; y = y'; % x,y 需要轉(zhuǎn)置才可以使用X = [ones(size(y)) x]; % 定義回歸方程[B,BINT,R,RINT,STATS] = regress(y,X); % 回歸 % [B,BINT,R,RINT,STATS] = REGRESS(Y,X,alpha) % B 為回歸系數(shù)向量 % BINT 回歸系數(shù)的估計區(qū)間 % R 殘差 % RINT 置信區(qū)間 % STATS 用于檢驗回歸模型的統(tǒng)計量。有4個數(shù)值：判定系數(shù)r2，F統(tǒng)計量觀測值，檢驗的p的值，誤差方差的估計 % r2越接近1，回歸方程越顯著；F>F1?α(k,n?k?1)時拒絕H0，F越大，回歸方程越顯著；p<α?xí)r拒絕H0 % ALPHA 顯著性水平（缺少時默認(rèn)0.05） y2 = B(2).* x + B(1); %預(yù)測值 plot(x', y' ,'+', x, y2); %回歸效果圖 disp('回歸分析統(tǒng)計量為'); disp([ 'R^2 : ',num2str(STATS(1))]); % 模型由 這么多原始數(shù)據(jù) 解釋，最優(yōu)為 1 disp([ 'F : ',num2str(STATS(2))]); disp([ '檢驗 p : ',num2str(STATS(3))]); % 日常小于 0.05 最好 disp([ '誤差方差 : ',num2str(STATS(4))]); % 越小越好7.4. 交叉檢驗
% 交叉檢驗，渡邊筆記 % 提高數(shù)據(jù)的利用率 % 有效的避免過學(xué)習(xí)以及欠學(xué)習(xí)狀態(tài)的發(fā)生,最后得到的結(jié)果也比較具有說服性 % 常用的 K 值有 3,6,10 等 clear,clc;close all;% 導(dǎo)入數(shù)據(jù) data = reshape(1:144,12,12); % 每行的是這個樣本屬性的數(shù)據(jù)，最后一個數(shù)據(jù)是標(biāo)簽 [data_r, data_c] = size(data); % 將數(shù)據(jù)樣本隨機分割為 n 部分，做 n 折交叉檢驗 n = 5; indices = crossvalind('Kfold', data_r, n); for i = 1 : n% 獲取第i份測試數(shù)據(jù)的索引邏輯值test = (indices == i);% 取反，獲取第i份訓(xùn)練數(shù)據(jù)的索引邏輯值train = ~test;% 測試用的數(shù)據(jù)test_data = data(test, 1 : data_c - 1);test_label = data(test, data_c);% 訓(xùn)練用的數(shù)據(jù)train_data = data(train, 1 : data_c - 1);train_label = data(train, data_c);% 下面放入檢驗數(shù)據(jù)的代碼end8. 差分
% 數(shù)據(jù)差分，渡邊筆記 % 差分就是后一個減去前一個，使得數(shù)據(jù)平穩(wěn)化 % n 階就是做 n 次差分，具體直接看例子就好 clc,clear;close all; Y = [0 5 15 30 50 75 105]; disp(['未差分時為 : ',num2str(Y)]) Y_1= diff(Y, 1); disp(['一階差分結(jié)果為 : ',num2str(Y_1)]) Y_2 = diff(Y, 2); disp(['二階差分結(jié)果為 : ',num2str(Y_2)]) Y_3 = diff(Y, 3); disp(['三階差分結(jié)果為 : ',num2str(Y_3)])9. 立法白噪音
% 立法數(shù)白噪聲檢驗，渡邊筆記 % 如果易知原始數(shù)據(jù)不是平穩(wěn)的，則不能做隨機性檢驗 % 接下來要求差分,目的: 變成平穩(wěn)的數(shù)據(jù) % p 如果比 0.05 小就不是白噪聲序列，可以使用時間序列 % 某種現(xiàn)象的指標(biāo)數(shù)值按照時間順序排列而成的數(shù)值序列 % 因為時間序列是某個指標(biāo)數(shù)值長期變化的數(shù)值表現(xiàn) % 所以時間序列數(shù)值變化背后必然蘊含著數(shù)值變換的規(guī)律性 % 這些規(guī)律性就是時間序列分析的切入點 % 如果原數(shù)據(jù)平穩(wěn)，但是沒通過，可以直接差分 clc,clear;close all;x = []; % 時間，一般做題就是順序時間排列 yanchi=[6 12 18]; % 通常是做6 12 18 24步延遲，這個數(shù)據(jù)的選擇上限請根據(jù)報錯來調(diào)整 y = [970279 1259308 1127571 1163959 1169540 1076938 991350 953275 951508 904434 889381 864015 836236 ]';; [h1] = adftest(y); %檢驗是否平穩(wěn) if h1 == 1disp('數(shù)據(jù)是平穩(wěn)的');y_1 = y; elsedisp('數(shù)據(jù)是不平穩(wěn)的');i = 1;while 1y_1 = diff(y,i); % 在這里對數(shù)據(jù)進行 n 階差分[h1,p1,adf,ljz] = adftest(y_1); %檢驗是否平穩(wěn)if h1 == 1disp(['差分后是平穩(wěn)的，做了 ',num2str(i),' 階差分']);subplot(2,2,3)plot(y_1); % 一階差分?jǐn)?shù)據(jù)時序圖title('一階差分?jǐn)?shù)據(jù)時序圖')subplot(2,2,4)autocorr(y_1); % 一階差分?jǐn)?shù)據(jù)的自相關(guān)系數(shù)圖title('一階差分?jǐn)?shù)據(jù)自相關(guān)系數(shù)圖');breakendi = i + 1;end end % 隨時間的變化值 subplot(2,2,1) plot(y); % 原始數(shù)據(jù)時序圖 title('原始數(shù)據(jù)時序圖') subplot(2,2,2) autocorr(y); % 原始數(shù)據(jù)的自相關(guān)系數(shù)圖 title('原始數(shù)據(jù)自相關(guān)系圖像')[H,pValue,Qstat,CriticalValue]=lbqtest(y_1,'lags',yanchi); % H.結(jié)果，pValue.p值, Qstat.卡方統(tǒng)計量 fprintf('%15s%15s%15s','延遲階數(shù)','卡方統(tǒng)計量','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); % 在這里對數(shù)據(jù)進行 n 階差分[h1,p1,adf,ljz] = adftest(y_1); %檢驗是否平穩(wěn)if h1 == 1disp(['差分后是平穩(wěn)的，做了 ',num2str(i),' 階差分']);subplot(2,2,3)plot(y_1); % 一階差分?jǐn)?shù)據(jù)時序圖title('一階差分?jǐn)?shù)據(jù)時序圖')subplot(2,2,4)autocorr(y_1); % 一階差分?jǐn)?shù)據(jù)的自相關(guān)系數(shù)圖title('一階差分?jǐn)?shù)據(jù)自相關(guān)系數(shù)圖');breakendi = i + 1;end end% 再來一次[H,pValue,Qstat,CriticalValue]=lbqtest(y_1,'lags',yanchi); % H.結(jié)果，pValue.p值, Qstat.卡方統(tǒng)計量 fprintf('%15s%15s%15s','延遲階數(shù)','卡方統(tǒng)計量','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('數(shù)據(jù)通過立法白噪聲檢驗'); elsedisp('數(shù)據(jù)沒通過立法白噪聲檢驗'); end10. 熵權(quán)法
% 渡邊筆記 熵權(quán)法 clc;clear;close all; % 熵權(quán)法的思想是：變量數(shù)值變化越大，變異程度越大，則其權(quán)重應(yīng)該更大；反之權(quán)重則越小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); % 人為修權(quán)，1代表不修改計算后的指標(biāo)權(quán)重 for i=1:nx(:,i)=(x(:,i)-min(x(:,i)))/(max(x(:,i))-min(x(:,i)))+1; % 對原始數(shù)據(jù)進行非負(fù)數(shù)化、歸一化處理，值介于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);% 指標(biāo)權(quán)重計算 end for j=1:nw(j)=w(j)*lamda(j)/sum(w.*lamda);% 修改指標(biāo)權(quán)重 end for i=1:mscore(i,1)=sum(x(i,:).*w);% 計算綜合分?jǐn)?shù) end disp('各指標(biāo)權(quán)重為：') disp(w) disp('各項綜合分?jǐn)?shù)為：') disp(score) Out = mean (score,1)11. 數(shù)據(jù)修復(fù)
% 判斷缺失值和異常值并修復(fù)，順便光滑噪音，渡邊筆記 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('原始數(shù)據(jù)');%% 判斷數(shù)據(jù)中是否存在缺失值 if sum(isnan(testdata(:)))disp('數(shù)據(jù)存在缺失值'); elsedisp('數(shù)據(jù)不存在缺失值'); end%% 判斷數(shù)據(jù)中是否存在異常值 % 1.mean 三倍標(biāo)準(zhǔn)差法 2.median 離群值法 3.quartiles 非正態(tài)的離群值法 % 4.grubbs 正態(tài)的離群值法 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('數(shù)據(jù)存在異常值'); elsedisp('數(shù)據(jù)不存在異常值'); end%% 對異常值賦空值 F = find(yi_chang == 1); y(F) = NaN; % 令數(shù)據(jù)點缺失 testdata = [x' y'];subplot(2,2,2); plot(testdata(:,1),testdata(:,2)); title('去除差異值');%% 對數(shù)據(jù)進行補全 % 數(shù)據(jù)補全方法選擇 % 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('數(shù)據(jù)補全結(jié)果');%% 進行數(shù)據(jù)平滑處理 % 濾波器選擇 1.Savitzky-golay 2.rlowess 3.rloess choice_3 = 2; lvboqi = char('Savitzky-golay', 'rlowess', 'pchip', 'rloess'); % 通過求 n 元素移動窗口的中位數(shù)，來對數(shù)據(jù)進行平滑處理 windows = 8; testdata_2 = smoothdata(testdata_1(:,2),strtrim(lvboqi(choice_3,:)),windows) ;subplot(2,2,4); plot(x,testdata_2) title('數(shù)據(jù)平滑結(jié)果');12. 相關(guān)性分析
clc,clear,close all % 相關(guān)性分析，渡邊筆記 % 列為指標(biāo)，行為數(shù)據(jù)data = [1 22 93 4 ]; % Pearson相關(guān)系數(shù) r1 = corr(data,'type','Pearson'); disp(r1); % Kendall相關(guān)系數(shù) r2 = corr(data,'type','Kendall'); disp(r2); % Spearman相關(guān)系數(shù) r3 = corr(data,'type','Spearman'); disp(r3);13. 馬氏鏈
13.1. 一般馬氏鏈
% 馬氏鏈模型，渡邊筆記 % 系統(tǒng)在每個時期所處的狀態(tài)是隨機的 % 從一時期到下時期的狀態(tài)按一定概率轉(zhuǎn)移 % 下時期狀態(tài)只取決于本時期狀態(tài)和轉(zhuǎn)移概率。即已知現(xiàn)在，將來與過去無關(guān)（無后效性） % 過去不影響未來的判斷，是馬氏鏈的本質(zhì) % P_0 是初始分布，P_n 是轉(zhuǎ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]; % 輸入統(tǒng)計矩陣 P_0 = [0.4 0.3 0.4]; % 初始概率% 求狀態(tài)轉(zhuǎn)移矩陣 [m,~] = size(a); ni=sum(a,2); %計算矩陣f的行和 for i = 1:mfor j = 1:mZT(i,j) = a(i,j)./ni(i); %求狀態(tài)轉(zhuǎn)移的頻率end end disp('轉(zhuǎn)移矩陣為 : '); 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 ]; % 輸入統(tǒng)計矩陣 R = [9 33 -7]; % 輸入利潤矩陣% 對統(tǒng)計矩陣，求狀態(tài)轉(zhuǎn)移矩陣 [m,~] = size(a); ni=sum(a,2); %計算矩陣f的行和 for i = 1:mfor j = 1:mZT(i,j) = a(i,j)./ni(i); %求狀態(tài)轉(zhuǎn)移的頻率end end disp('轉(zhuǎn)移矩陣為 : '); 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(['處于狀態(tài) ',num2str(i),' 的 ',num2str(n),' 期利潤為 ',num2str(V_Sum(i,:))])elsedisp(['處于狀態(tài) ',num2str(i),' 的 ',num2str(n),' 期利潤為 ',num2str(Sum(i,:))])end end% 帶利潤的馬氏鏈一般不存在平衡 % 因為錢是賺不完的14. 蒙特卡羅
% 渡邊筆記 % Monte_Carlo 是一種方法的概述，不是什么特定的算法 % 主要思想是利用大量隨機數(shù)來實現(xiàn)真實值的擬合 %% 例1 計算 0-3 上 x^2 的積分，正解應(yīng)該是 9 N=500; x=unifrnd(0,3,N,1); % 設(shè)定隨機數(shù)及其終止上限 M = mean(3*x.^2);disp(['M = ',num2str(M)]); % 將隨機數(shù)代入原式求平均值，特別的，蒙特卡洛需要式乘以積分區(qū)間,下同 %% 例2 求相交面積 %給定曲線y =2 – x2 和曲線y3 = x2，曲線的交點為：P1(–1，1)、P2(1，1),函數(shù)圖像如下 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') %曲線圍成平面有限區(qū)域，用蒙特卡羅方法計算區(qū)域面積。 P=rand(1000,2); %隨機產(chǎn)生1000個點 %定義x y 的范圍 x=2*P(:,1)-1; % 這一步是巧妙設(shè)置在(-1,2) y=2*P(:,2); II=find(y<=2-x.^2&y.^3>=x.^2); %找出在函數(shù)范圍的數(shù) %計算索引的長度 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; %總的實驗次數(shù) m = 0; %落在圓中點的次數(shù) %循環(huán)實驗 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) %加權(quán)平均型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('模型賦值不當(dāng)');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. 運用數(shù)學(xué)方法
clc,clear % 模擬退火(求最小) temperature = 100; % 初始溫度 final_temperature = 1; % 最低溫度 time_temperature = 1; % 用于計算步數(shù) time_tuihuo = 10; % 退火步長 cooling_rate = 0.95; % 降溫步數(shù) % 初始化隨機數(shù)生成器，以使結(jié)果具備可重復(fù)性 % rng(0,'twister'); % 生成范圍a-b的隨機數(shù) a = -10; b = 10; % rang_math = (b-a).*rand(1000,1) + a; rang_math = (b-a).*rand + a; f = @(x)(...x.^2 ... % 定義函數(shù)，求最小值); current_old = f(rang_math);while final_temperature < temperaturerang_math = (b-a).*rand + a; % 隨機數(shù)current_new = f(rang_math); % 生成新解diff = current_new - current_old;% Metropolis 準(zhǔn)則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])); % 原函數(shù)圖像 hold on; plot(route,current_old,'ko', 'linewidth', 2); % 作解的圖像16.2. 運用 matlab 自帶工具
######## 2.1. 一元示例
% 渡邊自己手打的 % 模擬退火解一元函數(shù) clc,clear;close all;% 設(shè)置區(qū)間 x = -10:1:10;f = @(x)(... x.^4 ... % 定義函數(shù)，求最小值 );% f = @(x)f_xy(x(1),x(2)); x0 = rand(1,1); % 退火開始點 lb = []; ub = []; % 退火實施范圍，可以不設(shè)置 % 實時觀測 options = saoptimset('MaxIter',20,... % 迭代次數(shù) '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('最優(yōu)解為 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. 二元示例
% 渡邊自己手打的 % 模擬退火解二元函數(shù) clc,clear;close all;% 設(shè)置區(qū)間 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 ... % 定義函數(shù)，求最小值 ); 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 = []; % 退火實施范圍，可以不設(shè)置 % 實時觀測 options = saoptimset('MaxIter',20,... % 迭代次數(shù) '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('最優(yōu)解為 x = %.10f, y = %.10f\n', jie(1),jie(2) ); fprintf('最優(yōu)值為 z = %.10f\n',fval); figure; surf(x,y,fun); hold on; plot3(jie(1),jie(2),fval,'ko', 'linewidth', 3);17. 神經(jīng)網(wǎng)絡(luò)
17.1. BP
% 渡辺筆記，BP 神經(jīng)網(wǎng)絡(luò) clc,clear,close all; %輸入數(shù)據(jù)矩陣,每行為對應(yīng)的屬性 p = [1 1 12 2 23 3 3]; %目標(biāo)（輸出）數(shù)據(jù)矩陣 t = [6 6 6];%對訓(xùn)練集中的輸入數(shù)據(jù)矩陣和目標(biāo)數(shù)據(jù)矩陣進行歸一化處理 [pn, inputStr] = mapminmax(p); [tn, outputStr] = mapminmax(t);%建立BP神經(jīng)網(wǎng)絡(luò) net = newff(pn, tn, [3 7 2], {'purelin', 'logsig', 'purelin'}); %每10輪回顯示一次結(jié)果 net.trainParam.show = 10; %最大訓(xùn)練次數(shù) net.trainParam.epochs = 5000; %網(wǎng)絡(luò)的學(xué)習(xí)速率 net.trainParam.lr = 0.05; %訓(xùn)練網(wǎng)絡(luò)所要達(dá)到的目標(biāo)誤差 net.trainParam.goal = 0.65 * 10^(-3); %網(wǎng)絡(luò)誤差如果連續(xù)6次迭代都沒變化，則matlab會默認(rèn)終止訓(xùn)練。為了讓程序繼續(xù)運行，用以下命令取消這條設(shè)置 net.divideFcn = ''; %開始訓(xùn)練網(wǎng)絡(luò) net = train(net, pn, tn); %使用訓(xùn)練好的網(wǎng)絡(luò)，基于訓(xùn)練集的數(shù)據(jù)對BP網(wǎng)絡(luò)進行仿真得到網(wǎng)絡(luò)輸出結(jié)果 %(因為輸入樣本（訓(xùn)練集）容量較少，否則一般必須用新鮮數(shù)據(jù)進行仿真測試) answer = sim(net, pn); %反歸一化 answer1 = mapminmax('reverse', answer, outputStr);17.2. 灰色
clc,clear,close all; % 渡邊筆記 灰色神經(jīng)網(wǎng)絡(luò)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);%灰色關(guān)聯(lián)分析，調(diào)整網(wǎng)絡(luò)隱層節(jié)點p=0.5;e=0.5;%此兩個系數(shù)的設(shè)定是根據(jù)一些論文，已經(jīng)實驗的嘗試得出的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('灰色神經(jīng)網(wǎng)絡(luò)','生態(tài)用水量'); xlabel('年份'); RES0=YC-YY; res0=RES0./YY; figure bar(years,res0); xlabel('訓(xùn)練次數(shù)'); ylabel('訓(xùn)練誤差');17.3. 擬合
擬合神經(jīng)網(wǎng)絡(luò)是什么我忘記了，反正能用就行了
shenjin.m
function [net,tr] = shenjin(choice,P,T) net=newff(minmax(P),[20,1],... % 隱單元的神經(jīng)元數(shù) 20，輸出 1，可調(diào) {'tansig','purelin'}); if(choice==1)% 采用 L-M 優(yōu)化算法 TRAINLMnet.trainFcn='trainlm';% 設(shè)置訓(xùn)練參數(shù)net.trainParam.epochs = 500;net.trainParam.goal = 1e-6;net=init(net);% 重新初始化 elseif(choice==2)% 采用貝葉斯正則化算法 TRAINBRnet.trainFcn='trainbr';% 設(shè)置訓(xùn)練參數(shù)net.trainParam.epochs = 500;randn('seed',192736547);net = init(net);% 重新初始化 elseif(choice==3)% 采用動量梯度下降算法 TRAINGDM% 當(dāng)前輸入層權(quán)值和閾值inputWeights=net.IW{1,1};inputbias=net.b{1};% 當(dāng)前網(wǎng)絡(luò)層權(quán)值和閾值layerWeights=net.LW{2,1};layerbias=net.b{2};% 設(shè)置訓(xùn)練參數(shù)net.trainParam.show = 50;net.trainParam.lr = 0.05;net.trainParam.mc = 0.9;net.trainParam.epochs = 1000;net.trainParam.goal = 1e-3; end % 調(diào)用相應(yīng)算法訓(xùn)練 BP 網(wǎng)絡(luò) [net,tr]=train(net,P,T); end調(diào)用方法
% L-M 優(yōu)化算法（trainlm）和貝葉斯正則化算法（trainbr），用以訓(xùn)練 BP 網(wǎng)絡(luò) % 使其能夠擬合某一附加有白噪聲的正弦樣本數(shù)據(jù)，渡邊筆記 % 神經(jīng)網(wǎng)絡(luò)的初始和常數(shù)設(shè)置，請在shenjin.m修改 clc,clear;close all;% 0.全選，如果數(shù)據(jù)不多時可以操作 % 1.L-M 優(yōu)化算法 TRAINLM % 2.貝葉斯正則化算法 TRAINBR % 3.動量梯度下降算法 TRAINGDM % trainbfg---BFGS算法（擬牛頓反向傳播算法）訓(xùn)練函數(shù)； % trainbr---貝葉斯歸一化法訓(xùn)練函數(shù)； % traincgb---Powell-Beale共軛梯度反向傳播算法訓(xùn)練函數(shù)； % traincgp---Polak-Ribiere變梯度反向傳播算法訓(xùn)練函數(shù)； % traingd---梯度下降反向傳播算法訓(xùn)練函數(shù)； % traingda---自適應(yīng)調(diào)整學(xué)習(xí)率的梯度下降反向傳播算法訓(xùn)練函數(shù)； % traingdm---附加動量因子的梯度下降反向傳播算法訓(xùn)練函數(shù)； % traingdx---自適應(yīng)調(diào)整學(xué)習(xí)率并附加動量因子的梯度下降反向傳播算法訓(xùn)練函數(shù)； % trainlm---Levenberg-Marquardt反向傳播算法訓(xùn)練函數(shù)； % trainoss---OSS（one step secant）反向傳播算法訓(xùn)練函數(shù)； % trainrp---RPROP（彈性BP算法）反向傳播算法訓(xùn)練函數(shù)； % trainscg---SCG（scaled conjugate gradient）反向傳播算法訓(xùn)練函數(shù)； % trainb---以權(quán)值/閾值的學(xué)習(xí)規(guī)則采用批處理的方式進行訓(xùn)練的函數(shù)； % trainc---以學(xué)習(xí)函數(shù)依次對輸入樣本進行訓(xùn)練的函數(shù)； % trainr---以學(xué)習(xí)函數(shù)隨機對輸入樣本進行訓(xùn)練的函數(shù)。 % 當(dāng)調(diào)用train函數(shù)時，上述訓(xùn)練函數(shù)被用于訓(xùn)練網(wǎng)絡(luò) L = 3; %共有三種方法,我覺得夠用了 choice = 0;% P 為輸入矢量 P = [1 2 3 4 5]; % T 為目標(biāo)矢量 T = [1 2 3 4 5];% 繪制樣本數(shù)據(jù)點 % 3種全選，比較誰更合適 % choice = 0 if choice == 0for i = 1:L% 神經(jīng)網(wǎng)絡(luò)[net,tr] = shenjin(i,P,T);% 對 BP 網(wǎng)絡(luò)進行仿真A = sim(net,P);title(['神經(jīng)網(wǎng)絡(luò) 方法',num2str(i),' 計算結(jié)果']);% 計算仿真誤差E = T - A;MSE(i) = mse(E); % mse 函數(shù)用于計算均方誤差% 均方誤差（mean-square error, MSE）是反映估計量與被估計量之間差異程度的一種度量endfor i = 1:Ldisp(['神經(jīng)網(wǎng)絡(luò) 方法',num2str(i),' 的_均方誤差_為 : ',num2str(MSE(i))])end[~,Min] = min(MSE);disp(['神經(jīng)網(wǎng)絡(luò) 方法',num2str(Min),' 比較優(yōu)']) elsesubplot(1,2,1)plot(P,T,'+');title('原始數(shù)據(jù)點');% 神經(jīng)網(wǎng)絡(luò)[net,tr] = shenjin(i,P,T);% 對 BP 網(wǎng)絡(luò)進行仿真A = sim(net,P);title(['神經(jīng)網(wǎng)絡(luò) 方法',num2str(choice),' 計算結(jié)果']);% 計算仿真誤差E = T - A;MSE = mse(E); % mse 函數(shù)用于計算均方誤差% 均方誤差（mean-square error, MSE）是反映估計量與被估計量之間差異程度的一種度量disp(['神經(jīng)網(wǎng)絡(luò) 方法',num2str(choice),' 的_均方誤差_為 : ',num2str(MSE)]) end% 這個用于分析是否可以很好的進行線性回歸 % [m,b,r] = postreg(A, P); % 神經(jīng)網(wǎng)絡(luò)的簡單使用方法 % P = [1 2 3]; % [Y] = sim(net,P); % 輸入即可輸出，維度應(yīng)與元數(shù)據(jù)保存一致18. 貪心算法
% 渡邊筆記 % 貪心算法是一種思想，即尋找當(dāng)前最優(yōu)解，至使結(jié)果最優(yōu)解 % 貪心算法必須經(jīng)過證明才可以使用 % 其證明圍繞著：整個問題的最優(yōu)解一定由在貪心策略中存在的子問題的最優(yōu)解得來的 % 一旦貪心算法得以證明，它就是最快最優(yōu)的算法 % 實際上，貪心算法適用的情況很少 % 判斷一個問題是否適合用貪心法求解，目前還沒有一個通用的方法，在競賽中，需要憑個人的經(jīng)驗來判斷 % 貪心算法畢竟是貪心算法，只能緩一時，而不是長久之計， % 問題的模型、參數(shù)對貪心算法求解結(jié)果具有決定性作用，這在某種程度上是不能接受的 % 于是聰明的人類就發(fā)明了各種智能算法（也叫啟發(fā)式算法） % 但在我看來所謂的智能算法本質(zhì)上就是貪心算法和隨機化算法結(jié)合 % 例如傳統(tǒng)遺傳算法用的選擇策略就是典型的貪心選擇，正是這些貪心算法和隨機算法的結(jié)合，我們才看到今天各種各樣的智能算法 %% 例1 設(shè)有n個正整數(shù)，將它們連接成一排，組成一個最大的多位整數(shù)。 % n=3時，3個整數(shù)13，312，343，連成的最大整數(shù)為34331213。 % n=4時，4個整數(shù)7，13，4，246，連成的最大整數(shù)為7424613。 % 此題很容易想到使用貪心法，比如把整數(shù)按從大到小的順序連接起來，例子符合，但測試結(jié)果不全對 % 反例：12，121應(yīng)該組成12121而非12112 % 正確的標(biāo)準(zhǔn)是：先把整數(shù)轉(zhuǎn)換成字符串，然后比較a+b和b+a，如果a+b>=b+a，就把a排在b的前面，反之則把a排在b的后面 %% 例2 最短路徑問題 clear,clc;close all; n=4; % 設(shè)置點數(shù) x = zeros(1,n); % 產(chǎn)生一個行向量存儲坐標(biāo) y = zeros(1,n); luxian = 1:1:n; % 生成1到n的矩陣，存儲路線 paixu = 1:1:n;% 這一步是作出隨機點的圖 for i = 1 : nx(i) = randn * n ; %生成正態(tài)分布的隨機坐標(biāo)點y(i) = randn * n ;hold ontext(x(i)+0.3,y(i)+0.3,num2str(i)) %用text做好標(biāo)記， 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; %設(shè)置起點，路線的起點是1 num = 1; for a = 1:(n-2) %去除起點和終點需要n-2次判斷 paixu(:,1)=[]; %刪除上一個最優(yōu)點 dis = zeros(1,(n-a)); %用來存剩下各個點的距離for b = 1:(n-a) %用來獲取剩下各個點的距離dis(b) = d (num, paixu(b));end num1 = find( dis == min(dis) ); %根據(jù)距離最小得到最優(yōu)點位置 t = paixu(1); %通過中間變量互換位置paixu(1) = paixu(num1); %最優(yōu)點位置換到第一個paixu(num1) = t;num = paixu(1); %獲取下次進行操作的數(shù)luxian(a+1) = paixu(1); %將最優(yōu)點存入luxian數(shù)組（空出開頭） end luxian(n) = paixu(2); %補上最后一個點plot(x(luxian),y(luxian),'--o') ; %標(biāo)出點用虛線相連得到路徑 grid on;title('貪心算法計算最優(yōu)路徑');19. 小波分析
%%初始化程序 clear,clc t1=clock;%% 載入噪聲信號數(shù)據(jù)，數(shù)據(jù)為.mat格式，并且和程序放置在同一個文件夾下 x = 0:0.06:10; YSJ = sin(x)+0.2*rand(size(x));plot(YSJ); title('原始數(shù)據(jù)');%% 數(shù)據(jù)預(yù)處理，數(shù)據(jù)可能是存儲在矩陣或者是EXCEL中的二維數(shù)據(jù)，銜接為一維的，如果數(shù)據(jù)是一維數(shù)據(jù)，此步驟也不會影響數(shù)據(jù) [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('橫坐標(biāo)'); ylabel('縱坐標(biāo)'); title('原始信號');%% 硬閾值處理 lev=3; xd=wden(Y,'heursure','h','one',lev,'db4');%硬閾值去噪處理后的信號序列 figure(2) plot(X,xd) xlabel('橫坐標(biāo)'); ylabel('縱坐標(biāo)'); title('硬閾值去噪處理') set(gcf,'Color',[1 1 1])%% 軟閾值處理 lev=3; xs=wden(Y,'heursure','s','one',lev,'db4');%軟閾值去噪處理后的信號序列 figure(3) plot(X,xs) xlabel('橫坐標(biāo)'); ylabel('縱坐標(biāo)'); title('軟閾值去噪處理') set(gcf,'Color',[1 1 1]) %% 固定閾值后的去噪處理 lev=3; xz=wden(Y,'sqtwolog','s','sln',lev,'db4');%固定閾值去噪處理后的信號序列 figure(4) plot(X,xz); xlabel('橫坐標(biāo)'); ylabel('縱坐標(biāo)'); 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); %% 輸出結(jié)果 disp('-------------三種閾值設(shè)定方式的降噪處理結(jié)果---------------'); 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. 遺傳算法
我是非常不建議使用這個的，太復(fù)雜太麻煩了，每次換算遺傳因子都要想半天，而且代碼冗余，又長又臭
20.1. TSP 示例
cross.m
%交叉操作函數(shù) 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為需要交叉的位數(shù) p=(L-W+1);%隨機產(chǎn)生一個交叉位置 %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)); endendexchange.m
%對調(diào)函數(shù) exchange.mfunction [x,y]=exchange(x,y) temp=x; x=y; y=temp;endfit.m
%適應(yīng)度函數(shù)fit.m，每次迭代都要計算每個染色體在本種群內(nèi)部的優(yōu)先級別，類似歸一化參數(shù)。越大約好！ 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; endMutation.m
%變異函數(shù) 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;endmyLength.m
%染色體的路程代價函數(shù) 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)); endendplot_route.m
%連點畫圖函數(shù) 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調(diào)用方法
% 渡邊筆記 % 多種群遺傳算法解決 TSP 問題 % TSP問題即旅行商問題 % 經(jīng)典的TSP可以描述為 % 一個商品推銷員要去若干個城市推銷商品，該推銷員從一個城市出發(fā)，需要經(jīng)過所有城市后，回到出發(fā)地 % 應(yīng)如何選擇行進路線，以使總的行程最短。從圖論的角度來看，該問題實質(zhì)是在一個帶權(quán)完全無向圖中，找一個權(quán)值最小的哈密爾頓回路clear,clc;close all; % 輸入?yún)?shù) M=100; %種群的個數(shù) ITER=200; %迭代次數(shù) m=2; %適應(yīng)值歸一化淘汰加速指數(shù) Pc=0.8; %交叉概率 Pmutation=0.05; %變異概率 % 城市坐標(biāo) pos=[0 01 1-1 -12 3] ;%生成城市之間距離矩陣 [N,~] = size(pos); %%城市的個數(shù) 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,:);%初始化種群及其適應(yīng)函數(shù) fitness = zeros(M,1); len = zeros(M,1); for i=1:M % 計算每個染色體對應(yīng)的總長度len(i,1)=myLength(D,popm(i,:)); end maxlen=max(len);%最大回路 minlen=min(len);%最小回路 fitness=fit(len,m,maxlen,minlen); rr=find(len==minlen);%找到最小值的下標(biāo)，賦值為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);%求當(dāng)代最佳個體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%%求適應(yīng)度函數(shù)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%畫出所有城市坐標(biāo) 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); % 最短距離隨迭代次數(shù)的變化 title('最短距離的變化');21. 蟻群算法
21.1. TSP 示例
yiyu%% 蟻群算法求解TSP問題%% 清空環(huán)境變量 clear clc%% 導(dǎo)入數(shù)據(jù) 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%% 初始化參數(shù) m = 50; % 螞蟻數(shù)量 alpha = 1; % 信息素重要程度因子 beta = 5; % 啟發(fā)函數(shù)重要程度因子 rho = 0.1; % 信息素?fù)]發(fā)因子 Q = 1; % 常系數(shù) Eta = 1./D; % 啟發(fā)函數(shù) Tau = ones(n,n); % 信息素矩陣 Table = zeros(m,n); % 路徑記錄表 iter = 1; % 迭代次數(shù)初值 iter_max = 200; % 最大迭代次數(shù) Route_best = zeros(iter_max,n); % 各代最佳路徑 Length_best = zeros(iter_max,1); % 各代最佳路徑的長度 Length_ave = zeros(iter_max,1); % 各代路徑的平均長度%% 迭代尋找最佳路徑 while iter <= iter_max% 隨機產(chǎn)生各個螞蟻的起點城市start = zeros(m,1);for i = 1:mtemp = randperm(n);start(i) = temp(1);endTable(:,1) = start;% 構(gòu)建解空間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;% 計算城市間轉(zhuǎn)移概率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;% 迭代次數(shù)加1，清空路徑記錄表iter = iter + 1;Table = zeros(m,n); end%% 結(jié)果顯示 [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('城市位置橫坐標(biāo)') ylabel('城市位置縱坐標(biāo)') title(['蟻群算法優(yōu)化路徑(最短距離:' num2str(Shortest_Length) ')']) figure(2) plot(1:iter_max,Length_best,'b',1:iter_max,Length_ave,'r:') legend('最短距離','平均距離') xlabel('迭代次數(shù)') 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 297522. 元胞自動機
22.1. 森林火災(zāi)示例
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 end22.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; endns.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)調(diào)用方法
clc; clear; [flux,vmean] = ns(0.2,1e-9,25,30);23. 主成分分析
% 渡邊筆記 主成分分析 % 當(dāng)研究的問題涉及多個變量，并且變量間相關(guān)性明顯，即包含的信息有所重疊時 % 可以考慮用主成分分析的方法，這樣更容易抓住事物的主要矛盾，使問題簡化 clc,clear;close all; % 與主成分分析相關(guān)的 MATLAB 函數(shù)有pcacov、princomp函數(shù)% pcacov函數(shù)用來根據(jù) 協(xié)方差矩陣 或 相關(guān)系數(shù)矩陣 進行主成分分析 % [COEFF,latent,explained]=pcacov（V） % V 是總體或樣本的協(xié)方差矩陣或相關(guān)系數(shù)矩陣，對于p維總體，V是pxp的矩陣 % 輸出參數(shù) COEFF 是p個主成分的系數(shù)矩陣，它是pxp的矩陣，它的第i列是第i個主成分的系數(shù)向量 % 輸出參數(shù) latent 是p個主成分的方差構(gòu)成的向量，即V的p個特征值的大小(從大到小)構(gòu)成的向量 % 輸出參數(shù) explained 是p個主成分的貢獻率向量，已經(jīng)轉(zhuǎn)化為百分比 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函數(shù)式用于求導(dǎo)數(shù)和差分的. result_1(2: a ,2) = num2cell(-diff(latent)); % cumsum函數(shù)通常用于計算一個數(shù)組各行的累加值。 result_1(2: a+1 ,3:4) = num2cell([explained, cumsum(explained)]); % disp(result_1);% princomp函數(shù)用來根據(jù)樣本觀測值矩陣進行主成分分析，其調(diào)用格式如下： % [COEFF,SCORE,latent,tsquare] = princomp(X,'econ') % 根據(jù)樣本觀測值矩陣X進行主成分分析。輸入?yún)?shù)X是n行p列的矩陣，每一行對應(yīng)一個觀測(樣品)，每一列對應(yīng)一個變量 % 輸出參數(shù)COEFF是p個主成分析的系數(shù)矩陣，是pxp的矩陣，它的第i列對應(yīng)第i個主成分的系數(shù)向量 % 輸出參數(shù)SCORE是n個樣品的p個主成分得分矩陣，它是n行p列的矩陣 % 每一行對應(yīng)一個觀測，每一列對應(yīng)一個主成分，第i行第j列元素表示第i個樣品的第j個主成分得分 % SCORE與X是一一對應(yīng)的關(guān)系，是X在新坐標(biāo)系中的數(shù)據(jù)，可以通過X*系數(shù)矩陣得到 % 返回樣本協(xié)方差矩陣的特征值向量latent，它是由p個特征值構(gòu)成的列向量，其中特征值按降序排列。 % 返回一個包含n個元素的列向量tsquare，它的第i個元素的第i個觀測對應(yīng)的霍特林T^2統(tǒng)計量 % 表述了第i個觀測與數(shù)據(jù)集(樣本觀測矩陣)的中心之間的距離，可用來尋找遠(yuǎn)離中心的極端數(shù)據(jù) % 通過設(shè)置參數(shù)‘econ’參數(shù)，使得當(dāng) n <= p 時 % 只返回latent中的前n-1個元素(去掉不必要的0元素)及COEFF和SCORE矩陣中相應(yīng)的列 X = [2000 1000 3000 900 951900 990 3000 700 662000 900 3400 800 56]; % 行為樣本，列為特征量 XZ = zscore(X); % 數(shù)據(jù)標(biāo)準(zhǔn)化 [COEFF_2,SCORE,latent,tsquare] = pca(XZ); [m, n] = size(X); % 求X的行數(shù)和列數(shù) result_2 = cell(m, 4); % 定義一個n+1行，4列的元胞數(shù)組 result_2(1,:) = {'特征值', '差值', '貢獻率', '累積貢獻率'}; explained = 100*latent/sum(latent);% 計算貢獻率 %result_2中第1列的第2行到最后一行存放的數(shù)據(jù)(latent)特征值 result_2(2:end,1) = num2cell(latent); %result1中第2列的第2行到倒數(shù)第2行存放的數(shù)據(jù)(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); % 數(shù)據(jù)標(biāo)準(zhǔn)化 [COEFF_2,SCORE,latent,tsquare] = pca(XZ); [m, n] = size(X); % 求X的行數(shù)和列數(shù) result_2 = cell(m, 4); % 定義一個n+1行，4列的元胞數(shù)組 result_2(1,:) = {'特征值', '差值', '貢獻率', '累積貢獻率'}; explained = 100*latent/sum(latent);% 計算貢獻率 %result_2中第1列的第2行到最后一行存放的數(shù)據(jù)(latent)特征值 result_2(2:end,1) = num2cell(latent); %result1中第2列的第2行到倒數(shù)第2行存放的數(shù)據(jù)(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 需要轉(zhuǎn)置才可以使用 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 為回歸系數(shù)向量 % BINT 回歸系數(shù)的估計區(qū)間 % R 殘差 % RINT 置信區(qū)間 % STATS 用于檢驗回歸模型的統(tǒng)計量。有4個數(shù)值：判定系數(shù)r2，F統(tǒng)計量觀測值，檢驗的p的值，誤差方差的估計 % r2越接近1，回歸方程越顯著；F>F1?α(k,n?k?1)時拒絕H0，F越大，回歸方程越顯著；p<α?xí)r拒絕H0 % ALPHA 顯著性水平（缺少時默認(rèn)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)'; %預(yù)測值 plot([2004:2019],num_3','ro-',years,y2,'b*'); legend('生態(tài)用水量', '主成分回歸'); xlabel('時間/年'); ylabel('城市生態(tài)用水量/億立方米'); disp('回歸分析統(tǒng)計量為'); disp([ 'R^2 : ',num2str(STATS(1))]); % 模型由 這么多原始數(shù)據(jù) 解釋，最優(yōu)為 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)'; %預(yù)測值 plot([2017:2019],y3,'m*'); line([2016.5,2016.5],[0,17],'linestyle','--'); axis([2004 2019 0 17]); legend('城市生態(tài)用水量-原始數(shù)據(jù)', '主成分回歸在數(shù)據(jù)集的結(jié)果','主成分回歸在測試集的結(jié)果'); 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);%選定的這個數(shù)挖掉，相當(dāng)于挖坑while i<jwhile (i<j)&&(a(j)>=temp)%從右往左，找到第一個小于設(shè)定的數(shù)，j=j-1;enda(i)=a(j);%將找到的第一個小于設(shè)定的數(shù)填坑到最開始挖的坑里面去while (i<j)&&(a(i)<=temp)%從做到由，找到第一個大于選定的數(shù)i=i+1;enda(j)=a(i);%將找到的第一個大于選定的數(shù)填入上一步挖的坑里面去 % if i==j % a(j)=temp; % endenda(j)=temp;%最后，i=j，將選定的數(shù)再填到上一步挖的坑里面去a=QuickSort(a,leftIndex,j-1);%對左邊序列進行遞歸a=QuickSort(a,i+1,rightIndex);%對右邊序列進行遞歸 end end調(diào)用方法
% 分治算法，渡邊筆記 % 分治算法是一種思路，把大問題歸納為小問題來解決 % 使用分治法，問題應(yīng)該滿足 % 該問題的規(guī)?？s小到一定的程度就可以容易地解決 % 該問題可以分解為若干個規(guī)模較小的相同問題，即該問題具有最優(yōu)子結(jié)構(gòu)性質(zhì) % 利用該問題分解出的子問題的解可以合并為該問題的解 % 該問題所分解出的各個子問題是相互獨立的，即子問題之間不包含公共的子子問題 % 第一條特征是絕大多數(shù)問題都可以滿足的，因為問題的計算復(fù)雜性一般是隨著問題規(guī)模的增加而增加 % 第二條特征是應(yīng)用分治法的前提它也是大多數(shù)問題可以滿足的，此特征反映了遞歸思想的應(yīng)用 % 第三條特征是關(guān)鍵，能否利用分治法完全取決于問題是否具有第三條特征，如果具備了第一條和第二條特征，而不具備第三條特征，則可以考慮用貪心法或動態(tài)規(guī)劃法。 % 第四條特征涉及到分治法的效率，如果各子問題是不獨立的則分治法要做許多不必要的工作， % 重復(fù)地解公共的子問題，此時雖然可用分治法，但一般用動態(tài)規(guī)劃法較好 clc,clear;close all; a=[72 57 88 60 42 83 73 48 85]; % 輸入排序數(shù)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調(diào)用方法
% 分治法求最近點，渡邊筆記 clear;clc n = 20; % 隨機生成20個點 A=rand(n,2)*10; % 將20個點按橫坐標(biāo)升序排列 A = sortrows(A,1); % 分治法求隨機點的最近點對 [mindist1,y1,y2] = cloest(A,1,n); % 使用紅色的點標(biāo)記隨機點 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 %使用綠色十字標(biāo)記分治法的最近點對 plot(y1(1),y1(2),'go', 'linewidth', 2);hold on plot(y2(1),y2(2),'go', 'linewidth', 2);hold on26. 雜項
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. 數(shù)據(jù)擬合
懶得修改了，反正方法也就那樣
clc,clear % 產(chǎn)生數(shù)據(jù) 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; end26.3. 標(biāo)準(zhǔn)化
% 標(biāo)準(zhǔn)化處理，渡邊筆記 % 歸一化的依據(jù)非常簡單，不同變量往往量綱不同 % 歸一化可以消除量綱對最終結(jié)果的影響，使不同變量具有可比性 % 如兩個人體重差10KG，身高差0.02M % 在衡量兩個人的差別時體重的差距會把身高的差距完全掩蓋，歸一化之后就不會有這樣的問題 % 數(shù)據(jù)歸一化應(yīng)該針對對于的 y，而不是針對每條數(shù)據(jù)，針對每條數(shù)據(jù)是完全沒有意義的，因為只是等比例縮放，對之后的分類沒有任何作用 clc,clear;close all; x = [1 2 31 2 3]; % Min-max 極大極小標(biāo)準(zhǔn)化 % 當(dāng)有新數(shù)據(jù)加入時，可能導(dǎo)致xmax和xmin的變化，需要重新計算 X_min_max = mapminmax(x, 0, 1); % 默認(rèn)標(biāo)準(zhǔn)化區(qū)間為 0-1 disp('極大極小標(biāo)準(zhǔn)化為 :') disp(X_min_max);% z-score 標(biāo)準(zhǔn)化 % 其結(jié)果沒有任何意義，只能用于比較 % X_z_score = zscore(x); % disp('z-score 標(biāo)準(zhǔn)化 :') % disp(X_z_score);26.4. excel 數(shù)據(jù)轉(zhuǎn)置
clc,clear % 省會城市 順序表 % 北京 天津 石家莊 太原 呼和浩特 沈陽 長春 哈爾濱 上海 南京 % 杭州 合肥 福州 南昌 濟南 鄭州 武漢 長沙 廣州 南寧 % ?？?重慶 成都 貴陽 昆明 拉薩 西安 蘭州 西寧 銀川 % 烏魯木齊 filename_lujin_1 = 'C:\Users\99791\Desktop\生態(tài)用水\數(shù)據(jù)\國家統(tǒng)計局\分省\'; % 文件的前面的路徑名稱 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\生態(tài)用水\數(shù)據(jù)\'; 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. 線性規(guī)劃
lingo 還是比 matlab 強，下面的方法好久不用了
% function 函數(shù)解多元方程 最小值 % function 基本語法 % [x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) % x的返回值是決策向量x的取值，fval的返回值是目標(biāo)函數(shù)f(x)的取值 % fun是用 M 文件定義的函數(shù)f(x),代表了(非)線性目標(biāo)函數(shù) % x0是x的初始值 % A,b,Aeq,beq 定義了線性約束 ,如果沒有線性約束，則A=[],b=[],Aeq=[],beq=[] % lb和ub是變量x的下界和上界，如果下界和上界沒有約束，則lb=[],ub=[],也可以寫成lb的各分量都為 -inf,ub的各分量都為inf % nonlcon是用 M 文件定義的非線性向量函數(shù)約束 % options 定義了優(yōu)化參數(shù)，不填寫表示使用 Matla b默認(rèn)的參數(shù)設(shè)置 clc,clear;close all; options = optimset('Algorithm','active-set'); % 可以不設(shè)置 [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]... % 設(shè)置范圍 ,[],options); disp(['函數(shù)的滿足最小值的 X 解為 : ',num2str(X)]); disp(['函數(shù)的最小值的 Y 解為 : ',num2str(Y)]);總結(jié)
以上是生活随笔為你收集整理的数学建模方法总结(matlab)的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： ebs克隆oracle not,Orac
- 下一篇： 直扩同步的跟踪 matlab,基于FPG