(71B) Shammas Polynomials for curve fitting (#1)

[Namir]

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

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. To give you a heads up, the HP-71B BASIC program presented below runs best on the Windows emulator for the HP-71B.

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 BASIC code listing below performs the calculations for the 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))/SQA(X(I)). You need to edit the assignment statement in line 440 for Y(I) to reflect the function you want to fit. Also edit line 430 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 430 with:

430 READ X(I),Y(I)

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

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 also the best values for a and b. 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 for parameter a:

A = .735050947921

2.2) Press [f][+] to resume the program. The program displays the value for parameter b:

B = 7.28164004367E-2

2.3) Press [f][+] to resume the program. The program displays the value of the adjusted coefficient of determination:

R2ADJ= .999863251812

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

C0=-.576521994341

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

C1= 1.94936466499

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

*X^.807867348358

2.7) Repeat steps 2.5 and 2.6 for the other terms to obtain the following LCD output:

C2=-.288735431868
*X^1.54291829628

C3= 1.89382407804E-2
*X^2.2779692442

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

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

R2ADJ= .999408874877

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

C0= .143492804145

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

C1= 1.14843234603

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

*X^1

2.13) Repeat steps 2.11 and 2.12 for the other terms to obtain the following LCD output:

C2=-9.94877382567E-2
^X^2

C3= 3.77333142368E-3
*X^3

2.14) The program stops after displaying the last output.

The adjusted coefficient of determination is slightly higher than that of the regular polynomial. That is good news!

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