Quantum Shammas Polynomial curve fitting (#2)

[Namir]

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

The Quantum 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 website 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 Quantum Shammas Polynomial has a general form of:

Qn(x) = c0 + c1*x^a(1) + c2*x^a(2) + ... + cn*x^a(n)

Where a(i) = i-1+r(i), and r(i) is a random value in the range around 1, like (0.6, 1.3) for all values of i. Each value of a(i) is distinct. Choosing values for parameters a(i) 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 for each a(i).

The code below has three functions:

1) QmakeData 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:=QmakeData(1,10,0.1)

The global matrix M1 contains two columns of x and y data, ready to be used by function QUANTPOLY. 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 QmakeData 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 QmakeData, with the expressions for other functions you want to fit with polynomials.

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

3) QUANTPOLY 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 Quantum 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:

M3:=QUANTPOLY(M1,5,1000)

The resulting matrix M3 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 Quantum 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 Quantum Shammas Polynomial. The values in columns 3 and 4 belong to the regular polynomial.

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

Code:

M3|    1    |   2     |   3     |    4    |
--+---------+---------+---------+---------+
1 |-1.768553|0.9984757|-0.934081|0.9963811|
2 |2.2946718|1.6586765|1.3646969|        1|
  |                ...                    |
6 |3.2090E-4|5.6005136|1.9313E-4|        5| 
--+---------+---------+---------+---------+

So the adjusted R-square for the Quantum Shammas Polynomial is 0.9984757. 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-M3(1,2) and 1-M3(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 -1.768553. 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:=QmakeData(1,10,0.1) and click Enter.
2) Enter M3:=QUANTPOLY(M1,5,1000) and click Enter.
3) Select the Matrix option (use [Shift][4] keys) from the keyboard and inspect matrix M3.

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

Code:

EXPORT QmakeData(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:=x+xStep;
  END;
  RETURN MatX;
END;

EXPORT QUANTPOLY(DSMat,order,maxIters)
BEGIN
  LOCAL i,iter,m; 
  LOCAL L1,Rsqr,RegCoeff; 
  LOCAL bestR2, bestA, bestB;
  LOCAL resMat,a,bestA, bestCoeff;
  
  m:=order+1;
  a:=MAKEMAT(0,order,1);
  bestA:=MAKEMAT(0,order,1);
  resMat:=MAKEMAT(0,m,4);
  bestR2:=0;
  
  FOR iter FROM 1 TO maxIters DO 
    FOR i FROM 1 TO order DO
      a(i,1) := RANDOM(0.6,1.3);
    END;
    L1:=QPolyReg(DSMat,order,a);
    Rsqr:=L1(1);
    IF Rsqr > bestR2 THEN
      bestR2:=Rsqr;
      bestCoeff:=L1(2);
      bestA:=a;
    END;
  END;
  a:=bestA;
  RegCoeff:=bestCoeff;
  Rsqr:=bestR2;
  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):=i-1+a(i-1,1);
  END;
  
  // now perform calculations for regular polynomial
  FOR i FROM 1 TO order DO
    a(i,1) := 1;
  END;
  L1:=QPolyReg(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 QPolyReg(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^(j-1+a(j,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^(j-1+a(j,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;