ANYFIT BASJLINFIT BASMMOSTFIT BASjD >POLYFIT BAS,,CURVFIT DOCJ MTEST1 DAT$CURTEST DAT%"LINTEST2DAT&dMTEST2 DAT'&POLYTESTDAT(LINTEST1DAT)U1000 ' *** PROGRAM ANYFIT *** 1010 ' 1020 ' A DRIVER PROGRAM FOR THE CURFIT SUBROUTINES THAT WILL FIND 1030 ' THE PARAMETERS IN A GIVEN FUNCTION THAT GIVE THE BEST FIT 1040 ' (IN A LEAST SQUARES SENSE) TO A SET OF INPUT DATA. X & Y ARE 1050 ' THE INDEPENDENT AND DEPENDENT VARIABLES RESPECTIVELY. THE 1060 ' FUNCTION THAT DETERMINES Y GIVEN A VALUE OF X MUST BE SPECIFIED 1070 ' IN THE SUBROUTINE FUNCTN (SEE END OF SOURCE CODE). SEE COMMENTS 1080 ' IN CURFIT SUBROUTINE FOR EXPLANATION OF VARIABLES. 1090 ' 1100 ' Programmer: G. M. Resch, Feb. 79 1110 ' 1120 ' DIMENSION ALL ARRAYS HERE 1130 ' THIS VERSION IS SET UP TO HANDLE UP TO 100 DATA POINTS AND 1140 ' AS MANY AS 10 PARAMETERS. 1150 DIM X(100), Y(100), SIGY(100), YFIT(100), WGHT(100) 1160 DIM A(10), DELA(10), IK(10), JK(10), BETA(10), B(10), DERIV(10) 1170 DIM ALPA#(10,10), ARY#(10,10) 1180 ' 1190 MAXN= 100 : MAXP=10 'SET MAXIMUM # OF POINT AND PARMS 1200 ' 1210 ' OPEN & READ THE INPUT DATA SET 1215 PRINT" Program ANYFIT, G. M. Resch, Feb. 79" 1220 INPUT "Name of Input data set=";N$ 1230 OPEN"I",2,N$ 1240 FOR I=1 TO MAXN 1250 IF EOF(2) GOTO 1280 1260 INPUT#2,X(I),Y(I),SIGY(I) 1270 NEXT I 1280 NPTS=I-1 : CLOSE#2 1290 PRINT"THERE ARE ";NPTS;" DATA POINTS" 1300 ' 1310 ' INPUT THE APRIORI VALUES FOR THE PARAMETERS 1320 PRINT" MODE=+1, WGHT(I)=1/SIGY(I)^2" 1330 PRINT" = 0, WGHT(I)=1" 1340 PRINT" =-1, WGHT(I)=1/Y(I)" 1350 PRINT" =-2, WGHT(I)=SIGY(I)" 1360 INPUT"SPECIFY THE WEIGHTING MODE (+1,0,-1, OR -2) -";MODE 1370 IF (MODE<=+1 AND MODE>=-2) GOTO 1390 1380 PRINT"ERROR IN MODE SPECIFICATION" : GOTO 1320 1390 PRINT 1400 INPUT"HOW MANY PARAMETERS? NTRM=";NTRM 1410 IF (NTERM<=10 AND NTRM>0) GOTO 1430 1420 PRINT"ILLEGAL NUMBER OF PARAMETERS" : GOTO 1390 1430 PRINT 1440 FOR I=1 TO NTRM 1450 PRINT"SPECIFY THE STARTING VALUE AND DELTA FOR NTRM=";I 1460 INPUT" NTRM=";A(I) 1470 INPUT" DELTA=";DELA(I) 1480 NEXT I 1490 ' 1500 ' SETUP THE CALL TO CURFIT 1510 FIRST$="TRUE" 'THIS IS THE FIRST CALL 1520 LSTCHI=1E+08 'SET TO REDICULOUSLY HIGH VALUE FOR START 1530 FLMD=.01 'GOOD VALUE FOR START 1540 LOOP=0 'INITIALIZE LOOP COUNT 1550 LOOP=LOOP+1 1560 GOSUB 2040 'CALL CURFIT 1570 PRINT"LOOP #";LOOP,"CHISQ=";CHISQ,"FLAMDA=",FLMD 1580 IF (LSTCHI-CHISQ)<.001 GOTO 1640 'EXIT IF NO CHANGE IN CHISQ 1590 LSTCHI=CHISQ 1600 IF FLMD<.000001 THEN FLMD=.00001 1610 IF FLMD<1E+06 GOTO 1550 'KEEP LOOPING IF FLMD IS NOT TOO LARGE 1620 ' 1630 ' WE FOUND AN ACCEPTABLE FIT, SUMMARIZE THE RESULTS 1640 PRINT 1650 FOR J=1 TO NTRM 1660 PRINT"A(";J;")=";A(J),"SIGMAA(";J;")=";SIGA(J) 1670 NEXT J 1680 PRINT"CHI SQUARE=";CHISQ 1690 PRINT 1700 PRINT"X(I)","Y(I)","YFIT","Y-YFIT" 1710 FOR I=1 TO NPTS 1720 PRINT X(I),Y(I),YFIT(I),(Y(I)-YFIT(I)) 1730 NEXT I 1740 END 1750 REM*************************************************************** 1760 ' SUBROUTINE CURFIT(X, Y, SIGY, NPTS, NTRM, MODE, A, DELA, 1770 ' SIGA, FLMD, YFIT, CHISQ, $FIRST) 1780 ' 1790 ' A ROUTINE THAT WILL FIND THE LEAST-SQUARES FIT TO A NON-LINEAR 1800 ' FUNCTION WITH A LINEARIZATION OF THE FITTING FUNCTION. 1830 ' 1840 ' PARAMETERS: 1850 ' X=ARRAY OF INDEPEND. VARIABLE, Y=ARRAY OF DEPENDENT VARIABLE 1860 ' SIGY=STND DEV. OF Y, NPTS=# OF PAIRS OF DATA POINTS 1870 ' NTRM=# OF PARMS, A=ARRAY OF PARMS 1880 ' MODE=+1, WGHT(I)=1/SIGY(I)^2 1890 ' = 0, WGHT(I)=1 1900 ' =-1, WGHT(I)=1/Y(I) 1910 ' =-2, WGHT(I)=SIGY(I) 1920 ' DELA=ARRAY OF INCREMENTS OF A, SIGAA=STND DEV OF A 1930 ' FLMD=PROPORTION OF GRADIENT SEARCH INCLUDED 1940 ' CHISQ=REDUCED CHI SQUARE 1950 ' FIRST$= A FLAG THAT IS TRUE WHEN CURFIT IS FIRST CALLED 1960 ' 1970 ' SUBROUTINES CALLED; 1980 ' FUNCTN(X,I,A), FCHISQ(Y,SIGY,NPTS,NFRE,MODE,YFIT) 1990 ' FDERIV(X,I,A,DELA,NTRM,DERIV), MATINV(ARY#,NTRM,DET) 2000 ' 2010 '--- START OF PROGRAM CODE --- 2020 ' 2030 ' ROUTINE IS FIRST CALLED 2040 IF FIRST$="FALSE" GOTO 2150 'SKIP IF NOT FIRST CALL 2050 FIRST$="FALSE" 2060 NFRE=NPTS-NTRM : CHISQ=0 : IF NFRE<1 THEN RETURN 2070 ' EVALUATE WEIGHTS 2080 FOR I=1 TO NPTS 2090 IF MODE=1 THEN WGHT(I)=1/(SIGY(I)*SIGY(I)) : GOTO 2130 2100 IF MODE=0 THEN WGHT(I)=1 : GOTO 2130 2110 IF MODE=-1 THEN WGHT(I)=1/ABS(Y(I)) : GOTO 2130 2120 IF MODE=-2 THEN WGHT(I)=SIGY(I) 2130 NEXT I 2140 ' EVAL ALPHA & BETA MATRICES 2150 FOR J=1 TO NTRM 2160 BETA(J)=0 2170 FOR K=1 TO J 2180 ALPA#(J,K)=0 2190 NEXT K,J 2200 ' 2210 FOR I=1 TO NPTS 2220 GOSUB 2850 :REM* CALL FDERIV(X,I,A,DELA,NTRM,DERIV) 2230 GOSUB 3780 :REM* CALL FUNCTN(X,I,A) 2240 FOR J=1 TO NTRM 2250 BETA(J)=BETA(J) +WGHT(I)*(Y(I)-FUNCTN)*DERIV(J) 2260 FOR K=1 TO J 2270 ALPA#(J,K)=ALPA#(J,K) +WGHT(I)*DERIV(J)*DERIV(K) 2280 NEXT K,J 2290 NEXT I 2300 ' 2310 FOR J=1 TO NTRM 2320 FOR K=1 TO J 2330 ALPA#(K,J)=ALPA#(J,K) 2340 NEXT K,J 2350 ' EVAL CHI SQUARE AT START POINT 2360 FOR I=1 TO NPTS 2370 GOSUB 3780 :REM* CALL FUNCTN(X,I,A) 2380 YFIT(I)=FUNCTN 2390 GOSUB 2990 :REM* CALL FCHISQ(Y,SIGY,NPTS,NFRE,MODE,YFIT) 2400 CHI=FCHISQ 2410 NEXT I 2420 ' INVERT MODIFIED CURVATURE MATRIX TO FIND NEW PARMS 2430 FOR J=1 TO NTRM 2440 FOR K=1 TO NTRM 2450 ARY#(J,K)=ALPA#(J,K)/SQR(ALPA#(J,J)*ALPA#(K,K)) 2460 NEXT K 2470 ARY#(J,J)=1+FLMD 2480 NEXT J 2490 ' 2500 NOR=NTRM : GOSUB 3140 :REM* CALL MATINV(ARY#,NTRM,DET) 2510 FOR J=1 TO NTRM 2520 B(J)=A(J) 2530 FOR K=1 TO NTRM 2540 B(J)=B(J) +BETA(K)*ARY#(J,K)/SQR(ALPA#(J,J)*ALPA#(K,K)) 2550 NEXT K,J 2560 ' IF CHI SQUARE INCREASED, INCREASE FLMD AND TRY AGAIN 2570 FOR J=1 TO NTRM 2580 TEMP=A(J) : A(J)=B(J) : B(J)=TEMP 2590 NEXT J 2600 ' 2610 FOR I=1 TO NPTS 2620 GOSUB 3780 :REM* CALL FUNCTN(X,I,A) 2630 YFIT(I)=FUNCTN 2640 NEXT I 2650 FOR J=1 TO NTRM 2660 TEMP=B(J) : B(J)=A(J) : A(J)=TEMP 2670 NEXT J 2680 ' 2690 GOSUB 2990 :REM* CALL FCHISQ(Y,SIGY,NPTS,NFRE,MODE,YFIT) 2700 CHISQ=FCHISQ 2710 IF (CHI-CHISQ)>0 GOTO 2740 2720 FLMD=10*FLMD : GOTO 2430 2730 ' EVAL PARMS & UNCERTAINTIES 2740 FOR J=1 TO NTRM 2750 A(J)=B(J) : SIGA(J)= SQR(ARY#(J,J)/ALPA#(J,J)) 2760 NEXT J 2770 FLMD=FLMD/10 2780 RETURN 2790 REM************************************************************** 2800 ' SUBROUTINE FDERIV(X, I, A, DELA NTRM, DERIV) 2810 ' 2820 ' A ROUTINE TO EVALUATE THE DERIVATIVE OF A FUNCTION WITH 2830 ' RESPECT TO THE PARAMETERS A(J) 2840 ' 2850 FOR JJ=1 TO NTRM 2860 AJ=A(JJ) : DEL=DELA(JJ) : A(JJ)=AJ+DEL 2870 GOSUB 3780 :REM* CALL FUNCTN(X,I,A) 2880 YF=FUNCTN : A(JJ)=AJ-DEL 2890 GOSUB 3780 2900 DERIV(JJ)=(YF-FUNCTN)/(2*DEL) 2910 A(JJ)=AJ 2920 NEXT JJ 2930 RETURN 2940 REM************************************************************** 2950 ' SUBROUTINE FCHISQ(Y, WGHT, NPTS, NFRE, YFIT) 2960 ' 2970 ' A ROUTINE TO CALCULATE THE REDUCED CHI SQUARE FIT TO DATA 2980 ' 2990 CHS=0 : IF NFRE>0 GOTO 3020 3000 FCHISQ=0 : RETURN 3010 ' ACCUMULATE CHI SQUARE 3020 FOR II=1 TO NPTS 3030 CHS=CHS +WGHT(II)*(Y(II)-YFIT(II))*(Y(II)-YFIT(II)) 3040 NEXT II 3050 FCHISQ=CHS/NFRE 3060 RETURN 3070 REM************************************************************** 3080 ' SUBROUTINE MATINV(ARY,NOR,DET) 3090 ' 3100 ' A ROUTINE TO INVERT A SYMMETRIC MATRIX AND CALCULATE ITS 3110 ' DETERMINANT. ARY(I,J) IS THE INPUT MATRIX OF ORDER=NOR 3120 ' AND IS REPLACED BY ITS INVERSE 3130 ' 3140 DET=1! 3150 ' START THE K% LOOP 3160 FOR K%=1 TO NOR 3170 AMX#=0! 3180 FOR I%=K% TO NOR 3190 FOR J%=K% TO NOR 3200 IF ABS(AMX#)>ABS(ARY#(I%,J%)) GOTO 3220 3210 AMX#=ARY#(I%,J%) : IK(K%)=I% : JK(K%)=J% 3220 NEXT J%,I% 3230 IF AMX#<>0! GOTO 3250 3240 DET=0! : RETURN 'RETURN IF DETERMINANT IS ZERO 3250 I%=IK(K%) 3260 IF I%+1 GOTO 1930 1920 WGHT(I)=1/(SIGY(I)*SIGY(I)) : GOTO 2010 1930 IF MODE<>0 GOTO 1950 1940 WGHT(I)=1! : GOTO 2010 1950 IF MODE<>-1 GOTO 1980 1960 IF Y(I)<0 THEN WGHT(I)=1/(-Y(I)) ELSE WGHT(I)=1/Y(I) 1970 GOTO 2010 1980 IF MODE=-2 GOTO 2000 1990 PRINT"MODE SPECIFICATION ERROR" : RETURN 2000 WGHT(I)=SIGY(I) 2010 SUM#=SUM#+WGHT(I) 2020 YMN#=YMN#+WGHT(I)*Y(I) 2030 FOR J=1 TO NTRM 2040 GOSUB 3440 :REM* CALC FCTN 2050 XMN#(J)=XMN#(J)+WGHT(I)*FCTN 2060 NEXT J,I 2070 YMN#=YMN#/SUM# 2080 ' 2090 FOR J=1 TO NTRM : XMN#(J)=XMN#(J)/SUM# : NEXT J 2100 WMN=SUM#/NPTS 2110 FOR I=1 TO NPTS : WGHT(I)=WGHT(I)/WMN : NEXT I 2120 ' *** ACCUMULATE MATRICES R & ARY# 2130 FOR I=1 TO NPTS 2140 SIG#=SIG#+WGHT(I)*(Y(I)-YMN#)*(Y(I)-YMN#) 2150 FOR J=1 TO NTRM 2160 GOSUB 3440 :REM* CALC FCTN(X,I,J,M) 2170 SIGX#(J)=SIGX#(J) +WGHT(I)*(FCTN -XMN#(J))*(FCTN-XMN#(J)) 2180 R(J)=R(J) +WGHT(I)*(FCTN -XMN#(J))*(Y(I)-YMN#) 2190 FJ=FCTN 2200 FOR K=1 TO J 2210 JJ=J : J=K : GOSUB 3440 2220 FK=FCTN : J=JJ 2230 ARY#(J,K)=ARY#(J,K) +WGHT(I)*(FJ-XMN#(J))*(FK-XMN#(K)) 2240 NEXT K,J,I 2250 ' 2260 FRE1=NPTS-1 : SIG#=SQR(SIG#/FRE1) 2270 FOR J=1 TO NTRM 2280 SIGX#(J)=SQR(SIGX#(J)/FRE1) 2290 R(J)=R(J)/(FRE1*SIGX#(J)*SIG#) 2300 FOR K=1 TO J 2310 ARY#(J,K)= ARY#(J,K)/(FRE1*SIGX#(J)*SIGX#(K)) 2320 ARY#(K,J)=ARY#(J,K) 2330 NEXT K 2340 NEXT J 2350 ' *** INVERT SYMMETRIC MATRIX 2360 NOR=NTRM : GOSUB 2830 :REM* CALL MATINV 2370 IF DET<>0 GOTO 2410 2380 A0=0 : SIGA0=0 : RMUL=0 : CHSQ=0 : FTEST=0 2390 RETURN 2400 ' *** CALC COEFF, FIT, & CHI SQUARE 2410 A0=YMN# 2420 FOR J=1 TO NTRM 2430 FOR K=1 TO NTRM 2440 A(J)=A(J) +R(K)*ARY#(J,K) 2450 NEXT K 2460 A(J)=A(J)*SIG#/SIGX#(J) 2470 A0=A0-A(J)*XMN#(J) 2480 FOR I=1 TO NPTS 2490 GOSUB 3390 :REM* CALC FCTN(X,I,J,M) 2500 YFIT(I)=YFIT(I) +A(J)*FCTN 2510 NEXT I,J 2520 ' 2530 FOR I=1 TO NPTS 2540 YFIT(I)=YFIT(I) +A0 2550 CHI#=CHI# +WGHT(I)*(Y(I)-YFIT(I))*(Y(I)-YFIT(I)) 2560 NEXT I 2570 ' 2580 FREN=NPTS-NTRM-1 : CHSQ=CHI#*WMN/FREN 2590 ' *** CALC UNCERTAINTIES 2600 CHSQ=CHI# 2610 IF MODE=0 GOTO 2630 2620 VARN=1/WMN : GOTO 2640 2630 VARN=CHSQ 2640 FOR J=1 TO NTRM 2650 SIGA(J)=ARY#(J,J)*VARN/(FRE1*SIGX#(J)*SIGX#(J)) 2660 SIGA(J)=SQR(SIGA(J)) 2670 RMUL=RMUL +A(J)*R(J)*SIGX#(J)/SIG# 2680 NEXT J 2690 FREJ=NTRM : FTST=(RMUL/FREJ)/((1-RMUL)/FREN) 2700 RMUL=SQR(RMUL) : SIGA0=VARN/NPTS 2710 FOR J=1 TO NTRM 2720 FOR K=1 TO NTRM 2730 SIGA0=SIGA0 +VARN*XMN#(J)*XMN#(K)*ARY#(J,K)/(FRE1*SIGX#(J)*SIGX#(K)) 2740 NEXT K,J 2750 SIGA0=SQR(SIGA0) 2760 RETURN 2770 REM---------------------------------------------------------------------- 2780 ' SUBROUTINE MATINV(ARY#,NOR,DET) 2790 ' A ROUTINE TO INVERT A SYMMETRIC MATRIX AND CALCULATE ITS 2800 ' DETERMINANT. ARY#(I,J) IS THE INPUT MATRIX OF ORDER=NOR 2810 ' AND IS REPLACED BY ITS INVERSE 2820 ' 2830 DET=1! 2840 ' *** FIND LARGEST ELEMENT 2850 FOR KK%=1 TO NOR 2860 AMX#=0! 2870 FOR II%=KK% TO NOR 2880 FOR JJ%=KK% TO NOR 2890 IF ABS(AMX#)>ABS(ARY#(II%,JJ%)) GOTO 2910 2900 AMX#=ARY#(II%,JJ%) : IK(KK%)=II% : JK(KK%)=JJ% 2910 NEXT JJ%,II% 2920 IF AMX#<>0! GOTO 2950 2930 ' *** INTERCHANGE ROWS & COLUMNS 2940 DET=0! : RETURN 2950 II%=IK(KK%) 2960 IF II%+1 GOTO 1850 1840 WGHT(I)=1/(SIGY(I)*SIGY(I)) : GOTO 1930 1850 IF MODE<>0 GOTO 1870 1860 WGHT(I)=1 : GOTO 1930 1870 IF MODE<>-1 GOTO 1900 1880 IF Y(I)<0 THEN WGHT(I)= 1/(-Y1) ELSE WGHT(I)= 1/Y1 1890 GOTO 1930 1900 IF MODE=-2 GOTO 1920 1910 PRINT" ERROR IN MODE SPECIFICATION" : RETURN 1920 WGHT(I)= SIGY(I) 1930 XTERM#= WGHT(I) 1940 FOR N=1 TO NMAX 1950 SUMX#(N)= SUMX#(N)+XTERM# 1960 XTERM#= XTERM#*X1 1970 NEXT N 1980 YTERM#= WGHT(I)*Y1 1990 FOR N=1 TO NTRM 2000 SUMY#(N)= SUMY#(N)+YTERM# 2010 YTERM#= YTERM#*X1 2020 NEXT N 2030 CHISQ#= CHISQ# +WGHT(I)*Y1*Y1 2040 NEXT I 2050 ' *** CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS 2060 FOR J=1 TO NTRM 2070 FOR K=1 TO NTRM 2080 N= J+K-1 2090 ARRAY#(J,K)= SUMX#(N) 2100 NEXT K,J 2110 GOSUB 2480 :' CALL DETERM 2115 DELTA= DETERM 2120 IF DELTA<>0 GOTO 2160 2130 CHISQ#=0! 2140 FOR J=1 TO NTRM : A(J)=0! : NEXT J : RETURN 2150 ' 2160 FOR L=1 TO NTRM 2170 FOR J=1 TO NTRM 2180 FOR K=1 TO NTRM 2190 N= J+K-1 2200 ARRAY#(J,K)= SUMX#(N) 2210 NEXT K 2220 ARRAY#(J,L)= SUMY#(J) 2230 NEXT J 2240 GOSUB 2480 : A(L)= DETERM/DELTA 2250 NEXT L 2260 ' 2270 ' *** CALC CHI SQUARE 2280 FOR J=1 TO NTRM 2290 CHISQ#= CHISQ# -2*A(J)*SUMY#(J) 2300 FOR K=1 TO NTRM 2310 N=J+K-1 2320 CHISQ#= CHISQ# +A(J)*A(K)*SUMX#(N) 2330 NEXT K,J 2340 ' 2350 FREE=NPTS-NTRM : CHISQR=CHISQ#/FREE 2360 RETURN 2370 REM---------------------------------------------------------------------- 2380 ' SUBROUTINE DETERM( ARRAY,NORDER ) 2390 ' 2400 ' A ROUTINE THAT CALCULATES THE DETERMINANT 2410 ' OF A SQUARE MATRIX. 2420 ' 2430 ' ARRAY= MATRIX 2440 ' NORDER= ORDER OF DETERMINANT (DEGREE OF MATRIX) 2450 ' 2460 ' NOTE: THIS ROUTINE DESTROYS THE INPUT MATRIX 2470 ' 2480 DETERM=1! : NORDER=NTRM 2490 ' *** INTERCHANGE COLUMNS IF DIAGONAL ELEMENT IS ZERO 2500 FOR KK=1 TO NORDER 2510 IF ARRAY#(KK,KK)<>0 GOTO 2640 2520 FOR JJ=KK TO NORDER 2530 IF ARRAY#(KK,JJ)<>0 GOTO 2560 2540 NEXT JJ 2550 DETERM=0! : RETURN 2560 FOR II=KK TO NORDER 2570 SAVE#= ARRAY#(II,JJ) 2580 ARRAY#(II,JJ)= ARRAY#(II,KK) 2590 ARRAY#(II,KK)= SAVE# 2600 NEXT II 2610 DETERM= -DETERM 2620 ' 2630 ' ***SUBTRACT ROW KK FROM LOWER ROWS TO GET DIAGONAL MATRIX 2640 DETERM= DETERM*ARRAY#(KK,KK) 2650 IF (KK-NORDER)=>0 GOTO 2710 2660 K1= KK+1 2670 FOR II=K1 TO NORDER 2680 FOR JJ=K1 TO NORDER 2690 ARRAY#(II,JJ)= ARRAY#(II,JJ)-ARRAY#(II,KK)*ARRAY#(KK,JJ)/ARRAY#(KK,KK) 2700 NEXT JJ,II 2710 NEXT KK 2720 RETURN  2690 ARRAY#(II,JJ)= ARRAY#(II,JJ)-ARRAY#(II,KK) CURVFIT.LIB I yo ar eve i th positio o havin t determin ho wel rea dat conform t som mathematica functio the yo wil fin thes program useful Th program i thi librar provid th basi tool t d thi job a leas fo case involvin singl independen variable Th rational behin th algorithm tha ar use wil b foun i th book, "Dat Reductio an Erro Analysi fo th Physica Sciences" b Phili R Bevington McGraw-Hil Boo Co. Ne York N.Y. 1969 Th boo i availabl i paperbac an i yo inten t mak us o thes routine urg yo t ge copy. Th basi concep behin curv fittin i ver simple I th cours o a experimen on ma collec se o dat point X(I) an Y(I) wher X(I i th independen variabl an Y(I i th dependen variable Usuall thes dat point contai som "noise bu on ma suspec tha i theory i relate t b a analyti function Y f(X) Hence fo ever poin X(I ther i correspondin valu Y(I tha i give b th formul f(X) Th functio f(X generall contain parameter tha w wil denot A(J) sinc ther ma b severa o them Th bes fi betwee th dat an f(X occur whe yo fin th value o A(J tha minimiz th mea squar deviatio betwee th calculate value YFIT(I) f(A,X an th dat Y(I) Thi i know a th techniqu o leas squares Wit th program i thi packag yo ca solv fo u t 1 parameter A(J) i an analyti functio wit u t 10 dat points Th numbe o parameter an dat point ca b extende i yo hav th memor an ar willin t wai fo th progra t d its thing Th fou program i th librar provid th tool t perfor curv fittin a fou differen level o complexity Thes program ar named: LINFIT.BAS POLYFIT.BAS MOSTFIT.BAS ANYFIT.BAS an ar accompanie b .DA file tha ar t b use fo simpl testing The ar writte i Microsof BASIC Versio 5.0 Th routine ar reasonabl wel commente an explai th meanin o th interna variables Yo ma recogniz generi similarit betwee thes routine an th som o th function availabl i th IB Scientifi Subroutin Package Bevington' boo provide mor in-dept documentation. *** LINFIT.BAS Th progra LINFI wil determin th parameter o straigh lin tha bes describe you data Th onl decisio yo hav t mak wit LINFI i ho yo wan t weigh th data Th dat se labele LINTEST1.DA wa generate fro th function, Y(I) 2. 5.7*X(I) Whe yo ru LINFIT th firs thin i wil d i reques th inpu dat se name I yo tel i LINTEST1.DAT i wil rea th dat fil fro dis (th currentl logge disk) Th dat fil wa simpl create wit a edito an consist o thre number pe lin separate b comas Th firs numbe i th independen variabl X(I) th secon numbe i th dependen variabl Y(I) an th thir numbe i th erro associate wit th measuremen o Y(I an i calle SIGY(I) Nex th progra wil as yo ho yo wan t weigh th dat b specifyin th mode Ente zer t specif MODE= whic mean tha al dat weight ar se t unit (equivalen t n weighting) I jus secon o th tw th progra the wil lis th intercep valu o 2. an th slop valu o 5. o th bes fittin straigh line Als liste ar th statistica error o eac paramete an th correlatio coefficien betwee th dat an th bes fi straigh line Lastly th progra summarize i columna format th inpu data th bes fi dat (i.e th dat predicte b usin th recentl determine value o slop an intercept) an th differenc betwee th dat an th bes fi (i.e th residual) I thi cas ther i n dat nois s th progra recover th slop an intercep exactly Not tha th onl non-zer residua i a th roundof level. A secon tes create th fil LINTEST2.DA i whic "twiddled th value o th Y(I slightl t simulat nois an assigne sigm accordingly Th equatio tha wa use t generate th dat was, Y(I) 12. +7.3*X(I) Tr runnin th progra unde differen mode an not ho th "bes fit wil vary Yo wil no b abl t recove th value 12. an 7. exactl bu yo shoul com prett clos wit MODE=1 Lis th dat se an not whic dat poin ha th larges sigm an whic dat poin ha th larges residua i th fi summary *** POLYFIT.BAS POLYFI wil fi you dat t polynomia o th form, Y A +A2* A3*X* A4*X*X* +..... wher yo mus specif th orde o th polynomial I yo specif NTRM= th progra reduce t LINFIT Tr runnin i wit LINTEST2.DA an NTRM=2 Not tha i take bi longe t ru an instea o tellin yo abou th correlatio coefficien i no list ch squar - th goodnes o th fit A i an fittin procedur yo mus hav mor dat point tha parameter t determine Th dat se containe i th fil POLYTEST.DA wa generate wit th function, Y 3. +1.7* +0.5*X* -0.782*X*X*X Yo wil se tha rounde of th value o whe constructe th tes data A resul th progra doe no comput th exac parameter give i th precedin function I i shor objec lesso regardin th problem o limite precisio (don' forge th compute truncate i th 7t decima plac i norma singl precision). *** MOSTFIT.BAS MOSTFI i progra tha wil fi t functio tha yo mus specif i th subroutin labele FUNCTN Th formulatio mus b linea i th coefficient A(J) i.e i mus b o th form; y= A(1) +A(2)*f(x) +A(3)*g(x) +..... wher f(x) g(x) etc. d no contai th parameter A(J). I th tes dat se MTEST1.DA use th function, Y(I)= 3.2 + 0.75 *SIN( X(I) ) to generate the data and coded FUNCTN to be, FUNCTN= A(1) +A(2)*SIN(X(I)) Fo MTEST2.DAT use th sam functio bu se A(1)=0. an A(2) 2.0 Not tha NTRM= i thi cas an th residual ar a th leve o fe part pe thousand Thi represent th quantizatio leve o th dat - th poin a whic truncate th inpu dat set. ** CURFIT.BAS CURFI i th mos genera fittin progra i th librar an allow yo t specif th functio f(X i th subroutin calle FUNCTN Th functio doe no hav t b continuou bu i doe hav t b analytic i.e i mus b define everywher no doe i hav t b linea i th coefficient (a i doe i MOSTFIT) B ver cautiou abou acceptin th statistic o th parameter determine b CURFIT particularl i yo ar usin weightin scheme Yo woul b advise t us Mont Carl schem t determin a estimat o th variatio i thos parameters A a exampl hav include tes cas tha i take fro Byt Magazin an wa describe i a articl title "Fittin Curve t Data" b Marc S Cacec an Willia P Cacheri (Ma 1984) Thi articl describe anothe algorith fo curv fittin an i anothe usefu too i yo hav th urg t "curv fit" Th functio use t generat th tes dat se CURTEST.DA i o th form; Y(I)= A(1)*X(I)/(X(I)+A(2)). Whe yo ru th progra tr NTRM= an MODE=0 Not tha th progra wil as yo fo th " priori" o initia value o th parameter an th delt value t use Thi i becaus th progra find th "bes fit b searchin fo th smalles valu o ch squared Thi i apparen durin executio a th progra loop an list th curren value o ch square an th interna paramete FLAMD (whic yo ca thin o a a invers gai parameter). hav use thes program fo th las 1 year o IB 360's 370's VA 11/780's an variou micro includin m homebre system The wor fo m bu sugges tha yo tes the thoroughl befor usin them Whe us the typicall us skeleto for o th progra an re-writ th drive containin th I/ par o th progra fo eac ne application als generat on o mor tes dat set (usuall wit know amount o noise i orde t tes th program. Th program enclose ar i BASI (Microsof Versio 5.2 whic i ver usefu fo smal dat sets I yo hav lo o dat t munc sugges yo ge th BASI compile - hav see speedup b factor o t 1 wit compile BASIC Bevington' boo list th program i FORTRA II woul appreciat hearin abou an bugs/improvement o translation t othe languages plac thes program int th publi domai an gran permissio t SIG/ t distribut th file fo non- commercia us i appreciatio o th man fin tool tha have taken from the public domain over the past several years. George M. Resch 1514 Stephora Ave. Glendora, CA 91740 Telenet mailbox - RADMETRIC/DIV420/JETPRO May 1, 19850.0,3.2,1.0 2.0,3.882,1.0 4.0,2.632,1.0 6.0,2.9904,1.0 1.68,0.172,0.01 3.33,0.250,0.01 5.00,0.286,0.01 6.67,0.303,0.03 10.0,0.334,0.02 20.0,0.384,0.02 -2.5,-6.1,1.3 1.0,19.5,1.0 5.0,48.1,1.2 12.0,84.8,10.5 0,0.3,.001 1,1.98,.001 2,2.12,.001 3,0.58,.001 4,-1.21,.001 5,-1.62,.001 6,-0.26,.001 7,1.61,.001 -2.5,14.29,1.7 1.0,4.618,1.0 5.0,-73.55,1.8 7.0,-228.626,2.0 10.0,-711.8,2.5 12.0,-1255.7,3.0  -2, -9.1, 0.0 1, 8.0, 0.0 4, 25.1, 0.0