PCA设计

（初学PCA，不足之处还望指出，共同进步！）
一 实现平台：
本次设计的三个部分全部基于MATLAB2017a来实现。

二 要求：
共有三个部分：
第一部分是 给定的数据：PointNormals绿萝.txt 对绿萝植株的数据进行降维，要求将数据降低到两维和三维，由第1,2主元形成的二维点分布图；由第1,2,4主元形成的三维点分布图。生成的两幅分布图都应该是单色的。要求对这两幅图进行多角度展示。
第二部分是 自己构建一组存在噪声的三维点构成的片状集合，使用IPCA对点集进行平面拟合。并可视化的方法展示出平面拟合的过程与结果。
第三部分是 如果不使用以上的IPCA算法，设计出其他的办法使用迭代的思路从片状三维点集中拟合出平面方程，
三 PCA算法分析：
PCA算法主要用于降维，就是将样本数据从高维空间投影到低维空间中，用一种较少数量的主特征对样本进行描述以达到降低特征空间维数，并尽可能的在低维空间中表示原始数据的方法。
PCA算法计算步骤：
1: 给定n组数据{x1,x2,x3,…xn}，其中每组数据具有d维
2: 计算数据向量均值，进行中心化操作：（MATLAB中均值用mean()函数）
3: 构造数据矩阵X，随后计算协方差矩阵CX：（Matlab可使用cov()函数求样本的协方差矩阵)
4: 对CX进行特征值分解，CX=PTDP；其中D为对角阵，对角线上的元素按递减顺序排列，对应的特征向量在P矩阵中按行排列（Matlab可使用eigs() 函数求矩阵的特征值D和特征向量V,算出的结果是按照递增的顺序排列）
5: 计算Y=PX。那么: Y=[y1,y2,y3,…yn] ，其中每个向量是PCA得到的新数据。(如果第四步是CX=PDPT ，则P矩阵特征向量按列排列，第五步中变为Y=PTX)
6: 若要降维，只需要保留P中的前几行再进行第5步的计算。

Part1：
设计思路： 首先对已知的绿萝图数据进行处理，然后求解数据的均值得到均值m，这样可以对数据进行去中心化操作，求得去中心化后的新数据a参数的值。
数据处理后开始进行降维处理得到特征向量，首先利用cov()指令求解新数据的协方差，再利用eig()指令到它的特征向量。因为是从小到大的排列，与老师要求的正好相反，所以选择第6，5列为要求的1，2列来处理得到绿萝的二维图形，第6，5，3列为要求的1，2，4列来处理得到绿萝的三维图形。
图一为生成的绿萝的二维图像，图二与图三是对二维图进行旋转从不同的角度进行观察。
图四为生成的绿萝的三维图像，图五与图六是对绿萝的三维图的不同角度的观察。

具体程序如下：
%加载绿萝的数据
d=load(‘C:\Users\Administrator\Desktop\PointsNormals绿萝.txt’);
%pca降维处理
m=mean(d);%求数据的均值
m=repmat(m,[size(d,1),1]);
a=d-m;%绿萝的数据减去均值来对数据进行去中心化
c=cov(a);%求新数据的协方差
[V,D]=eig©;%V是特征向量，D是特征值构成的对角矩阵
V_6=V(:,6);%因为从小到大排列，所以取第6列，也就是最大列
V_5=V(:,5);%取第5列
V_3=V(:,3);%取第3列
V_56=[V_6,V_5]; %取12特征向量
V_356=[V_6,V_5,V_3];%取124特征向量
A=aV_56; %降维
B=aV_356;
X1=A(:,1); Y1=A(:,2);
X2=B(:,1); Y2=B(:,2);
Z2=B(:,3);
figure(1);%绘制二维图像
plot(X1,Y1,’.’,‘MarkerSize’,3);
title(‘绿萝的二维图像’);
figure(2);%绘制三维图像
plot3(X2,Y2,Z2,’.’,‘MarkerSize’,3);
title(‘绿萝的三维图像’);





Part2：
设计思路：运用PCA降维拟合平面
首先利用rand()指令随机生成20组空间的点（x,y,z）坐标。取其中的第一组的数据（x(1),y(1),z(1)）点作为迭代点用红圈出来。
然后对这20组的数据进行处理，求解数据的均值得到均值m，这样可以对数据进行去中心化操作，求得去中心化后的新数据a参数的值。数据处理后开始进行降维处理得到特征向量，首先利用cov()指令求解新数据的协方差，再利用eig()指令到它的特征向量。因为是从小到大的排列，所以选择第1列作为法向量(a,b,c)，根据公式d=-(ax(1)+by(1)+c*(1))来求出d参数，这样就可以运用点法式来求得平面。
运用点到平面的距离公式L=abs(ax(i)+by(i)+cz(i)+d)/sqrt(aa+bb+cc)来进行计算每个点到平面的距离，设置阈值为0.5。除了迭代点以外的其余19个点进行迭代计算，每次删除距离大于阈值的点，然后返回重新对剩余的三维点进行平面的拟合。直到所有的点都符合要求。
如图七是生成三维图，图八与图九是不同角度的观察图，图八可以看处三维点基本在最终的拟合平面附近，图九可以看出所有的平面都是以迭代点进行拟合的。
具体程序如下：

%随机生成一组空间三维的点
x=1+rand(20,1)2;
y=1+rand(20,1)2;
z=1+rand(20,1)2;
scatter3(x,y,z,‘filled’);hold on;
scatter3(x(1),y(1),z(1),‘rd’);hold on;
d=[x,y,z];
%pca降维处理
m=mean(d);%求数据的均值
a=d-m;%数据减去均值来对数据进行去中心化
c=cov(a);%求新数据的协方差
[V,D]=eig©;%V是特征向量，D是特征值构成的对角矩阵
%因为从小到大排列，所以取第1列，也就是最小列来作为法向量
a=V(1,1);%分别取出法向量的三个参数a,b,c的值
b=V(2,1);
c=V(3,1);
d=-(ax(1)+by(1)+c(1));%根据公式求出平面方程的d参数
%绘制三维拟合平面的图像
xfit = min(x):0.1:max(x);%取x轴的坐标参数
yfit = min(y):0.1:max(y);%取x轴的坐标参数
[XFIT,YFIT]= meshgrid (xfit,yfit);%将向量x和y定义的区域转换成矩阵X和Y
ZFIT = -( a * (XFIT-x(1)) + b ( YFIT-y(1)))/c+z(1);%根据已知的点和平面方程公式来计算对应的z坐标
mesh(XFIT,YFIT,ZFIT);%绘制三维拟合平面的图像
title(‘拟合平面图’);hold on;
for i=2:20;%循环检验点到平面的距离，阈值的取值为0.5
if abs(ax(i)+by(i)+cz(i)+d)/sqrt(aa+bb+cc)>0.5%点到平面的距离公式
%去除不符合要求的点
x(i)=[]
y(i)=[]
z(i)=[]
%重新降维，取法向量，拟合平面的步骤，直到全部符合要求为止
d=[x,y,z];
m=mean(d,1);
a=d-m;
c=cov(a);
[V,D]=eig©;
a=V(1,1);
b=V(2,1);
c=V(3,1);
d=-(ax(1)+by(1)+c(1));
xfit = min(x):0.1:max(x);
yfit = min(y):0.1:max(y);
[XFIT,YFIT]= meshgrid (xfit,yfit);
ZFIT = -( a * (XFIT-x(1)) + b *( YFIT-y(1)))/c+z(1);
mesh(XFIT,YFIT,ZFIT);
title(‘拟合平面图’);hold on;
end
end

Part3：

设计思路：运用最小二乘法的线性回归来拟合平面
首先利用rand()指令随机生成20组空间的点（x,y,z）坐标。
然后对这20组的数据进行处理，使用X=[ones(size(x,1),1) x y]方法，用x参数与y参数来构造一个新的X矩阵。使用命令regress()实现多元线性回归拟合平面求出b参数。运用平面方程z =b(1)+b(2)x+b(3)y便可以拟合平面。
迭代部分依旧采用点到平面的距离方公式L=abs(ax(i)+by(i)+cz(i)+d)/sqrt(aa+bb+cc)判断每次删除距离大于阈值的点，然后返回重新对剩余的三维点进行平面的拟合。直到所有的点都符合要求。
如图十是生成新的三维图，图十一是不同角度的观察图，可以看处三维点基本在最终的拟合平面附近，但是拟合的效果不如Part2的好，可能是没有固定的迭代点导致。

具体程序如下：

%随机生成一组空间三维的点
x=1+rand(20,1)*2;
y=1+rand(20,1)*2;
z=1+rand(20,1)*2;
scatter3(x,y,z,‘filled’);hold on;
X=[ones(size(x,1),1) x y];%构造一个新的矩阵X
b=regress(z,X);%使用命令regress()实现多元线性回归拟合平面
xfit = min(x):0.1:max(x);%取x轴的坐标参数
yfit = min(y):0.1:max(y);%取y轴的坐标参数
[XFIT,YFIT]= meshgrid (xfit,yfit);%将向量x和y定义的区域转换成矩阵X和Y
ZFIT =b(1)+b(2)*XFIT+b(3)*YFIT;%根据已知b的参数值来构造平面方程进而求解z的值
mesh(XFIT,YFIT,ZFIT);title(‘第二种拟合平面图’);hold on;
for i=1:20;
if abs(b(2)*x(i)+b(3)*y(i)-z(i)+b(1))/sqrt(b(2)*b(2)+b(3)*b(3)+1)>0.5%点到平面的距离公式
%去除不符合要求的点
x(i)=[]
y(i)=[]
z(i)=[]
%重复拟合平面的步骤
X=[ones(size(x,1),1) x y];
b=regress(z,X);
xfit = min(x):0.1:max(x);
yfit = min(y):0.1:max(y);
[XFIT,YFIT]= meshgrid (xfit,yfit);
ZFIT =b(1)+b(2)*XFIT+b(3)*YFIT;
mesh(XFIT,YFIT,ZFIT);title(‘第二种拟合平面图’);hold on;
end
end