Multivariable Polynomial Regression
10-09-2014, 08:27 PM (This post was last modified: 03-10-2017 09:22 PM by Han.)
Post: #1
 Han Senior Member Posts: 1,843 Joined: Dec 2013
Multivariable Polynomial Regression
Please see the second post and onward. This program below is obsoleted by the polynomial_regression() command. However, if you wish to do 3D surface regression, the program in post #2 and onward are relevant to you.

Edit #1: apparently there is also the built-in command: polynomial_regression() which I just now noticed
Edit #2: for a polynomial regression in 3D, see the second post in this thread

Here is a program that will create a polynomial regression based on an input matrix m (of dimension 2xn, for n data points) and d, the degree of the polynomial. The result is a list containing
1. a vector corresponding to the polynomial $$p(x) = a_0 + a_1 x + \dotsm + a_d x^d$$ that best-fits (via least squares) the data in the matrix supplied
2. a string representing the polynomial that can be used to define a function in the Function App

$\mathrm{polyfit}\left(\left[ \begin{array}{cccc} x_1 & x_2 & \dotsm & x_n \\ y_1 & y_2 & \dotsm & y_n \end{array} \right], d\right) \rightarrow \left\{ \left[ \begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_d \end{array} \right] , \text{}a_0*X^0 + a_1*X^1 + a_2*X^2 + \dotsm + a_d*X^d\text{''}\right\}$

If you have data from the the Stat2Var app, you can easily create the data matrix via:
Code:
L1:=concat(C1,C2);
M1:=list2mat(L1,SIZE(C1));
at the Home screen -- assuming your data are in C1 and C2 for $$x_i$$ and $$y_i$$ respectively.

To use the string representation of the polynomial, simply type in the Home screen:
Code:
mylist:=polyfit(M1,d);
F1:=mylist(2);
immediately after using the polyfit command. You can use this in the Stat2Var Symb view by selecting User Defined in the Type input selection box.

Below is the short source code for polynomial fit:

PHP Code:
pferror:={"polyfit(m,d):
m: matrix of points (size 2xn)
d: degree of polynomial regression"
,
"matrix argument has wrong dimensions",
"insufficient data",
"singular"}
;

export polyfit(m,d)
begin
local i
;

local s=size(m);

local p="";

if
type(m)<>4 then msgbox(pferror(1)); killend;
if
type(s)<>6 then msgbox(pferror(1)); killend;
if
s(1)<>2 then msgbox(pferror(2)); killend;
if
s(2)<d+1 then msgbox(pferror(3)); killend;

local row=[1]; row[s(2)]:=1;

local tvm=[[0]]; tvm(d+1,s(2)):=0;
for
i from 1 to s(2) do

row(i):=1;

end;

tvm(1):=row;
for
i from 2 to d+do

row:=row .* m(1);

tvm(i):=row

end;

local y=transpose(m(2));

local vm=transpose(tvm);

iferr
m
:=(tvm*vm)^(-1);

then
msgbox
(pferror(4));
return({[[
0]],""});
else

m:=m*tvm*y;
for
i from 0 to d do
if
i then p:=p+"+"end;

p:=p+m(i+1,1)+"*X^" i;

end;
return({
m,p});

end;

end

Graph 3D | QPI | SolveSys
10-13-2014, 10:28 PM (This post was last modified: 11-10-2014 08:18 PM by Han.)
Post: #2
 Han Senior Member Posts: 1,843 Joined: Dec 2013
RE: Polynomial Regression
Here is a program that is best used in conjunction with the Stat2Var app. Eventually this will evolve into a full-blown app that will enable users to create 3D surface regressions (using polynomials of a degree specified by the user).
Usage as follows:
$\mathrm{multireg}(\text{list_of_data},\text{degree}) \rightarrow \left\{ a_0, a_1, \dotsm, a_{k(d)-1} \right\}$
where $$k(d) = (d+1)(d+2)/2$$ and the 2-variable polynomial is of the form
$f(x,y) = a_0 + a_1 x + a_2 y + a_3 x^2 + a_4 xy + a_5 y^2 + \dotsm + a_{k(d)-2} xy^{d-1} + a_{k(d)-1} y^d$
and the list_of_data parameter is a list of three lists (for the $$x_i$$, $$y_i$$, and $$z_i$$ values). Typically the data would be pulled from C0 through C9 in the Stat2Var app. So list_of_data may be something like: { C1, C2, C3 }.

The algorithm is essentially the same as that for 1-variable polynomial regression. Both cases use an appropriate Vandermonde matrix to solve for the critical point in the least sum-of-squared-residuals equation.

(Please note that error checking for valid input data is very minimal.)

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

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;

Graph 3D | QPI | SolveSys
11-07-2014, 09:43 PM (This post was last modified: 11-10-2014 08:22 PM by Han.)
Post: #3
 Han Senior Member Posts: 1,843 Joined: Dec 2013
RE: Polynomial Regression
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;

Attached File(s) Thumbnail(s)

Graph 3D | QPI | SolveSys
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)