急！急！小生是MATLAB新手，跪求高手——偏最小二乘回归问题
本帖最后由 程成RQ 于 2012-12-29 16:33 编辑偏最小二乘回归问题（plsr）：
运行程序出现："??? Subscripted assignment dimension mismatch."错误提醒…………
我的plsr的MATLAB实施程序主要参考书上的，按照书上例子给的数据，运行程序过程没有错误提醒，结果和书上一样！
可是，换为我的实验数据，就出现错误提醒，搞了好就久都搞不明天，也在各个论坛搜索了相关帖子，还在没解决！！
——————————下面给出我的M文件说明（传文件时才发现不支持m文件，所以我都转为文本文件了）
1.PLSR_.m——函数初始化（使用的书本例子数据，其中调用pressf函数）；2.PLSR_yushengrong.m——函数初始化（使用自己的实验数据，其中调用pressf函数）；3.pressf.m——定义pressf函数，其中调用到函数pls_test和函数autoscaling；4.autoscaling.m——定义autoscaling函数；5.pls_test.m——定义pls_test.mf函数，其中调用到MATLAB工具箱自带函数plsr。6.顺便也上传MATLAB自带函数plsr的m文件——plsr.m。
以下是函数初始化加载到的数据说明：
书本例子数据_shujuX.mat和shujuY.mat;自己的实验数据：shuju_X.mat和shuju_Y.mat。每个mat文件下包含有与mat同名的变量。
——————————下面是遇到的问题————————
运行PLSR_.m文件，出结果，和书本一样；运行PLSR_yushengrong.m文件，出现错误提醒，如下：
??? Subscripted assignment dimension mismatch.
Error in ==> plsr at 54
w(:,i) = x'*y;
Error in ==> pressf at 31
=plsr(xreg,yreg,nox+1,lv);
Error in ==> PLSR_yushengrong at 11
=pressf(X,Y,xpre,alpha) %调用函数pressf.
————————————求高手指教——————
小生是MATLAB入门新手，我这问题，相信高手看来真不值一提，小菜一碟摆了！但是如果您们帮我解决了此难题，那对我的帮助可不小，小生必会万分感激！今后也向你们学习，热心帮助入门生！
感谢看帖的人，更感谢那些为他人伸出援助之手的人，有你们，【中国振动联盟】必定更强大！预祝祝大家元旦快乐，生活幸福美满，工作顺顺利利，2013更上一层楼！
补充内容 (2012-12-31 10:27):
__________PLSR_yushengrong.m___________
clear
clc
load shuju_X; %加载两个M文件，shuju_X.m文件里面包含有shuju_X变量，
load shuju_Y; %shuju_Y.m文件里面包含有shuju_Y变量.
%%%%%%%******以下对函数pressf输入变量进行初始化*****%%%%%
X=shuju_X; %把变量shuju_X赋予函数pressf的输入变量X
Y=shuju_Y; %把变量shuju_Y赋予函数pressf的输入变量Y
xpre=X; %xpre为预测集自变量矩阵，这里设置和变量X相同
alpha=0.05; %alpha为置信水平，这里取0.05.
=pressf(X,Y,xpre,alpha) %调用函数pressf.
_________________pressf.m___________
function=pressf(X,Y,xpre,alpha)
% pressf(留一法确定潜变量个数)
% 输入变量：
% X,Y: 自变量与因变量矩阵；xpre为预测集自变量矩阵
% 输出变量：
% lv: 最佳PLS潜变量个数，theta：PLSR回归系数向量，含义同plsr.m
% ycal-根据优化的lv所确定的plsr模型给出的因变量值；
% t为对theta给出的各plsr回归系数进行显著性检验的统计量t（不含常数项），其大小=nox（X的列数）；
% Tcrit为t双边分布函数在nox - lv - 1 自由度下的临界值（置信水平=alpha）
% STATUS存储最佳lv值下所得plsr模型的评价指标：其第1个元素r为建模样本的模型值与实测值之间的相关系数；
% 第2个元素为RMSEC（建模样本的均方根误差）；第3个元素e为模型给出的因变量值与实际值的绝对平方残差；
% 第4个元素R为plsr回归模型的复相关系数；
% 第5个元素F为plsr回归模型检验的统计量F比，F大于第6个元素Fcrit(Fd的理论临界值)时模型可以通过显著性检验。
=size(X);
=size(Y);
ax=autoscaling(X);
XX=;
YY=;
if nrx~= nry
error('自变量样本个数必须等于因变量的样本个数！');
return
end
for j=1:nrx
for i=j:(nrx+(j-2))
xreg(i,:)=XX(i,:);
yreg(i,:)=YY(i,:);
end
for lv=1:nox
=plsr(xreg,yreg,nox+1,lv);
ypredict(j,lv)=XX(nrx+(j-1),:)*theta'; %计算每个LV下留一检验样本的y值
end
j=j+1;
end
for v =1:nox
press(:,v)=sum((ypredict(:,v)-Y).^2);
% 计算对应的潜变量下的留一检验样本的因变量残差平方和press
end
plot(1:nox,press,'-o') %画出press随着潜变量个数的变化图
=min(press); %a=press最小时的潜变量个数
lv=a;
=pls_test(X,Y,lv,xpre,alpha);
%调用pls_test自编函数计算最佳潜变量下的plsr回归系数并对回归模型和回归系数进行统计检验
ycal=*theta';
RY=corrcoef(Y,ycal);
r=RY(1,2); %计算最佳LV下建模集模型值与实际值之的相关系数
RMSEC=sqrt(sum((Y-ycal).^2)/nry);%计算建模集均方残差
e=abs(yres); %计算建模集因变量的残差绝对值
rmse=mean(e); %计算建模集因变量的平均绝对残差
STATUS=;
%将最佳plsr回归模型的评价指标赋给STAUS数组
%end of pressf
______pls_test.m________
function =pls_test(xreg,yreg,LV,xpre,alpha)
% 偏小二乘回归建模、预测并对回归模型和回归系数进行检验的程序
% xreg为建模集自变量；yreg为建模集因变量；alpha为置信水平，取0.05或0.01
% LV为潜变量个数（LV﹦﹤自变量个数）
% xpre为预测集的自量矩阵（如果不做预测，xpre输xreg即可）
% 输出变量：theta为plsr回归系数组成的行向量，其个数=自变量个数+1
% ssqdif:PLSR解释的自变量和因变量方差百分率；yres：因变量残差.
% t为按照bi/(c^1/2*(S/n-m-1)^1/2)计算的LV个潜变量的回归系数统计检验量，Tcrit为t分布的临界值
% R为plsr回归模型的复相关系数；
% F为plsr回归所得回归模型的方差比值统计量，Fcrit为F分布的临界值
% yu为plsr回归模型给出的未知样本的因变量值
=size(xreg);
=size(yreg);
if nx ~= ny
error('建模集的自变量与因变量的样本数必须相等！');
return
end
=autoscaling(xreg);
% 对建模集自变量矩阵（每一列为一个自变量）进行自标度化处理并将结果赋给矩阵ax，返回参数mx与s分别为各变量的均值与标准方差
% 自编函数autoscaling
=plsr(,yreg,px+1,LV);
% 以含常数列的ax为自变量，yreg为因变量进行偏最小二乘回归
Yc=*theta';%求建模集因变量并赋给Yc
ym=mean(yreg);%求因变量的均值
Sc=sum(yres.^2);%求建模集残差平方和
Sh=sum((Yc-repmat(ym,ny,1)).^2);%求回归平方和
R=sqrt(Sh/sum((yreg-repmat(ym,ny,1)).^2));
%求plsr回归的模型的复相关系数R
F=(nx-LV-1)*Sh/(LV*Sc);%求plsr回归模型的统计检验量F比值
Fcrit=finv(1-alpha,LV,nx-LV-1);%求F分布的临界值
R,F,Fcrit %屏幕显示R,F,Fcrit的值
%%%%%%%%%%%%*****以下计算plsr回归系数（不含常数项）显著性检验的统计参数t值*****%%%%%%%
cc=inv(xreg'*xreg);
B=theta(:,2:px+1)';
t=B./(sqrt(diag(cc))*sqrt(sum(yres.^2))/(nx-px-1));
%根据式bi/(c^1/2*(S/n-m-1)^1/2)计算实际t所组成的列向量（不含常数项）
Tcrit=tinv((1-alpha/2),nx-px-1);%求t的临界值
t,Tcrit%屏幕显示t,Tcrit的值
%%%%%%%%%%%%%*****下面根据建立的plsr回归模型预测未知样本的y值*****%%%%%%%%%%%%%%%%%%%%
=size(xpre);
%将矩阵xpre的行数和列数分别赋给变量np，pr
axp=(xpre-repmat(mx,np,1))./repmat(s,np,1);
%将预测集样的自变量减建模集对应变量的均值并除以该变量的方差（对预测样本进行自标度化预处理）
yu=*theta';
%求未知样本的因变量yu
%end of pls_test
______autoscaling.m________
function=autoscaling(X) %定义一个自标度化函数，自标度化后的变量赋给矩阵autoX，变量的均值与标准方差分别返回给mx与s
=size(X);%确定矩阵X的行数（nh）和列数（nl）
mx=mean(X);%求X中每个变量的均值mx（是一个列数为nl的行向量）
s=std(X);%求X中每个变量的标准方差stdx（是一个列数为nl的行向量）
autoX=(X-repmat(mx,nh,1))./repmat(s,nh,1);%对X中每个变量进行自标度化预处理
_____matlab自带函数plsr.m__
function = plsr(xreg,yreg,ninput,lv,plotopt);
%PLSR Determine impulse response coefficients via Partial Least Squares.
% = plsr(xreg,yreg,ninput,lv,plotopt)
% Required Inputs:
% xreg, yreg: input & output data.
% ninput: number of input.
% lv: number of latent vectors, ie., number of directions chosen
% to maximize the cross variance between the input and output
% data in PLS routine.
% Optional Inputs:
% plotopt = 0 (default), no plot is produced;
% = 1, plot of actual output (yreg) and predicted output;
% = 2, plot of yreg and yfit, and plot of residues (yres).
% Outputs:
% w: orthogonal matrix of dimension n by lv consisting of
% principle directions maximizing the cross variance between
%input and output.n is number of impulse response coeff.
% cw: a vector containing coefficients associated with each column
% of w in determining theta.
% theta: impulse response coefficient matrix of dimension n by nu.
% nu is number of inputs. Column i corresponds to the impulse
%response coefficients for input i.
% ssqdif: percentage variances of input & output captured by PLS.
% yres: residue output.
% See also MLR, VALIDMOD, WRTREG.
% Copyright 1994-2003 The MathWorks, Inc.
% $Revision: 1.1.6.2 $
if nargin == 0
disp('Usage: = plsr(xreg,yreg,ninput,lv,plotopt)')
return
end
if nargin < 4
error('Minimum number of input arguments is 4!');
return
end
if nargin == 4
plotopt = 0;
end
= size(xreg);
= size(yreg);
if nrx ~= nry
error('Number of rows for input and output must be same!');
return
end
ssqx = sum(sum(xreg.^2)');
ssqy = sum(sum(yreg.^2)');
x = xreg;
y = yreg;
I = eye(ncx);
w = zeros(ncx,lv);
for i = 1:lv
w(:,i) = x'*y;
w(:,i) = w(:,i)/norm(w(:,i));
r(:,i) = x * w(:,i);
x = x -r(:,i) * w(:,i)';
y = y - r*(r\y);
ssq(i,1) = (sum(sum(x.^2)'))*100/ssqx;
ssq(i,2) = (sum(sum(y.^2)'))*100/ssqy;
end
cw = r\yreg;
thetam = w * cw;
ssqdif = zeros(lv,2);
ssqdif(1,1) = 100 - ssq(1,1);
ssqdif(1,2) = 100 - ssq(1,2);
for i = 2:lv
for j = 1:2
ssqdif(i,j) = -ssq(i,j) + ssq(i-1,j);
end
end
disp('')
disp(' Percent Variance Captured by PLS Model')
disp('')
disp(' ----X-Block------ ----Y-Block------')
disp(' LV# This LV Total This LV Total ')
disp([(1:lv)' ssqdif(:,1) cumsum(ssqdif(:,1)) ssqdif(:,2) cumsum(ssqdif(:,2))])
ypred = xreg * thetam;
yres = yreg - ypred;
rms = norm(yres);
if plotopt == 1 | plotopt == 2
clf
if plotopt == 2
subplot(211);
plot(1:nrx,yreg,1:nrx,yreg,'or',1:nrx,ypred,1:nrx,ypred,'+g');
title('Actual value (o) versus Predicted Value (+)');
xlabel('Sample Number');
ylabel('Output');
subplot(212);
plot(yres);
title('Output Residual or Prediction Error');
xlabel('Sample Number');
ylabel('Residual');
else
plot(1:nrx,yreg,1:nrx,yreg,'or',1:nrx,ypred,1:nrx,ypred,'+g');
title('Actual value (o) versus Predicted Value (+)');
xlabel('Sample Number');
ylabel('Output');
end
end
n = ncx / ninput;
for i = 1:ninput
theta(:,i) = thetam(n*(i-1)+1:n*i);
end
%End of function plsr
页:
[1]