Fortran Program for Least Squares Fit of Data

[Bravus]

(The usual format doesn't really fit this question, so I hope you won't mind if I set it up slightly differently)

The task is to make a Fortran program that will read in a given set of data from a file, stopping when it reaches the end, and calculate the gradient and intercept of a best fit line.

Here's my source code:

! Least Squares Fit program
!
! David Geelan, 2012 - Code free to use with acknowledgment
!
!23456
PROGRAM Least
IMPLICIT NONE
!
REAL :: E,a,a1,a2,b,x,y ! Same variables as worksheet
REAL :: sigmax, sigmay, sigmaxy, sigmax2 ! Sum variables
INTEGER :: i,n,IO ! Counting variables
!
! Open input file for reading
!
n=0 ! Initialise n
! Initialise all sum variables
sigmax=0
sigmay=0
sigmaxy=0
sigmax2=0
!
OPEN(unit=8,file="ws08_leastsquaresfit.dat",action='read')
!
100 FORMAT (1ES8.1,1X,1ES8.1)
DO
READ (8,100,IOSTAT=IO) x,y
!
IF (IO<0) EXIT
IF (IO>0) THEN
STOP "Something's up with the data"
!
END IF
!
n=n+1
sigmax=sigmax+x
sigmay=sigmay+y
sigmaxy=sigmaxy+(x*y)
sigmax2=sigmax2+x**2
WRITE (6,*) n,x,y,sigmax,sigmay,sigmaxy,sigmax2
!
END DO
!
a1=(n*sigmaxy)-(sigmax*sigmay)
a2=(n*sigmax2)-(sigmax)**2
a=a1/a2
b=(sigmay-(a*sigmax))/n
WRITE (6,*) a,b
!
REWIND(unit=8) ! Reset data file to start for reading
!
E=0 ! Initialise E
!
DO
READ (8,100,IOSTAT=IO) x,y
IF (IO<0) EXIT
IF (IO>0) THEN
STOP "Something's up with the data"
END IF
!
E=E+(y-(a*x)+b)**2
!
END DO
!
WRITE (6,*) 'For this model, a=',a,', b=',b
WRITE (6,*) 'Average square error=',E/n
!
CLOSE(unit=8) ! Close input file
!
END PROGRAM Least

It compiles and runs, and even outputs values... but they're not the right size! No idea what is going wrong. I set it to output the variables at each iteraction, and here's the full output when it's run:

bravus@bravus-VirtualBox:~/ws08_leastsquaresfit$ Least.exe
1 20.000000 3.0000000 20.000000 3.0000000 60.000000 400.00000
2 22.000000 4.0000000 42.000000 7.0000000 148.00000 884.00000
3 24.000000 4.0000000 66.000000 11.000000 244.00000 1460.0000
4 26.000000 5.0000000 92.000000 16.000000 374.00000 2136.0000
5 28.000000 6.0000000 120.00000 22.000000 542.00000 2920.0000
6 31.000000 8.0000000 151.00000 30.000000 790.00000 3881.0000
7 33.000000 7.0000000 184.00000 37.000000 1021.0000 4970.0000
8 34.000000 9.0000000 218.00000 46.000000 1327.0000 6126.0000
9 37.000000 8.0000000 255.00000 54.000000 1623.0000 7495.0000
0.34444445 -3.7592595
For this model, a= 0.34444445 , b= -3.7592595
Average square error= 56.968864

The instructor said the error should be at or less than one, so clearly something is badly wrong.

Should the formulae for a (which I broke into 2 steps for clarity) and b be inside the loop and iterated up, rather than working with the sums after the loop finishes?

Any help and hints gratefully received!

 

[Bravus]

http://www.bravus.com/formulae.jpg

 

[Bravus]

Thanks but never mind: a and b were being calculated correctly (I was misled by something else), and E just needed an extra set of brackets to sort the calculation order.