Multivariable Polynomial Regression

[Han]

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)); kill; end;
  if type(s)<>6 then msgbox(pferror(1)); kill; end;
  if s(1)<>2 then msgbox(pferror(2)); kill; end;
  if s(2)<d+1 then msgbox(pferror(3)); kill; end;

  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+1 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;