DECLARE SUB qwlinfit (n%, x!(), y!(), m!, sm!, b!, sb!)
n% = 7

DIM x!(1 TO n%)
DIM y!(1 TO n%)

x!(1) = -3!
x!(2) = -2!
x!(3) = -1!
x!(4) = 0!
x!(5) = 1!
x!(6) = 2!
x!(7) = 3!

y!(1) = -.0018278
y!(2) = -.0012146
y!(3) = -.0006021
y!(4) = .0000113
y!(5) = .0006262
y!(6) = .0012412
y!(7) = .001854

CALL qwlinfit(n%, x!(), y!(), m!, sm!, b!, sb!)

PRINT "m! = "; m!; " +- "; sm!
PRINT "b! = "; b!; " +- "; sb!


SUB qwlinfit (n%, x!(), y!(), m!, sm!, b!, sb!)

'****************************************************************************
'This subroutine does a least squares linear fit to the (x,y) data points.
'It returns the slope and intercept of the line with errors.
'See Bevington, pp. 104, 114
'
'Input:
' n%, the number of (x,y) points
' x!(1 to n%), the x values
' y!(1 to n%), the y values
'
'Output:
' m!, the fitted line's slope
' sm!, error estimate for m
' b!, the fitted line's intercept
' sb!, error estimate for b
'
'Zachary Wolf
'5/28/95
'****************************************************************************

'Compute various sums
SUMxy! = 0!
SUMxx! = 0!
SUMx! = 0!
SUMy! = 0!
FOR i% = 1 TO n%
 SUMxy! = SUMxy! + x!(i%) * y!(i%)
 SUMxx! = SUMxx! + x!(i%) * x!(i%)
 SUMx! = SUMx! + x!(i%)
 SUMy! = SUMy! + y!(i%)
NEXT i%

'Compute the determinant
d! = n% * SUMxx! - SUMx! * SUMx!

'Compute the slope and intercept
m! = (n% * SUMxy! - SUMx! * SUMy!) / d!
b! = (SUMxx! * SUMy! - SUMxy! * SUMx!) / d!

'Compute the errors

'Estimate the error on each y value
sigsq! = 0!
FOR i% = 1 TO n%
 sigsq! = sigsq! + (y!(i%) - (m! * x!(i%) + b!)) ^ 2
NEXT i%
sigsq! = sigsq! / (n% - 2)

'Use the estimated y error to find sb! and sm!
sb! = SQR(sigsq! * SUMxx / d!)
sm! = SQR(n% * sigsq! / d!)

END SUB