跪求非负独立成分分析的matlab程序？

您好， 这样的：

一、FastICA算法的基本步骤：

1. 对悄桥观测数据进行中心化，使它的均值为0；

2. 对数据进行白化，。

3. 选择需要估计的分量的个数，设迭代次数

4. 选择一个初始权矢量（随机的）。

5. 令,非线性函数的选取见前文。

6. 。

7. 令。

8. 假如不收敛的话，返回第5步。

9．令，禅运仿如贺纤果，返回第4步。

二．MATLAB源程序及说明：

%下程序为ICA的调用函数，输入为观察的信号，输出为解混后的信号

function Z=ICA(X)

%-----------去均值---------

[M,T] = size(X)%获取输入矩阵的行/列数，行数为观测数据的数目，列数为采样点数

average= mean(X')' %均值

for i=1:M

X(i,:)=X(i,:)-average(i)*ones(1,T)

end

%---------白化/球化------

Cx = cov(X',1) %计算协方差矩阵Cx

[eigvector,eigvalue] = eig(Cx)%计算Cx的特征值和特征向量

W=eigvalue^(-1/2)*eigvector' %白化矩阵

Z=W*X %正交矩阵

%----------迭代-------

Maxcount=10000 %最大迭代次数

Critical=0.00001 %判断是否收敛

m=M %需要估计的分量的个数

W=rand(m)

for n=1:m

WP=W(:,n) %初始权矢量（任意）

% Y=WP'*Z

% G=Y.^3%G为非线性函数，可取y^3等

% GG=3*Y.^2 %G的导数

count=0

LastWP=zeros(m,1)

W(:,n)=W(:,n)/norm(W(:,n))

while abs(WP-LastWP)&abs(WP+LastWP)>Critical

count=count+1 %迭代次数

LastWP=WP %上次迭代的值

% WP=1/T*Z*((LastWP'*Z).^3)'-3*LastWP

for i=1:m

WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i)

end

WPP=zeros(m,1)

for j=1:n-1

WPP=WPP+(WP'*W(:,j))*W(:,j)

end

WP=WP-WPP

WP=WP/(norm(WP))

if count==Maxcount

fprintf('未找到相应的信号)

return

end

end

W(:,n)=WP

end

Z=W'*Z

%以下为主程序，主要为原始信号的产生，观察信号和解混信号的作图

clear allclc

N=200n=1:N%N为采样点数

s1=2*sin(0.02*pi*n)%正弦信号

t=1:Ns2=2*square(100*t,50)%方波信号

a=linspace(1,-1,25)s3=2*[a,a,a,a,a,a,a,a]%锯齿信号

s4=rand(1,N)%随机噪声

S=[s1s2s3s4]%信号组成4*N

A=rand(4,4)

X=A*S%观察信号

%源信号波形图

figure(1)subplot(4,1,1)plot(s1)axis([0 N -5,5])title('源信号')

subplot(4,1,2)plot(s2)axis([0 N -5,5])

subplot(4,1,3)plot(s3)axis([0 N -5,5])

subplot(4,1,4)plot(s4)xlabel('Time/ms')

%观察信号（混合信号）波形图

figure(2)subplot(4,1,1)plot(X(1,:))title('观察信号（混合信号）')

subplot(4,1,2)plot(X(2,:))

subplot(4,1,3)plot(X(3,:))subplot(4,1,4)plot(X(4,:))

Z=ICA(X)

figure(3)subplot(4,1,1)plot(Z(1,:))title('解混后的信号')

subplot(4,1,2)plot(Z(2,:))

subplot(4,1,3)plot(Z(3,:))

subplot(4,1,4)plot(Z(4,:))xlabel('Time/ms')

在软桥猜件Matlab中实现主成分分析可以采取两种方式实现：一是通过编程来实现；二是直接调用Matlab中自带程序实现。

通过直接调用Matlab中的程序可以实现主成分分析：

式中：X为输入数据矩阵

（一般要求n>m）

输出变量：

①pc 主分量fi的系数，也叫因子系数；注意：pcTpc=单位阵

②score是主分量下的得分值；得分矩阵与数据矩阵X的阶数是一致的；

③variance是score对应列的方差向量，即A的特征值；容易计算方差所占的百分比

percent-v = 100*variance/sum(variance)

④t2表示检验的t2-统计量（方差分析要用）

计算过程中应用到计算模型:

(要求p<m)

例：表1为某地区农业生态经济戚消明系统各区域单元相关指标数据，运用主成分分析方法可以用更少的指标信息较为精确地描述该地区农业生态经济的发展状况。

表1 某农业生态经济系统各区域单元的有关数据

样本序号 x1：人口密度(人/km2) x 2：人均耕地面积(ha) x 3：森林覆盖率(%) x 4：农民人均纯收入(元/人) x 5：人均粮食产量 (kg/人) x 6：经济作物占高告农作物播面比例(％) x 7：耕地占土地面积比率(％) x 8：果园与林地面积之比(％) x 9：灌溉田占耕地面积之比(％)

1 363.912 0.352 16.101 192.11 295.34 26.724 18.492 2.231 26.262

2 141.503 1.684 24.301 1 752.35 452.26 32.314 14.464 1.455 27.066

3 100.695 1.067 65.601 1 181.54 270.12 18.266 0.162 7.474 12.489

4 143.739 1.336 33.205 1 436.12 354.26 17.486 11.805 1.892 17.534

5 131.412 1.623 16.607 1 405.09 586.59 40.683 14.401 0.303 22.932

6 68.337 2.032 76.204 1 540.29 216.39 8.128 4.065 0.011 4.861

7 95.416 0.801 71.106 926.35 291.52 8.135 4.063 0.012 4.862

8 62.901 1.652 73.307 1 501.24 225.25 18.352 2.645 0.034 3.201

9 86.624 0.841 68.904 897.36 196.37 16.861 5.176 0.055 6.167

10 91.394 0.812 66.502 911.24 226.51 18.279 5.643 0.076 4.477

11 76.912 0.858 50.302 103.52 217.09 19.793 4.881 0.001 6.165

12 51.274 1.041 64.609 968.33 181.38 4.005 4.066 0.015 5.402

13 68.831 0.836 62.804 957.14 194.04 9.110 4.484 0.002 5.790

14 77.301 0.623 60.102 824.37 188.09 19.409 5.721 5.055 8.413

15 76.948 1.022 68.001 1 255.42 211.55 11.102 3.133 0.010 3.425

16 99.265 0.654 60.702 1 251.03 220.91 4.383 4.615 0.011 5.593

17 118.505 0.661 63.304 1 246.47 242.16 10.706 6.053 0.154 8.701

18 141.473 0.737 54.206 814.21 193.46 11.419 6.442 0.012 12.945

19 137.761 0.598 55.901 1 124.05 228.44 9.521 7.881 0.069 12.654

20 117.612 1.245 54.503 805.67 175.23 18.106 5.789 0.048 8.461

21 122.781 0.731 49.102 1 313.11 236.29 26.724 7.162 0.092 10.078

对于上述例子，Matlab进行主成分分析，可以得到如下结果。

① 以及每一个主成分的贡献率和累计贡献率，如表2和图1。

表2. 特征根及主成分贡献率

主成分 特征值 贡献率％累积贡献率％

1 4.661 51.791 51.791

2 2.089 23.216 75.007

3 1.043 11.589 86.596

4 0.507 5.638 92.234

5 0.315 3.502 95.736

6 0.193 2.140 97.876

7 0.114 1.271 99.147

8 4.533E-02 0.504 99.650

9 3.147E-02 0.350 100.000

图1 特征根

② 前3几个主成分的载荷系数如表3所示。

表3 前三个主成分在原变量上的载荷

前三个主成分

变量 1 2 3

X1 0.158 -0.255 -0.059

X2 0.026 0.424 -0.027

X3 -0.207 0.046 0.091

X4 0.009 0.415 0.036

X5 0.174 0.212 -0.011

X6 0.176 0.086 0.120

X7 0.200 -0.064 -0.241

X8 0.042 -0.048 0.930

X9 0.207 -0.012 0.088

clc,clear

load x.txt %原始的x组的数据氏桐保存在纯文本文件x.txt中

load y.txt %原始的y组的数据保存在纯文本文件y.txt中

n1=size(x,2)n2=size(y,2)

x=zscore(x)y=zscore(y) %标准化数据

n=size(x,1)

%a,b返回的是典型变量的系数，r返回的是典型相关系数

%u,v返回的是典型变量的值，stats返回的是假设检歼晌坦验的一些统计量的值

[a,b,r,u,v,stats]=canoncorr(x,y)

x_u_r=x'*u/(n-1) %计算x,u的相关系数

y_v_r=y'*v/(n-1) %计算y,v的相关系谨简数

x_v_r=x'*v/(n-1) %计算x,v的相关系数

y_u_r=y'*u/(n-1) %计算y,u的相关系数

mu=sum(x_u_r.^2)/n1 %x组原始变量被u_i解释的方差比例

mv=sum(x_v_r.^2)/n1 %x组原始变量被v_i解释的方差比例

nu=sum(y_u_r.^2)/n2 %y组原始变量被u_i解释的方差比例

nv=sum(y_v_r.^2)/n2 %y组原始变量被v_i解释的方差比例

val=r.^2 %典型系数的平方

---