10   '
20   '    LEAST SQUARES FIT OF A DATA SET TO A STRAIGHT LINE
30   '
40   ' ---------------------- 1/5/88 -----------------------
50   '
60   '          --> LIST OF PROCEDURES <--
70   '
80        GOSUB  150  ' INITIALIZATIONS
90        GOSUB  370  ' STORE DATA IN ARRAYS
100       GOSUB  850  ' CHISQ CALCULATION
110       GOSUB 1170  ' OUTPUT
120  '
130       END
140  '------------------------------------------------------
150  'INITIALIZATIONS
160  '
170  DIM X(100),Y(100),DY(100),R(100),Z$(101)
180  DATA  " ",".","*","O"
190  READ BLANK$,DOT$,ASTER$,OH$
200  '
210  PRINT
220  PRINT"BEFORE RUNNING THIS PROGRAM,"
230  PRINT"   YOU MUST PUT ALL OF YOUR DATA IN"
240  PRINT"   DATA STATEMENTS AFTER THE LAST PROGRAM LINE"
250  PRINT
260  PRINT"DATA MUST BE IN THE FORM:"
270  PRINT"     N"
280  PRINT"     X(1),Y(1),DY(1)"
290  PRINT"      .     .    ."
300  PRINT"      .     .    ."
310  PRINT"     X(N),Y(N),DY(N)"
320  PRINT
330  NR = 61
340  '
350  RETURN
360  '------------------------------------------------------
370  'STORE DATA IN ARRAYS
380  '
390  READ N
400  PRINT
410  PRINT N;" DATA POINTS, CURRENTLY."
420  PRINT
430  '
440  ' -- put in array: initial coordinates 
450  '    [X(1), Y(1)] & std. dev. DY(1)
460  READ X(1),Y(1),DY(1)
470  '
480  ' --  store initial values in sums
490  '
500  SIG = DY(1) ^ 2
510  IF DY(1) < = 0 THEN  LET SIG = 1!
520  S1 = 1 / SIG
530  XS = X(1) / SIG
540  XXS = X(1) ^ 2 / SIG
550  YS = Y(1) / SIG
560  XYS = X(1) * Y(1) / SIG
570  XMAX = X(1)
580  XLOW = X(1)
590  YLOW = Y(1)
600  YMAX = Y(1)
610  N =  INT (N)
620  '
630  ' -- put remaining points in arrays and do sume
640  '
650  FOR J = 2 TO N
660    READ X(J),Y(J),DY(J)
670    SIG = DY(J) ^ 2
680    IF DY(J) <  = 0 THEN  LET SIG = 1!
690    S1 = S1 + 1 / SIG
700    XS = XS + X(J) / SIG
710    YS = YS + Y(J) / SIG
720    XXS = XXS + X(J) ^ 2 / SIG
730    XYS = XYS + X(J) * Y(J) / SIG
740    '
750    ' -- find maximum and minimum values of X and Y
760    '
770    IF X(J) < XLOW THEN  LET XLOW = X(J)
780    IF X(J) > XMAX THEN  LET XMAX = X(J)
790    IF Y(J) < YLOW THEN  LET YLOW = Y(J)
800    IF Y(J) > YMAX THEN  LET YMAX = Y(J)
810    NEXT J
820    '
830  RETURN
840  '------------------------------------------------------
850  'CHISQ CALCULATION
860  '
870  ' -- compute intercept A1 and slope A2
880  '
890  A1 = (XXS * YS - XS * XYS) / (S1 * XXS - XS * XS)
900  A2 = (XS * YS - S1 * XYS) / (XS * XS - XXS * S1)
910  '
920  ' -- compute residuals R and chi-square CHISQ
930  '
940  CHISQ = 0
950  '
960  FOR J = 1 TO N
970    Y1 = A1 + A2 * X(J)
980    R(J) = Y(J) - Y1
990    SIG = DY(J) ^ 2
1000   IF DY(J) < = 0 THEN SIG = 1!
1010   CHISQ = CHISQ + R(J) ^ 2 / SIG
1020   NEXT J
1030   '
1040 ' -- Compute uncertainties D1A and D2A:
1050 '
1060 D1A =  SQR (XXS / (S1 * XXS - XS * XS))
1070 D2A =  SQR (S1 / (S1 * XXS - XS * XS))
1080 '
1090 ' -- if DY(1) = 0 then calculate 
1100 '    the uncertainties from chi-square
1110 IF DY(1) <> 0 THEN RAV = 1 ELSE RAV = SQR (CHISQ) / N
1120 D1A = D1A * RAV
1130 D2A = D2A * RAV
1140 '
1150 RETURN
1160 '------------------------------------------------------
1170 'OUTPUT
1180 '
1190 INPUT "PRINTER (Y/N)? ", ANS$
1200 IF ANS$ = "y" THEN ANS$ = "Y"
1210 IF ANS$ = "Y" THEN OPEN "LPT1:" FOR OUTPUT AS #1
1220 IF ANS$ <>"Y" THEN OPEN "SCRN:" FOR OUTPUT AS #1
1230 '
1240 PRINT #1,
1250 PRINT #1,"  X*     Y*    DY*     Y    RESIDUALS"
1260 FOR J = 1 TO N
1270   PRINT #1, USING "####  "; X(J);
1280   PRINT #1, USING "##.##  "; Y(J);
1290   PRINT #1, USING " #.##  "; DY(J);
1300   YP = A1 + A2 * X(J)
1310   PRINT #1, USING "##.##  "; YP;
1320   PRINT #1, USING " ##.#####"; R(J)
1330   NEXT J
1340   '
1350 PRINT #1,
1360 '
1370 ' --   DX and DY1 are step sizes (used for making plot)
1380 '
1390 DX = 1.2 * (XMAX - XLOW) / 50
1400 DY1 = 1.2 * (YMAX - YLOW) / (NR - 1)
1410 '
1420 ' the next two lines are an attempt to make the *'s
     ' form a straight line: remove the 2 lines if they give
     ' trouble
1440 '
1450 II= INT (A2 * DX / DY1 + .49999)
1460 DX = II * DY1 / A2
1470 '
1480 ' -- set first X value for plot
1490 '
1500 X1 = XLOW - 5 * DX
1510 FOR J = 1 TO NR
1520   Z$(J) = BLANK$
1530   NEXT J
1540 PRINT #1,"LEAST SQUARES FIT PARAMETERS ( Y = A1 + A2 @+
         +@* X ):"
1550 PRINT #1,USING "  A1 = ##.#####"; A1;
1560 PRINT #1,"  +/- ";
1570 PRINT #1,USING "##.#####"; D1A
1580 PRINT #1,USING "  A2 = ##.#####"; A2;
1590 PRINT #1,"  +/- ";
1600 PRINT #1,USING "##.#####"; D2A
1610 PRINT #1,USING "  CHISQ =####.#####"; CHISQ
1620 PRINT #1,
1630 PRINT #1,"   X      Y    ";
1640 FOR J = 1 TO NR
1650   PRINT #1,DOT$;
1660   NEXT J
1670 PRINT #1,
1680 '
1690 FOR J = 1 TO 51
1700   Z$(1) = DOT$
1710   Z$(NR) = DOT$
1720   '
1730   ' -- Y1 is the fitted Y value
1740   '
1750   Y1 = A1 + A2 * X1
1760   M1 =  INT ((Y1 - YLOW) / DY1) + NR / 10
1770   '
1780   ' -- if M1 is between 1 and NR, store an asterisk in
       '    Z$(M1)
1790   '
1800   IF M1 > 0 AND M1 < NR + 1 THEN  LET Z$(M1) = ASTER$
1810   FOR K = 1 TO N
1820     '
1830     ' -- if this X1 is within DX of any X(K),
         '    K=1,...,N, then store the letter O in Z$(M2)
1850     '
1860     IF  ABS (X1 - X(K)) > DX / 2 THEN 1890
1870     M2 =  INT ((Y(K) - YLOW) / DY1) + NR / 10
1880     Z$(M2) = OH$
1890     NEXT K
1900   PRINT #1,USING "##.###"; X1;
1910   PRINT #1,USING "###.###"; Y1;
1920   PRINT #1, "  ";
1930   FOR K = 1 TO NR
1940     PRINT #1,Z$(K);
1950     NEXT K
1960   PRINT #1,
1970   '
1980   ' -- restore Z$(M1) and Z$(M2) to blanks
1990   '
2000   Z$(M2) = BLANK$
2010   IF M1 > 0 AND M1 < NR + 1 THEN  LET Z$(M1) = BLANK$
2020   '
2030   ' -- do next fitted X value
2040   '
2050   X1 = X1 + DX
2060   NEXT J
2070 PRINT #1,TAB(16);STRING$(NR,DOT$)
2080 CLOSE #1
2090 RETURN
2100 '------------------- END OF PROGRAM LINES--------------
2110 'DATA TO BE FIT
2120 '
2130 DATA 15
2140 DATA 1,.7,.1,2,1.3,.1,3,1.9,.1,4,2.5,.1,5,3.1,.1
2150 DATA 6,3.7,.1,7,4.1,.1,8,4.9,.1,9,5.6,.1,10,6.1,.1,11@+
          +@,6.7,.1
2160 DATA 12,7.5,.1,13,7.9,.1,14,8.5,.1,15,9.1,.1
2170 '---------------------- END OF DATA -------------------