1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Least Squares Fit In Fortran

  1. Mar 25, 2012 #1
    (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 acknowledgement
    !
    !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!
     
  2. jcsd
  3. Mar 25, 2012 #2
    http://www.bravus.com/formulae.jpg [Broken]
     
    Last edited by a moderator: May 5, 2017
  4. Mar 25, 2012 #3
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook