数据说明:采用的数据源是从别人那里拷的2012年全年的Sea Surface Temperature(海标温度,SST)数据,一直想找一份比较好的主成分分析数据,也没找到。
Matlab自身有主成分分析的函数princomp,其中返回的第二个数据就是样本经过K-L变换后的各个成分数据,第三个参数就是特征值大小。(第一个参数貌似是协方差矩阵,我还没看)
y=Tx,
式中X为待变换图像数据矩阵,Y为变换后的数据矩阵,T为实现这一线性变换的变换矩阵。如果变换矩阵T是正交矩阵,并且它是由原始图像数据矩阵X的协方差矩阵S的特征向量所组成,则该线性变换成为主成分分析,并且成Y矩阵的每一行矢量为变换后的一个主成分。
具体的算法步骤为:
马上要走了,偷个懒,直接把握之前做的一个PPT放出来了(话说也真是懒啊,这n多年没更新一点内容能不懒嘛!!!(⊙o⊙)…)
以下为代码实现部分:
close all; clear all; clc;%% PCA_Testpath = 'data\sst\2012\'; filelist = dir('data\sst\2012\*.sst.txt'); len = length(filelist); % savepath = 'data\sst_pic\';figure; % 读取文件 for i = 1:lenfilename = [path filelist(i).name];date = filename(15:20);sst_data(i,:) = dlmread(filename);%读取min_val = min(sst_data(i,:));max_val = max(sst_data(i,:));pic = reshape(sst_data(i,:),360,180);pic = rot90(pic);subplot(4,3,i),imshow(pic,[min_val max_val]),title(date);colormap; % % 保存数据图像 % pic = reshape(sst_data(i,:),360,180); % pic = rot90(pic); % strI = int2str(i); % save_path = [savepath strI '.tif']; % imwrite(pic,save_path);end% 取得样本数据 X = sst_data'; tempX = X;% 取样本大小:维度m,样本数n [m n] = size(X);% 求各样本平均值 meanVal = mean(X);% 样本矩阵中心化 X = X - repmat(meanVal,64800,1);% 计算协方差 S = X' * X ./ (m - 1); %".*"表示矩阵元素对应相乘% 计算特征值eg和特征向量Ev [Ev eg] = eig(S);% 由大到小排列对应 Ev = fliplr(Ev); eg = fliplr(eg);% % 特征向量转置 % EvT = Ev';% 得到新的成分 Y = tempX * Ev;% 将各个成分进行输出 figure; for i = 1:noutpic = Y(:,i);min_num = min(outpic);max_num = max(outpic);outpic = reshape(outpic,360,180);outpic = rot90(outpic); %outtitle = ['特征值=' num2str(latent(i))];subplot(4,3,i),imshow(outpic,[min_num max_num]),title(outtitle);%显示,并输出特征值colorbar; end
处理的结果对比显示:
主成分结果显示(按照特征值由大到小):
这个结果是和它内部的函数运算出来的结果是一样的:
% PCA [coeff score latent] = princomp(sst_data'); [m,n] = size(score); % figure; for i = 1:lenoutpic = score(:,i);min_num = min(outpic);max_num = max(outpic);outpic = reshape(outpic,360,180);outpic = rot90(outpic); outtitle = ['特征值=' num2str(latent(i))];subplot(4,3,i),imshow(outpic,[min_num max_num]),title(outtitle);%显示,并输出特征值colorbar; end
唉哟,这个博客园居然对Matlab着色显示的时候都不怎么想要支持,看来还是我程序界都不怎么用Matlab哦
正在写Python的实现方式,过两天有时间了再放上来看看吧~~~
有点仓促,就此搁笔。
2016年1月11日22:38
原文链接:http://www.cnblogs.com/leonwen/p/5122811.html