HƒCURVEFITBASFäŸCURVEFITDOCG «÷ÿ10 REM 20 REM POLYNOMIAL CURVEFIT 30 REM 40 REM-------------------MAIN PROGRAM------------------------ 50 GOSUB 110 REM DIMENSION ARRAYS AND INTRODUCTION 60 GOSUB 220 REM INITIALIZE ARRAY VALUES/INPUT DATA 70 GOSUB 520 REM GENERATE VALUES FOR ARRAY 80 GOSUB 1150 REM SUM OF ERRORS SQUARED 90 GOSUB 1270 REM OUTPUT DATA 100 GOSUB 1460 REM MAIN OPTIONS MENU 110 REM-----------------DIMENSION ARRAYS--------------------- 120 DIM A(20,20),B(20),C(20),X(20),Y(20) 130 DIM XX(120),YY(120),Q(45),YZ(20) 140 REM------------------INTRODUCTION----------------------- 150 SCREEN 0,0,0:CLS 160 PRINT " *************** POLYNOMIAL CURVE FIT ****************" 170 PRINT:PRINT 180 PRINT"THIS PROGRAM PERFORMS A LEAST SQUARED CURVE FIT FOR A GENERAL POLYNOMIAL" 190 PRINT"IN X AND Y. IT WILL ALSO EVALUATE AND PLOT THE FUNCTION IF YOU WISH." 200 PRINT 210 RETURN 220 REM 230 REM--------------SET ARRAY VALUES TO ZERO--------------- 240 FOR I%=0 TO 20 250 B(I%)=0 260 Q(I%)=0 270 NEXT I% 280 IF Z% = 2 THEN GOTO 420 290 IF Z% = 5 THEN MYANSWER$ = "Y" 300 IF Z% = 5 THEN GOTO 400 310 REM--------------------INPUT DATA----------------------- 320 INPUT"How many DATA POINTS to be entered";DATAPTS% 330 CLS:SCREEN 1,0,0 340 PRINT"NOW INPUT THE DATA." 350 FOR I%=1 TO DATAPTS% 360 PRINT USING "X(##),Y(##)";I%,I%; 370 INPUT X(I%),Y(I%) 380 NEXT I% 390 INPUT"Any CORRECTIONS (Y/N)";MYANSWER$ 400 IF (MYANSWER$ = "Y") OR (MYANSWER$ = "y") THEN GOTO 2380 410 CLS 420 INPUT"What ORDER fit:";ORDER% 430 N%=ORDER%+1 440 REM 450 REM--------------VALIDITY CHECK------------------------- 460 REM 470 IF ORDER%>=DATAPTS% THEN PRINT"NOT ENOUGH DATA FOR THAT ORDER FIT": GOTO 420 480 CLS 490 RETURN 500 REM--------------MATRIX LOAD---------------------------- 510 REM--------------GENERATE VALUES FOR ARRAY-------------- 520 FOR I%=0 TO ORDER%*2 530 FOR J%=1 TO DATAPTS% 540 Q(I%)=X(J%)^I%+Q(I%) 550 NEXT J% 560 Q(I%)=2*Q(I%) 570 NEXT I% 580 REM LOAD 590 FOR I%=0 TO ORDER% 600 FOR J%=0 TO ORDER% 610 A(I%+1,J%+1)=Q(I%+J%) 620 NEXT J% 630 NEXT I% 640 FOR I%=0 TO ORDER% 650 FOR J%=1 TO DATAPTS% 660 B(I%+1)=2*X(J%)^I%*Y(J%)+B(I%+1) 670 NEXT J% 680 NEXT I% 690 REM------------SIMULTANEOUS EQU SOLN------------------ 700 FLAG%=0 710 IF N%<>0 THEN GOTO 750 720 IF A(1,1)=0 THEN FLAG%=1: RETURN 730 X(1)=B(1)/A(1,1) 740 RETURN 750 M%=N%-1 760 FOR I%=1 TO M% 770 H=ABS(A(I%,I%)) 780 L%=I% 790 I1%=I%+1 800 FOR J%=I1% TO N% 810 IF ABS(A(J%,I%)) < H THEN 840 820 H=ABS(A(J%,I%)) 830 L%=J% 840 NEXT J% 850 IF H=0 THEN FLAG%=1: RETURN 860 IF L%=I% THEN 950 870 FOR J%=1 TO N% 880 G=A(L%,J%) 890 A(L%,J%)=A(I%,J%) 900 A(I%,J%)=G 910 NEXT J% 920 G=B(L%) 930 B(L%)=B(I%) 940 B(I%)=G 950 FOR J%=I1% TO N% 960 T=A(J%,I%)/A(I%,I%) 970 FOR K%=I1% TO N% 980 A(J%,K%)=A(J%,K%)-T*A(I%,K%) 990 NEXT K% 1000 B(J%)=B(J%)-T*B(I%) 1010 NEXT J% 1020 NEXT I% 1030 IF A(N%,N%)=0 THEN FLAG%=1: RETURN 1040 C(N%)=B(N%)/A(N%,N%) 1050 I%=N%-1 1060 S=0 1070 I1%=I%+1 1080 FOR J%=I1% TO N% 1090 S=S+A(I%,J%)*C(J%) 1100 NEXT J% 1110 C(I%)=(B(I%)-S)/A(I%,I%) 1120 I%=I%-1 1130 IF I%>0 THEN 1060 1140 RETURN 1150 REM--------------SUM OF ERRORS SQRD------------------- 1160 FOR I%=1 TO DATAPTS% 1170 YZ(I%)=0 1180 NEXT I% 1190 E=0 1200 FOR I%=1 TO DATAPTS% 1210 FOR J%=1 TO N% 1220 YZ(I%)=C(J%)*X(I%)^(J%-1)+YZ(I%) 1230 NEXT J% 1240 E=E+(Y(I%)-YZ(I%))^2 1250 NEXT I% 1260 RETURN 1270 REM--------------OUTPUT DATA---------------------------- 1280 IF FLAG%=1 THEN PRINT "AMBIGUOUS DATA. NO SOLUTION FOUND":RETURN 1290 REM 1300 PRINT" POWER OF X"," COEFFICIENT":PRINT 1310 FOR I%=N% TO 1 STEP -1 1320 PRINT USING " X^#";I%-1, 1330 PRINT USING " ###.###";C(I%) 1340 NEXT I% 1350 PRINT USING " SUM OF ERRORS SQRD= ###.####";E 1360 IF N%=2 THEN GOTO 2510 1370 IF (N%=2) AND ((DATAPTS%-2)>0) THEN SEE = SQR(E/(DATAPTS%-2)) 1380 IF (N%=2) AND ((DATAPTS%-2)>0) THEN PRINT USING " STD ERR OF EST OF Y ON X= ###.####";SEE 1390 IF N%=2 THEN PRINT USING " COEFF OF DETN = #.#####";COD 1400 PRINT 1410 FOR I%=1 TO DATAPTS% 1420 PRINT USING " X=###.###";X(I%), 1430 PRINT USING " Y=###.###";YZ(I%), 1440 PRINT USING " INPUT Y=###.###";Y(I%) 1450 NEXT I% 1460 REM--------------OPTION MENU------------------------ 1470 REM 1480 INPUT"Press RETURN to continue";A$ 1490 SCREEN 1,0,0 1500 CLS 1510 PRINT"Enter the desired option" 1520 PRINT"1)CALCULATE VALUES OF Y" 1530 PRINT"2)TRY A DIFFERENT ORDER FIT" 1540 PRINT"3)TRY A NEW DATA SET" 1550 PRINT"4)PLOT THE FUNCTION" 1560 PRINT"5)CHANGE DATA PAIR VALUES" 1570 PRINT"6)EXIT THE PROGRAM" 1580 PRINT:INPUT"Enter your choice";Z% 1590 ON Z% GOTO 1630,60,60,1850,60,1600 1600 REM---------------EXIT THE PROGRAM------------------- 1610 SCREEN 0,0,0 1620 CLS:PRINT "EXITING THE CURVE FITTING PROGRAM":END 1630 REM------------CALCULATE FITTED DATA------------- 1640 CLS 1650 SCREEN 0,0,0 1660 INPUT "ENTER NUMBER OF X VALUES:";NUMXS% 1670 FOR I%=1 TO NUMXS% 1680 PRINT USING "ENTER VALUE OF X(##):";I%; 1690 INPUT XX(I%) 1700 NEXT I% 1710 FOR I%=1 TO NUMXS% 1720 YY(I%)=0 1730 FOR J%=1 TO N% 1740 YY(I%)=C(J%)*XX(I%)^(J%-1)+YY(I%) 1750 NEXT J% 1760 NEXT I% 1770 CLS 1780 FOR I%=1 TO NUMXS% 1790 PRINT USING "X = ###.###";XX(I%), 1800 PRINT USING " Y = ###.###";YY(I%) 1810 IF I%=7 OR I%=14 OR I%=21 THEN INPUT "Hit return to continue........";DUMMY$ 1820 NEXT I% 1830 INPUT "Hit return to continue......";DUMMY 1840 GOTO 1510 1850 REM--------------PLOT FUNCTION-------------------- 1851 CLS 1852 PRINT" PLOT MENU" 1853 PRINT:PRINT"Choose one of the following" 1854 PRINT"1 - Plot normalized on X AXIS" 1855 PRINT"2 - Plot normalized on Y AXIS" 1856 PRINT"3 - Plot not normalized" 1857 PRINT:INPUT"Enter your choice: ";PLOT% 1870 XPTMAX = X(1):XMAXI=1 1880 YPTMAX = Y(1):YMAXI=1 1890 FOR I=1 TO DATAPTS%-1 1900 IF XPTMAX <= X(I+1) THEN XPTMAX = X(I+1):XMAXI=I+1 1910 ELSE XPTMAX = XPTMAX 1920 IF YPTMAX <= Y(I+1) THEN YPTMAX = Y(I+1):YMAXI=I+1 1930 ELSE YPTMAX = YPTMAX 1940 NEXT I 1945 ON PLOT% GOTO 1955,1956,1950 1946 GOTO 1851 1950 SFACTOR = 1:SFACTORX = (Y(YMAXI)-C(1))/X(YMAXI):GOTO 1970 1955 SFACTOR = (Y(XMAXI)-C(1))/XPTMAX:SFACTORX = (Y(XMAXI)-C(1))/XPTMAX:GOTO 1970 1956 SFACTOR = (YPTMAX-C(1))/X(YMAXI):SFACTORX = (YPTMAX-C(1))/X(YMAXI):GOTO 1970 1970 CLS 1980 SCREEN 3,0,0 1990 LINE (0,60)-(250,60):LINE (2,0)-(2,63) 2000 SCALE=0:YMAX=0 2010 FOR I%=1 TO N% 2020 DUMMY = (C(I%)*(XPTMAX)^(I%-1)) 2030 YMAX = YMAX + DUMMY 2040 NEXT I% 2050 SCALE = YMAX/50:XINT=YMAX/250 2060 FOR XPLOT=0 TO 250/SFACTOR STEP 1/SFACTOR 2070 YPLOT=0 2080 FOR J%=1 TO N% 2090 DUMMY1=(C(J%)*((XINT*XPLOT)^(J%-1)))/SCALE 2100 YPLOT = YPLOT+DUMMY1 2110 NEXT J% 2120 YPLOT = 60-YPLOT 2130 IF YPLOT < 0 THEN GOTO 2160 2140 PSET (XPLOT*SFACTOR+2,YPLOT),1 2150 NEXT XPLOT 2160 FOR I%=1 TO DATAPTS% 2170 YPLOT1 = 0 2180 XPLOT1 = (X(I%)/XINT) 2190 FOR J%=1 TO N% 2200 DUMMY1=(C(J%)*(X(I%)^(J%-1)))/SCALE 2210 YPLOT1 = YPLOT1+DUMMY1 2220 NEXT J% 2230 YPLOT1X = 60-(YPLOT1*Y(I%)/(YZ(I%)+.001))*SFACTORX 2240 YPLOT1 = 60-YPLOT1*Y(I%)/(YZ(I%)+.001) 2250 IF YPLOT1 < 0 THEN GOTO 2350 2260 IF (SFACTORX<1) AND (SFACTOR=1) THEN GOTO 2290 2270 LINE (XPLOT1*SFACTOR+1,YPLOT1-1)-(XPLOT1*SFACTOR+3,YPLOT1+1),,BF 2280 GOTO 2310 2290 REM SCALING FOR NONNORMALIZED PLOT SPECIAL CASE SLOPE<1 2300 LINE (XPLOT1*SFACTORX+1,YPLOT1X-1)-(XPLOT1*SFACTORX+3,YPLOT1X+1),,BF 2310 NEXT I% 2320 IF PLOT% = 3 THEN LOCATE 55,3,0:PRINT"UNNORMALIZED PLOT OF":GOTO 2340 2330 LOCATE 55,3,0:PRINT" NORMALIZED PLOT OF" 2340 LOCATE 55,4,0:PRINT USING"FUNCTION OF ORDER = ##";N%-1 2350 LOCATE 55,8,0:INPUT"HIT TO CONTINUE";DUMMY$ 2360 SCREEN 1,0,0 2370 GOTO 1500 2380 REM---------DATA CORRECTION ROUTINE------------------ 2390 PRINT "Which DATA PAIR do you want to CHANGE"; 2400 INPUT CHANGES% 2410 SCREEN 0,0,0 2420 PRINT USING " VALUE X(##),Y(##)=###.###,###.###";CHANGES%,CHANGES%,X(CHANGES%),Y(CHANGES%) 2430 PRINT " CHANGE TO (X,Y): "; 2440 INPUT X(CHANGES%),Y(CHANGES%) 2450 CLS:SCREEN 1,0,0 2460 FOR I%=1 TO DATAPTS% 2470 PRINT USING "X(##),Y(##)=###.###,###.###";I%,I%,X(I%),Y(I%) 2480 NEXT I% 2490 INPUT"Any Corrections (Y/N)";MYANSWER$ 2500 GOTO 400 2510 REM--------COEFFICIENT OF DETERMINATION-------------- 2520 SUMXY=0:SUMX=0:SUMY=0:SUMX2=0:SUMY2=0 2530 FOR I%=1 TO DATAPTS% 2540 SUMXY=X(I%)*Y(I%)+SUMXY 2550 SUMX=X(I%)+SUMX 2560 SUMX2=X(I%)^2+SUMX2 2570 SUMY=Y(I%)+SUMY 2580 SUMY2=Y(I%)^2+SUMY2 2590 NEXT I% 2600 COD1=(SUMXY-(SUMX*SUMY/DATAPTS%))^2 2610 COD2=(SUMX2-SUMX^2/DATAPTS%)*(SUMY2-SUMY^2/DATAPTS%) 2620 COD = COD1/COD2 2630 GOTO 1370 This program performs a least squared curve fit for a general polynomial in x and y. It will also evaluate and plot the function on the screen of the PX-8. I have tested it for linear and parabolic fits and the results agree with my HP-65, so I know it works for first and second order fits. This program was taken from the Capital Osborne Users Group Library and "extensively" modified to clean up the formatting, provide functions not originally included, and in general make it directly applicable to the PX-8. I made no attempt to rename variables used by the original programmer, so you will find some parts of it rather obtuse. I tried to make my own additions and corrections reasonably clear. In general most prompts are self-explanatory. Order fits are linear (1), polynomial (2), etc. The plotting routine is scaled for about 40 columns for a slope of 1 and the onscreen portion of the curve will shrink in the horizontal direction as the slope increases. A scaling factor could be included here and perhaps I will do that later. The plotting routine plots the function and then plots the original x,y values that were fitted. After the plot is completed and the ? appears in the left hand corner, hit to get out of the plot routine. If you have any questions or problems, leave me a note (73157,1527). I hope to add some other useful programs to the portable library as I work with my PX-8 over the next few months. This is a THANKS in return for all the software I got from the library, which got me up and running over the last month. Mike Epstein