% Author: Gerhard Nieuwenhuiys
% Date of last revision : 2009/02/04
% Thank you for downloading / using mregres V1.
% mregres can only handle explicit functions of the form z=f(x,y). Should you
% like to obtain more information on its operating instructions or you would like
% to give comments for improvement, feel free to contact me at gerhardn654@yahoo.com
% mregres need the following inputs:
% fstring = can have up to 10 independent parameters. z=f(x,y)
% v = starting estimates for parameters. Number of elements in v must equal
% number of parameters in fstring
% xexp = Experimental values of the first dimension
% yexp = Experimental values of the second dimension
% zexp = Experimental values of the third dimension
%
% Example 1:
% >>xexp=[-3 -2 -1 1 2 3 -3 -2 -1 1 2 3 -3 -2 -1 1 2 3];
% >>yexp=[-3 -3 -3 -3 -3 -3 -2 -2 -2 -2 -2 -2 -1 -1 -1 -1 -1 -1];
% >>zexp=[17.5 13.1 10.1 10.2 13.2 17.2 13.1 7.9 4.8 4.7 7.8 12.7 9.8 5.2 2.1 2.3 4.8 10.2];
% >>mregres('p1.*x.^2+p2.*y.^2',[1 1],xexp,yexp,zexp);
% p1=0.97953 p2=0.99681 R2=0.99649
%
% Notice how the vectors are put together. The second experimental point is
% [-2,-3,13.1] as [x,y,z]. You would then construct the inputs as xexp=[-2 ...]
% yexp=[-3 ...] and zexp=[13.1 ...]
function mregres(fstring,v,xexp,yexp,zexp);
[w1,w2,w3,w4,w5,w6,w7,w8,w9,w10]=deal(strfind(fstring,'p1'),strfind(fstring,'p2'),strfind(fstring,'p3'),...
strfind(fstring,'p4'),strfind(fstring,'p5'),strfind(fstring,'p6'),strfind(fstring,'p7'),strfind(fstring,'p8'),...
strfind(fstring,'p9'),strfind(fstring,'p10'));
[k1,k2,k3,k4,k5,k6,k7,k8,k9,k10]=deal(length(w1),length(w2),length(w3),length(w4),length(w5),length(w6),length(w7),...
length(w8),length(w9),length(w10));
szx=size(xexp);
szy=size(yexp);
[x,fval]=fminsearch(@errfun,v,[],fstring,xexp,yexp,zexp,szy,szx);
[xi,yi]=meshgrid(min(xexp):(max(xexp)-min(xexp))/30:max(xexp),min(yexp):(max(yexp)-min(yexp))/30:max(yexp));
[yi2,xi2]=meshgrid(min(yexp):(max(yexp)-min(yexp))/30:max(yexp),min(xexp):(max(xexp)-min(xexp))/30:max(xexp));
sse=fval;
zs=sum(sum(zexp))/length(zexp+yexp+xexp);
sst=sum(sum((zexp-zs).^2));
n=length(zexp+yexp+xexp);
if k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1 & k9>=1 & k10>=1;
p=10;
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1 & k9>=1;
p=9;
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1;
p=8;
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1;
p=7;
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1;
p=6;
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1;
p=5;
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1;
p=4;
elseif k1>=1 & k2>=1 & k3>=1;
p=3;
elseif k1>=1 & k2>=1;
p=2;
elseif k1>=1;
p=1;
end;
R2=1-((sse/(n-p))/(sst/(n-1)));
if k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1 & k9>=1 & k10>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' p3=',num2str(x(3)),' p4=',num2str(x(4)),' p5=',num2str(x(5)),...
' p6=',num2str(x(6)),' p7=',num2str(x(7)),' p8=',num2str(x(8)),' p9=',num2str(x(9)),' p10=',num2str(x(10)),...
' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8','p9','p10');
zi=f(xi,yi,x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10));
f2=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8','p9','p10');
zi2=f2(xi2,yi2,x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10));
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1 & k9>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' p3=',num2str(x(3)),' p4=',num2str(x(4)),' p5=',num2str(x(5)),...
' p6=',num2str(x(6)),' p7=',num2str(x(7)),' p8=',num2str(x(8)),' p9=',num2str(x(9)),' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8','p9');
zi=f(xi,yi,x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9));
f2=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8','p9');
zi2=f2(xi2,yi2,x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9));
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' p3=',num2str(x(3)),' p4=',num2str(x(4)),' p5=',num2str(x(5)),...
' p6=',num2str(x(6)),' p7=',num2str(x(7)),' p8=',num2str(x(8)),' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8');
zi=f(xi,yi,x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8));
f2=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8');
zi2=f2(xi2,yi2,x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8));
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' p3=',num2str(x(3)),' p4=',num2str(x(4)),' p5=',num2str(x(5)),...
' p6=',num2str(x(6)),' p7=',num2str(x(7)),' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7');
zi=f(xi,yi,x(1),x(2),x(3),x(4),x(5),x(6),x(7));
f2=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7');
zi2=f2(xi2,yi2,x(1),x(2),x(3),x(4),x(5),x(6),x(7));
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' p3=',num2str(x(3)),' p4=',num2str(x(4)),' p5=',num2str(x(5)),...
' p6=',num2str(x(6)),' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6');
zi=f(xi,yi,x(1),x(2),x(3),x(4),x(5),x(6));
f2=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6');
zi2=f2(xi2,yi2,x(1),x(2),x(3),x(4),x(5),x(6));
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' p3=',num2str(x(3)),' p4=',num2str(x(4)),' p5=',num2str(x(5)),...
' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2','p3','p4','p5');
zi=f(xi,yi,x(1),x(2),x(3),x(4),x(5));
f2=inline(fstring,'x','y','p1','p2','p3','p4','p5');
zi2=f2(xi2,yi2,x(1),x(2),x(3),x(4),x(5));
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' p3=',num2str(x(3)),' p4=',num2str(x(4)),' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2','p3','p4');
zi=f(xi,yi,x(1),x(2),x(3),x(4));
f2=inline(fstring,'x','y','p1','p2','p3','p4');
zi2=f2(xi2,yi2,x(1),x(2),x(3),x(4));
elseif k1>=1 & k2>=1 & k3>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' p3=',num2str(x(3)),' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2','p3');
zi=f(xi,yi,x(1),x(2),x(3));
f2=inline(fstring,'x','y','p1','p2','p3');
zi2=f2(xi2,yi2,x(1),x(2),x(3));
elseif k1>=1 & k2>=1;
disp([' p1=',num2str(x(1)),' p2=',num2str(x(2)),' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1','p2');
zi=f(xi,yi,x(1),x(2));
f2=inline(fstring,'x','y','p1','p2');
zi2=f2(xi2,yi2,x(1),x(2));
elseif k1>=1;
disp([' p1=',num2str(x(1)),' R2=',num2str(R2)]);
f=inline(fstring,'x','y','p1');
zi=f(xi,yi,x(1));
f2=inline(fstring,'x','y','p1');
zi2=f2(xi2,yi2,x(1));
end;
plot3(xexp,yexp,zexp,'r.',xi,yi,zi,'b-',xi2,yi2,zi2,'b-'),ylabel('y'),xlabel('x'),zlabel('z');
function err=errfun(x,fstring,xexp,yexp,zexp,szy,szx);
[w1,w2,w3,w4,w5,w6,w7,w8,w9,w10]=deal(strfind(fstring,'p1'),strfind(fstring,'p2'),strfind(fstring,'p3'),...
strfind(fstring,'p4'),strfind(fstring,'p5'),strfind(fstring,'p6'),strfind(fstring,'p7'),strfind(fstring,'p8'),...
strfind(fstring,'p9'),strfind(fstring,'p10'));
[k1,k2,k3,k4,k5,k6,k7,k8,k9,k10]=deal(length(w1),length(w2),length(w3),length(w4),length(w5),length(w6),length(w7),...
length(w8),length(w9),length(w10));
if k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1 & k9>=1 & k10>=1;
[p1,p2,p3,p4,p5,p6,p7,p8,p9,p10]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)),x(3)*ones(szy(1),szx(2)),...
x(4)*ones(szy(1),szx(2)),x(5)*ones(szy(1),szx(2)),x(6)*ones(szy(1),szx(2)),x(7)*ones(szy(1),szx(2)),...
x(8)*ones(szy(1),szx(2)),x(9)*ones(szy(1),szx(2)),x(10)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8','p9','p10');
g=f(xexp,yexp,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10);
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1 & k9>=1;
[p1,p2,p3,p4,p5,p6,p7,p8,p9]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)),x(3)*ones(szy(1),szx(2)),...
x(4)*ones(szy(1),szx(2)),x(5)*ones(szy(1),szx(2)),x(6)*ones(szy(1),szx(2)),x(7)*ones(szy(1),szx(2)),...
x(8)*ones(szy(1),szx(2)),x(9)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8','p9');
g=f(xexp,yexp,p1,p2,p3,p4,p5,p6,p7,p8,p9);
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1 & k8>=1;
[p1,p2,p3,p4,p5,p6,p7,p8]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)),x(3)*ones(szy(1),szx(2)),...
x(4)*ones(szy(1),szx(2)),x(5)*ones(szy(1),szx(2)),x(6)*ones(szy(1),szx(2)),x(7)*ones(szy(1),szx(2)),...
x(8)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7','p8');
g=f(xexp,yexp,p1,p2,p3,p4,p5,p6,p7,p8);
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1 & k7>=1;
[p1,p2,p3,p4,p5,p6,p7]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)),x(3)*ones(szy(1),szx(2)),...
x(4)*ones(szy(1),szx(2)),x(5)*ones(szy(1),szx(2)),x(6)*ones(szy(1),szx(2)),x(7)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6','p7');
g=f(xexp,yexp,p1,p2,p3,p4,p5,p6,p7);
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1 & k6>=1;
[p1,p2,p3,p4,p5,p6]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)),x(3)*ones(szy(1),szx(2)),...
x(4)*ones(szy(1),szx(2)),x(5)*ones(szy(1),szx(2)),x(6)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2','p3','p4','p5','p6');
g=f(xexp,yexp,p1,p2,p3,p4,p5,p6);
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1 & k5>=1;
[p1,p2,p3,p4,p5]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)),x(3)*ones(szy(1),szx(2)),...
x(4)*ones(szy(1),szx(2)),x(5)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2','p3','p4','p5');
g=f(xexp,yexp,p1,p2,p3,p4,p5);
elseif k1>=1 & k2>=1 & k3>=1 & k4>=1;
[p1,p2,p3,p4]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)),x(3)*ones(szy(1),szx(2)),...
x(4)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2','p3','p4');
g=f(xexp,yexp,p1,p2,p3,p4);
elseif k1>=1 & k2>=1 & k3>=1;
[p1,p2,p3]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)),x(3)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2','p3');
g=f(xexp,yexp,p1,p2,p3);
elseif k1>=1 & k2>=1;
[p1,p2]=deal(x(1)*ones(szy(1),szx(2)),x(2)*ones(szy(1),szx(2)));
f=inline(fstring,'x','y','p1','p2');
g=f(xexp,yexp,p1,p2);
elseif k1>=1;
p1=x(1)*ones(szy(1),szx(2));
f=inline(fstring,'x','y','p1');
g=f(xexp,yexp,p1);
end;
z=double(g);
err=sum(sum((z-zexp).^2));