(71B) Quantum Shammas Polynomials for curve fitting (#2)

[Namir]

The Quantum Shammas Polynomial Curve Fitting for the HP-71B (#2)

Note: This presentation is part of a small series of related topics about curve fitting with quasi-polynomials with the HP-71B calculator.

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-71B program for one of these polyterms.

A Quantum 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) = i-1+r(i), with r(i) is a random number in the range of (0.6, 1.3). This range is not set in stone and can be fine-tuned by the user. Each r(i) is different. Choosing values for each r(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 of for each r(i).

The BASIC code listing below performs the calculations for the Quantum Shammas Polynomials. I coded the program so that it uses labels instead of line numbers in branching. The program uses a MATH ROM. The program works much faster with the Windows version of the HP-71B simulator (for which the MATH ROM is easy found and installed) with the Authentic Calculator Speed option checked off (found in the File => Settings menu).

Notice the following program segments:

1) Line 30 has a READ statement that reads the number of data points (N), the order of the polynomials (M), and the maximum number of iterations for the random search (K9). Line 40 has a DATA statement that supplies the values for the three inputs. You can edit the values to suite your own case that you are studying, or alter the operational parameters for the example in the listing.

2) The label 'GETDATA' calculates the values for arrays X and Y. The current code is set to calculate each Y(I) as LOG(X(I))/SQR(X(I)). You need to edit the assignment statement in line 370 for Y(I) to reflect the function you want to fit. Also edit line 360 to set the sampling of array X, Alternatively, you can replace the statements after label 'GETDATA' with READ and DATA statements to supply values for arrays X and Y. Replace line 370 with:

370 READ X(I),Y(I)

And delete line 380. You need to insert an adequate number of DATA statements, say after line 650.

3) The label 'CALC' calculates the regression coefficients and adjusted coefficient of determination for the Shammas Polynomial and the regular polynomial. You need not change the code in this segment.

The BASIC code shown below is ready to run as a demo program. Remember that since it uses random numbers the results will vary each time you run the program. Once you have performed any edits, you run the program by pressing the [RUN] key. The program performs the following:

1) Displays intermediate results which shows the iteration number and current adjusted coefficient of determination and a message that a better point was found. The program uses the WAIT command to auto-resume excution.

2) The final results start to appear when the program displays "SHAMMAS POLYNOMIAL". The sequence of output is:

2.1) The program displays the value of the adjusted coefficient of determination:

R2ADJ= .974727628685

2.2) Press [f][+] to resume the program. The program displays the constant term:

C0=-.608423568928

2.3) Press [f][+] to resume the program. The program displays the coefficient for the first term:

C1= .856937996154

2.4) Press [f][+] to resume the program. The program displays the power for the first term:

*X^.60894732893

2.5) Repeat steps 2.3 and 2.4 for the other terms to obtain the following LCD output:

C2=-8.35985156185E-2
*X^1.65694731315
C3= 3.87567991994E-3
*X^2.62744974863

2.6) Press [f][+] to resume the program. The program displays "REGULAR POLYNOMIAL" before displaying the results for the regular polynomial regression.

2.7) The program displays the value of the adjusted coefficient of determination:

R2ADJ= .950314751473

2.8) Press [f][+] to resume the program. The program displays the constant term:

C0=-9.10950919331E-2

2.9 Press [f][+] to resume the program. The program displays the coefficient for the first term:

C1= .363826983255

2.10) Press [f][+] to resume the program. The program displays the power for the first term:

*X^1

2.11) Repeat steps 2.09 and 2.10 for the other terms to obtain the following LCD output:

C2=-5.05940470771E-2
*X^2
C3= 2.23138815699E-3
*X^3

2.12) The program stops after displaying the last output.

The adjusted coefficient of determination of the Quantum Shammas Polynomial is higher than that of the regular polynomial.

Keep in mind since the program uses random search optimization, the results will be different each time you run the program.

Here is the program listing.

Code:

10 REM QUANTUM SHAMMAS POLYNOMIALS VER 1
20 DESTROY ALL
30 READ N,M,K9
40 DATA 50,3,200
50 DIM X(N),Y(N),M1(N,M+1),M2(M+1,M+1),M3(N),C(M+1)
60 DIM A(M),A8(M)
70 GOSUB 'GETDATA'
80 R8=0
90 FOR K=1 TO K9
100 FOR I=1 TO M @ A(I)=0.6+0.7*RND @ NEXT I
110 GOSUB 'CALC'
120 IF R8>=R2 THEN 'END1'
130 MAT A8=A
140 R8=R2
150 DISP 'FOUND BETTER POINT" @ WAIT 1
160 'END1':
170 NEXT K
180 BEEP
190 MAT A=A8
200 FOR K=1 TO 2
210 IF K=1 THEN DISP "QUANTUM SHAMMAS POLYNOMIAL" @ WAIT 1
220 IF K=2 THEN  DISP "REGULAR POLYNOMIAL" @ WAIT 1
230 GOSUB 'CALC'
240 DISP "R2ADJ=";R2 @ PAUSE
250 DISP "C0=";C(1) @ PAUSE
260 FOR I=2 TO M+1
270 DISP "C";STR$(I-1);"=";C(I) @ PAUSE
280 DISP "*X^";STR$(I-2+A(I-1)) @ PAUSE
290 NEXT I
300 REM SET A AND B FOR REGULAR POLYNOMIALS
310 FOR I=1 TO M @ A(I) = 1 @ NEXT I
320 NEXT K
330 END
340 'GETDATA':
350 Y5 = 0
360 FOR I=1 TO N
370 X(I) = 1 + 0.2*I
380 Y(I) = LOG(X(I))/SQR(X(I))
390 Y5 = Y5 + Y(I)
400 NEXT I
410 Y5 = Y5/N
420 RETURN
430 'CALC':
440 FOR I=1 TO N
450 M1(I,1)=1
460 FOR J=2 TO M+1
470 M1(I,J) = X(I)^(J-2+A(J-1))
480 NEXT J
490 NEXT I
500 MAT M2 = TRN(M1)*M1
510 MAT M3 = TRN(M1)*Y
520 MAT C = SYS(M2,M3)
530 S1 = 0 @ S2 = 0
540 FOR I=1 TO N
550 Y2=C(1)
560 FOR J=2 TO M+1
570 Y2=Y2+C(J)*X(I)^(J-2+A(J-1))
580 NEXT J
590 S1=S1+(Y2-Y(I))^2
600 S2=S2+(Y(I)-Y5)^2
610 NEXT I
620 R2=1-S1/S2
630 R2=1-(1-R2)*(N-1)/(N-M-1)
640 DISP K;":";R2 @ WAIT 1
650 RETURN