No BSD License  

Highlights from
advstat

advstat

by

 

10 Nov 2000 (Updated )

Adds some additional capabilities to the Stats TB

[b0,b,B]=leasrsm(x,type)
function  [b0,b,B]=leasrsm(x,type)
%function  [b0,b,B]=leasrsm(x{,type})
% x=data file [x1 x2 ..xn y]
% Uses LSCurvefit FROM OPT PACKAGE
% type= 'q' for full BOX composite model with quad terms,else Factorial
% model with only cross terms is fit {}=optional
% This program sets up the data for a call to modelrsm.m
% for Least Squares calc of parameters for General Quadratic
% Model ===> used in Response Surfaces	y=b0 + x*b +x*B*x'
% Uses Subroutines: LSQCURVEFIT ,MODELRSM,SYMUNPCK
%	    A.Jutan. UWO '98
% modified for 5.3 4-28-99,5-19-99,for6.1 8-23-01

[l,w]=size(x);
y=x(:,w);
xx=x;xx(:,w)=[];  % remove last col store x's in xx
w=w-1; % Number of x's
% INITIAL GUESS'S for params in b0 + x*b +x*B*x'
b0=1.0;
b=2.0*ones(w,1);
B=ones(w,w); % guess ones matrix no zeros
%**********************************************************
% FOR PURE FACTORIAL DATA I.E. NO BOX COMPOSITE (STAR) PTS
%**********************************************************
% SET DIAGONAL OF GUESSED B TO ZEROS
if( nargin == 1)
 B=B -diag(diag(B));
end;
%***********************************************************
%
% PACK GUESS INTO ONE LONG VECTOR
%
B=tril(B); % take lower half only since B symmetric
p1=B(:); % form long vector--includes zeros
i=find(p1==0);% find the zero positions
p1(i)=[]; % eliminate the zero's ==> reduces p vector
% pack into single column vector: length (n^2 +n)/2
% pin=[b0 b1 b2 b3 B11 B12 B13 B22 B23 B33]' for 3 x's (n=3)
% pin=[b0 b1 b2 b3  B12 B13  B23 ] for 3 x s no SQUARE terms
pin=[b0;b(:);p1(:)]; % set up as long vector incl b0 for init param guess
% CHECK FOR SUFFICIENT DATA TO EST PARAMS
lpin=length(pin);
if (l<lpin),error('NOT ENOUGH DATA TO ESTIMATE ALL PARAMETERS'),end;
options=optimset('LevenbergMarquardt','on');%,'TolX',1e-4,'MaxFunEvals',1000)%added aj 5-18-99
[p,resn,error,exf,out,ll,jac]=lsqcurvefit(@modelrsm,pin,xx,y,[],[],options);
ymodel=y+error;
[std,varresid,r2,cor,vcv,varinf]=regdata(p,ymodel,y,jac);
% reshape parameter matrices
disp(' parameter matrices and 95% Confidence Limits +/- CL95')
b0=p(1)
stdb0=std(1)	 ;CL95_b0=2*stdb0
b=p(2:w+1)';
stdb=std(2:w+1)' ;CL95_b=2*stdb
lp=length(p);
pp=p(2+w:lp); % select only Bij terms in p
stdpp=std(2+w:lp); % select only Bij terms in std
lpp=length(pp);
% check for diag B=0 ie NO SQUARE TERMS
if lpp~=(w^2+w)/2 ,w1=w-1;else,w1=w;,end; % Unpack B without diag(quad) PARAMS
BF=symunpck(pp,w1);% recreates symmetric matrix of B from vector form in pp
B=tril(BF);    % take lower half
stdB=tril( symunpck(stdpp,w1));
% Add Back Zero Diagonals for no quadratic terms
      if lpp~=(w^2+w)/2
      B=[zeros(1,w);B,zeros(w1,1)]; % advanced fiddle see p.2-39 MATLAB
      stdB=[zeros(1,w);stdB,zeros(w1,1)];
      end
B  % show B matrix and 95% CL
CL95_B=2* stdB
disp('p=[b0 b1 b2 b3 B11 B12 B13 B22 B23 B33] for 3 x s (n=3)')
disp('p=[b0 b1 b2 b3  0  B12 B13  0  B23  0 ] for 3 x s no SQUARE terms')
disp('b12=2*B12 in y= b0 + b1.x1 +b2.x2 +b12.x1.x2 + b11.x1^2+b22.x2^2')
disp(' non cross terms b11 b22 equal to B11 B22,etc')

% Display results
% MAKE SURE YOU USE FULL B=BF (not B from leasqr1 =lower B)


%================
function y=yrsm(x,b0,b,BF)
% MAKE SURE YOU USE FULL B (not B from leasqr1 =lower B)
xc=x';

if nargin > 1 , b0_=b0;b_=b;B_=BF; end
y=b0_ + xc*b_' +xc*B_*xc'; % calc y for rsm model

%=====================================
function y=modelrsm(p,x)
% To calc general Quadratic Form y=b0 + x*b +x*B*x'
% Includes all main effects , 2 cross products x(i)*x(j) also x(i)*x(i) (quads)
% used in Response Surface Methodology
% If no Quadratic terms {length(p) is < (w^2+w)/2}  then diag(B) is zeroed
% p=vector(b0,b,B)===> need to reshape p from long vector
% B was further packed since it is symmetric and zeros removed
% p=[b0 b1 b2 b3 B11 B12 B13 B22 B23 B33] for 3 x's (w=3)
% ALL CROSS PARAM'S E.G. B12 (NOT B11) ==> USE 2*B12 IN y=b0 +b1.x +b12.x1.x2
% ie b12=2*B12

% A Jutan May 98

[l,w]=size(x);% order of B =w*w
lp=length(p);
b0=p(1);
b=p(2:w+1);
pp=p(2+w:lp); % select only Bij terms in p
% Check for zero diag in B i.e. no Square terms ***PURE FACTORIAL DATA****
if length(pp) ~=(w^2 +w)/2
   % diags on B=0 reconstruct one order lower B then add diag of zeros
   a1=symunpck(pp,w-1);
   a2=tril(a1); % take lower triangle
   a3=[zeros(1,w-1);a2]; %add top row of zeros
   a4=[a3 , zeros(w,1)]; %add RHS row of zeros
   B=a4+a4';% form symmetric B with zeros on diag
   else
   B=symunpck(pp,w); % recreates symmetric matrix of B from vector form in pp
end
for i=1:l
y(i)=b0 + x(i,:)*b + x(i,:)*B*x(i,:)';
end
y=y';	% make sure col vector returned to leasqr
return
%===============
function a=symunpck(p,n)
%function a=symunpck(p,[n])
% reforms first lower diag matrix a from vector	in p
% then full symmetric matrix
% p contains cols of original a without zeros on upper triangle
% n is order of a, p must contain full lower triangle as a vector
% If n not given it is calculated from length of p=n(n+1)/2

if nargin==1,n=(-1 + sqrt(1+ 8.* (length(p))))/2;end
p=p(:); %form column vector
for i=1:n
a(:,i)=[zeros(i-1,1);p(1:n-i+1)];
p(1:n-i+1)=[]; %  cut of retreived part
end
% a is now lower diag. Can form full symmetric a as:
c=a+a';% c is almost full symm but diag valuels 2* true values
d=diag(diag(c)); % gives diag matrix with 2*diagonal values
a=c-d/2;

% In some applications we need to match the  quadratic form with the 
% polynomial form:
% in order for x'ax to match cross terms in the polynomial form
% the off diagonal elements in 'a'
% must be divided by 2. Diagonal terms are not affected
%a=(a-diag(diag(a)))/2 + diag(diag(a)); % sub diag divide cross terms by
				       % 2 then add back diag matrix

Contact us