Shammas Polynomial curve fitting (#1)

[Namir]

This is article #1 in a series of similar articles that show how to work with quasi-polynomials.

The Shammas Polynomial curve fitting for the HP Prime.
===================================

Note: I have often posted in the "Not remotely HP Calculators" comments section about using various types of quasi-polynomials with non-integer powers. Recently I introduced the term "polyterms" to define these polynomials more formally. I use MATLAB code in my web site when I present results for using these polyterms. Well, I am glad to share here an HP-Prime program for one of these polyterms.

A Shammas Polynomial has a general form of:

Pn(x) = a0 + a1*x^p(1) + a2*x^p(2) + ... + an*x^p(n)

Where p(i) = a*f(i)+b, with a and be as parameters and f(i) is a function of term i. Each term uses the same values for parameters a and b. I am using f(i)= i in the accompanying program. So, p(i) = a*i+b for the program below. My aim is the have values of a and b such that p(i) < i. Choosing values for parameters a and b in an arbitrary way and expecting excellent results is a pipe dream! Thus, fitting data with the above equation requires optimization to find the best values of a and b.

The code below has three functions:

1) makeData creates and returns a two-column matrix for x and y. The parameters xMin, xMax, and xStep define the range and steps used to calculate the values for x. The function also calculates the corresponding values of y. The function is currently coded to calculate values of y using the expression LN(x)*SQRT(x). Here is a sample call which generates values of x and y for x in the range of (1, 10) in steps of 0.1:

M1:=makeData(1,10,0.1)

The global matrix M1 contains two columns of x and y data, ready to be used by function SHAMPOLY. Here is a partial listing of matrix M1:

Code:

M1|    1    |   2     |
--+---------+---------+
1 |        1|        0|
2 |      1.1|9.0875E-2|
3 |      1.2|0.1664360|
  |        ...        |
91|       10|0.7281413|
--+---------+---------+

You can bypass using the function makeData and enter the x and y data manually into one of the global matrices. Remember to store the x values in column 1 and the y values in column 2.

Of course you need to replace the expression LN(x)/SQRT(x), in function makeData, with the expressions for other functions you want to fit with polynomials.

2) PolyReg is a function that performs the regression for Shammas Polynomials and for regular polynomials. It is called by function SHAMPOLY and returns a list containing the adjusted R-square value and the array of regression coefficients.

3) SHAMPOLY performs the optimization and polynomial regression to return a matrix of results. The parameter DSMat specifies the two-column matrix that contains x and y data. The parameter order specifies the order for the Shammas Polynomial (and comparable regular polynomial). The parameter maxIters specifies the maximum number of iterations used for random search optimization. This method is very simple to code and is relatively effective because the search ranges are small. Here is a sample call that fits the data in matrix M1 using 5th order polynomials, using 1000 iterations for the random search:

M2:=SHAMPOLY(M1,5,1000)

The resulting matrix M2 has four columns containing the following information:

1) The first row has two pairs of constant terms and Adjusted R-square values. The ones in columns 1 and 2 belong to the Shammas Polynomial. The values in columns 3 and 4 belong to the regular polynomial.
2) Rows 2 to 6 contain pairs of regression coefficients and powers for the corresponding terms. The ones in columns 1 and 2 belong to the Shammas Polynomial. The values in columns 3 and 4 belong to the regular polynomial.

Here is a partial listing of the contents of matrix M2 (your results will be different every time you run the program since it uses random numbers):

Code:

M2|  1      |   2     |   3     |    4    |
--+---------+---------+---------+---------+
1 |-2.352815|0.9994502|-0.934081|0.9963811|
2 |4.1410055|0.7296846|1.3646969|        1|
  |                ...                    |
6 |5.4076E-3|3.5439532|-1.9313E-4|       5|
--+---------+---------+---------+---------+

So the adjusted R-square for the Shammas Polynomial is 0.9994502. The adjusted R-square for the regular polynomial is 0.9963811. IF either or both adjusted R-square values are rounded to 1, you can calculate 1-M2(1,2) and 1-M2(1,4) to calculate the values of 1 - the adjusted R-square for each polynomial, and then compare the results. The smaller of the two values indicates which polynomial offers a slightly better fit.

The constant term for the Shammas Polynomial is -2.352815. The constant term for the regular polynomial is -0.934081. Again, the results in the above table will vary each time you run the program.

I am using the adjusted R-square values so that you can compare results obtained from using different polynomials orders.

I used the HP-Prime emulator and was very pleased with the speed of calculations.

So to recap, to test the program perform the following:

1) Enter M1:=makeData(1,10,0.1) and click Enter.
2) Enter M2:=SHAMPOLY(M1,5,1000) and click Enter.
3) Select the Matrix option (use [Shift][4] keys) from the keyboard and inspect matrix M2.

Here is the listing. You can cut-and-paste it into the HP-Prime emulator.

Code:

EXPORT makeData(xMin,xMax,xStep)
BEGIN
  LOCAL i,x; 
  LOCAL lstDimX,NumRows; 
  LOCAL MatX; 
   
  x:=xMin;
  NumRows:=0;
  WHILE x <= xMax DO
    x:=x+xStep;
    NumRows:=NumRows+1;
  END;
  // initialize matrix and vector 
  MatX:=MAKEMAT(0,NumRows,2); 
  x:=xMin;
  FOR i FROM 1 to NumRows DO
    MatX(i,1):=x;
    MatX(i,2):=LN(x)/SQRT(x);
    x:=xMin+xStep
  END;
  RETURN MatX;
END;

EXPORT SHAMPOLY(DSMat,order,maxIters)
BEGIN
  LOCAL i,iter,m; 
  LOCAL L1,Rsqr,RegCoeff; 
  LOCAL bestR2, bestA, bestB, bestCoeff;
  LOCAL resMat,a;
  
  m:=order+1;
  a:=MAKEMAT(0,2,1);
  resMat:=MAKEMAT(0,m,4);
  bestR2:=0;
  
  FOR iter FROM 1 TO maxIters DO 
    a(1,1):=RANDOM(0.7,1.3);
    a(2,1):=RANDOM(0,0.3);
    L1:=PolyReg(DSMat,order,a);
    Rsqr:=L1(1);
    IF Rsqr > bestR2 THEN
      bestR2:=Rsqr;
      bestCoeff:=L1(2);
      bestA := a(1,1);
      bestB := a(2,1);
    END;
  END;
  a(1,1):=bestA;
  a(2,1):=bestB;
  Rsqr:=bestRsr;
  RegCoeff:=bestCoeff;
  resMat(1,1):=RegCoeff(1,1);
  resMat(1,2):=Rsqr;
  FOR i FROM 2 TO m DO
    resMat(i,1):=RegCoeff(i,1);
    resMat(i,2):=a(1,1)*(i-1)+a(2,1);
  END;
  
  // now perform calculations for regular polynomial
  a(1,1):=1;
  a(2,1):=0;
  L1:=PolyReg(DSMat,order,a);
  Rsqr:=L1(1);
  RegCoeff:=L1(2);
  resMat(1,3):=RegCoeff(1,1);
  resMat(1,4):=Rsqr;
  FOR i FROM 2 TO m DO
    resMat(i,3):=RegCoeff(i,1);
    resMat(i,4):=i-1;
  END;
  RETURN resMat;
END;

EXPORT PolyReg(DSMat,order,a) 
BEGIN 
  LOCAL i,j,x,y; 
  LOCAL lstDimX,NumRows; 
  LOCAL MatX, VectY, RegCoeff; 
  LOCAL YMean,Sum1,Sum2,Yhat,Rsqr; 
   
  lstDimX:=SIZE(DSMat); 
  NumRows:=lstDimX(1); 
  // initialize matrix and vector 
  MatX:=MAKEMAT(1,NumRows,order+1); 
  VectY:=MAKEMAT(1,NumRows,1); 
  Sum1:=0; 
  // populate matrix X and vector y with data 
  FOR i FROM 1 TO NumRows DO 
    y:=DSMat(i,2); 
    x:=DSMat(i,1); 
    VectY(i,1):=y; 
    Sum1:=Sum1+y; 
    FOR j FROM 1 TO order DO 
      MatX(i,j+1):=x^(a(1,1)*j+a(2,1)); 
    END;   
  END;  
  // calculate ymean 
  YMean:=Sum1/NumRows; 
  // calculate regression coefficients 
  RegCoeff:=LSQ(MatX,VectY); 
   
  // calculate the sums of squares of  
  // (yhat - ymean) and (y - ymean) 
  Sum1:=0; 
  Sum2:=0; 
  FOR i FROM 1 TO NumRows DO 
    y:=DSMat(i,2); 
    x:=DSMat(i,1); 
    Yhat:=RegCoeff(1,1); 
    FOR j FROM 1 TO order DO 
     Yhat:=Yhat+RegCoeff(j+1,1)*x^(a(1,1)*j+a(2,1)); 
    END; 
    Sum1:=Sum1+(Yhat-YMean)^2; 
    Sum2:=Sum2+(y-YMean)^2; 
  END; 
  // calculate coefficient of determination 
  Rsqr:=Sum1/Sum2;   
  Rsqr:=1-(1-Rsqr)*(NumRows-1)/(NumRows-order-1);
  RETURN {Rsqr,RegCoeff}; 
END;