九、預(yù)測模型

        預(yù)測模型主要涉及灰色預(yù)測模型、神經(jīng)網(wǎng)絡(luò)預(yù)測模型。

1.灰色預(yù)測模型

Matlab代碼：

%%  輸入原始數(shù)據(jù)并做出時間序列圖
clear;clc
year =[1995:1:2004]';  % 橫坐標(biāo)表示年份，寫成列向量的形式（加'就表示轉(zhuǎn)置）
x0 = [174,179,183,189,207,234,220.5,256,270,285]';  %原始數(shù)據(jù)序列，寫成列向量的形式（加'就表示轉(zhuǎn)置）
% year = [2009:2015]; % 其實本程序?qū)懗闪诵邢蛄恳部梢裕驗槲遗履銈冋娴倪@么寫了，所以在后面會有判斷。
% x0 = [730, 679, 632, 599, 589, 532, 511];
% year = [2010:2017]';   % 該數(shù)據(jù)很特殊，可以通過準(zhǔn)指數(shù)規(guī)律檢驗，但是預(yù)測效果卻很差
% x0 = [1.321,0.387,0.651,0.985,1.235,0.987,0.854,1.021]';
% year = [2014:2017]';
% x0 = [2.874,3.278,3.337,3.390]';
% 畫出原始數(shù)據(jù)的時間序列圖
figure(1); % 因為我們的圖形不止一個，因此要設(shè)置編號
plot(year,x0,'o-'); grid on;  % 原式數(shù)據(jù)的時間序列圖
set(gca,'xtick',year(1:1:end))  % 設(shè)置x軸橫坐標(biāo)的間隔為1
xlabel('年份');  ylabel('排污總量');  % 給坐標(biāo)軸加上標(biāo)簽
%% 因為我們要使用GM(1,1)模型，其適用于數(shù)據(jù)期數(shù)較短的非負(fù)時間序列
ERROR = 0;  % 建立一個錯誤指標(biāo)，一旦出錯就指定為1
% 判斷是否有負(fù)數(shù)元素
if sum(x0<0) > 0  % x0<0返回一個邏輯數(shù)組(0-1組成)，如果有數(shù)據(jù)小于0，則所在位置為1，如果原始數(shù)據(jù)均為非負(fù)數(shù)，那么這個邏輯數(shù)組中全為0，求和后也是0~
    disp('親，灰色預(yù)測的時間序列中不能有負(fù)數(shù)哦')
    ERROR = 1;
end
% 判斷數(shù)據(jù)量是否太少
n = length(x0);  % 計算原始數(shù)據(jù)的長度
disp(strcat('原始數(shù)據(jù)的長度為',num2str(n)))    % strcat()是連接字符串的函數(shù)，第一講學(xué)了，可別忘了哦
if n<=3
    disp('親，數(shù)據(jù)量太小，我無能為力哦')
    ERROR = 1;
end
% 數(shù)據(jù)太多時提示可考慮使用其他方法（不報錯）
if n>10
    disp('親，這么多數(shù)據(jù)量，一定要考慮使用其他的方法哦，例如ARIMA，指數(shù)平滑等')
end
% 判斷數(shù)據(jù)是否為列向量，如果輸入的是行向量則轉(zhuǎn)置為列向量
if size(x0,1) == 1
    x0 = x0';
end
if size(year,1) == 1
    year = year';
end
%% 對一次累加后的數(shù)據(jù)進(jìn)行準(zhǔn)指數(shù)規(guī)律的檢驗(注意，這個檢驗有時候即使能通過，也不一定能保證預(yù)測結(jié)果非常好，例如上面的第三組數(shù)據(jù))
if ERROR == 0   % 如果上述錯誤均沒有發(fā)生時，才能執(zhí)行下面的操作步驟
    disp('------------------------------------------------------------')
    disp('準(zhǔn)指數(shù)規(guī)律檢驗')
    x1 = cumsum(x0);   % 生成1-AGO序列，cumsum是累加函數(shù)哦~    注意：1.0e+03 *0.1740的意思是科學(xué)計數(shù)法,即10^3*0.1740 = 174
    rho = x0(2:end) ./ x1(1:end-1) ;   % 計算光滑度rho(k) = x0(k)/x1(k-1)
    % 畫出光滑度的圖形，并畫上0.5的直線，表示臨界值
    figure(2)
    plot(year(2:end),rho,'o-',[year(2),year(end)],[0.5,0.5],'-'); grid on;
    text(year(end-1)+0.2,0.55,'臨界線')   % 在坐標(biāo)(year(end-1)+0.2,0.55)上添加文本
    set(gca,'xtick',year(2:1:end))  % 設(shè)置x軸橫坐標(biāo)的間隔為1
    xlabel('年份');  ylabel('原始數(shù)據(jù)的光滑度');  % 給坐標(biāo)軸加上標(biāo)簽
    disp(strcat('指標(biāo)1：光滑比小于0.5的數(shù)據(jù)占比為',num2str(100*sum(rho<0.5)/(n-1)),'%'))
    disp(strcat('指標(biāo)2：除去前兩個時期外，光滑比小于0.5的數(shù)據(jù)占比為',num2str(100*sum(rho(3:end)<0.5)/(n-3)),'%'))
    disp('參考標(biāo)準(zhǔn)：指標(biāo)1一般要大于60%, 指標(biāo)2要大于90%，你認(rèn)為本例數(shù)據(jù)可以通過檢驗嗎？')
    Judge = input('你認(rèn)為可以通過準(zhǔn)指數(shù)規(guī)律的檢驗嗎？可以通過請輸入1，不能請輸入0：');
    if Judge == 0
        disp('親，灰色預(yù)測模型不適合你的數(shù)據(jù)哦~ 請考慮其他方法吧 例如ARIMA，指數(shù)平滑等')
        ERROR = 1;
    end
    disp('------------------------------------------------------------')
end
%% 當(dāng)數(shù)據(jù)量大于4時，我們利用試驗組來選擇使用傳統(tǒng)的GM(1,1)模型、新信息GM(1,1)模型還是新陳代謝GM(1,1)模型； 如果數(shù)據(jù)量等于4，那么我們直接對三種方法求一個平均來進(jìn)行預(yù)測
if ERROR == 0   % 如果上述錯誤均沒有發(fā)生時，才能執(zhí)行下面的操作步驟
    if  n > 4  % 數(shù)據(jù)量大于4時，將數(shù)據(jù)分為訓(xùn)練組和試驗組(根據(jù)原數(shù)據(jù)量大小n來取，n為5-7個則取最后兩年為試驗組，n大于7則取最后三年為試驗組)
        disp('因為原數(shù)據(jù)的期數(shù)大于4，所以我們可以將數(shù)據(jù)組分為訓(xùn)練組和試驗組')   % 注意，如果試驗組的個數(shù)只有1個，那么三種模型的結(jié)果完全相同，因此至少要取2個試驗組
        if n > 7
            test_num = 3;
        else
            test_num = 2;
        end
        train_x0 = x0(1:end-test_num);  % 訓(xùn)練數(shù)據(jù)
        disp('訓(xùn)練數(shù)據(jù)是: ')
        disp(mat2str(train_x0'))  % mat2str可以將矩陣或者向量轉(zhuǎn)換為字符串顯示, 這里加一撇表示轉(zhuǎn)置，把列向量變成行向量方便觀看
        test_x0 =  x0(end-test_num+1:end); % 試驗數(shù)據(jù)
        disp('試驗數(shù)據(jù)是: ')
        disp(mat2str(test_x0'))  % mat2str可以將矩陣或者向量轉(zhuǎn)換為字符串顯示
        disp('------------------------------------------------------------')
        % 使用三種模型對訓(xùn)練數(shù)據(jù)進(jìn)行訓(xùn)練，返回的result就是往后預(yù)測test_num期的數(shù)據(jù)
        disp(' ')
        disp('***下面是傳統(tǒng)的GM(1,1)模型預(yù)測的詳細(xì)過程***')
        result1 = gm11(train_x0, test_num); %使用傳統(tǒng)的GM(1,1)模型對訓(xùn)練數(shù)據(jù)，并預(yù)測后test_num期的結(jié)果
        disp(' ')
        disp('***下面是進(jìn)行新信息的GM(1,1)模型預(yù)測的詳細(xì)過程***')
        result2 = new_gm11(train_x0, test_num); %使用新信息GM(1,1)模型對訓(xùn)練數(shù)據(jù)，并預(yù)測后test_num期的結(jié)果
        disp(' ')
        disp('***下面是進(jìn)行新陳代謝的GM(1,1)模型預(yù)測的詳細(xì)過程***')
        result3 = metabolism_gm11(train_x0, test_num); %使用新陳代謝GM(1,1)模型對訓(xùn)練數(shù)據(jù)，并預(yù)測后test_num期的結(jié)果
        % 現(xiàn)在比較三種模型對于試驗數(shù)據(jù)的預(yù)測結(jié)果
        disp(' ')
        disp('------------------------------------------------------------')
        % 繪制對試驗數(shù)據(jù)進(jìn)行預(yù)測的圖形（對于部分?jǐn)?shù)據(jù)，可能三條直線預(yù)測的結(jié)果非常接近）
        test_year = year(end-test_num+1:end);  % 試驗組對應(yīng)的年份
        figure(3)
        plot(test_year,test_x0,'o-',test_year,result1,'*-',test_year,result2,'+-',test_year,result3,'x-'); grid on;
        set(gca,'xtick',year(end-test_num+1): 1 :year(end))  % 設(shè)置x軸橫坐標(biāo)的間隔為1
        legend('試驗組的真實數(shù)據(jù)','傳統(tǒng)GM(1,1)預(yù)測結(jié)果','新信息GM(1,1)預(yù)測結(jié)果','新陳代謝GM(1,1)預(yù)測結(jié)果')  % 注意：如果lengend擋著了圖形中的直線，那么lengend的位置可以自己手動拖動
        xlabel('年份');  ylabel('排污總量');  % 給坐標(biāo)軸加上標(biāo)簽
        % 計算誤差平方和SSE
        SSE1 = sum((test_x0-result1).^2);
        SSE2 = sum((test_x0-result2).^2);
        SSE3 = sum((test_x0-result3).^2);
        disp(strcat('傳統(tǒng)GM(1,1)對于試驗組預(yù)測的誤差平方和為',num2str(SSE1)))
        disp(strcat('新信息GM(1,1)對于試驗組預(yù)測的誤差平方和為',num2str(SSE2)))
        disp(strcat('新陳代謝GM(1,1)對于試驗組預(yù)測的誤差平方和為',num2str(SSE3)))
        if SSE1<SSE2
            if SSE1<SSE3
                choose = 1;  % SSE1最小，選擇傳統(tǒng)GM(1,1)模型
            else
                choose = 3;  % SSE3最小，選擇新陳代謝GM(1,1)模型
            end
        elseif SSE2<SSE3
            choose = 2;  % SSE2最小，選擇新信息GM(1,1)模型
        else
            choose = 3;  % SSE3最小，選擇新陳代謝GM(1,1)模型
        end
        Model = {'傳統(tǒng)GM(1,1)模型','新信息GM(1,1)模型','新陳代謝GM(1,1)模型'};
        disp(strcat('因為',Model(choose),'的誤差平方和最小，所以我們應(yīng)該選擇其進(jìn)行預(yù)測'))
        disp('------------------------------------------------------------')
        %% 選用誤差最小的那個模型進(jìn)行預(yù)測
        predict_num = input('請輸入你要往后面預(yù)測的期數(shù)： ');
        % 計算使用傳統(tǒng)GM模型的結(jié)果，用來得到另外的返回變量：x0_hat, 相對殘差relative_residuals和級比偏差eta
        [result, x0_hat, relative_residuals, eta] = gm11(x0, predict_num);  % 先利用gm11函數(shù)得到對原數(shù)據(jù)擬合的詳細(xì)結(jié)果
        % % 判斷我們選擇的是哪個模型，如果是2或3，則更新剛剛由模型1計算出來的預(yù)測結(jié)果
        if choose == 2
            result = new_gm11(x0, predict_num);
        end
        if choose == 3
            result = metabolism_gm11(x0, predict_num);
        end
        %% 輸出使用最佳的模型預(yù)測出來的結(jié)果
        disp('------------------------------------------------------------')
        disp('對原始數(shù)據(jù)的擬合結(jié)果：')
        for i = 1:n
            disp(strcat(num2str(year(i)), ' ： ',num2str(x0_hat(i))))
        end
        disp(strcat('往后預(yù)測',num2str(predict_num),'期的得到的結(jié)果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result(i))))
        end
        %% 如果只有四期數(shù)據(jù)，那么我們就沒必要選擇何種模型進(jìn)行預(yù)測，直接對三種模型預(yù)測的結(jié)果求一個平均值~
    else
        disp('因為數(shù)據(jù)只有4期，因此我們直接將三種方法的結(jié)果求平均即可~')
        predict_num = input('請輸入你要往后面預(yù)測的期數(shù)： ');
        disp(' ')
        disp('***下面是傳統(tǒng)的GM(1,1)模型預(yù)測的詳細(xì)過程***')
        [result1, x0_hat, relative_residuals, eta] = gm11(x0, predict_num);
        disp(' ')
        disp('***下面是進(jìn)行新信息的GM(1,1)模型預(yù)測的詳細(xì)過程***')
        result2 = new_gm11(x0, predict_num);
        disp(' ')
        disp('***下面是進(jìn)行新陳代謝的GM(1,1)模型預(yù)測的詳細(xì)過程***')
        result3 = metabolism_gm11(x0, predict_num);
        result = (result1+result2+result3)/3;
        disp('對原始數(shù)據(jù)的擬合結(jié)果：')
        for i = 1:n
            disp(strcat(num2str(year(i)), ' ： ',num2str(x0_hat(i))))
        end
        disp(strcat('傳統(tǒng)GM(1,1)往后預(yù)測',num2str(predict_num),'期的得到的結(jié)果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result1(i))))
        end
        disp(strcat('新信息GM(1,1)往后預(yù)測',num2str(predict_num),'期的得到的結(jié)果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result2(i))))
        end
        disp(strcat('新陳代謝GM(1,1)往后預(yù)測',num2str(predict_num),'期的得到的結(jié)果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result3(i))))
        end
        disp(strcat('三種方法求平均得到的往后預(yù)測',num2str(predict_num),'期的得到的結(jié)果：'))
        for i = 1:predict_num
            disp(strcat(num2str(year(end)+i), ' ： ',num2str(result(i))))
        end
    end
    %% 繪制相對殘差和級比偏差的圖形(注意：因為是對原始數(shù)據(jù)的擬合效果評估，所以三個模型都是一樣的哦~~~)
    figure(4)
    subplot(2,1,1)  % 繪制子圖（將圖分塊）
    plot(year(2:end), relative_residuals,'*-'); grid on;   % 原數(shù)據(jù)中的各時期和相對殘差
    legend('相對殘差'); xlabel('年份');
    set(gca,'xtick',year(2:1:end))  % 設(shè)置x軸橫坐標(biāo)的間隔為1
    subplot(2,1,2)
    plot(year(2:end), eta,'o-'); grid on;   % 原數(shù)據(jù)中的各時期和級比偏差
    legend('級比偏差'); xlabel('年份');
    set(gca,'xtick',year(2:1:end))  % 設(shè)置x軸橫坐標(biāo)的間隔為1
    disp(' ')
    disp('****下面將輸出對原數(shù)據(jù)擬合的評價結(jié)果***')
    %% 殘差檢驗
    average_relative_residuals = mean(relative_residuals);  % 計算平均相對殘差 mean函數(shù)用來均值
    disp(strcat('平均相對殘差為',num2str(average_relative_residuals)))
    if average_relative_residuals<0.1
        disp('殘差檢驗的結(jié)果表明：該模型對原數(shù)據(jù)的擬合程度非常不錯')
    elseif average_relative_residuals<0.2
        disp('殘差檢驗的結(jié)果表明：該模型對原數(shù)據(jù)的擬合程度達(dá)到一般要求')
    else
        disp('殘差檢驗的結(jié)果表明：該模型對原數(shù)據(jù)的擬合程度不太好，建議使用其他模型預(yù)測')
    end
    %% 級比偏差檢驗
    average_eta = mean(eta);   % 計算平均級比偏差
    disp(strcat('平均級比偏差為',num2str(average_eta)))
    if average_eta<0.1
        disp('級比偏差檢驗的結(jié)果表明：該模型對原數(shù)據(jù)的擬合程度非常不錯')
    elseif average_eta<0.2
        disp('級比偏差檢驗的結(jié)果表明：該模型對原數(shù)據(jù)的擬合程度達(dá)到一般要求')
    else
        disp('級比偏差檢驗的結(jié)果表明：該模型對原數(shù)據(jù)的擬合程度不太好，建議使用其他模型預(yù)測')
    end
    disp(' ')
    disp('------------------------------------------------------------')
    %% 繪制最終的預(yù)測效果圖
    figure(5)  % 下面繪圖中的符號m:洋紅色 b:藍(lán)色
    plot(year,x0,'-o',  year,x0_hat,'-*m',  year(end)+1:year(end)+predict_num,result,'-*b' );   grid on;
    hold on;
    plot([year(end),year(end)+1],[x0(end),result(1)],'-*b')
    legend('原始數(shù)據(jù)','擬合數(shù)據(jù)','預(yù)測數(shù)據(jù)')  % 注意：如果lengend擋著了圖形中的直線，那么lengend的位置可以自己手動拖動
    set(gca,'xtick',[year(1):1:year(end)+predict_num])  % 設(shè)置x軸橫坐標(biāo)的間隔為1
    xlabel('年份');  ylabel('排污總量');  % 給坐標(biāo)軸加上標(biāo)簽
end

---

2.BP神經(jīng)網(wǎng)絡(luò)

---

十、降維模型

1.奇異值分解（Singular Value Decomposition,SVD）

        奇異值分解（Singular Value Decomposition）是線性代數(shù)中一種重要的矩陣分解，其在圖形學(xué)、統(tǒng)計學(xué)、推薦系統(tǒng)、信號處理等領(lǐng)域有重要應(yīng)用。此處介紹奇異值分解在圖形壓縮中的運用，并將介紹Matlab對于圖形和視頻的處理。

 

Matlab代碼：

A = [4 0 1 6;0 0 5 1;2 1 3 2]  % A : 3*4
% 注意：U*S*(V的轉(zhuǎn)置) == A    U：3*3    S：3*4    V：4*4
[U,S,V] = svd(A)  
U(:,1:2)*S(1:2,1:2)*V(:,1:2)'  
U(:,1:3)*S(1:3,1:3)*V(:,1:3)' % 就是A

 <mysvd.m>

function [compress_A] = mysvd(A, ratio)
    % 函數(shù)作用：使用奇異值分解將矩陣A壓縮到指定的特征比例
    % 輸入變量
    %     A：要壓縮的m*n維的矩陣
    %     ratio：要保留原矩陣的特征比例（100%表示不壓縮）
    % 輸出變量
    %     compress_A: 壓縮后的矩陣
    [U,S,V] = svd(A);   % U:m*m         S：m*n      V ： n*n
    eigs = diag(S);  % diag函數(shù)可以返回S的主對角線元素，即矩陣A的奇異值，并將其保存到列向量中
    SUM = sum(eigs);  % 計算所有奇異值的總和
    temp = 0;   % 新建臨時變量，用于下面的循環(huán)
    for  i = 1: length(eigs)  % 循環(huán)
        temp =temp + eigs(i);  % 每循環(huán)一次，就更新temp的值為原來的temp值+接下來的一個奇異值
        if (temp/SUM) > ratio  % 如果現(xiàn)在的比例超過了ratio,就退出循環(huán)
            break
        end
    end
    disp(['壓縮后保留原矩陣的比例特征為：',num2str(roundn(100*temp/SUM,-2)),'%'])  %roundn(x,-2)可將x四舍五入到小數(shù)點后兩位
    compress_A= U(:,1:i)*S(1:i,1:i)*V(:,1:i)';
end

 <photo_compress.m>

function []= photo_compress(photo_address, save_address, ratio, greycompress)
    % 函數(shù)作用：利用SVD對圖形進(jìn)行壓縮
    % 輸入變量
    %     photo_address：要壓縮的圖片存放的位置（建議輸入完整的路徑）
    %     save_address：將壓縮后的圖片保存的位置（建議輸入完整的路徑）
    %     ratio：要保留原矩陣的特征比例（100%表示不壓縮）
    %     greycompress: 如果該值等于1，則會彩色的原圖片轉(zhuǎn)換為灰色圖片后再壓縮；默認(rèn)值為0，表示不進(jìn)行轉(zhuǎn)換
    % 輸出變量
    %     無（不需要輸出，因為函數(shù)運行過程中已經(jīng)將圖片保存了~）
    if nargin == 3  % 判斷用戶輸入的參數(shù)，如果只輸入了前三個參數(shù)，則默認(rèn)最后的參數(shù)greycompress=0
        greycompress = 0;
    end
    img = double(imread(photo_address));
    % 圖片保存的對象是 'uint8' 類型，需要將其轉(zhuǎn)換為double類型才能進(jìn)行奇異值分解的操作
    % 注意:  img是圖形的像素矩陣，如果是彩色圖片則是三維矩陣，如果是灰色圖片(R=G=B)則是二維矩陣
    % '赫本.jpg'是灰色的圖片，得到的img類型是[914×1200]double
    % '千與千尋.jpg'是彩色的圖片，得到的img類型是[768×1024×3]double
    % 因此我們可利用第三個維度的大小來判斷圖片是否為灰色的
    % 灰色圖片的只有兩個維度，所以size(img ,3) == 1
    if (greycompress == 1) && (size(img ,3) == 3)  % 如果圖片為彩色，且greycompress的值等于1，則會彩色的原圖片轉(zhuǎn)換為灰色圖片后再壓縮
        img = double(rgb2gray(imread(photo_address)));    % rgb2gray函數(shù)可以將彩色圖片轉(zhuǎn)換為灰色圖片, 注意：輸入的變量要為默認(rèn)的'uint8' 類型的圖片對象
    end  % 注意： grey(英)和gray(美)都表示灰色
    if size(img ,3) == 3   % 判斷圖片是否為彩色的
        R=img(:,:,1);       % RGB色彩模式三要素：紅色
        G=img(:,:,2);       % RGB色彩模式三要素：綠色
        B=img(:,:,3);       % RGB色彩模式三要素：藍(lán)色
        disp(['正在壓縮:  ',photo_address,'的紅色要素'])
        r = mysvd(R, ratio);  % 調(diào)用自定義函數(shù)將R矩陣壓縮成r
        disp(['正在壓縮:  ',photo_address,'的綠色要素'])
        g = mysvd(G, ratio); % 調(diào)用自定義函數(shù)將G矩陣壓縮成g
        disp(['正在壓縮:  ',photo_address,'的藍(lán)色要素'])
        b = mysvd(B, ratio); % 調(diào)用自定義函數(shù)將B矩陣壓縮成b
        compress_img=cat(3,r,g,b);  % 根據(jù)三個RGB矩陣（壓縮后的r、g、b）生成圖片對象
    else  % 如果圖片是灰色的要執(zhí)行的步驟
        disp(['正在壓縮灰色圖片:  ',photo_address])
        compress_img = mysvd(img, ratio);  %如果是灰色圖片的話，直接壓縮img矩陣就好了
    end
    % 將壓縮后的圖片保存
    imwrite(uint8(compress_img), save_address);   % 如果你的矩陣是double格式的，導(dǎo)出時會自動將范圍認(rèn)為是[0 1]，需要重新轉(zhuǎn)換為uint8類型
end

 

---

2.主成分分析（Principal Component Analysis,PCA）

        主成分分析(Principal Component Analysis)是一種降維算法，它能將多個指標(biāo)轉(zhuǎn)換為少數(shù)幾個主成分，這些主成分是原始變量的線性組合，且彼此之間互不相關(guān)，其能反映出原始數(shù)據(jù)的大部分信息。一般來說， 當(dāng)研究的問題涉及到多變量且變量之間存在很強的相關(guān)性時，我們可考慮使用主成分分析的方法來對數(shù)據(jù)進(jìn)行簡化。

 <PCA.m>

load data1.mat   % 主成分聚類
%  load data2.mat   % 主成分回歸
% 注意，這里可以對數(shù)據(jù)先進(jìn)行描述性統(tǒng)計
% 描述性統(tǒng)計的內(nèi)容見第5講.相關(guān)系數(shù)
[n,p] = size(x);  % n是樣本個數(shù)，p是指標(biāo)個數(shù)
%% 第一步：對數(shù)據(jù)x標(biāo)準(zhǔn)化為X
X=zscore(x);   % matlab內(nèi)置的標(biāo)準(zhǔn)化函數(shù)（x-mean(x)）/std(x)
%% 第二步：計算樣本協(xié)方差矩陣
R = cov(X);
%% 注意：以上兩步可合并為下面一步：直接計算樣本相關(guān)系數(shù)矩陣
R = corrcoef(x);
disp('樣本相關(guān)系數(shù)矩陣為：')
disp(R)
%% 第三步：計算R的特征值和特征向量
% 注意：R是半正定矩陣，所以其特征值不為負(fù)數(shù)
% R同時是對稱矩陣，Matlab計算對稱矩陣時，會將特征值按照從小到大排列哦
% eig函數(shù)的詳解見第一講層次分析法的視頻
[V,D] = eig(R);  % V 特征向量矩陣  D 特征值構(gòu)成的對角矩陣
%% 第四步：計算主成分貢獻(xiàn)率和累計貢獻(xiàn)率
lambda = diag(D);  % diag函數(shù)用于得到一個矩陣的主對角線元素值(返回的是列向量)
lambda = lambda(end:-1:1);  % 因為lambda向量是從小大到排序的，我們將其調(diào)個頭
contribution_rate = lambda / sum(lambda);  % 計算貢獻(xiàn)率
cum_contribution_rate = cumsum(lambda)/ sum(lambda);   % 計算累計貢獻(xiàn)率  cumsum是求累加值的函數(shù)
disp('特征值為：')
disp(lambda')  % 轉(zhuǎn)置為行向量，方便展示
disp('貢獻(xiàn)率為：')
disp(contribution_rate')
disp('累計貢獻(xiàn)率為：')
disp(cum_contribution_rate')
disp('與特征值對應(yīng)的特征向量矩陣為：')
% 注意：這里的特征向量要和特征值一一對應(yīng)，之前特征值相當(dāng)于顛倒過來了，因此特征向量的各列需要顛倒過來
%  rot90函數(shù)可以使一個矩陣逆時針旋轉(zhuǎn)90度，然后再轉(zhuǎn)置，就可以實現(xiàn)將矩陣的列顛倒的效果
V=rot90(V)';
disp(V)
%% 計算我們所需要的主成分的值
m =input('請輸入需要保存的主成分的個數(shù):  ');
F = zeros(n,m);  %初始化保存主成分的矩陣（每一列是一個主成分）
for i = 1:m
    ai = V(:,i)';   % 將第i個特征向量取出，并轉(zhuǎn)置為行向量
    Ai = repmat(ai,n,1);   % 將這個行向量重復(fù)n次，構(gòu)成一個n*p的矩陣
    F(:, i) = sum(Ai .* X, 2);  % 注意，對標(biāo)準(zhǔn)化的數(shù)據(jù)求了權(quán)重后要計算每一行的和
end
%% (1)主成分聚類 ： 將主成分指標(biāo)所在的F矩陣復(fù)制到Excel表格，然后再用Spss進(jìn)行聚類
% 在Excel第一行輸入指標(biāo)名稱（F1,F2, ..., Fm）
% 雙擊Matlab工作區(qū)的F,進(jìn)入變量編輯中，然后復(fù)制里面的數(shù)據(jù)到Excel表格
% 導(dǎo)出數(shù)據(jù)之后，我們后續(xù)的分析就可以在Spss中進(jìn)行。
%%（2）主成分回歸：將x使用主成分得到主成分指標(biāo)，并將y標(biāo)準(zhǔn)化，接著導(dǎo)出到Excel，然后再使用Stata回歸
% Y = zscore(y);  % 一定要將y進(jìn)行標(biāo)準(zhǔn)化哦~
% 在Excel第一行輸入指標(biāo)名稱（Y,F1, F2, ..., Fm）
% 分別雙擊Matlab工作區(qū)的Y和F,進(jìn)入變量編輯中，然后復(fù)制里面的數(shù)據(jù)到Excel表格
% 導(dǎo)出數(shù)據(jù)之后，我們后續(xù)的分析就可以在Stata中進(jìn)行。

 

 

 

---

3.因子分析

        因子分析由斯皮爾曼在1904年首次提出，其在某種程度上可以被看成是主成分分析的推廣和擴(kuò)展。 因子分析法通過研究變量間的相關(guān)系數(shù)矩陣，把這些變量間錯綜復(fù)雜的關(guān)系歸結(jié)成少數(shù)幾個綜合因子，由于歸結(jié)出的因子個數(shù)少于原始變量的個數(shù)，但是它們又包含原始變量的信息，所以，這一分析過程也稱為降維。由于因子往往比主成分更易得到解釋，故因子分析比主成分分析更容易成功，從而有更廣泛的應(yīng)用。

 

---

4.典型相關(guān)分析（Canonical Correlation analysis）

         典型相關(guān)分析是研究兩組變量（每組變量中都可能有多個指標(biāo)）之間相關(guān)關(guān)系的一種多元統(tǒng)計方法，它能夠揭示出兩組變量之間的內(nèi)在聯(lián)系。

 

 

---

內(nèi)容原作者：數(shù)學(xué)建模清風(fēng)

學(xué)習(xí)用途，僅作參考。