RBBS in a Box Volume 1 #3.1/RBBSIABOX31.cdr/rock/regress.bas

BASIC Source File | 1984-06-23 | 6KB | 128 lines

10 ' Programme to calculate correlation coeffiecnt and best fit line
20 ' through a set of data points by Mark A. Herkommer
30 ' Programme is written for use on a IBM PC running Basica.
40 '
50 OPTION BASE 1:KEY OFF:DEFINT I-N:SCREEN 0,0,0
60 '
70 DIM X(1000),Y(1000)
80 '
90 CLS:PRINT STRING$(34,"=");" REGRESS.BAS ";STRING$(33,"=")
100 PRINT "Mark A. Herkommer":PRINT "4005 Burning Tree Lane":PRINT "Garland, TX
110 PRINT
120 PRINT "REGRESS.BAS is a program for the IBM PC written in BASICA language.
130 PRINT "The program calculates the best-fit line through a set of XY data
140 PRINT "points by the method of linear least-squares.  The best fitting
150 PRINT "line is the regression curve of Y on X (that is, Y is the dependent
160 PRINT "variable and X is the independent variable).  Standard Error,
170 PRINT "Standard Deviation, and the Correlation Coefficent are determined.
180 PRINT "The equation of the line in slope-intercept form is :":PRINT
190 PRINT TAB(19)"Y = slope * X + Y-intercept  (Y = m*X + b)":PRINT
200 PRINT "Up to 1000 XY data points can be used in the calculation, however,
210 PRINT "input data must be in a sequential disk file and in ASCII format.
220 PRINT "Use shift-Prtsc to have the statistics printed on your line printer.
230 PRINT "After the statistics are calculated, press any key to make a graph.
240 PRINT "When the graph is finished (indicated with a beep), press any key to
250 PRINT "end the program.":LOCATE 25,1
260 INPUT "Enter the name of the file that holds the X,Y data  : ",INFIL$
270 '
280 ' read input data and calculate sums
290 '
300 XMIN=1E+32:XMAX=-1E+32:YMIN=1E+32:YMAX=-1E+32
310 IF INFIL$="" THEN 90
320 OPEN INFIL$ FOR INPUT AS #1
330 IF EOF(1) THEN 420 ELSE I=I+1
340 INPUT #1,X(I),Y(I)
350 SUMX=SUMX+X(I):SUMY=SUMY+Y(I)
360 IF X(I)>XMAX THEN XMAX=X(I):INDXMAX=I
370 IF X(I)<XMIN THEN XMIN=X(I):INDXMIN=I
380 IF Y(I)>YMAX THEN YMAX=Y(I)
390 IF Y(I)<YMIN THEN YMIN=Y(I)
400 SUMXY=SUMXY+X(I)*Y(I):SUMXX=SUMXX+X(I)*X(I):SUMYY=SUMYY+Y(I)*Y(I)
410 GOTO 330
420 CLOSE #1:CLS:NP=I
430 XMEAN=SUMX/NP:YMEAN=SUMY/NP
440 '
450 PRINT STRING$(80,"=")
460 PRINT "Total number of points read from ";INFIL$;TAB(57)"="NP:PRINT
470 PRINT "X minimum        =";XMIN;TAB(40);"Y minimum        =";YMIN
480 PRINT "X maximum        =";XMAX;TAB(40);"Y maximum        =";YMAX
490 PRINT "Range  of X      =";XMAX-XMIN;TAB(40);"Range of Y       =";YMAX-YMIN
500 PRINT "Mean value of X  =";XMEAN;TAB(40);"Mean value of Y  =";YMEAN
510 '
520 ' Linear least-squares fit of the data
530 '
540 PRINT:PRINT STRING$(80,"="):PRINT "Equation of the best fit line =====>";
550 SLOPE=(NP*SUMXY-SUMX*SUMY)/(NP*SUMXX-SUMX*SUMX)      ' slope of line
560 B=YMEAN-SLOPE*XMEAN                                  ' y-intercept
570 PRINT TAB(40)"Y = ";SLOPE;"* X +";B
580 '
590 ' Calculate Standard deviation and Standard Error
600 '
610 FOR I=1 TO NP
620 YC=SLOPE*X(I)+B     ' calculated y-value
630 D=Y(I)-YC           ' deviation
640 SUMDD=SUMDD+D*D     ' sum of the deviation squared
650 DX=X(I)-XMEAN:SUMDXDX=SUMDXDX+DX*DX    ' X-deviation squared
660 DY=Y(I)-YMEAN:SUMDYDY=SUMDYDY+DY*DY    ' Y-deviation squared
670 SUMDXDY=SUMDXDY+DX*DY
680 NEXT
690 STDERR=SQR(SUMDD/(NP-1))                        ' standard error
700 SIGMAX=SQR((SUMXX-SUMX*SUMX/NP)/(NP-1))         ' standard deviation of X
710 SIGMAY=SQR((SUMYY-SUMY*SUMY/NP)/(NP-1))         ' standard deviation of Y
720 IF SIGMAY=0 THEN PRINT "Y values are invariate-R is undefined":GOTO 800
730 R=SUMDXDY/SQR(SUMDXDX*SUMDYDY)
740 PRINT
750 PRINT STRING$(80,"="):PRINT "Elementary statistics on the best fit line"
760 PRINT:PRINT "Standard Error              =";STDERR
770 PRINT "Standard Deviation of X     =";SIGMAX
780 PRINT "Standard Deviation of Y     =";SIGMAY
790 PRINT "Correlation Coefficent      =";R
800 PRINT:PRINT STRING$(80,"=");
810 LOCATE 25,1:PRINT "press any key to produce a graph...";
820 I$=INKEY$:IF I$="" THEN 820 ELSE GOSUB 860
830 SCREEN 0,0,0:CLS:END
840 '
850 ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
860 ' Subroutine plots a graph of any XY data
870 '
880 SCREEN 2:CLS
890 '
900 ' Calculate increment and axis
910 '
920 XLEN=530:YLEN=176:NTICS=5
930 XMAX=INT(XMAX+.9990001):XMIN=INT(XMIN-.9990001)
940 XINC=(XMAX-XMIN)/XLEN
950 YTEMP=SLOPE*X(INDXMIN)+B:IF YTEMP<YMIN THEN YMIN=YTEMP ELSE:IF YTEMP>YMAX       THEN YMAX=YTEMP
960 YTEMP=SLOPE*X(INDXMAX)+B:IF YTEMP<YMIN THEN YMIN=YTEMP ELSE:IF YTEMP>YMAX       THEN YMAX=YTEMP
970 YMAX=INT(YMAX+.9990001):YMIN=INT(YMIN-.9990001)
980 YINC=(YMAX-YMIN)/YLEN
990 '
1000 ' Plot axis
1010 '
1020 LINE (79,YLEN+3)-(79+XLEN,YLEN+3)                   ' x-axis
1030 LINE (79,YLEN+3)-(79,3):LINE (80,YLEN+3)-(80,3)     ' y-axis
1040 XLABINC=(XMAX-XMIN)/NTICS:YLABINC=(YMAX-YMIN)/NTICS
1050 FOR I=0 TO NTICS
1060 XLOC=79+I*(XLEN-1)/NTICS:LINE (XLOC,YLEN+4)-(XLOC,YLEN+6):XLOC=XLOC+1:          LINE (XLOC,YLEN+4)-(XLOC,YLEN+6)                    ' x-axis tic marks
1070 XLABLOC=6+I*(XLEN+1)/(8*NTICS):LOCATE 24,XLABLOC:                               PRINT USING "##.#^^^^";XMIN+XLABINC*I;
1080 YLOC=I*YLEN/NTICS:LINE (73,YLOC+3)-(79,YLOC+3)      ' y-axis tic marks
1090 YLABLOC=1+I*22/NTICS:LOCATE YLABLOC,1:PRINT USING "##.#^^^^";YMAX-YLABINC*I
1100 NEXT
1110 LOCATE 25,1:PRINT "Y / X";
1120 '
1130 ' Plot the points
1140 '
1150 FOR I=1 TO NP
1160 XLOC=((X(I)-XMIN)/(XMAX-XMIN))*XLEN+80
1170 YLOC=((YMAX-Y(I))/(YMAX-YMIN))*YLEN+3
1180 PSET (XLOC,YLOC)                            ' center of data point
1190 DRAW "U2 L1 D4 R1 U2 NL4 NR3"               ' data point symbol
1200 NEXT
1210 '
1220 ' Draw best fit line
1230 '
1240 PSET (((X(INDXMIN)-XMIN)/(XMAX-XMIN))*XLEN+80,((YMAX-(SLOPE*X(INDXMIN)+B))      /(YMAX-YMIN))*YLEN+4)
1250 LINE -(((X(INDXMAX)-XMIN)/(XMAX-XMIN))*XLEN+80,((YMAX-(SLOPE*X(INDXMAX)+B))     /(YMAX-YMIN))*YLEN+4)
1260 BEEP
1270 I$=INKEY$:IF I$="" THEN 1270 ELSE RETURN
IN))*XLEN+80,((YMAX-(SL