数学建模方法总结(matlab)
生活随笔
收集整理的這篇文章主要介紹了
数学建模方法总结(matlab)
小編覺(jué)得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.
1. 前言
在這篇文章中，我會(huì)把我所接觸的數(shù)學(xué)建模的知識(shí)的代碼分享給大家，有的是我自己寫(xiě)的，更多的是網(wǎng)上借鑒并修改為可執(zhí)行可用的代碼
需要說(shuō)明的是，我不太了解其中的數(shù)學(xué)實(shí)現(xiàn)的具體方法和原理，我只分享源碼和作用條件以及使用的經(jīng)驗(yàn)說(shuō)明（更詳細(xì)見(jiàn)注釋），網(wǎng)上有不少文章對(duì)某些方法講得非常透徹，因此我沒(méi)必要贅述
各位只需要忘記代碼或使用條件時(shí)從我這里 copy 就好了
2. 灰色預(yù)測(cè)模型
這是用一種在數(shù)學(xué)上對(duì)數(shù)據(jù)進(jìn)行累加和累減以取巧性地削弱預(yù)測(cè)數(shù)據(jù)與原始數(shù)據(jù)的關(guān)聯(lián)性
最好用于短期預(yù)測(cè)，只能用于遞增序列
2.1. GM（1，1）
我也不懂為什么叫 GM(1,1)，但是用得最多的就是它，又稱一階灰色方程
clc,clear;close all y = [0.6 1.1 1.6] yuceshu = 3; % 本程序主要用來(lái)計(jì)算根據(jù)灰色理論建立的模型的預(yù)測(cè)值 % 應(yīng)用的數(shù)學(xué)模型是 GM(1,1) % 原始數(shù)據(jù)的處理方法是一次累加 % 渡辺曜改進(jìn)了這個(gè)代碼,通常GM預(yù)測(cè)第一個(gè)數(shù)是不夠優(yōu)秀的 % y 是應(yīng)該輸入的數(shù)組，yuceshu是你想預(yù)測(cè)的后幾個(gè)數(shù)據(jù) 我們一般就預(yù)測(cè)兩個(gè) % 示例 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擬合的第一個(gè)數(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(['后驗(yàn)差比值為：',num2str(c)]); if c<0.35disp('系統(tǒng)預(yù)測(cè)精度好') else if c<0.5disp('系統(tǒng)預(yù)測(cè)精度合格')else if c<0.65disp('系統(tǒng)預(yù)測(cè)精度勉強(qiáng)')elsedisp('系統(tǒng)預(yù)測(cè)精度不合格')endend enddisp(['下個(gè)擬合值為 ',num2str(ys(n+1))]); for i=2:yuceshudisp(['再下個(gè)擬合值為',num2str(ys(n+i))]); end2.2. 分?jǐn)?shù)階 GM(1,1)
分?jǐn)?shù)階灰色方程的預(yù)測(cè)效果優(yōu)于原始灰色模型，我只找到了 python 版本的可運(yùn)行代碼
首先，將 GM(1,1)打包成一個(gè)類
GM.py
##渡辺筆記 ##傳統(tǒng)的灰色模型 import numpy as np import xlrd import xlwt import datetime class Grey_model(object):def __init__(self,input_value):##初始化過(guò)程，先建立好變量空間，即累加序列，背景值序列，B和Y矩陣，擬合序列，預(yù)測(cè)值序列等。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()##計(jì)算累加序列def accumulation(self):for i in range(len(self.input_value)):self.accumulation_value[i]=np.sum(self.input_value[0:i+1])##計(jì)算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##計(jì)算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]##計(jì)算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]##計(jì)算參數(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()##計(jì)算預(yù)測(cè)值def predict(self,number_of_forecast):self.u_matrix()self.predicted_accumulation_value=[]##使用了float(%.2f% 來(lái)調(diào)整輸出值的小數(shù)的個(gè)數(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開(kāi)始算，故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])))##計(jì)算誤差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]+ '%'## 打開(kāi)xls，只能用這種格式了 work_book = xlrd.open_workbook('first_shuju.xls') ## 按sheet表名稱獲取sheet對(duì)象，名稱分大小寫(xiě) 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ù)測(cè)值')print(A.predicted_value)print('\nMAPE值')print(A.MAPE)## xlwt_data = A.predicted_valuedef save_excel(target_list, output_file_name):將數(shù)據(jù)寫(xiě)入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)在這個(gè)類的基礎(chǔ)上，運(yùn)用分?jǐn)?shù)階進(jìn)行優(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## 打開(kāi)xls，只能用這種格式了 work_book = xlrd.open_workbook('first_shuju.xls') ## 按sheet表名稱獲取sheet對(duì)象，名稱分大小寫(xiě) 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模型的計(jì)算模塊 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ù)測(cè)值的還原方式 A.test_error()##誤差用MAPE值來(lái)衡量 print('預(yù)測(cè)值') print(A.predicted_value) print('\nMAPE值') print(A.MAPE)def save_excel(target_list, output_file_name):將數(shù)據(jù)寫(xiě)入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
感覺(jué)大多都能用 floyd 解決
% 這是floyd算法，以下注釋為渡邊所加 % floyd不能解決多條最短路徑和負(fù)權(quán)回路的問(wèn)題 % 解稠密圖效果最佳，邊權(quán)可正可負(fù)，是雙源解法 % 對(duì)于稠密圖，效率要高于執(zhí)行|V|次Dijkstra算法，也要高于執(zhí)行|V|次SPFA算法。 % 缺點(diǎn)：時(shí)間復(fù)雜度比較高，不適合計(jì)算大量數(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é)點(diǎn)，stop為指定終止結(jié)點(diǎn) D = W; %最短距離矩陣賦初值 n = length(D); %n為結(jié)點(diǎn)個(gè)數(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ù)個(gè)數(shù)輸入有誤！', 'Warning!') elsedis = D(start, stop); %指定兩結(jié)點(diǎn)間的最短距離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é)點(diǎn)之間的最短路徑 end3.2. Dijkstra
Dijkstra 運(yùn)算速度快，但是用不了負(fù)權(quán)圖
% 這是 dijkstra 算法，以下注釋為渡邊所加 % 不能用于負(fù)權(quán)圖，但是由于時(shí)間復(fù)雜度低，適合在大數(shù)據(jù)中運(yùn)算 % 是單源解法 % 輸入示例 % 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 搜索的起點(diǎn) e 搜索的終點(diǎn) n=length(W);%節(jié)點(diǎn)數(shù) D = W(st,:); visit= ones(1:n); visit(st)=0; parent = zeros(1,n);%記錄每個(gè)節(jié)點(diǎn)的上一個(gè)節(jié)點(diǎn)path =[];for i=1:n-1temp = [];%從起點(diǎn)出發(fā)，找最短距離的下一個(gè)點(diǎn)，每次不會(huì)重復(fù)原來(lái)的軌跡，設(shè)置visit判斷節(jié)點(diǎn)是否訪問(wèn)for j=1:nif visit(j)temp =[temp D(j)];elsetemp =[temp inf];endend[value,index] = min(temp);visit(index) = 0;%更新 如果經(jīng)過(guò)index節(jié)點(diǎn)，從起點(diǎn)到每個(gè)節(jié)點(diǎn)的路徑長(zhǎng)度更小，則更新，記錄前趨節(jié)點(diǎn)，方便后面回溯循跡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)問(wè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. 層次分析法
感覺(jué)這個(gè)方法就像我們土木工程所謂的“拍腦袋”，先拍腦袋得出一個(gè)比較矩陣，待反復(fù)測(cè)試符合規(guī)范后，就可以拿去比較了
% 本代碼只需要在提示輸入處，輸入構(gòu)成的比較矩陣 就會(huì)輸出各指標(biāo)的權(quán)重值。 % 例子： 選拔干部考慮5個(gè)條件：品德x1，才能x2，資歷x3，年齡x4，群眾關(guān)系x5。某決策人用成對(duì)比較法，得到成對(duì)比較陣如下： % [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; 怎么來(lái)的呢？ % 其實(shí)這些x的值就是由決策人決定了，對(duì)應(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)個(gè)數(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); %對(duì)應(yīng)特征向量 CI=(B-n)/(n-1); %計(jì)算一致性檢驗(yàn)指標(biāo)CI CR=CI/RI(1,n); if CR<0.10disp('CI=');disp(CI);disp('CR=');disp(CR);disp('對(duì)比矩陣A通過(guò)一致性檢驗(yàn)，各向量權(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('對(duì)比矩陣A未通過(guò)一致性檢驗(yàn)，需對(duì)對(duì)比矩陣A重新構(gòu)造'); end sc = Q';5. 聚類分析
這類方法對(duì)于我的專業(yè)毫無(wú)實(shí)際用處，我并不需要知道管道應(yīng)該分成幾類，因?yàn)橐?guī)范早就規(guī)定好了的
聚類分析方法豐富，要根據(jù)實(shí)際情況謹(jǐn)慎選取
5.1. k-mean
% k_mean 算法，渡邊筆記 % 對(duì)處理大型數(shù)據(jù)集，該算法保持可收縮性，高效性； % 當(dāng)簇接近正態(tài)分布時(shí)，效果最好； % 若簇中含有異常點(diǎn)，將導(dǎo)致均值偏離嚴(yán)重；（即對(duì)噪聲、孤立點(diǎn)敏感）； % 不適用于發(fā)現(xiàn)非凸形狀的簇，或大小差別大的簇 % 中k是事先給定的，這個(gè)k值難以估計(jì)，很多時(shí)候并不知道數(shù)據(jù)集應(yīng)該分成多少個(gè)類別最合適 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. 層次聚類
% 層次聚類，渡邊筆記 % 看圖說(shuō)話，分類清晰明了 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); % 計(jì)算矩陣X中樣本兩兩之間的距離，但得到的Y是個(gè)行向量 D = squareform(Y); % 將行向量的Y轉(zhuǎn)換成方陣，方便觀察兩點(diǎn)距離(方陣的對(duì)角線元素都是0) Z = linkage(Y); % 產(chǎn)生層次聚類樹(shù)，默認(rèn)采用最近距離作為類間距離的計(jì)算公式 dendrogram(Z); % 圖示層次聚類樹(shù) title('層次聚類法');ylabel('標(biāo)準(zhǔn)距離');5.3. 直接聚類
% 直接聚類，渡筆記 % 不管機(jī)理，不管距離，比較愚蠢但是有用且方便 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時(shí)，分到不同類(簇)中 T2=clusterdata(X,3); % 如果cutoff是一個(gè)≥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. 最小生成樹(shù)
好像是個(gè)示例，太久遠(yuǎn)忘了干什么用的了
zhuixiao_shengchengshu.m
clc,clear,close all; % 使用 Kruskal 的算法來(lái)計(jì)算 無(wú)向圖的最小生成樹(shù) 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)建并繪制一個(gè)立方體圖 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);% 計(jì)算并在圖上方繪制圖的最小生成樹(shù)。T 包含的節(jié)點(diǎn)與 G 相同，但包含的邊僅為后者的子集 % figure; [T,pred] = minspantree(G); % p = plot(G,'EdgeLabel',G.Edges.Weight); highlight(p,T) % 高亮顯示7. 誤差檢驗(yàn)
7.1. 一般檢驗(yàn)
看著用吧，我覺(jué)得沒(méi)什么意義，反正計(jì)算機(jī)計(jì)算出的誤差都挺小的，最多能在論文中湊字?jǐn)?shù)
% 渡辺筆記 % 各類檢驗(yàn)，除決定系數(shù)是1最好，都是0最好 YReal = [1 2 3 4 5]; YPred = [1 2 3 4 5.1]; % 平均絕對(duì)百分比誤差（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); % 平均絕對(duì)誤差（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 檢驗(yàn)
% 渡邊筆記 t檢驗(yàn) % 待檢驗(yàn)的數(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è)置一個(gè)值作為t檢驗(yàn)執(zhí)行的顯著性水平 % 在不設(shè)置ALPHA的情況下默認(rèn)ALPHA為0.05，即計(jì)算x和y在5%的顯著性水平下是否來(lái)自同一分布(假設(shè)是否被接受） % H = 0，P > 0.05，表明零假設(shè)被接受，即x，y在統(tǒng)計(jì)上可看做來(lái)自同一分布的數(shù)據(jù) % H = 1，P < 0.05，表明零假設(shè)被拒絕，即x，y在統(tǒng)計(jì)上認(rèn)為是來(lái)自不同分布的數(shù)據(jù)，即有區(qū)分度 disp(['H = ',num2str(H)]); disp(['P = ',num2str(P)]); if H == 0disp('x，y在統(tǒng)計(jì)上可看做來(lái)自同一分布的數(shù)據(jù)') elsedisp('x，y在統(tǒng)計(jì)上可看做來(lái)自不同分布的數(shù)據(jù)') end7.3. 回歸檢驗(yàn)
clc,clear;close all; % 回歸檢驗(yàn)；渡邊筆記 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ù)的估計(jì)區(qū)間 % R 殘差 % RINT 置信區(qū)間 % STATS 用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量。有4個(gè)數(shù)值：判定系數(shù)r2，F統(tǒng)計(jì)量觀測(cè)值，檢驗(yàn)的p的值，誤差方差的估計(jì) % r2越接近1，回歸方程越顯著；F>F1?α(k,n?k?1)時(shí)拒絕H0，F越大，回歸方程越顯著；p<α?xí)r拒絕H0 % ALPHA 顯著性水平（缺少時(shí)默認(rèn)0.05） y2 = B(2).* x + B(1); %預(yù)測(cè)值 plot(x', y' ,'+', x, y2); %回歸效果圖 disp('回歸分析統(tǒng)計(jì)量為'); disp([ 'R^2 : ',num2str(STATS(1))]); % 模型由 這么多原始數(shù)據(jù) 解釋，最優(yōu)為 1 disp([ 'F : ',num2str(STATS(2))]); disp([ '檢驗(yàn) p : ',num2str(STATS(3))]); % 日常小于 0.05 最好 disp([ '誤差方差 : ',num2str(STATS(4))]); % 越小越好7.4. 交叉檢驗(yàn)
% 交叉檢驗(yàn)，渡邊筆記 % 提高數(shù)據(jù)的利用率 % 有效的避免過(guò)學(xué)習(xí)以及欠學(xué)習(xí)狀態(tài)的發(fā)生,最后得到的結(jié)果也比較具有說(shuō)服性 % 常用的 K 值有 3,6,10 等 clear,clc;close all;% 導(dǎo)入數(shù)據(jù) data = reshape(1:144,12,12); % 每行的是這個(gè)樣本屬性的數(shù)據(jù)，最后一個(gè)數(shù)據(jù)是標(biāo)簽 [data_r, data_c] = size(data); % 將數(shù)據(jù)樣本隨機(jī)分割為 n 部分，做 n 折交叉檢驗(yàn) n = 5; indices = crossvalind('Kfold', data_r, n); for i = 1 : n% 獲取第i份測(cè)試數(shù)據(jù)的索引邏輯值test = (indices == i);% 取反，獲取第i份訓(xùn)練數(shù)據(jù)的索引邏輯值train = ~test;% 測(cè)試用的數(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);% 下面放入檢驗(yàn)數(shù)據(jù)的代碼end8. 差分
% 數(shù)據(jù)差分，渡邊筆記 % 差分就是后一個(gè)減去前一個(gè)，使得數(shù)據(jù)平穩(wěn)化 % n 階就是做 n 次差分，具體直接看例子就好 clc,clear;close all; Y = [0 5 15 30 50 75 105]; disp(['未差分時(shí)為 : ',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ù)白噪聲檢驗(yàn)，渡邊筆記 % 如果易知原始數(shù)據(jù)不是平穩(wěn)的，則不能做隨機(jī)性檢驗(yàn) % 接下來(lái)要求差分,目的: 變成平穩(wěn)的數(shù)據(jù) % p 如果比 0.05 小就不是白噪聲序列，可以使用時(shí)間序列 % 某種現(xiàn)象的指標(biāo)數(shù)值按照時(shí)間順序排列而成的數(shù)值序列 % 因?yàn)闀r(shí)間序列是某個(gè)指標(biāo)數(shù)值長(zhǎng)期變化的數(shù)值表現(xiàn) % 所以時(shí)間序列數(shù)值變化背后必然蘊(yùn)含著數(shù)值變換的規(guī)律性 % 這些規(guī)律性就是時(shí)間序列分析的切入點(diǎn) % 如果原數(shù)據(jù)平穩(wěn)，但是沒(méi)通過(guò)，可以直接差分 clc,clear;close all;x = []; % 時(shí)間，一般做題就是順序時(shí)間排列 yanchi=[6 12 18]; % 通常是做6 12 18 24步延遲，這個(gè)數(shù)據(jù)的選擇上限請(qǐng)根據(jù)報(bào)錯(cuò)來(lái)調(diào)整 y = [970279 1259308 1127571 1163959 1169540 1076938 991350 953275 951508 904434 889381 864015 836236 ]';; [h1] = adftest(y); %檢驗(yàn)是否平穩(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); % 在這里對(duì)數(shù)據(jù)進(jìn)行 n 階差分[h1,p1,adf,ljz] = adftest(y_1); %檢驗(yàn)是否平穩(wěn)if h1 == 1disp(['差分后是平穩(wěn)的，做了 ',num2str(i),' 階差分']);subplot(2,2,3)plot(y_1); % 一階差分?jǐn)?shù)據(jù)時(shí)序圖title('一階差分?jǐn)?shù)據(jù)時(shí)序圖')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 % 隨時(shí)間的變化值 subplot(2,2,1) plot(y); % 原始數(shù)據(jù)時(shí)序圖 title('原始數(shù)據(jù)時(shí)序圖') 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)計(jì)量 fprintf('%15s%15s%15s','延遲階數(shù)','卡方統(tǒng)計(jì)量','p值'); fprintf('\n'); for i=1:length(yanchi) % i = 1,時(shí)候?yàn)?，i = 2時(shí)候?yàn)?2fprintf('%18f%19f%19f',yanchi(i),Qstat(i),pValue(i));fprintf('\n'); endif sum(find(pValue > 0.05))disp('但是沒(méi)通過(guò)立法白噪聲檢驗(yàn)');i = 1;while 1y_1 = diff(y,i); % 在這里對(duì)數(shù)據(jù)進(jìn)行 n 階差分[h1,p1,adf,ljz] = adftest(y_1); %檢驗(yàn)是否平穩(wěn)if h1 == 1disp(['差分后是平穩(wěn)的，做了 ',num2str(i),' 階差分']);subplot(2,2,3)plot(y_1); % 一階差分?jǐn)?shù)據(jù)時(shí)序圖title('一階差分?jǐn)?shù)據(jù)時(shí)序圖')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% 再來(lái)一次[H,pValue,Qstat,CriticalValue]=lbqtest(y_1,'lags',yanchi); % H.結(jié)果，pValue.p值, Qstat.卡方統(tǒng)計(jì)量 fprintf('%15s%15s%15s','延遲階數(shù)','卡方統(tǒng)計(jì)量','p值'); fprintf('\n'); for i=1:length(yanchi) % i = 1,時(shí)候?yàn)?，i = 2時(shí)候?yàn)?2fprintf('%18f%19f%19f',yanchi(i),Qstat(i),pValue(i));fprintf('\n'); end if sum(find(pValue < 0.05))disp('數(shù)據(jù)通過(guò)立法白噪聲檢驗(yàn)'); elsedisp('數(shù)據(jù)沒(méi)通過(guò)立法白噪聲檢驗(yàn)'); 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代表不修改計(jì)算后的指標(biāo)權(quán)重 for i=1:nx(:,i)=(x(:,i)-min(x(:,i)))/(max(x(:,i))-min(x(:,i)))+1; % 對(duì)原始數(shù)據(jù)進(jìn)行非負(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)重計(jì)算 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ì)算綜合分?jǐn)?shù) end disp('各指標(biāo)權(quán)重為：') disp(w) disp('各項(xiàng)綜合分?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%% 對(duì)異常值賦空值 F = find(yi_chang == 1); y(F) = NaN; % 令數(shù)據(jù)點(diǎn)缺失 testdata = [x' y'];subplot(2,2,2); plot(testdata(:,1),testdata(:,2)); title('去除差異值');%% 對(duì)數(shù)據(jù)進(jìn)行補(bǔ)全 % 數(shù)據(jù)補(bǔ)全方法選擇 % 1.線性插值 linear 2.分段三次樣條插值 spline 3.保形分段三次樣條插值 pchip % 4.移動(dòng)滑窗插補(bǔ) 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); % 窗口長(zhǎng)度為 10 的移動(dòng)均值 endsubplot(2,2,3); plot(testdata_1(:,1),testdata_1(:,2)); title('數(shù)據(jù)補(bǔ)全結(jié)果');%% 進(jìn)行數(shù)據(jù)平滑處理 % 濾波器選擇 1.Savitzky-golay 2.rlowess 3.rloess choice_3 = 2; lvboqi = char('Savitzky-golay', 'rlowess', 'pchip', 'rloess'); % 通過(guò)求 n 元素移動(dòng)窗口的中位數(shù)，來(lái)對(duì)數(shù)據(jù)進(jìn)行平滑處理 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)在每個(gè)時(shí)期所處的狀態(tài)是隨機(jī)的 % 從一時(shí)期到下時(shí)期的狀態(tài)按一定概率轉(zhuǎn)移 % 下時(shí)期狀態(tài)只取決于本時(shí)期狀態(tài)和轉(zhuǎn)移概率。即已知現(xiàn)在，將來(lái)與過(guò)去無(wú)關(guān)（無(wú)后效性） % 過(guò)去不影響未來(lái)的判斷，是馬氏鏈的本質(zhì) % P_0 是初始分布，P_n 是轉(zhuǎn)移矩陣，則在 n 未來(lái)的概率為 % P_0 * ( P_n ^ n ) clc,clear;close all;% 第一次買(mǎi)鹽，概率為 P_0 % 已知 買(mǎi)鹽傾向 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)計(jì)矩陣 P_0 = [0.4 0.3 0.4]; % 初始概率% 求狀態(tài)轉(zhuǎn)移矩陣 [m,~] = size(a); ni=sum(a,2); %計(jì)算矩陣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)]);% 概率平衡時(shí)，有 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 趨向于無(wú)窮的 X 的解為 : ',num2str(x)]);13.2. 帶利潤(rùn)的馬氏鏈
% 帶利潤(rùn)的馬氏鏈模型，渡邊筆記 clc,clear;close all;a=[5 54 6 ]; % 輸入統(tǒng)計(jì)矩陣 R = [9 33 -7]; % 輸入利潤(rùn)矩陣% 對(duì)統(tǒng)計(jì)矩陣，求狀態(tài)轉(zhuǎn)移矩陣 [m,~] = size(a); ni=sum(a,2); %計(jì)算矩陣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 期利潤(rùn) n = 3;% n = 1 時(shí) for i = 1:mfor j = 1:mv(i,j) = R(i,j) * ZT(i,j);endSum = sum(v,2); end % n = n 時(shí)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),' 期利潤(rùn)為 ',num2str(V_Sum(i,:))])elsedisp(['處于狀態(tài) ',num2str(i),' 的 ',num2str(n),' 期利潤(rùn)為 ',num2str(Sum(i,:))])end end% 帶利潤(rùn)的馬氏鏈一般不存在平衡 % 因?yàn)殄X(qián)是賺不完的14. 蒙特卡羅
% 渡邊筆記 % Monte_Carlo 是一種方法的概述，不是什么特定的算法 % 主要思想是利用大量隨機(jī)數(shù)來(lái)實(shí)現(xiàn)真實(shí)值的擬合 %% 例1 計(jì)算 0-3 上 x^2 的積分，正解應(yīng)該是 9 N=500; x=unifrnd(0,3,N,1); % 設(shè)定隨機(jī)數(shù)及其終止上限 M = mean(3*x.^2);disp(['M = ',num2str(M)]); % 將隨機(jī)數(shù)代入原式求平均值，特別的，蒙特卡洛需要式乘以積分區(qū)間,下同 %% 例2 求相交面積 %給定曲線y =2 – x2 和曲線y3 = x2，曲線的交點(diǎn)為：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ū)域，用蒙特卡羅方法計(jì)算區(qū)域面積。 P=rand(1000,2); %隨機(jī)產(chǎn)生1000個(gè)點(diǎn) %定義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ù) %計(jì)算索引的長(zhǎng)度 M=length(II); %計(jì)算面積 S = 4*M/1000;disp(['S = ',num2str(S)]); subplot(1,2,2); plot(x(II),y(II),'g.'); close all % 刪除這個(gè)就可以顯示圖像 %% 例3 求圓周率 n = 1000; %總的實(shí)驗(yàn)次數(shù) m = 0; %落在圓中點(diǎn)的次數(shù) %循環(huán)實(shí)驗(yà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. 模糊評(píng)價(jià)
先打包成一個(gè) function
fuzzy_zhpj.m
function[B]=fuzzy_zhpj(model,A,R) %模糊綜合評(píng)判 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. 運(yùn)用數(shù)學(xué)方法
clc,clear % 模擬退火(求最小) temperature = 100; % 初始溫度 final_temperature = 1; % 最低溫度 time_temperature = 1; % 用于計(jì)算步數(shù) time_tuihuo = 10; % 退火步長(zhǎng) cooling_rate = 0.95; % 降溫步數(shù) % 初始化隨機(jī)數(shù)生成器，以使結(jié)果具備可重復(fù)性 % rng(0,'twister'); % 生成范圍a-b的隨機(jī)數(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; % 隨機(jī)數(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. 運(yùn)用 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); % 退火開(kāi)始點(diǎn) lb = []; ub = []; % 退火實(shí)施范圍，可以不設(shè)置 % 實(shí)時(shí)觀測(cè) 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); % 退火開(kāi)始點(diǎn) lb = []; ub = []; % 退火實(shí)施范圍，可以不設(shè)置 % 實(shí)時(shí)觀測(cè) 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ù)矩陣,每行為對(duì)應(yīng)的屬性 p = [1 1 12 2 23 3 3]; %目標(biāo)（輸出）數(shù)據(jù)矩陣 t = [6 6 6];%對(duì)訓(xùn)練集中的輸入數(shù)據(jù)矩陣和目標(biāo)數(shù)據(jù)矩陣進(jìn)行歸一化處理 [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次迭代都沒(méi)變化，則matlab會(huì)默認(rèn)終止訓(xùn)練。為了讓程序繼續(xù)運(yùn)行，用以下命令取消這條設(shè)置 net.divideFcn = ''; %開(kāi)始訓(xùn)練網(wǎng)絡(luò) net = train(net, pn, tn); %使用訓(xùn)練好的網(wǎng)絡(luò)，基于訓(xùn)練集的數(shù)據(jù)對(duì)BP網(wǎng)絡(luò)進(jìn)行仿真得到網(wǎng)絡(luò)輸出結(jié)果 %(因?yàn)檩斎霕颖?#xff08;訓(xùn)練集）容量較少，否則一般必須用新鮮數(shù)據(jù)進(jìn)行仿真測(cè)試) 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é)點(diǎn)p=0.5;e=0.5;%此兩個(gè)系數(shù)的設(shè)定是根據(jù)一些論文，已經(jīng)實(shí)驗(yàn)的嘗試得出的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;%進(jìn)行測(cè)試 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)% 采用動(dòng)量梯度下降算法 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è)置，請(qǐng)?jiān)趕henjin.m修改 clc,clear;close all;% 0.全選，如果數(shù)據(jù)不多時(shí)可以操作 % 1.L-M 優(yōu)化算法 TRAINLM % 2.貝葉斯正則化算法 TRAINBR % 3.動(dòng)量梯度下降算法 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---附加動(dòng)量因子的梯度下降反向傳播算法訓(xùn)練函數(shù)； % traingdx---自適應(yīng)調(diào)整學(xué)習(xí)率并附加動(dòng)量因子的梯度下降反向傳播算法訓(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ī)則采用批處理的方式進(jìn)行訓(xùn)練的函數(shù)； % trainc---以學(xué)習(xí)函數(shù)依次對(duì)輸入樣本進(jìn)行訓(xùn)練的函數(shù)； % trainr---以學(xué)習(xí)函數(shù)隨機(jī)對(duì)輸入樣本進(jìn)行訓(xùn)練的函數(shù)。 % 當(dāng)調(diào)用train函數(shù)時(shí)，上述訓(xùn)練函數(shù)被用于訓(xùn)練網(wǎng)絡(luò) L = 3; %共有三種方法,我覺(jué)得夠用了 choice = 0;% P 為輸入矢量 P = [1 2 3 4 5]; % T 為目標(biāo)矢量 T = [1 2 3 4 5];% 繪制樣本數(shù)據(jù)點(diǎn) % 3種全選，比較誰(shuí)更合適 % choice = 0 if choice == 0for i = 1:L% 神經(jīng)網(wǎng)絡(luò)[net,tr] = shenjin(i,P,T);% 對(duì) BP 網(wǎng)絡(luò)進(jìn)行仿真A = sim(net,P);title(['神經(jīng)網(wǎng)絡(luò) 方法',num2str(i),' 計(jì)算結(jié)果']);% 計(jì)算仿真誤差E = T - A;MSE(i) = mse(E); % mse 函數(shù)用于計(jì)算均方誤差% 均方誤差（mean-square error, MSE）是反映估計(jì)量與被估計(jì)量之間差異程度的一種度量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ù)點(diǎn)');% 神經(jīng)網(wǎng)絡(luò)[net,tr] = shenjin(i,P,T);% 對(duì) BP 網(wǎng)絡(luò)進(jìn)行仿真A = sim(net,P);title(['神經(jīng)網(wǎng)絡(luò) 方法',num2str(choice),' 計(jì)算結(jié)果']);% 計(jì)算仿真誤差E = T - A;MSE = mse(E); % mse 函數(shù)用于計(jì)算均方誤差% 均方誤差（mean-square error, MSE）是反映估計(jì)量與被估計(jì)量之間差異程度的一種度量disp(['神經(jīng)網(wǎng)絡(luò) 方法',num2str(choice),' 的_均方誤差_為 : ',num2str(MSE)]) end% 這個(gè)用于分析是否可以很好的進(jìn)行線性回歸 % [m,b,r] = postreg(A, P); % 神經(jīng)網(wǎng)絡(luò)的簡(jiǎn)單使用方法 % P = [1 2 3]; % [Y] = sim(net,P); % 輸入即可輸出，維度應(yīng)與元數(shù)據(jù)保存一致18. 貪心算法
% 渡邊筆記 % 貪心算法是一種思想，即尋找當(dāng)前最優(yōu)解，至使結(jié)果最優(yōu)解 % 貪心算法必須經(jīng)過(guò)證明才可以使用 % 其證明圍繞著：整個(gè)問(wèn)題的最優(yōu)解一定由在貪心策略中存在的子問(wèn)題的最優(yōu)解得來(lái)的 % 一旦貪心算法得以證明，它就是最快最優(yōu)的算法 % 實(shí)際上，貪心算法適用的情況很少 % 判斷一個(gè)問(wèn)題是否適合用貪心法求解，目前還沒(méi)有一個(gè)通用的方法，在競(jìng)賽中，需要憑個(gè)人的經(jīng)驗(yàn)來(lái)判斷 % 貪心算法畢竟是貪心算法，只能緩一時(shí)，而不是長(zhǎng)久之計(jì)， % 問(wèn)題的模型、參數(shù)對(duì)貪心算法求解結(jié)果具有決定性作用，這在某種程度上是不能接受的 % 于是聰明的人類就發(fā)明了各種智能算法（也叫啟發(fā)式算法） % 但在我看來(lái)所謂的智能算法本質(zhì)上就是貪心算法和隨機(jī)化算法結(jié)合 % 例如傳統(tǒng)遺傳算法用的選擇策略就是典型的貪心選擇，正是這些貪心算法和隨機(jī)算法的結(jié)合，我們才看到今天各種各樣的智能算法 %% 例1 設(shè)有n個(gè)正整數(shù)，將它們連接成一排，組成一個(gè)最大的多位整數(shù)。 % n=3時(shí)，3個(gè)整數(shù)13，312，343，連成的最大整數(shù)為34331213。 % n=4時(shí)，4個(gè)整數(shù)7，13，4，246，連成的最大整數(shù)為7424613。 % 此題很容易想到使用貪心法，比如把整數(shù)按從大到小的順序連接起來(lái)，例子符合，但測(cè)試結(jié)果不全對(duì) % 反例：12，121應(yīng)該組成12121而非12112 % 正確的標(biāo)準(zhǔn)是：先把整數(shù)轉(zhuǎn)換成字符串，然后比較a+b和b+a，如果a+b>=b+a，就把a(bǔ)排在b的前面，反之則把a(bǔ)排在b的后面 %% 例2 最短路徑問(wèn)題 clear,clc;close all; n=4; % 設(shè)置點(diǎn)數(shù) x = zeros(1,n); % 產(chǎn)生一個(gè)行向量存儲(chǔ)坐標(biāo) y = zeros(1,n); luxian = 1:1:n; % 生成1到n的矩陣，存儲(chǔ)路線 paixu = 1:1:n;% 這一步是作出隨機(jī)點(diǎn)的圖 for i = 1 : nx(i) = randn * n ; %生成正態(tài)分布的隨機(jī)坐標(biāo)點(diǎn)y(i) = randn * n ;hold ontext(x(i)+0.3,y(i)+0.3,num2str(i)) %用text做好標(biāo)記， end% 這一步是計(jì)算點(diǎn)與點(diǎn)之間的距離 d = zeros(n) ; % 生成一個(gè)n*n的矩陣用來(lái)存儲(chǔ)點(diǎn)與點(diǎ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è)置起點(diǎn)，路線的起點(diǎn)是1 num = 1; for a = 1:(n-2) %去除起點(diǎn)和終點(diǎn)需要n-2次判斷 paixu(:,1)=[]; %刪除上一個(gè)最優(yōu)點(diǎn) dis = zeros(1,(n-a)); %用來(lái)存剩下各個(gè)點(diǎn)的距離for b = 1:(n-a) %用來(lái)獲取剩下各個(gè)點(diǎn)的距離dis(b) = d (num, paixu(b));end num1 = find( dis == min(dis) ); %根據(jù)距離最小得到最優(yōu)點(diǎn)位置 t = paixu(1); %通過(guò)中間變量互換位置paixu(1) = paixu(num1); %最優(yōu)點(diǎn)位置換到第一個(gè)paixu(num1) = t;num = paixu(1); %獲取下次進(jìn)行操作的數(shù)luxian(a+1) = paixu(1); %將最優(yōu)點(diǎn)存入luxian數(shù)組（空出開(kāi)頭） end luxian(n) = paixu(2); %補(bǔ)上最后一個(gè)點(diǎn)plot(x(luxian),y(luxian),'--o') ; %標(biāo)出點(diǎn)用虛線相連得到路徑 grid on;title('貪心算法計(jì)算最優(yōu)路徑');19. 小波分析
%%初始化程序 clear,clc t1=clock;%% 載入噪聲信號(hào)數(shù)據(jù)，數(shù)據(jù)為.mat格式，并且和程序放置在同一個(gè)文件夾下 x = 0:0.06:10; YSJ = sin(x)+0.2*rand(size(x));plot(YSJ); title('原始數(shù)據(jù)');%% 數(shù)據(jù)預(yù)處理，數(shù)據(jù)可能是存儲(chǔ)在矩陣或者是EXCEL中的二維數(shù)據(jù)，銜接為一維的，如果數(shù)據(jù)是一維數(shù)據(jù)，此步驟也不會(huì)影響數(shù)據(jù) [c,l]=size(YSJ); Y=[]; for i=1:cY=[Y,YSJ(i,:)]; end [c1,l1]=size(Y); X=[1:l1];%% 繪制噪聲信號(hào)圖像 figure(1); plot(X,Y); xlabel('橫坐標(biāo)'); ylabel('縱坐標(biāo)'); title('原始信號(hào)');%% 硬閾值處理 lev=3; xd=wden(Y,'heursure','h','one',lev,'db4');%硬閾值去噪處理后的信號(hào)序列 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');%軟閾值去噪處理后的信號(hào)序列 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');%固定閾值去噪處理后的信號(hào)序列 figure(4) plot(X,xz); xlabel('橫坐標(biāo)'); ylabel('縱坐標(biāo)'); title('固定閾值后的去噪處理') set(gcf,'Color',[1 1 1]) %% 計(jì)算信噪比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); %% 計(jì)算均方根誤差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(['------------------運(yùn)行耗時(shí)',num2str(tim),'秒-------------------'])20. 遺傳算法
我是非常不建議使用這個(gè)的，太復(fù)雜太麻煩了，每次換算遺傳因子都要想半天，而且代碼冗余，又長(zhǎng)又臭
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);%隨機(jī)產(chǎn)生一個(gè)交叉位置 %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
%對(duì)調(diào)函數(shù) exchange.mfunction [x,y]=exchange(x,y) temp=x; x=y; y=temp;endfit.m
%適應(yīng)度函數(shù)fit.m，每次迭代都要計(jì)算每個(gè)染色體在本種群內(nèi)部的優(yōu)先級(jí)別，類似歸一化參數(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
%染色體的路程代價(jià)函數(shù) mylength.m function len=myLength(D,p)%p是一個(gè)排列 [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
%連點(diǎn)畫(huà)圖函數(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 問(wèn)題 % TSP問(wèn)題即旅行商問(wèn)題 % 經(jīng)典的TSP可以描述為 % 一個(gè)商品推銷員要去若干個(gè)城市推銷商品，該推銷員從一個(gè)城市出發(fā)，需要經(jīng)過(guò)所有城市后，回到出發(fā)地 % 應(yīng)如何選擇行進(jìn)路線，以使總的行程最短。從圖論的角度來(lái)看，該問(wèn)題實(shí)質(zhì)是在一個(gè)帶權(quán)完全無(wú)向圖中，找一個(gè)權(quán)值最小的哈密爾頓回路clear,clc;close all; % 輸入?yún)?shù) M=100; %種群的個(gè)數(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); %%城市的個(gè)數(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);%隨機(jī)排列，比如[2 4 5 6 1 3] end %隨機(jī)選擇一個(gè)種群 R=popm(1,:);%初始化種群及其適應(yīng)函數(shù) fitness = zeros(M,1); len = zeros(M,1); for i=1:M % 計(jì)算每個(gè)染色體對(duì)應(yīng)的總長(zhǎ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); %%各次迭代的最小的種群的路徑總長(zhǎng) 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)代最佳個(gè)體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%畫(huà)出所有城市坐標(biāo) subplot(1,3,1); scatter(pos(:,1),pos(:,2),'rx');fprintf('最短路程是 '); for i=1:Nfprintf('%d ',R(i));%把R順序打印出來(lái) end fprintf('\n最短距離是 %d\n',minlen); title('點(diǎn)圖'); subplot(1,3,2); plot_route(pos,R); % 路程圖 title('路程圖'); subplot(1,3,3); plot(distance_min); % 最短距離隨迭代次數(shù)的變化 title('最短距離的變化');21. 蟻群算法
21.1. TSP 示例
yiyu%% 蟻群算法求解TSP問(wèn)題%% 清空環(huán)境變量 clear clc%% 導(dǎo)入數(shù)據(jù) load citys_data.mat%% 計(jì)算城市間相互距離 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); % 各代最佳路徑的長(zhǎng)度 Length_ave = zeros(iter_max,1); % 各代路徑的平均長(zhǎng)度%% 迭代尋找最佳路徑 while iter <= iter_max% 隨機(jī)產(chǎn)生各個(gè)螞蟻的起點(diǎn)城市start = zeros(m,1);for i = 1:mtemp = randperm(n);start(i) = temp(1);endTable(:,1) = start;% 構(gòu)建解空間citys_index = 1:n;% 逐個(gè)螞蟻路徑選擇for i = 1:m% 逐個(gè)城市路徑選擇for j = 2:ntabu = Table(i,1:(j - 1)); % 已訪問(wèn)的城市集合(禁忌表)allow_index = ~ismember(citys_index,tabu);allow = citys_index(allow_index); % 待訪問(wèn)的城市集合P = allow;% 計(jì)算城市間轉(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);% 輪盤(pán)賭法選擇下一個(gè)訪問(wèn)城市Pc = cumsum(P);target_index = find(Pc >= rand);target = allow(target_index(1));Table(i,j) = target;endend% 計(jì)算各個(gè)螞蟻的路徑距離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% 計(jì)算最短路徑距離及平均距離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);% 逐個(gè)螞蟻計(jì)算for i = 1:m% 逐個(gè)城市計(jì)算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),' 起點(diǎn)'); text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),' 終點(diǎn)'); 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('各代最短距離與平均距離對(duì)比')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. 元胞自動(dòng)機(jī)
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)研究的問(wèn)題涉及多個(gè)變量，并且變量間相關(guān)性明顯，即包含的信息有所重疊時(shí) % 可以考慮用主成分分析的方法，這樣更容易抓住事物的主要矛盾，使問(wèn)題簡(jiǎn)化 clc,clear;close all; % 與主成分分析相關(guān)的 MATLAB 函數(shù)有pcacov、princomp函數(shù)% pcacov函數(shù)用來(lái)根據(jù) 協(xié)方差矩陣 或 相關(guān)系數(shù)矩陣 進(jìn)行主成分分析 % [COEFF,latent,explained]=pcacov（V） % V 是總體或樣本的協(xié)方差矩陣或相關(guān)系數(shù)矩陣，對(duì)于p維總體，V是pxp的矩陣 % 輸出參數(shù) COEFF 是p個(gè)主成分的系數(shù)矩陣，它是pxp的矩陣，它的第i列是第i個(gè)主成分的系數(shù)向量 % 輸出參數(shù) latent 是p個(gè)主成分的方差構(gòu)成的向量，即V的p個(gè)特征值的大小(從大到小)構(gòu)成的向量 % 輸出參數(shù) explained 是p個(gè)主成分的貢獻(xiàn)率向量，已經(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,:) = {'特征值', '差值', '貢獻(xiàn)率', '累積貢獻(xiàn)率'}; result_1(2: a+1 ,1) = num2cell(latent); % diff函數(shù)式用于求導(dǎo)數(shù)和差分的. result_1(2: a ,2) = num2cell(-diff(latent)); % cumsum函數(shù)通常用于計(jì)算一個(gè)數(shù)組各行的累加值。 result_1(2: a+1 ,3:4) = num2cell([explained, cumsum(explained)]); % disp(result_1);% princomp函數(shù)用來(lái)根據(jù)樣本觀測(cè)值矩陣進(jìn)行主成分分析，其調(diào)用格式如下： % [COEFF,SCORE,latent,tsquare] = princomp(X,'econ') % 根據(jù)樣本觀測(cè)值矩陣X進(jìn)行主成分分析。輸入?yún)?shù)X是n行p列的矩陣，每一行對(duì)應(yīng)一個(gè)觀測(cè)(樣品)，每一列對(duì)應(yīng)一個(gè)變量 % 輸出參數(shù)COEFF是p個(gè)主成分析的系數(shù)矩陣，是pxp的矩陣，它的第i列對(duì)應(yīng)第i個(gè)主成分的系數(shù)向量 % 輸出參數(shù)SCORE是n個(gè)樣品的p個(gè)主成分得分矩陣，它是n行p列的矩陣 % 每一行對(duì)應(yīng)一個(gè)觀測(cè)，每一列對(duì)應(yīng)一個(gè)主成分，第i行第j列元素表示第i個(gè)樣品的第j個(gè)主成分得分 % SCORE與X是一一對(duì)應(yīng)的關(guān)系，是X在新坐標(biāo)系中的數(shù)據(jù)，可以通過(guò)X*系數(shù)矩陣得到 % 返回樣本協(xié)方差矩陣的特征值向量latent，它是由p個(gè)特征值構(gòu)成的列向量，其中特征值按降序排列。 % 返回一個(gè)包含n個(gè)元素的列向量tsquare，它的第i個(gè)元素的第i個(gè)觀測(cè)對(duì)應(yīng)的霍特林T^2統(tǒng)計(jì)量 % 表述了第i個(gè)觀測(cè)與數(shù)據(jù)集(樣本觀測(cè)矩陣)的中心之間的距離，可用來(lái)尋找遠(yuǎn)離中心的極端數(shù)據(jù) % 通過(guò)設(shè)置參數(shù)‘econ’參數(shù)，使得當(dāng) n <= p 時(shí) % 只返回latent中的前n-1個(gè)元素(去掉不必要的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); % 定義一個(gè)n+1行，4列的元胞數(shù)組 result_2(1,:) = {'特征值', '差值', '貢獻(xiàn)率', '累積貢獻(xiàn)率'}; explained = 100*latent/sum(latent);% 計(jì)算貢獻(xiàn)率 %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行到最后一行分別存放主成分的貢獻(xiàn)率和累積貢獻(xiàn)率 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); % 定義一個(gè)n+1行，4列的元胞數(shù)組 result_2(1,:) = {'特征值', '差值', '貢獻(xiàn)率', '累積貢獻(xiàn)率'}; explained = 100*latent/sum(latent);% 計(jì)算貢獻(xiàn)率 %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行到最后一行分別存放主成分的貢獻(xiàn)率和累積貢獻(xiàn)率 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ù)的估計(jì)區(qū)間 % R 殘差 % RINT 置信區(qū)間 % STATS 用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量。有4個(gè)數(shù)值：判定系數(shù)r2，F統(tǒng)計(jì)量觀測(cè)值，檢驗(yàn)的p的值，誤差方差的估計(jì) % r2越接近1，回歸方程越顯著；F>F1?α(k,n?k?1)時(shí)拒絕H0，F越大，回歸方程越顯著；p<α?xí)r拒絕H0 % ALPHA 顯著性水平（缺少時(shí)默認(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ù)測(cè)值 plot([2004:2019],num_3','ro-',years,y2,'b*'); legend('生態(tài)用水量', '主成分回歸'); xlabel('時(shí)間/年'); ylabel('城市生態(tài)用水量/億立方米'); disp('回歸分析統(tǒng)計(jì)量為'); disp([ 'R^2 : ',num2str(STATS(1))]); % 模型由 這么多原始數(shù)據(jù) 解釋，最優(yōu)為 1 disp([ 'F : ',num2str(STATS(2))]); disp([ '檢驗(yàn) 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ù)測(cè)值 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é)果','主成分回歸在測(cè)試集的結(jié)果'); xlswrite('save_4.xlsx',y3);25. 分治算法
25.1. 快速排序
QuickSort.m
function a=QuickSort(a,leftIndex,rightIndex) % a是待排序序列 %leftIndex和rightIndex是開(kāi)始的左右兩個(gè)邊界 %%----------------------------------------------------------%% % if leftIndex>rightIndex % break; % end if leftIndex<rightIndexi=leftIndex;j=rightIndex;temp=a(i);%選定的這個(gè)數(shù)挖掉，相當(dāng)于挖坑while i<jwhile (i<j)&&(a(j)>=temp)%從右往左，找到第一個(gè)小于設(shè)定的數(shù)，j=j-1;enda(i)=a(j);%將找到的第一個(gè)小于設(shè)定的數(shù)填坑到最開(kāi)始挖的坑里面去while (i<j)&&(a(i)<=temp)%從做到由，找到第一個(gè)大于選定的數(shù)i=i+1;enda(j)=a(i);%將找到的第一個(gè)大于選定的數(shù)填入上一步挖的坑里面去 % if i==j % a(j)=temp; % endenda(j)=temp;%最后，i=j，將選定的數(shù)再填到上一步挖的坑里面去a=QuickSort(a,leftIndex,j-1);%對(duì)左邊序列進(jìn)行遞歸a=QuickSort(a,i+1,rightIndex);%對(duì)右邊序列進(jìn)行遞歸 end end調(diào)用方法
% 分治算法，渡邊筆記 % 分治算法是一種思路，把大問(wèn)題歸納為小問(wèn)題來(lái)解決 % 使用分治法，問(wèn)題應(yīng)該滿足 % 該問(wèn)題的規(guī)模縮小到一定的程度就可以容易地解決 % 該問(wèn)題可以分解為若干個(gè)規(guī)模較小的相同問(wèn)題，即該問(wèn)題具有最優(yōu)子結(jié)構(gòu)性質(zhì) % 利用該問(wèn)題分解出的子問(wèn)題的解可以合并為該問(wèn)題的解 % 該問(wèn)題所分解出的各個(gè)子問(wèn)題是相互獨(dú)立的，即子問(wèn)題之間不包含公共的子子問(wèn)題 % 第一條特征是絕大多數(shù)問(wèn)題都可以滿足的，因?yàn)閱?wèn)題的計(jì)算復(fù)雜性一般是隨著問(wèn)題規(guī)模的增加而增加 % 第二條特征是應(yīng)用分治法的前提它也是大多數(shù)問(wèn)題可以滿足的，此特征反映了遞歸思想的應(yīng)用 % 第三條特征是關(guān)鍵，能否利用分治法完全取決于問(wèn)題是否具有第三條特征，如果具備了第一條和第二條特征，而不具備第三條特征，則可以考慮用貪心法或動(dòng)態(tài)規(guī)劃法。 % 第四條特征涉及到分治法的效率，如果各子問(wèn)題是不獨(dú)立的則分治法要做許多不必要的工作， % 重復(fù)地解公共的子問(wèn)題，此時(shí)雖然可用分治法，但一般用動(dòng)態(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. 最近點(diǎn)對(duì)
cloest.m
%分治求最近點(diǎn)對(duì) 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)用方法
% 分治法求最近點(diǎn)，渡邊筆記 clear;clc n = 20; % 隨機(jī)生成20個(gè)點(diǎn) A=rand(n,2)*10; % 將20個(gè)點(diǎn)按橫坐標(biāo)升序排列 A = sortrows(A,1); % 分治法求隨機(jī)點(diǎn)的最近點(diǎn)對(duì) [mindist1,y1,y2] = cloest(A,1,n); % 使用紅色的點(diǎn)標(biāo)記隨機(jī)點(diǎ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 %使用綠色十字標(biāo)記分治法的最近點(diǎn)對(duì) plot(y1(1),y1(2),'go', 'linewidth', 2);hold on plot(y2(1),y2(2),'go', 'linewidth', 2);hold on26. 雜項(xiàng)
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ǎn)單，不同變量往往量綱不同 % 歸一化可以消除量綱對(duì)最終結(jié)果的影響，使不同變量具有可比性 % 如兩個(gè)人體重差10KG，身高差0.02M % 在衡量?jī)蓚€(gè)人的差別時(shí)體重的差距會(huì)把身高的差距完全掩蓋，歸一化之后就不會(huì)有這樣的問(wèn)題 % 數(shù)據(jù)歸一化應(yīng)該針對(duì)對(duì)于的 y，而不是針對(duì)每條數(shù)據(jù)，針對(duì)每條數(shù)據(jù)是完全沒(méi)有意義的，因?yàn)橹皇堑缺壤s放，對(duì)之后的分類沒(méi)有任何作用 clc,clear;close all; x = [1 2 31 2 3]; % Min-max 極大極小標(biāo)準(zhǔn)化 % 當(dāng)有新數(shù)據(jù)加入時(shí)，可能導(dǎo)致xmax和xmin的變化，需要重新計(jì)算 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é)果沒(méi)有任何意義，只能用于比較 % 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 % 省會(huì)城市 順序表 % 北京 天津 石家莊 太原 呼和浩特 沈陽(yáng) 長(zhǎng)春 哈爾濱 上海 南京 % 杭州 合肥 福州 南昌 濟(jì)南 鄭州 武漢 長(zhǎng)沙 廣州 南寧 % 海口 重慶 成都 貴陽(yáng) 昆明 拉薩 西安 蘭州 西寧 銀川 % 烏魯木齊 filename_lujin_1 = 'C:\Users\99791\Desktop\生態(tài)用水\數(shù)據(jù)\國(guó)家統(tǒng)計(jì)局\分省\'; % 文件的前面的路徑名稱 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'; % 改變寫(xiě)入 列 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 強(qiáng)，下面的方法好久不用了
% function 函數(shù)解多元方程 最小值 % function 基本語(yǔ)法 % [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 定義了線性約束 ,如果沒(méi)有線性約束，則A=[],b=[],Aeq=[],beq=[] % lb和ub是變量x的下界和上界，如果下界和上界沒(méi)有約束，則lb=[],ub=[],也可以寫(xiě)成lb的各分量都為 -inf,ub的各分量都為inf % nonlcon是用 M 文件定義的非線性向量函數(shù)約束 % options 定義了優(yōu)化參數(shù)，不填寫(xiě)表示使用 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)的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： ebs克隆oracle not,Orac
- 下一篇： python分析北京租房现状，最后的价格