急！急！小生是MATLAB新手，跪求高手——偏最小二乘回归问题

程成RQ

PLSR.rar (11.68 KB, 下载次数: 1)

2012-12-29 16:28 上传

点击文件名下载附件
里面有多个文件

                                                                       偏最小二乘回归问题（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
          [theta,w,cw,ssqdif,yres]=plsr(xreg,yreg,nox+1,lv);

Error in ==> PLSR_yushengrong at 11
[lv,theta,ycal,t,Tcrit,STATUS,ypre]=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.
[lv,theta,ycal,t,Tcrit,STATUS,ypre]=pressf(X,Y,xpre,alpha)   %调用函数pressf.

_________________pressf.m___________

function[lv,theta,ycal,t,Tcrit,STATUS,ypre]=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的理论临界值)时模型可以通过显著性检验。
[nrx,nox]=size(X);
[nry,noy]=size(Y);
ax=autoscaling(X);
XX=[ones(nrx,1) ax;ones(nrx,1) ax];
YY=[Y;Y];
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
            [theta,w,cw,ssqdif,yres]=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随着潜变量个数的变化图
[pressvalue,a]=min(press);         %a=press最小时的潜变量个数
lv=a;
[theta,ssqdif,yres,t,Tcrit,R,F,Fcrit,ypre]=pls_test(X,Y,lv,xpre,alpha);
                            %调用pls_test自编函数计算最佳潜变量下的plsr回归系数并对回归模型和回归系数进行统计检验
ycal=[ones(nrx,1) ax]*theta';
RY=corrcoef(Y,ycal);
r=RY(1,2);                         %计算最佳LV下建模集模型值与实际值之的相关系数
RMSEC=sqrt(sum((Y-ycal).^2)/nry);  %计算建模集均方残差
e=abs(yres);                       %计算建模集因变量的残差绝对值
rmse=mean(e);                      %计算建模集因变量的平均绝对残差
STATUS=[r RMSEC rmse R F Fcrit];
                                   %将最佳plsr回归模型的评价指标赋给STAUS数组
%end of pressf

______pls_test.m________
function [theta,ssqdif,yres,t,Tcrit,R,F,Fcrit,yu]=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回归模型给出的未知样本的因变量值

[nx,px]=size(xreg);
[ny,py]=size(yreg);
if nx ~= ny
    error('建模集的自变量与因变量的样本数必须相等！');
    return
end
[ax,mx,s]=autoscaling(xreg);
% 对建模集自变量矩阵（每一列为一个自变量）进行自标度化处理并将结果赋给矩阵ax，返回参数mx与s分别为各变量的均值与标准方差
% 自编函数autoscaling
[theta,w,cw,ssqdif,yres]=plsr([ones(nx,1) ax],yreg,px+1,LV);
% 以含常数列的ax为自变量，yreg为因变量进行偏最小二乘回归
Yc=[ones(nx,1) ax]*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值*****%%%%%%%%%%%%%%%%%%%%
[np,pr]=size(xpre);
%将矩阵xpre的行数和列数分别赋给变量np，pr
axp=(xpre-repmat(mx,np,1))./repmat(s,np,1);
%将预测集样的自变量减建模集对应变量的均值并除以该变量的方差（对预测样本进行自标度化预处理）
yu=[ones(np,1) axp]*theta';
%求未知样本的因变量yu
%end of pls_test

______autoscaling.m________
function[autoX,mx,s]=autoscaling(X) %定义一个自标度化函数，自标度化后的变量赋给矩阵autoX，变量的均值与标准方差分别返回给mx与s
[nh,nl]=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 [theta,w,cw,ssqdif,yres] = plsr(xreg,yreg,ninput,lv,plotopt);
%PLSR Determine impulse response coefficients via Partial Least Squares.
%    [theta,w,cw,ssqdif,yres] = 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: [theta,w,cw,ssqdif,yres] = 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
[nrx,ncx] = size(xreg);
[nry,ncy] = 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