偏最小二乘(plsr)回归问题 ::"??? Subscripted assignment dimension mismatch."
运行程序出现:"??? Subscripted assignment dimension mismatch."错误提醒…………我的plsr的MATLAB实施程序主要参考书上的,按照书上例子给的数据,运行程序过程没有错误提醒,结果和书上一样!
可是,换为我的实验数据,就出现错误提醒,搞了好就久都搞不明天,也在各个论坛搜索了相关帖子,还在没解决!!
——————————下面给出我的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入门新手,我这问题,相信高手看来真不值一提,小菜一碟摆了!但是如果您们帮我解决了此难题,那对我的帮助可不小,小生必会万分感激!今后也向你们学习,热心帮助入门生!
感谢看帖的人,更感谢那些为他人伸出援助之手的人,有你们,【BCCN】必将更强大!预祝祝大家元旦快乐,生活幸福美满,工作顺顺利利,2013更上一层楼!!
PLSR.rar
(11.93 KB)
————补充——————
__________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.
[ 本帖最后由 程成RQ 于 2012-12-31 10:30 编辑 ]