Multivariable Polynomial Regression

[Han]

Updated: Nov. 10, 2014

Here's a preview of the next release:

Points \( (x_i, y_i, z_i) \) were created in the Stat2Var app in columns C1, C2, and C3. The black lines show the distance from the regression surface to the points above the surface, and the red lines from the surface to the points below it. I've included a beta release of the program below.

How to use the programs

multireg( { list1, list2, list3 }, degree );
plotmreg( { list1, list2, list3 }, regression );

where list1, list2, and list3 are each lists (e.g. C1, C2, and C3 in the Stat2Var app), degree is the degree of the 2-variable polynomial, and regression is result from multireg().


In the Stat2Var app, create your columns of data. Let's say your data are in C1, C2, and C3. Run the multireg() program using:

L1:=multireg({C1,C2,C3},deg);

where deg is the degree of the polynomial \( f(X,Y) \). Then run the plotmreg() program. Currently, you will also need manually adjust the viewing window (top of the source code):

plotmreg({C1,C2,C3},L1);

Code:

// viewing window

grid:=14;
xmin:=0;
xmax:=11;
ymin:=0;
ymax:=7;
zmin:=-6;
zmax:=14;

// do not edit below this line

mrerr:={
"Invalid input",
"Insufficient data",
"Singular data",
"Invalid coefficient vector"
};

export ptdef1,ptdef2,ldef1,tridef;

xc:=(xmax+xmin)/2;
yc:=(ymax+ymin)/2;
zc:=(zmax+zmin)/2;
xwin:=-160; ywin:=-120;


// min/max values
mz1:=0; Mz1:=0;
mz2:=0; Mz2:=0;
mz:=0; Mz:=0;

mx:=0; Mx:=0;
my:=0; My:=0;
minD:=10;
dx:=0; dy:=0;

// rotations
rx:=0; ry:=0; rz:=0;
zoom:=10;
RotM;
zclip;
mrparam3d;

mr_rotate();
mr_3dparams();
mr_draw();
mr_keys();

export multireg(m,d)
begin
  local i,n;
  local s1,s2,s3;
  local xi,yi,zi;
  local l1:={},l2:={};

  xi:=m(1); yi:=m(2); zi:=m(3);
  s1:=size(xi); s2:=size(yi); s3:=size(zi);

  if (s1-s2) OR (s2-s3) OR (s1-s3) then
    msgbox(mrerr(1));
    return([[0]]);
  end;

  if (s1 < (d+1)*(d+2)/2) then
    msgbox(mrerr(2));
    return([[0]]);
  end;

  local tmat:=makelist(1,X,1,s1);
  local l:=tmat;

  for i from 1 to d do
    for n from 0 to i do
      if n then l2:=yi^n; else l2:=l; end;
      if (i-n) then l1:=xi^(i-n); else l1:=l; end;
      tmat:=concat(tmat,l1*l2);
    end; 
  end; 

  local tvm:=list2mat(tmat,s1);
  local vm:=transpose(tvm);

  iferr
    vm:=(tvm*vm)^(-1);
  then
    msgbox(mrerr(3));
    return([[0]]);
  else
    return(mat2list(vm*tvm*list2mat(zi,1)));
  end;

end;


export MULTIREG(m,d)
begin
  multireg(m,d);
end;


export plotmreg(m,v)
begin

  local n:=size(v);
  local d:=ceiling((2*n-1/4)^.5-3/2);
  local fxy:="", oldfxy:="";
  local i,j,z;
  local xi:=m(1), yi:=m(2), zi:=m(3);
  local s:=size(xi);
  local g1:=grid+1;

  if ((d+1)*(d+2)/2==n) then

    mx:=min(xi); Mx:=max(xi);
    my:=min(yi); My:=max(yi);
    mz1:=min(zi); Mz1:=max(zi);
    dx:=(xmax-xmin)/grid; dy:=(ymax-ymin)/grid;

    n:=1; fxy:=fxy+v(1);

    for i from 1 to d do
      for j from 0 to i do
        n:=n+1;
        fxy:=fxy+"+"+v(n);
        if (i-j) then fxy:=fxy+"*X^"+(i-j); end;
        if j then fxy:=fxy+"*Y^"+j; end;
      end;
    end;


    iferr oldfxy:=STRING(V0); then end;
    V0:=fxy;
    ptdef1:=concat(
      makelist([xi(X)-xc,yi(X)-yc,V0(xi(X),yi(X))-zc],X,1,s),
      makelist([xi(X)-xc,yi(X)-yc,zi(X)-zc],X,1,s)
      );
    ptdef2:=makelist(
      [
        xmin+((X-1) MOD g1)*dx-xc,
        ymin+IP((X-1)/g1)*dy-yc,
        V0(xmin+((X-1) MOD g1)*dx,ymin+IP((X-1)/g1)*dy)
-zc
      ],
      X,1,g1^2
    );
    V0:=oldfxy;

    ldef1:=makelist({X,X+s,IFTE(ptdef1(X,3)<ptdef1(X+s,3),#0h,#FF0000h)},X,1,s);
//    ldef1:=concat({0},ldef1);

    tridef:={};
    for i from 1 to grid do
      tridef:=concat(tridef,
        makelist({X,X+1,X+g1,#FFh,128},X,g1*(i-1)+1,g1*i-1)
      );
      tridef:=concat(tridef,
        makelist({X,X+1,X-grid,#FFh,128},X,g1*i+1,g1*(i+1)-1)
      );
    end; 

    for i from 1 to s do
      z:=ptdef1(i,3);
      if (z<mz2) then mz2:=z; end;
      if (z>Mz2) then Mz2:=z; end;
    end;

    mz:=min(mz1,mz2); Mz:=max(Mz1,Mz2);
    minD:=10+(max((zc-mz)^2,(mz-zc)^2)+(xmax-xc)^2+(ymax-yc)^2)^(.5);

    mr_keys();

    freeze;

  else
    msgbox(mrerr(4));
  end;

end;

mr_3dparams()
begin
  mr_rotate();

  mrparam3d:={
    RotM, "N", {xwin,ywin,minD*zoom},
    {xmin-xc,xmax-xc,ymin-yc,ymax-yc,zmin-zc,zmax-zc}
  };
end;


mr_draw()
begin

  mr_3dparams();
  dimgrob_p(G1,320,240);
  zclip:=triangle();
  local t:=triangle_p(G1,ptdef2,tridef,mrparam3d,
zclip);
  t:=line_p(G1,ptdef1,ldef1,mrparam3d);
  blit_p(G0,G1);

end;


mr_rotate()
begin
  
  local a=rx+90,b=ry,c=rz;

  if (AAngle==1) OR (NOT (AAngle+HAngle)) then
    a:=a*PI/180; b:=b*PI/180; c:=c*PI/180;
  end;

  local Rx:=[[1,0,0],[0,COS(a),-SIN(a)],[0,SIN(a),COS(a)]];
  local Ry:=[[COS(b),0,-SIN(b)],[0,1,0],[SIN(b),0,COS(b)]];
  local Rz:=[[COS(c),-SIN(c),0],[SIN(c),COS(c),0],[0,0,1]];
  RotM:=Rx*Ry*Rz;
  RotM(3,4):=minD; 
end;

mr_keys()
begin
  local key, run=1;

  mr_draw();

  while run do

    key:=wait(-1);

    case

    if key==-1 then run:=0; end;

    // default
    repeat
      case
        if key==2 then rx:=(rx-5) MOD 360; end;
        if key==12 then rx:=(rx+5) MOD 360; end;
        if key==4 then kill; end;
        if key==7 then rz:=(rz-5) MOD 360; end;
        if key==8 then rz:=(rz+5) MOD 360; end;
        if key==1 then ry:=(ry+5) MOD 360; end;
        if key==3 then ry:=(ry-5) MOD 360; end;
        if key==37 then xwin:=xwin+5; end;
        if key==39 then xwin:=xwin-5; end;
        if key==33 then ywin:=ywin+5; end;
        if key==43 then ywin:=ywin-5; end;
        if key==38 then xwin:=-160; ywin:=-120; end;
        if key==45 then zoom:=max(zoom/1.05,.01); end;
        if key==50 then zoom:=zoom*1.05; end;
      end;
      mr_draw();
    until NOT iskeydown(key);

    end;

  end; // end while
  freeze;

end;