C++与matlab混合编程基于主成份分析算法的数值分析
转自http://blog.csdn.net/codesofmonste/article/details/2467521
程序设计思路：
利用c++和matlab混合编程，完成基于pca（主成份分析算法）的数值分析程序。
关键字：c++、matlab、pca
function GRpca424(input)
%input为待分析的原始矩阵
%eigenValue是主成份分析过后所得到的原始矩阵的特征值（按照由大到小的顺序给出）
%eigenVector是主成分分析过后得到的对应特征值的特征向量矩阵（每列为对应）
x=importdata('C:/matlab/infile.txt');
input=x;
[n,p]=size(input);
%求每一列的平均值，并将列平均值保存在向量x2中
x=input;
for j=1:p
x1(j,1)=input(1,j);
for i=2:n
x1(j,1)=x1(j,1)+input(i,j);
end
x2(j,1)=x1(j,1)/n;
end
%对输入矩阵进行标准化计算过程
for j=1:p
s1(j,1)=(input(1,j)-x2(j,1))*(input(1,j)-x2(j,1));
for i=2:n
s1(j,1)=s1(j,1)+(input(i,j)-x2(j,1))*(input(i,j)-x2(j,1));
end
s(j,1)=sqrt(s1(j,1)/(n-1));
end
for i=1:n
for j=1:p
standard(i,j)=(input(i,j)-x2(j,1))/s(j,1);
end
end
%求相关系数矩阵cormatrix
for j=1:p
Ave1(j,1)=standard(1,j);
for i=2:n
Ave(j,1)=Ave1(j,1)+standard(i,j);
end
Ave(j,1)=Ave(j,1)/n;%standard列平均值
end
for i=1:n
ave1(i,1)=standard(i,1);
for j=2:p
ave(i,1)=ave1(i,1)+standard(i,j);
end
ave(i,1)=ave(i,1)/p;%standard行平均值
end
if n>p
for i=1:p
for j=1:p
Rsum(i,j)=(standard(1,i)-ave(i,1))*(standard(1,j)-Ave(j,1));
for l=2:n
Rsum(i,j)=Rsum(i,j)+(standard(l,i)-ave(i,1))*(standard(l,j)-
Ave(j,1));
end
end
end
for i=1:p
for j=1:p
Cov(i,j)=Rsum(i,j)/(n-1);
end
end
for i=1:p
for j=1:p
corMatrix(i,j)=Cov(i,j)/(sqrt(Cov(i,i))*sqrt(Cov(j,j)));
end
end
end
corMatrix
%计算相关系数矩阵的特征值、特征向量
disp('EigenVector为特征向量构成的标准化正交矩阵；Eigenvalue为特征值');
[EigenVector,Eigenvalue]=eig(corMatrix);
%对EigenVector和Eigenvalue进行排序，是特征值按照从大到小的顺序进行排列
if n>p
g=p;
end
if n<=p
g=n;
end
for i=1:(g-1)
if Eigenvalue(i,i)<Eigenvalue(i+1,i+1)
t=Eigenvalue(i,i);
Eigenvalue(i,i)=Eigenvalue(i+1,i+1);
Eigenvalue(i+1,i+1)=t;
for j=1:g
s=EigenVector(j,i);
EigenVector(j,i)=EigenVector(j,i+1);
EigenVector(j,i+1)=s;
end
end
end
Eigenvalue
EigenVector
for i=1:g
eigenValue(i,1)=Eigenvalue(i,i);
end
eigenValue
%计算分量方差贡献率
[m]=size(eigenValue);
%如果主成份的特征值矩阵有负数的话，就将其取相反数吧
for i=1:m
if eigenValue(i,1)<0
eigenValue(i,1)=-eigenValue(i,1)
end
end
%各分量的方差贡献率保存在main中
subM=0;
for i=1:m
subM=subM+eigenValue(i,1);
end
for i=1:m
main(i,1)=eigenValue(i,1)/subM;
end
disp('各分量的方差贡献率是：');
main
%将主成份方差贡献率输出到文本
dlmwrite('c:/matlab/方差贡献率.txt',main,'precision','%
10.5f','newline','pc');
%对主成份所占百分比进行比例分析，并根据用户设置的权重数进行最终的主成份判定
y=importdata('C:/matlab/var.txt');
n=y;%用户通过C++交互界面定义的达标的方差贡献率
k=main(1,1);
for i=1:m-1
if k>=n
disp('主成份向量个数为：');
m=i;
m
break;
else
k=k+main(i+1,1);
end
end
k=0;
for s=1:i
k=k+main(s,1);
end
%对特征向量矩阵进行标准化，将其转化为Vector
[p,p]=size(EigenVector);
x=EigenVector;
for j=1:p
x1(j,1)=EigenVector(1,j);
for i=2:p
x1(j,1)=x1(j,1)+EigenVector(i,j);
end
x2(j,1)=x1(j,1)/p;
end
for j=1:p
v1(j,1)=(EigenVector(1,j)-x2(j,1))*(EigenVector(1,j)-x2(j,1));
for i=2:p
v1(j,1)=s1(j,1)+(EigenVector(i,j)-x2(j,1))*(EigenVector(i,j)-x2
(j,1));
end
v(j,1)=sqrt(v1(j,1)/(p-1));
end
for i=1:p
for j=1:p
Vector(i,j)=(EigenVector(i,j)-x2(j,1))/v(j,1);
end
end
%样本的主成份值矩阵为
if m<=p
for i=1:m
for j=1:p
vector(j,i)=Vector(j,i);
end
end
end
F=standard*vector;
%F为输出矩阵，是经过变化之后的相关主成份分析值，具体用途可根据工程实际需要而
定
%将F输出到指定的文件中
dlmwrite('c:/matlab/主成份矩阵.txt',F,'precision','%10.5f','-
append','newline','pc');
前面的代码主要是一个可执行的M文件。我们可以通过C++与matlab的混合编程，即用C++调用matlab计算引擎，来完成，从可视化界面到高速计算实现主成分分析算法计算的全部过程。
用C++调用matlab的函数很简单：
程序执行必须的头文件是：
#include "engine.h"
调用matlab计算引擎的代码是：
Engine *ep=engOpen(NULL);--------（1）
engEvalString(ep,"GRpca424");-----（2）
engClose(ep);-----------------------------（3）
（1）定义个指针变量，可以看作是指向matlab的字符，通过它来定位matalb引擎位置；
（2）函数功能是调用matlab引擎“GRpca424”是可执行M文件的名字；
（3）关闭matlab引擎。
非常 乔木和小乔 网友 的学习总结，等十一之后准备自己也实现下！
