The Symmetric Shammas Polynomial curve fitting (#3)
01-07-2024, 02:02 PM
Post: #1
 Namir Senior Member Posts: 1,012 Joined: Dec 2013
The Symmetric Shammas Polynomial curve fitting (#3)
This is article #3 in a series of similar articles that show how to work with quasi-polynomials.

The Symmetric 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 Symmetric Shammas Polynomial is an even-ordered polyterm (half with positive power terms and half with negative ones) and has a general form of:

Pn(x) = bn/x^p(n) + ... + b1/x^p(1) + a0 + ... + + a2*x^p(1) + 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 four functions:

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

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

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

3) SRegPoly performs regular polynomial least-square calcualtions.

4) SYMSHAMPOLY 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:

M4:=SYMSHAMPOLY(M1,4,1000)

The resulting matrix M4 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 M4 (your results will be different every time you run the program since it uses random numbers):

Code:
M4|  1      |   2     |   3     |    4    | --+---------+---------+---------+---------+ 1 |0.9170903|0.9999994|-0.592398|0.9868934| 2 |0.1760767|-1.302522|0.8642195|        1|   |                ...                    | 5 |3.2186E-4|1.3025222|-8.336E-4|        4| --+---------+---------+---------+---------+

So the adjusted R-square for the Shammas Polynomial is 0.9999994. The adjusted R-square for the regular polynomial is 0.9868934. IF either or both adjusted R-square values are rounded to 1, you can calculate 1-M4(1,2) and 1-M4(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 0.9170903. The constant term for the regular polynomial is -0.592398. 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:=SymMakeData(1,10,0.1) and click Enter.
2) Enter M4:=SYMSHAMPOLY(M1,4,1000) and click Enter.
3) Select the Matrix option (use [Shift][4] keys) from the keyboard and inspect matrix M4.

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

Code:
EXPORT SymMakeData(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 SYMSHAMPOLY(DSMat,order,maxIters) BEGIN   LOCAL i,k,iter,m,order2;    LOCAL L1,Rsqr,RegCoeff;    LOCAL bestR2, bestA, bestB, bestCoeff;   LOCAL resMat,a;      IF 2*IP(order/2) < order THEN     order:=order+1;   END;      order2:=IP(order/2);   m:=order+1;   a:=MAKEMAT(0,1,2);   resMat:=MAKEMAT(0,m,4);   bestR2:=0;      FOR iter FROM 1 TO maxIters DO      a(1,1):=RANDOM(0.7,1.3);     a(1,2):=RANDOM(0,0.3);     L1:=SymPolyReg(DSMat,order,a);     Rsqr:=L1(1);     IF Rsqr > bestR2 THEN       bestR2:=Rsqr;       bestCoeff:=L1(2);       bestA := a(1,1);       bestB := a(1,2);     END;   END;   a(1,1):=bestA;   a(1,2):=bestB;   Rsqr:=bestR2;   RegCoeff:=bestCoeff;   resMat(1,1):=RegCoeff(1,1);   resMat(1,2):=Rsqr;   k:=2;   FOR i FROM order2 DOWNTO 1 DO    resMat(k,1):=RegCoeff(i+1,1);    resMat(k,2):=-((1,1)*(i-1)+a(1,2));    k:=k+1;   END;   FOR i FROM 1 TO order2 DO     resMat(k,1):=RegCoeff(i+order2+1,1);     resMat(k,2):=a(1,1)*(i-1)+a(1,2);     k:=k+1;   END;      // now perform calculations for regular polynomial   L1:=SPolyReg(DSMat,order);   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 SymPolyReg(DSMat,order,a)  BEGIN    LOCAL i,j,k,x,y,order2;    LOCAL lstDimX,NumRows;    LOCAL MatX, VectY, RegCoeff;    LOCAL YMean,Sum1,Sum2,Yhat,Rsqr;        order2:=IP(order/2);   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;      k:=2;     FOR j FROM order2 DOWNTO 1 DO       MatX(i,k):=1/x^(a(1,1)*j+a(1,2));        k:=k+1;     END;     FOR j FROM 1 TO order2 DO        MatX(i,k):=x^(a(1,1)*j+a(1,2));        k:=k+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);      k:=2;     FOR j FROM order2 DOWNTO 1 DO       Yhat:=Yhat+RegCoeff(k,1)/x^(a(1,1)*j+a(1,2));       k:=k+1;     END;          FOR j FROM 1 TO order2 DO       Yhat:=Yhat+RegCoeff(k,1)*x^(a(1,1)*j+a(1,2));       k:=k+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;  EXPORT SPolyReg(DSMat,order)  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;      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;      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;
 « Next Oldest | Next Newest »

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