數(shù)學(xué)建模常用方法及MATLAB代碼
二分法
我們通常使用二分法計算非線性方程或者超越方程近似根，MATLAB代碼為：
// 二分法求根
function [x0,n]=dichotomy(a,b,err,f_x)
%輸入?yún)?shù)a為根的區(qū)間左端點
%輸入?yún)?shù)b為根的區(qū)間右端點
%輸入?yún)?shù)err為誤差精度
% 輸入?yún)?shù)f_x為待求根函數(shù)
%輸出參數(shù)x0為滿足精度要求的根
% 輸出參數(shù)n為迭代的次數(shù)
n=ceil((log((b-a)./err)./log(2)-1));%求迭代次數(shù)
while(sign(f_x(a))==0)%如果a的函數(shù)值為0停止迭代輸出x0為a的值，迭代次數(shù)為0
x0=a;
n=0;
return;
end
while(sign(f_x(b))==0)%如果b的函數(shù)值為0停止迭代輸出x0為b的值，迭代次數(shù)為0
x0=b;
n=0;
return;
end
while(sign(f_x(a))~=sign(f_x(b)))%a、b符號不同進(jìn)行區(qū)間二分
x0=a/2+b/2; %區(qū)間中點
if(sign(f_x(a))~=sign(f_x(x0)))%判斷區(qū)間中點函數(shù)值與區(qū)間端點符號
b=x0; %將中點賦值給符號相同的端點
if(b-a<=err||abs(f_x(x0))<=err)%判斷是否滿足精度要求
x0=a/2+b/2;
return;
else %不滿足精度循環(huán)調(diào)用二分法
dichotomy(a,b,err,f_x);
end
else
a=x0; %將中點賦值給符號相同的端點
if(b-a<=err||abs(f_x(x0))<=err)%判斷是否滿足精度要求
x0=a/2+b/2;
return;
else %不滿足精度循環(huán)調(diào)用二分法
dichotomy(a,b,err,f_x);
end
end
end
量綱分析法
π定理的解題步驟:
(1)確定關(guān)系式：根據(jù)對所研究的現(xiàn)象的認(rèn)識，確定影響這個現(xiàn)象的各個物理量及其關(guān)系式:
(2)確定基本量：從n個物理量中選取所包含的m個基本物理量作為基本量綱的代表，一般取m=3。在管流中，一般選d，v，ρ三個作基本變量，而在明渠流中，則常選用H，v，ρ。
(3)確定π數(shù)的個數(shù)N(π)=(n-m)，并寫出其余物理量與基本物理量組成的π表達(dá)式
(4)確定無量綱π參數(shù):由量綱和諧原理解聯(lián)立指數(shù)方程，求出各π項的指數(shù)x，y，z，從而定出各無量綱π參數(shù)。π參數(shù)分子分母可以相互交換，也可以開方或乘方，而不改變其無因次的性質(zhì)。
(5)寫出描述現(xiàn)象的關(guān)系式或顯解一個π參數(shù)，或求得一個因變量的表達(dá)式。
選擇基本量時的注意原則:
1)基本變量與基本量綱相對應(yīng)。即若基本量綱(M，L，T)為三個，那么基本變量也選擇三個;倘若基本量綱只出現(xiàn)兩個，則基本變量同樣只須選擇兩個。
2)選擇基本變量時，應(yīng)選擇重要的變量。換句話說，不要選擇次要的變量作為基本變量，否則次要的變量在大多數(shù)項中出現(xiàn)，往往使問題復(fù)雜化，甚至要重新求解。
3)不能有任何兩個基本變量的因次是完全一樣的，換言之，基本變量應(yīng)在每組量綱中只能選擇一個。
圖論法
其MATLAB代碼為：
// 圖論算法
function P = dgraf( A )
%A為圖的鄰接矩陣
%P為圖的可達(dá)矩陣
n=size(A,1);
P=A;
for i=2:n
P=P+A^i;
end
P(P~=0)=1;
P;
end
應(yīng)用舉例：
A=[0 1 1 1;1 0 1 1;1 1 0 1;1 1 10];
P=dgraf(A);
P=[1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1];
差分法
背景差分法是采用圖像序列中的當(dāng)前幀和背景參考模型比較來檢測運動物體的一種方法,其性能依賴于所使用的背景建模技術(shù)。背景構(gòu)建的方法有多種，簡單的有均值法、中值法，復(fù)雜點的有卡爾曼濾波器模型法、單高斯分布模型法、雙高斯分布模型法等，這里我用的是均值法。以下是相應(yīng)的matlab代碼和輸出結(jié)果：
// 差分法
%背景差分法
clear all;
clc;
avi=MMREADER('I:\\電影\\test.avi');
VidFrames=read(avi,[10000,10050]);
N=4;
start=11;
threshold=15;
bg(start).cdata=0;
i=1;
for k=start-10:start-1 bg(k).cdata=rgb2gray(VidFrames(:,:,:,k)); bg(start).cdata=abs((bg(start).cdata+bg(k).cdata)/i); %均值法構(gòu)建背景
i=i+1;
end
for k=1+start:N+1+start mov(k).cdata=rgb2gray(VidFrames(:,:,:,k)); %轉(zhuǎn)化成灰度圖
end
[row,col]=size(mov(1+start).cdata);
alldiff=zeros(row,col,N);
bgpic=zeros(row,col,1);
bgdata=bg(start).cdata>threshold;
bgpic(:,:,1)=double(bgdata);
figure(1);
imshow(bgpic(:,:,1)) %輸出構(gòu)建的背景
for k=1+start:N+start
diff=abs(mov(k).cdata-bg(start).cdata);
idiff=diff>threshold;
alldiff(:,:,k)=double(idiff);
end
j=2;
for k=1+start:N+start
figure(j);
imshow(alldiff(:,:,k)) %輸出測試幀減去構(gòu)建背景的效果
title(strcat(num2str(k),'幀','-','背景'));
j=j+1;
end
變分法
//變分法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I=imread('toys.bmp'); % load image
I=double(I(20:120,10:105)); % cut a piece, convert to double
%%% Parameters
std_n=10; var_n=std_n^2; % Gaussian noise standard deviation
reduced_pw = 1.5*var_n; % power to reduce in first phase
sig_w = 5; ws=4*sig_w+1; % window size
%%%%%%%%%%%%%%
%%% Add noise
In = randn(size(I))*std_n;
I0 = I + In; % noisy input image
% show original and noisy images
figure(1); imshow(uint8(I)); title('Original')
figure(2); imshow(uint8(I0)); title('Noisy image')
snr_noisy=db(I,I0)
% run normal tv - strong denoising
tic;
J=I0;
ep_J=0.1; % minimum mean change in image J
J_old=0;
lam=0; iter=10; dt=0.2; ep=1;
while (mean(mean(abs(J - J_old))) > ep_J), % iterate until convergence
J_old=J;
J=tv(J,iter,dt,ep,lam,I0); % scalar lam
lam = calc_lam(J,I0,reduced_pw); % update lambda (fidelity term)
end % while
% figure(3); imshow(uint8(J)); title('residue TV')
% snr_residue= db(I,J)
Ir=I0-J; % Ir scalar
Pr = mean(mean(Ir.^2)); % power of residue
LV = loc_var(Ir,ws,sig_w^2); % local variance (local power of the residue )
% P=mean(mean(LV));
Pxy=1*(var_n^2)./LV; %% Sxy inverse proportional to the LV
%%% Varying Lambda
lamxy=zeros(size(I0));
J=I0; J_old=0;
ep_J=0.001;
%eps=0.01;
while (mean(mean(abs(J - J_old))) > ep_J), % iterate until convergence
J_old=J;
J=tv(J,iter,dt,ep,lamxy,I0); % adaptive lam
%J=tv(J,iter,dt,ep_J,lamxy,I0);
lamxy = calc_lamxy(J,I0,Pxy,sig_w); % update lambda (fidelity term)
end % while
figure(3); imshow(uint8(J)); title('Adaptive TV')
snr_adap= db(I,J)
toc
t1=toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run scalar TV denoising for comparision
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
J=I0;
% params
ep_J = 0.01; % minimum mean change in image J
lam=0; J_old=0;
i=0;
while (mean(mean(abs(J - J_old))) > ep_J), % iterate until convergence
J_old = J;
J=tv(J,iter,dt,ep,lam,I0); % scalar lam
lam = calc_lam(J,I0,var_n,ep); % update lambda (fidelity term)
end % for i
% Ir=I0-J; % Ir scalar
% Pr = mean(mean(Ir.^2)); % power of residue
function Ig=gauss(I,ks,sigma2)
%private function: gauss (by Guy Gilboa):
% Ig=gauss(I,ks,sigma2)
% ks - kernel size (odd number)
% sigma2 - variance of Gaussian
[Ny,Nx]=size(I);
hks=(ks-1)/2; % half kernel size
if (Ny<ks) % 1d convolutin
x=(-hks:hks); %x:1x ks
flt=exp(-(x.^2)/(2*sigma2)); % 1D gaussian
flt=flt/sum(sum(flt)); % normalize
% expand
x0=mean(I(:,1:hks)); xn=mean(I(:,Nx-hks+1:Nx));%x0,xn :1 x hks;
eI=[x0*ones(Ny,ks) I xn*ones(Ny,ks)]; % ???
Ig=conv(eI,flt);
Ig=Ig(:,ks+hks+1:Nx+ks+hks); % truncate tails of convolution
else
%% 2-d convolution
x=ones(ks,1)*(-hks:hks); y=x';
flt=exp(-(x.^2+y.^2)/(2*sigma2)); % 2D gaussian
flt=flt/sum(sum(flt)); % normalize
% expand
if (hks>1)
xL=mean(I(:,1:hks)')'; xR=mean(I(:,Nx-hks+1:Nx)')';% xL,xR :Ny x 1
else
xL=I(:,1); xR=I(:,Nx);
end
eI=[xL*ones(1,hks) I xR*ones(1,hks)]; % Ny x Nx+2hks
if (hks>1)
xU=mean(eI(1:hks,:)); xD=mean(eI(Ny-hks+1:Ny,:)); % xU,xD: 1 x Nx+2hks
else
xU=eI(1,:); xD=eI(Ny,:);
end
數(shù)據(jù)擬合法
以下以一個較為簡單的數(shù)據(jù)組作為示例：
//數(shù)據(jù)擬合法
clc;
clear;
% the first one;
R=[0.68 0.805 0.863 0.893 0.9122 0.916];
x=[1.07 2.24 3.86 5.91 8.5 11.5];
y=R./(1-R);
plot(x,y,'ro')
hold on
x=x'; % transfer Row matrix to Column matrix;
y=y';
p=fittype('(x/a)*exp(-x/b)') % Fitting function
f=fit(x,y,p)
plot(f,x,y);
回歸分析法
在 Matlab 中，可以直接調(diào)用命令實現(xiàn)回歸分析，
（ 1 ） [b,bint,r,rint,stats]=regress(y,x) ，其中 b 是回歸方程中的參數(shù)估計值， bint 是 b 的置信區(qū)間， r 和 rint 分別表示殘差及殘差對應(yīng)的置信區(qū)間。stats 包含三個數(shù)字，分別是相關(guān)系數(shù)， F 統(tǒng)計量及對應(yīng)的概率 p 值。
（ 2 ） recplot(r,rint) 作殘差分析圖。
（ 3 ） rstool(x,y) 一種交互式方式的句柄命令。
例：現(xiàn)有多個樣本的因變量和自變量的數(shù)據(jù)，下面我們利用 Matlab ，通過回歸分析建立兩者之間的回歸方程。
//一元線性回歸分析
% 輸入數(shù)據(jù)
x = [1, 2, 3, 4, 5]; % 自變量
y = [2, 3, 4, 3, 5]; % 因變量
% 計算回歸系數(shù)
n = length(x); % 數(shù)據(jù)點個數(shù)
xy = x .* y;
xx = x .* x;
b1 = (n * sum(xy) - sum(x) * sum(y)) / (n * sum(xx) - sum(x)^2);
b0 = mean(y) - b1 * mean(x);
% 計算預(yù)測值和殘差
y_pred = b0 + b1 * x; % 預(yù)測值
residuals = y - y_pred; % 殘差
% 計算置信區(qū)間
X = [ones(n, 1), x'];
[~, ~, ~, ~, stats] = regress(y', X); % 使用regress函數(shù)計算置信區(qū)間
CI = stats(1:2); % 置信區(qū)間上界和下界
% 繪制數(shù)據(jù)點、回歸線和置信區(qū)間
subplot(2,1,1);
scatter(x, y); % 繪制散點圖
hold on;
plot(x, y_pred, 'r'); % 繪制回歸線
fill([min(x), max(x), max(x), min(x)], [CI(1), CI(1), CI(2), CI(2)], 'g', 'FaceAlpha', 0.2); % 填充置信區(qū)間
xlabel('x');
ylabel('y');
legend('數(shù)據(jù)點', '回歸線', '置信區(qū)間');
% 繪制殘差圖
subplot(2,1,2);
scatter(y_pred, residuals); % 繪制殘差圖
hold on;
plot([min(y_pred), max(y_pred)], [0, 0], 'r--'); % 繪制零線
xlabel('預(yù)測值');
ylabel('殘差');
legend('殘差', '零線');
% 輸出回歸系數(shù)和置信區(qū)間
disp(['回歸系數(shù) b0 = ', num2str(b0)]);
disp(['回歸系數(shù) b1 = ', num2str(b1)]);
disp(['置信區(qū)間 = [', num2str(CI(1)), ', ', num2str(CI(2)), ']']);
其輸出結(jié)果為：
//輸出結(jié)果
回歸系數(shù) b0 = 1.6
回歸系數(shù) b1 = 0.6
置信區(qū)間 = [0.69231, 6.75]
多元回歸分析
x1=[5.5 2.5 8 3 3 2.9 8 9 4 6.5 5.5 5 6 5 3.5 8 6 4 7.5 7]';
x2=[31 55 67 50 38 71 30 56 42 73 60 44 50 39 55 70 40 50 62 59]';
x3=[10 8 12 7 8 12 12 5 8 5 11 12 6 10 10 6 11 11 9 9]';
x4=[8 6 9 16 15 17 8 10 4 16 7 12 6 4 4 14 6 8 13 11]';
% 輸入因變量數(shù)據(jù)
y=[79.3 200.1 163.1 200.1 146.0 177.7 30.9 291.9 160 339.4 159.6 86.3 237.5 107.2 155 201.4 100.2 135.8 223.3 195]';
X=[ones(size(x1)),x1,x2,x3,x4];
[b,bint,r,rint,stats]=regress(y,X)% 回歸分析
Q=r'*r
sigma=Q/18
rcoplot(r,rint);
結(jié)果為：
回歸系數(shù) b0 = 1.6
回歸系數(shù) b1 = 0.6
置信區(qū)間 = [0.69231, 6.75]
>> untitled
b =
191.9906
-0.7751
3.1718
-19.6849
-0.4494
bint =
103.1866 280.7946
-7.1474 5.5972
2.0634 4.2802
-25.1686 -14.2012
-3.7276 2.8288
r =
-6.3088
-4.2260
5.0635
-3.1701
0.0269
6.6185
-10.2294
32.1836
-2.8287
26.5217
1.2446
0.2377
12.3762
-5.9684
-10.0799
-82.0142
5.2187
8.4492
23.4779
3.4069
rint =
-56.5100 43.8924
-51.7978 43.3459
-40.5391 50.6662
-49.7192 43.3789
-44.3792 44.4331
-33.2139 46.4509
-53.1383 32.6796
-9.9179 74.2851
-52.7738 47.1163
-17.9441 70.9875
-50.2831 52.7722
-50.3838 50.8593
-37.1054 61.8578
-57.6997 45.7628
-58.2285 38.0687
-94.9997 -69.0286
-47.1632 57.6005
-44.4731 61.3716
-26.9733 73.9292
-50.3016 57.1154
stats =
0.9034 35.0555 0.0000 644.5831
Q =
9.6687e+03
sigma =
537.1526
逐步回歸
X1=[x1,x2,x3,x4];
stepwise(X1,y,[1,2,3])% 逐步回歸
% X2=[ones(size(x1)),x2,x3];
% X3=[ones(size(x1)),x1,x2,x3];
% X4=[ones(size(x1)),x2,x3,x4];
% [b1,b1int,r1,r1int,stats1]=regress(y,X2)
% [b2,b2int,r2,r2int,stats2]=regress(y,X3);
% [b3,b3int,r3,r3int,stats3]=regress(y,X4);
機理分析
- 機理方法概述
（1）自身存在的發(fā)展規(guī)律和理由–分析事物的內(nèi)在因素，研究其內(nèi)在關(guān)系，得到內(nèi)在規(guī)律—機理模型。
（2）如何從事物的內(nèi)在因素和條件中研究其內(nèi)在關(guān)系與規(guī)律？？------機理分析建模方法。
（3）利用機理分析方法所建立的模型有：代數(shù)方程、函數(shù)方程、微分方程、積分方程和一般的動力學(xué)模型等。
一般不能直接應(yīng)用某種現(xiàn)成的方法得到模型，或直接套用現(xiàn)成的模型得到結(jié)論。 - 機理建模的基本方法
類比分析法：根據(jù)一些物理定律，經(jīng)濟規(guī)律，數(shù)學(xué)原理等建立不同事物之間的類比關(guān)系，建立問題的數(shù)學(xué)模型。
量綱分析法：通過分析問題相關(guān)物理量的量綱，根據(jù)量綱一致性原則建立各物理量之間的關(guān)系。
幾何分析法：針對實際問題，利用平面幾何、立體幾何、解析幾何的原理等建立模型。
邏輯分析法：一句問題的客觀條件和實際情況，利用邏輯推理和邏輯運算建立模型
比較分析法：對照各個事物，確定事務(wù)間的共同點和差異點，通過文字描述、圖表等方式對事物特征進(jìn)行分析，建立模型。
推理分析法：在掌握一定的已知事實，數(shù)據(jù)信息或者因素相關(guān)性的基礎(chǔ)上，通過因果關(guān)系或其他相關(guān)關(guān)系順次，逐步地推論得出新結(jié)論，建立模型。 - 機理分析建模流程
針對實際問題-----了解問題背景----分析問題-----明確相關(guān)因素和參數(shù)-----分析其內(nèi)在關(guān)系—用適當(dāng)數(shù)學(xué)方法—建立關(guān)聯(lián)模型–選用實際數(shù)據(jù)–確定未知數(shù)據(jù)—求解模型—用結(jié)果解釋實際問題—用實際數(shù)據(jù)或模擬檢驗?zāi)Ｐ汀M(jìn)一步擴展模型。
例：人才吸引力評價（2018深圳杯A題）—中國大學(xué)生在線可看此論文
將定性分析轉(zhuǎn)換成定量分析，比如：發(fā)展前景、經(jīng)濟收入、社會環(huán)境
問題的關(guān)鍵：首先搞清楚：什么是人才吸引力？水平如何？優(yōu)勢與不足？如何提升？
（1）柯布–道格拉斯函數(shù)
人才吸引力：發(fā)展前景、經(jīng)濟收入、社會環(huán)境
發(fā)展前景：歷史數(shù)據(jù)、當(dāng)前數(shù)據(jù)
經(jīng)濟收入：人均收入、物價水平
社會環(huán)境：空氣環(huán)境指數(shù)、
（2）萬有引力定律
（3）類歐姆定律
例：出租車模型
（1）建立合理性模型，并分析不同時空出租車資源的供求匹配程度。
（2）分析相關(guān)公司的出租車補貼方案是否對緩解打車難問題有幫助？
（3）試創(chuàng)建一個打車軟件服務(wù)平臺，計一個補貼方案，并論證其合理性
深圳杯的題目----深圳一天的數(shù)據(jù) 都可以查到
關(guān)系的問題===核心問題
基例—情景—數(shù)據(jù)
1.目前出租車資源的配置是否合理？
2.如何實現(xiàn)出租車供求關(guān)系的良好匹配？
3.打車軟件是否能有效解決打車難問題？
解決問題：分析與出租車供求相關(guān)的因素和關(guān)系機理，建立問題的機理模型。
需要收集某區(qū)域或城市出租車的相關(guān)數(shù)據(jù)，通過對數(shù)據(jù)的分析研究，統(tǒng)計挖掘出相關(guān)規(guī)律，來支持所建立的數(shù)學(xué)哦行和模型的結(jié)論。
出租車運行軌跡數(shù)據(jù)----------------不需要
優(yōu)秀論文：-----------中國大學(xué)生在線可看到
（1）匹配指標(biāo)如何建立？出租車?yán)锍汤寐屎凸┣蟊嚷实慕?jīng)驗公式
（2）某城市的滴滴數(shù)據(jù)補貼的前后數(shù)據(jù)分析，客觀形成好單和壞單要分階段進(jìn)行討論考慮，而不是單單說有利還是無利。
通過機理分析建立分區(qū)動態(tài)實時補貼模型，給出符合實際的補貼方案：針對司機和乘客不同的設(shè)計。
在各區(qū)各時段的補貼保持平衡的前提下，建立平衡補貼的方程組模型。
排隊方法
這個就直接上代碼了。
//排隊論
% /M/PH/1(k);
% p=stationary_prob(4,beta,S,order,lambda);
function [p,p_minus,p_add,mean,time]=stationary_prob(k,beta,S,lambda) order=length(beta);
e=ones(order,1);
R=lambda*inv(lambda*eye(size(S))-lambda*e*beta-S);
R_k=eye(size(S));
for i=1:k
R_k = R_k+R^i;
end
p0=inv(beta*(R_k-lambda*R^k*inv(S))*e);
p_t=p0;
for i=1:k
p_t=[p_t,p0*beta*R^i];
end
p_t=[p_t,p0*beta*(R^k)*(-lambda*inv(S))];
% the queue length distribution in any time
% add every phase in the same level
p=p0;
for i=1:k+1
sum=0;
for j=1:order
sum =sum+p_t(1+(i-1)*order+j);
end
p=[p,sum];
end
p_minus = p;
% the queue length distribution at departure
% p_add = p_minus(i)/(1-p_minus(k+1))
p_add=[0];
for i=0:k
p_add=[p(k-i+1)/(1-p(k+2)),p_add];
end
% the mean queue length
mean=0;
for i=1:k
mean =mean+i*p(i+1);
end
層次分析法(AHP)
基本步驟：
1.建立層次分析結(jié)構(gòu)模型。常見的有目標(biāo)層-準(zhǔn)則層-方案層模型。
2.構(gòu)造成對比較矩陣。常用的有1-9尺度。尺度大小取決于下層的元素個數(shù)。
3.計算權(quán)向量并作一致性檢驗。引入一致性指標(biāo)CI，CI越大，不一致越嚴(yán)重。引入隨機一致性指標(biāo)RI。定義一致性比率CR=CI/RI。CR<0.1時，通過一致性檢驗。
4.計算組合權(quán)向量作組合一致性檢驗，組合權(quán)向量可作為決策的定量依據(jù)。方案層對目標(biāo)層的組合權(quán)向量為ww*w。選擇組合權(quán)向量大的元素作為輸出結(jié)果。
注意，該方法具有較強的主觀性
//AHP代碼
clc
a=[1,1,1,4,1,1/2
1,1,2,4,1,1/2
1,1/2,1,5,3,1/2
1/4,1/4,1/5,1,1/3,1/3
1,1,1/3,3,1,1
2,2,2,3,3,1];%一致矩陣
[x,y]=eig(a);eigenvalue=diag(y);lamda=max(eigenvalue);
ci1=(lamda-6)/5;cr1=ci1/1.24
w1=x(:,1)/sum(x(:,1))
b1=[1,1/4,1/2;4,1,3;2,1/3,1];
[x,y]=eig(b1);eigenvalue=diag(y);lamda=eigenvalue(1);
ci21=(lamda-3)/2;cr21=ci21/0.58
w21=x(:,1)/sum(x(:,1))
b2=[1 1/4 1/5;4 1 1/2;5 2 1];
[x,y]=eig(b2);eigenvalue=diag(y);lamda=eigenvalue(1);
ci22=(lamda-3)/2;cr22=ci22/0.58
w22=x(:,1)/sum(x(:,1))
b3=[1 3 1/3;1/3 1 1/7;3 7 1];
[x,y]=eig(b3);eigenvalue=diag(y);lamda=eigenvalue(1);
ci23=(lamda-3)/2;cr23=ci23/0.58
w23=x(:,1)/sum(x(:,1))
b4=[1 1/3 5;3 1 7;1/5 1/7 1];
[x,y]=eig(b4);eigenvalue=diag(y);lamda=eigenvalue(1);
ci24=(lamda-3)/2;cr24=ci24/0.58
w24=x(:,1)/sum(x(:,1))
b5=[1 1 7;1 1 7;1/7 1/7 1];
[x,y]=eig(b5);eigenvalue=diag(y);lamda=eigenvalue(1);
ci25=(lamda-3)/2;cr25=ci25/0.58
w25=x(:,1)/sum(x(:,1))
b6=[1 7 9;1/7 1 1 ;1/9 1 1];
[x,y]=eig(b6);eigenvalue=diag(y);lamda=eigenvalue(1);
ci26=(lamda-3)/2;cr26=ci26/0.58
w26=x(:,1)/sum(x(:,1))
w_sum=[w21,w22,w23,w24,w25,w26]*w1
ci=[ci21,ci22,ci23,ci24,ci25,ci26];
cr=ci*w1/sum(0.58*w1)
結(jié)果為：
cr1 =
0.0996
w1 =
0.1507
0.1792
0.1886
0.0472
0.1464
0.2879
cr21 =
0.0158
w21 =
0.1365
0.6250
0.2385
cr22 =
0.0212
w22 =
0.0974
0.3331
0.5695
cr23 =
0.0061
w23 =
0.2426
0.0879
0.6694
cr24 =
0.0559
w24 =
0.2790
0.6491
0.0719
cr25 =
-3.8284e-16
w25 =
0.4667
0.4667
0.0667
cr26 =
0.0061
w26 =
0.7986
0.1049
0.0965
w_sum =
0.3952
0.2996
0.3052
cr =
0.0117
主成分分析法
主成分分析法是利用降維的思想，把多指標(biāo)轉(zhuǎn)化為少數(shù)幾個綜合指標(biāo)（即主成分），其中每個主成分都能夠反映原始變量的大部分信息，且所含信息互不重復(fù)。這種方法在引進(jìn)多方面變量的同時將復(fù)雜因素歸結(jié)為幾個主成分，使問題簡單化，同時得到的結(jié)果更加科學(xué)有效的數(shù)據(jù)信息。
例如，做一件上衣需要測量很多尺寸，如身高，袖長，腰圍，胸圍，肩寬等十幾項指標(biāo)，但是廠家不可能把尺寸型號分這么多，而是從這十幾種指標(biāo)中綜合成幾個少數(shù)的綜合指標(biāo)作為分類型號，例如綜合成反映身高，反映胖瘦和反應(yīng)特體的三項指標(biāo)，這就是主成分的思想。主要的方法有特征值分解，SVD(奇異值分解)和NMF(非負(fù)矩陣分解)；
原始指標(biāo)數(shù)據(jù)的標(biāo)準(zhǔn)化采集p維隨機向量
x
=
(
X
i
1
,
X
i
2
,
?
?
,
X
i
p
)
T
x=(X_{i1},X_{i2},\cdots,X_{ip})^{T}
x=(Xi1?,Xi2?,?,Xip?)Tn個樣品
x
i
=
(
x
i
1
,
x
i
2
,
?
?
,
x
i
p
T
)
x_{i}=(x_{i1},x_{i2},\cdots,x_{ip}^{T})
xi?=(xi1?,xi2?,?,xipT?),i=1,2,…n,n>p,構(gòu)造樣本陣,對樣本陣元進(jìn)行如下變化標(biāo)準(zhǔn)化變換：
Z
i
j
=
x
i
j
?
x
￣
j
s
j
,
i
=
1
,
2
,
?
?
,
n
;
j
=
1
,
2
,
?
?
,
p
Z_{ij}=\frac{x_{ij}-\overline{x}_{j}}{s_{j}},i=1,2,\cdots,n;j=1,2,\cdots,p
Zij?=sj?xij??xj??,i=1,2,?,n;j=1,2,?,p
其中
x
￣
j
=
∑
i
=
1
n
x
i
j
n
,
s
j
2
=
∑
i
=
1
n
(
x
i
j
?
x
￣
j
)
2
n
?
1
\overline{x}_{j}=\frac{\sum_{i=1}^{n}x_{ij}}{n},s_{j}^{2}=\frac{\sum_{i=1}^{n}(x_{ij}-\overline{x}_{j})^{2}}{n-1}
xj?=n∑i=1n?xij??,sj2?=n?1∑i=1n?(xij??xj?)2?得標(biāo)準(zhǔn)化陣Z。
2、對標(biāo)準(zhǔn)化陣Z求相關(guān)系數(shù)矩陣
R
=
[
r
i
j
]
p
x
p
=
Z
T
Z
n
?
1
R=[r_{ij}]_{p}xp=\frac{Z^{T}Z}{n-1}
R=[rij?]p?xp=n?1ZTZ?
其中，
r
i
j
=
∑
z
k
j
?
z
k
j
n
?
1
,
i
,
j
=
1
,
2
,
?
?
,
p
r_{ij}=\frac{\sum{z_{kj}·z_{kj}}}{n-1},i,j=1,2,\cdots,p
rij?=n?1∑zkj??zkj??,i,j=1,2,?,p
3、解樣本相關(guān)矩陣R的特征方程
∣
R
?
λ
I
p
∣
=
0
|R-\lambda{I_{p}}|=0
∣R?λIp?∣=0得出p個特征根，確定主成分滿足
∑
j
=
1
m
λ
j
∑
j
=
1
p
λ
j
≥
0.85
\frac{\sum_{j=1}^{m}\lambda_{j}}{\sum_{j=1}^{p}\lambda_{j}}\ge0.85
∑j=1p?λj?∑j=1m?λj??≥0.85確定m的值，使信息的利用率達(dá)85%以上，對每個
λ
j
,
j
=
1
,
2
,
?
?
,
m
\lambda_{j},j=1,2,\cdots,m
λj?,j=1,2,?,m,解方程組
R
b
=
λ
j
b
Rb=\lambda_{j}b
Rb=λj?b得單位特征向量
b
j
a
b_{j}^{a}
bja?。
4、將標(biāo)準(zhǔn)化后的指標(biāo)變量轉(zhuǎn)化為主成分
U
i
j
=
z
i
T
b
j
a
,
j
=
1
,
2
,
?
?
,
m
U_{ij}=z_{i}^{T}b_{j}^{a},j=1,2,\cdots,m
Uij?=ziT?bja?,j=1,2,?,m
U
1
U_{1}
U1?稱為第一主成分，
U
2
U_{2}
U2?稱為第二主成分，
?
\cdots
?,
U
p
U_{p}
Up?稱為第p主成分。
5、對m個主成分進(jìn)行綜合評價
對m個主成分進(jìn)行加權(quán)求和，即得最終評價值，權(quán)數(shù)為每個主成分的方差貢獻(xiàn)率。
以下為MATLAB代碼：文章來源：http://www.zghlxwxcb.cn/news/detail-771656.html
/主成分分析法
x=[1.2,3,-1.1,17;
1.5,5,-3,22;
1.3,4.0,-2,19;
0.7,3,-2.3,11
1,4,-1.2,20.8];
stdr=std(x); %求各變量的標(biāo)準(zhǔn)差；
[n,m]=size(x); %矩陣的行與列
sddata=x./stdr(ones(n,1),:); %標(biāo)準(zhǔn)化變換
[p,princ,egenvalue]=princomp(sddata); %調(diào)用主成分
p=p(:,1:3); %輸出前3主成分系
sc=princ(:,1:3); %前3主成分量；
egenvalue; %相關(guān)系數(shù)矩陣的特征值，即各主成分所占比例；
per=100*egenvalue/sum(egenvalue); %各個主成分所占百分比；
輸出結(jié)果為：文章來源地址http://www.zghlxwxcb.cn/news/detail-771656.html
| 成分 | Value |
|---|---|
| 1 | 68.2742 |
| 2 | 25.1274 |
| 3 | 6.5776 |
| 4 | 0.0208 |
到了這里，關(guān)于數(shù)學(xué)建模常用方法及MATLAB代碼的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
