One of the most common tasks of a computer is to perform statistical analysis and data manipulative type computations. Many of the computations involved are very tedious to perform by hand since they sometimes require a large number of reiterative calculations. Statistical analysis problems that take hours for a human (equipped only with pencil and paper) are solved in a matter of seconds with a computer. This article presents a program which computes some simple statistical quantities commonly used by engineers and scientists.
Data collected in the laboratory or in the field often follows some logical pattern. The quantity of interest depends on or is a function of another quantity. By recording a sufficient amount of data, this relationship may be determined by fitting a function to the set of data points. It is rare that data collected in the real world corresponds exactly to some function; rather it is scattered about with a certain probability of corresponding to the function.
Consider a function which would describe the ambient temperature in Cove City, North Carolina, as a function of time for the last twenty years. If the function were to include all of the daily temperature fluctuations, it would indeed be a very complicated function. However, if only the average weekly or monthly temperatures were recorded, then a reasonable function could be fitted to the finite number of data points. The function would resemble a sinosoid with a period of 365 days with maxima and minima occurring in the late summer and winter respectively. There would be some uncertainty associated with correlating a particular day in the future with some specific temperature. In this case the uncertainty lies within some confidence limits determined by the climate. In other cases the uncertainty depends primarily on the method or apparatus used to measure the data.
Fitting a function to a set of data points is desirable because it allows one to make predictions or extrapolations from hypothetical or trial data.
The Method of Least Squares
In many situations the function which relates one quantity to another is linear. That is, if the data are plotted, the points lie on or about straight line. Examples include: the current through a resistor as a function of voltage, the conductance of pure water as a function of the amount of salt added or the length of a suspended spring as a function of the load attached.
An approximate method of determining the “best” straight line through a set of data points is by plotting the points on graph paper and drawing a line which uniformly divides all the points. Of course, this method is prone to error because determining the “best” straight line is somewhat subjective and two persons working with the same data will rarely arrive at identical functions.
The “best” straight line through a set of data points is defined as that line in which the sum of the squares of the deviations of all points from the line is a minimum. This method of “least squares” is far superior because it calculates, not approximates, the line. The method yields two quantities: the slope and the intercept of the line, thereby defining the linear function by the well-known relation
y = ax + b
where a is the slope and b is the y-intercept.
The slope and intercept are determined by the formulas in Figure 1. Another useful quantity is the correlation coefficient, given in Figure 2, which is the relative amount of scatter about the line. A correlation coefficient of 1 means that all data points lie exactly on the least squares line. A correlation coefficient of less than 0.9 implies that the data are either invalid or are not linearly related.
Figure 1. (The summation signs appear to have dropped out of the printed formulas; they are reproduced as printed.)
x y
xy - ---
N
slope = a = ------------
( x)²
x² - ----
N
y - m x
y-intercept = b = ---------
N
Figure 2.
( x)²
x² - -----
N
----------
N - 1
correlation coefficient = c = m ----------------
( y)²
y² - -----
N
----------
N - 1
The Program
Though the program as written will not fit into a 1K system, it can easily be broken down into independent modules so that 1K owners should still be able to benefit from this article. The DIM statements at lines 10 and 11 limit the maximum number of data points to 20. The FOR-NEXT loop at lines 140-180 repeatedly invokes the subroutine at line 1000 to enter the data points into the X and Y arrays. The main program is made from lines 200-400. These statements produce a menu which gives the user several paths to follow. The user enters a number from 1 to 5 which is stored in the variable G. The program calculates the appropriate line number from G and branches accordingly. (Note that branching to a line number specified by a variable is one feature of the Sinclair Basic not found in many others. This feature offers many advantages including increased readability; it is much easier to comprehend a statement that says GOSUB PLOTSUB than one that says GOSUB 130.) The options given to the user are: calculating the least squares line, printing a table of the data, plotting the data, starting over or ending the program.
The analytical calculations are done in the subroutine starting at line 3000. As you can see, the summations indicated in the equations in Figures 1 and 2 are carried out in a FOR-NEXT loop from 1 TO N, where N is the number of points. Line 3340 prints the equation of the line from the calculated slope and y-intercept. The correlation coefficient is calculated at lines 3360-3400 and is printed at line 3420. This module, along with some way of getting data into the X and Y arrays, may be used as the main program in a 1K system.
The plotting routine is useful as a qualitative assessment of the data. By plotting the data on the screen one can see at a glance if the data are well behaved: in this case linear. In some cases the user may or may not wish to see a plot which includes the origin. For instance, a straight-line plot with a high positive y-intercept would not utilize the full resolution of the 44 x 64 pixel screen. Instead, the line would look compressed against the top of the screen. For this reason the user is given the option of having the graph include the origin or the point (0,0) which is in the lower left corner of the screen. This is a rather unsophisticated plotting routine and is only valid for positive values. The user may wish to substitute another. A more sophisticated routine would allow negative values and position the origin accordingly. Labeled axes would also be a useful enhancement.
Note the instruction POKE 16437,255 at line 5581. This is necessary for ZX80 8K ROM users. It must follow every PAUSE instruction. Failure to include it can result in a system crash. The reason for this is not clear to me, but it is apparently the result of a bug in the ROM. It has something to do with a counter that both the PAUSE instruction and video circuitry use. Without the extra instruction the system may work fine for as long as ten minutes then suddenly produce a blank screen with no response to keyboard input. A power down or reset is the only road to recovery which, of course, results in losing the program which you had just spent the last two hours typing in. I had to find out the hard way, and, after several weeks of testing, probing, and endangering my sanity, I happened to stumble across the solution on page 127 of ZX81 BASIC Programming, a typical case of “when all else fails, read the instructions.” Be forewarned.
Example
Perhaps the best explanation of the method of “least squares” is an example. Let’s say that we are looking for the function which describes the forward voltage drop across a transistor as a function of temperature. In order to use the least squares method we must have a function which we will get back as a straight line. Therefore, we must assume that the relation of the voltage drop across a transistor is linear. Well, it is not quite linear, but close enough when operated between 0 and 100 degrees Celsius. The data in Figure 3 were collected using a 4-digit multimeter and a National Bureau of Standards mercury thermometer with absolute accuracy to 1/100 of a degree.
Figure 3: Voltage/Temperature Table
X Y Voltage Drop Temperature (volts) (degree Celcius) (1) 0.6853 0.00 (2) 0.6752 5.00 (3) 0.6647 10.00 (4) 0.6542 15.00 (5) 0.6440 20.00 (6) 0.6337 25.00 (7) 0.6234 30.00 (8) 0.6024 40.00 (9) 0.5811 50.00 (10) 0.5595 60.00 (11) 0.5377 70.00 (12) 0.5264 75.00 (13) 0.5163 80.00 (14) 0.5054 85.00 (15) 0.4942 90.00 (16) 0.4821 95.00
From this table and with the aid of the program in Figure 4 we should easily be able to determine the temperature of the transistor casing (or ambient temperature if the device is in the surrounding environment) by measuring its forward voltage drop. Conversely, we could determine the forward voltage drop if we knew the temperature. Suppose we measure a voltage drop of 0.6295 volts across the transistor and we wish to determine the temperature accurately. We know that it should be between 25 and 30 degrees just by looking at the table. But we want to be a little more accurate than that. Let’s run the program and see.
The program first asks for the number of data points. We have collected 16 known data points so enter 16. The program then asks for the data, X1 then Y1, X2 then Y2, etc., until the Y component of the last data point is entered. Then a menu is displayed so we can tell the computer what to do with the data: 1) perform least squares analysis, 2) list the numerical values, 3) plot the data, 4) start over, or 5) quit. The best thing to do here is list the data to double check all of the entries. If a mistake is found, the program forgivingly offers the option of correcting it without entering the entire collection again. After entry corrections, the program waits for another NEWLINE before going back to the menu. With the corrected data now in the computer we are ready to perform a least squares analysis. The program responds with the slope, y-intercept, and the equation of the best straight line through the set of data points. The unknown quantity can be immediately calculated by substituting the corresponding known quantity into the equation.
y = ax + b
or
(y-b)
x = -----
a
In this case we know x.
y = -469.02146 deg/V* 0.6295 V + 321.99771 deg = 26.75 deg
The result is valid to the number of digits of the input data. The same results would have been obtained had we switched the coordinates, that is, if we had let the voltage drop be the X value and the temperature be the Y value. Then the unknown temperature would have been determined by substitution of the transistor’s voltage drop into y. Finally, we can plot the data by hitting NEWLINE to get back to the menu and then entering 3. The sample plot reveals that the data are hardly scattered, as also predicted by the correlation coefficient of very near unity.
Hopefully, this article will be of help or interest to many readers. This article barely scratches the surface of statistical analysis and function fitting. Bear in mind that the least squares technique is not limited to linear functions. Variations on the method can be used to fit polynomials of any degree to a set of data points.
Note: Readers are referred to Jon Passler’s “Linear Regression” in SYNC 2:1 for more information on relating variables.
Figure 4: Program Listing
[inv …] marks inverse-video characters.
1 REM ... LEAST SQURES
2 REM ... ANALYSIS PROGRAM
3 REM ... BY
4 REM ... ALGER SALT
5 REM ... JANUARY 1982
10 DIM X(20)
11 DIM Y(20)
100 PRINT "LEAST SQUARES DATA ANALYSIS"
105 PRINT
106 PRINT
110 PRINT "[inv ENTER NUMBER OF POINTS]"
120 INPUT N
121 CLS
140 FOR I=1 TO N
160 GOSUB 1000
180 NEXT I
200 PRINT
220 PRINT "1- LEAST SQUARES"
240 PRINT "2- LIST DATA"
260 PRINT "3- PLOT DATA"
270 PRINT "4- START AGAIN"
280 PRINT "5- QUIT"
290 PRINT
300 PRINT "[inv ENTER 1 TO 5]"
320 INPUT G
321 CLS
340 GOSUB G*1000+2000
360 INPUT A$
361 CLS
400 GOTO 200
999 STOP
1000 PRINT AT 10,10;"[inv ENTER] X(";I;")"
1020 INPUT X(I)
1021 CLS
1040 PRINT AT 10,10;"[inv ENTER] Y(";I;")"
1060 INPUT Y(I)
1061 CLS
1099 RETURN
3000 REM ... LEAST SQUARES
3010 LET S0=0
3020 LET S1=0
3040 LET S2=0
3060 LET S3=0
3080 LET S4=0
3100 FOR I=1 TO N
3110 LET S0=S0+Y(I)**2
3120 LET S1=S1+X(I)
3140 LET S2=S2+Y(I)
3160 LET S3=S3+X(I)*Y(I)
3180 LET S4=S4+X(I)**2
3200 NEXT I
3220 LET S5=S1**2
3240 LET M=((S1*S2/N)-S3)/((S5/N)-S4)
3260 LET B=((-M*S1)+S2)/N
3270 PRINT
3280 PRINT "SLOPE","Y INTCP"
3290 PRINT
3300 PRINT M,B
3310 PRINT
3320 PRINT
3330 PRINT
3335 PRINT "EQUATION OF THE LINE"
3336 PRINT
3340 PRINT "Y = ";M;"*X + (";B;")"
3360 LET SY=SQR ((S0-(S2**2)/N)/(N-1))
3380 LET SX=SQR ((S4-(S1**2)/N)/(N-1))
3400 LET CC=M*SX/SY
3410 PRINT
3411 PRINT
3420 PRINT "COR. COEF. = ";CC
3999 RETURN
4000 REM ... LISTING ROUTINE
4005 PRINT AT 0,0;"[inv I] ";" [inv X]","[inv Y]"
4010 PRINT
4100 FOR I=1 TO N
4110 LET R=I+1
4120 PRINT AT R,0;I;
4121 PRINT AT R,3;"[graphic]";
4125 PRINT AT R,4;X(I);
4130 PRINT AT R,16;Y(I)
4140 NEXT I
4160 PRINT
4200 PRINT "[inv DO YOU WANT TO CHANGE ONE?]"
4220 INPUT A$
4240 IF A$<>"Y" THEN GOTO 4300
4260 PRINT "[inv WHICH ONE?]"
4280 INPUT I
4281 CLS
4290 GOSUB 1000
4295 GOTO 4000
4300 RETURN
5000 REM ... PLOT ROUTINE
5010 LET XL=X(1)
5020 LET YL=Y(1)
5040 LET XH=X(1)
5060 LET YH=Y(1)
5099 REM ...FIND MAX AND MIN X,Y
5400 FOR I=1 TO N
5420 IF X(I)<XL THEN LET XL=X(I)
5440 IF Y(I)<YL THEN LET YL=Y(I)
5460 IF X(I)>XH THEN LET XH=X(I)
5480 IF Y(I)>YH THEN LET YH=Y(I)
5500 NEXT I
5510 PRINT "[inv INCLUDE ORIGIN?]"
5511 INPUT A$
5512 CLS
5513 IF A$<>"Y" THEN GOTO 5528
5516 LET XL=0
5517 LET YL=0
5520 FOR I=0 TO 63
5521 PLOT I,0
5523 NEXT I
5524 FOR I=0 TO 41
5525 PLOT 0,I
5526 NEXT I
5528 FOR I=1 TO N
5530 LET XN=(X(I)-XL)/(XH-XL)*41
5540 LET YN=(Y(I)-YL)/(YH-YL)*41
5560 PLOT XN,YN
5580 PAUSE 30
5581 POKE 16437,255
5600 NEXT I
5999 RETURN
6000 GOTO 1
9999 STOPNote: Type-in program listings on this website use ZMAKEBAS notation for graphics characters.