%程序說明：PCA=?pca（X），程序中X為?n*m階混合數據矩陣，m為信號個數即變量的個數，n為采樣點數即樣本個數
%?T為主分量矩陣P為得分向量矩陣，Q和T2為統計量。

%function?PCA=?pca（X）
tic;
load（‘高維特征.mat‘）
X=features（:1:30）;
[m?n]=size（X）;
if?X?==?0
error（‘You?must?supply?the?mixed?data?as?input?argument.‘）;%判斷是否輸入數據
end
if?length（size（X））>2
error（‘Input?data?can?not?have?more?than?two?dimensions.?‘）;%判斷輸入數據維數
end
if?any（any（isnan（X）））
error（‘Input?data?contains?NaN‘‘s.‘）;%判斷輸入數據是否含有非數值的數據
end

%——————————————標準化數據（首先去均值，然后標定到單位方差）————————————
a=0.99;
meanValue?=?mean（X）;?%求采集數據的樣本均值
e=ones（size（X1）1）;
s=std（X）;
d=diag（s.^-10）;
fprintf（‘Normalizing?the?input?data?X：\n‘）;
normX?=?（X?-?e*meanValue）*d;?%將采集數據去均值并標準化
[NumofSamplDim]?=?size（X）;?%最初采集數據的維數及每個變量樣本的個數
oldDim?=?Dim;??%最初變量個數，即數據的維數
fprintf（‘Number?of?process?variables:?%d\n‘Dim）;
fprintf（‘Number?of?samples?of?every?variable:?%d\n‘NumofSampl）;
fprintf（‘Calculate?PCA：\n‘）;

%———計算協方差矩陣的特征值大于閾值的個數numEig（實際求出特征值中非零的特征值個數）———
firstEig=1;
numEig?=?0;
covarianceMatrix?=?（1/（NumofSampl-1））*（normX‘*normX）;?%計算normX的協方差矩陣


[VD]?=?eig（covarianceMatrix）?;?%計算協方差矩陣的特征值和特征向量全部特征值構成對角陣D，特征向量為V的列向量
threshold?=?1e-5;
for?n=1:Dim
???if?D（nn）>threshold
???????numEig=numEig+1;
???else?if?D（nn）<=threshold
???????numEig=numEig;???????%實際求出的特征值大于零的非零特征值的個數
???????end
????end
end
%——————————降序排列各特征值——————————
eigenvalues?=?flipud（sort（diag（D）））;?%按從大到小順序排列的所有的特征值
t=eigenvalues‘;
%fprintf（‘The?eigenvalues?are?as?follows:?%g?\n‘eigenvalues）;
%fprintf（‘Number?of?non-zero?eigenvalues:?%d\n‘numEig）;

%—————————去掉較小的特征值（去掉較小值的閾值是去掉最大和最小特征值的和求平均）——————————
if?numEig?lowerValue=（sum（eigenvalues）-eigenvalues（1）-eigenvalues（2））/numEig;
else
lowerValue=（sum（eigenvalues）-eigenvalues（1）-eigenvalues（numEig））/numEig;?%特征值個數等于最初數據維數的情況
end
columns?=?diag（D）?>=lowerValue?;%求得去掉較小特征值后的特征值向量（邏輯向量，大于設定值的特征值為1，小于設定值的為0）

%—————————合并選擇的特征值，確定所選擇則的特征值——————————
selectedColumns?=columns?;?
newDim=sum（selectedColumns）;
%—————————輸出處理的結果信息—————————
fprintf（‘Selected[?%d?]?dimensions.\n‘newDim）;
fprintf（‘Smallest?remaining?（non-zero）?eigenvalue:?%g?\n‘eigenvalues（numEig））;
fprintf（‘Largest?remaining?（non-zero）?eigenvalue?:?%g?\n‘eigenvalues（firstEig））;

%———————選擇相應的特征值和特征向量———————
?B=V*fliplr（diag（selectedColumns））;??%根據特征值大小，選擇相應的特征向量
?y=find（sum（B）~=0）;????????%計算經過換算過的矩陣的列向量是否為0
?selectedEvector=B（:y）;??%根據所選擇的特征值，確定所選擇的單位特征向量組成的新矩陣

%——————————計算負荷向量構成的矩陣———————————
disp（‘Calculate?the?loading?vectors:‘）;
P=selectedEvector?;?%計算負荷向量P

%——————————提取主分量T————————————
disp（‘Calculate?the?score?vectors:‘）;
T=?normX*P??;???%計算得分向量，即主分量矩陣

%——————————計算Q統計量，即主元模型的SPE———————————
disp（‘Calculate?the?model?prediction?values:‘）;
normXpre=T*P‘;??%計算歸一化后的主元模型預測值矩陣
error=normX-normXpre;??%計算歸一化的相對于主元模型的預測的誤差
squaredError=error.^2??;%計算歸一化后的平方預測誤差
disp（‘Calculate?the?Q?statistic:‘）;
Q=sum（squaredError2）?;?%計算Q統計量

%——————————計算T統計量———————————
selectedEigvalue=diag（eigenvalues（1:newDim））;?%計算所選擇的特征值矩陣
disp（‘Calculate?the?T2?statistic:‘）;
T2=diag（T*inv（selec

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????5732??2018-10-15?10:02??pca降維算法\PCA.m

?????文件?????117313??2018-06-17?15:39??pca降維算法\高維特征.mat

?????目錄??????????0??2018-11-19?10:34??pca降維算法

-----------?---------??----------?-----??----

???????????????123045????????????????????3