Comp Sci Fortran Program for Least Squares Fit of Data

AI Thread Summary
The discussion revolves around a Fortran program designed to perform a least squares fit of data from a file, calculating the gradient and intercept of the best fit line. The user initially faced issues with incorrect output values, particularly with the average square error, which was expected to be less than one. After troubleshooting, it was determined that the calculations for the coefficients a and b were correct, but the average square error calculation required additional parentheses to ensure proper order of operations. The program compiles and runs successfully, but the user emphasizes the importance of careful attention to mathematical expressions in the code. Ultimately, the issue was resolved by adjusting the calculation of the average square error.
Bravus
Messages
22
Reaction score
0
(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!
 
Physics news on Phys.org
http://www.bravus.com/formulae.jpg
 
Last edited by a moderator:
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.
 
Back
Top