Fortran77 woes- single parameter least squares minimisation stuck

Click For Summary
SUMMARY

The forum discussion centers on a Fortran77 program designed for single parameter least squares minimization, specifically targeting the chi-squared value associated with experimental data fitting. The user, Nathan, encounters an infinite loop while attempting to minimize the chi-squared value, which should converge to 6618.681 but instead fluctuates between 6921.866 and 6920.031. Key issues identified include improper array indexing and unclear variable naming conventions, which hinder code readability and functionality.

PREREQUISITES
  • Understanding of Fortran77 programming language
  • Familiarity with least squares minimization techniques
  • Knowledge of chi-squared statistical analysis
  • Experience with debugging array operations in programming
NEXT STEPS
  • Review Fortran77 array handling and indexing practices
  • Study least squares fitting methods and their implementation in Fortran
  • Learn about debugging techniques specific to Fortran77
  • Explore best practices for variable naming and code organization in programming
USEFUL FOR

This discussion is beneficial for Fortran developers, data analysts working with statistical models, and anyone involved in numerical methods for data fitting and optimization.

NathanM
Messages
14
Reaction score
0
Hey guys! Having issues with a code I'm writing. It doesn't fit the typical format and I may have already accidentally posted it in the wrong forum, so here's a copy of what was there with the newest version of the code. It'll (eventually) be a four-parameter least squares minimisation, but for the moment I'm just trying to get the single parameter version to work. I've got the first module done but, when I run it it gets caught up in a loop and doesn't actually minimise anything! Good start huh?

Here's my code, bless me with your insight!

PROGRAM assignment

! A program designed to fit experiemental data, using the method
! of least squares to minimise the associated chi-squared and
! obtain the four control parameters A,B,h1 and h2.
!*****************************************************************

IMPLICIT NONE
INTEGER i,t0
DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
DOUBLE PRECISION deltah,Cs(17)

OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
DO i=1,17
READ(21,*)t(i),Ct(i),eCt(i)
END DO
CLOSE(21)

!Read in data.txt as three one dimensional arrays.
!*****************************************************************
!OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
!DO i=1,17
! WRITE(21,*)t(i),Ct(i),eCt(i)
!END DO
!CLOSE(21)
!
!Just to check input file is being read correctly.
!*****************************************************************
!****************************Parameters***************************

Aa= 0
Ba= 0
h1a= 0
h2a= 0

!*****************************************************************
!**********************Minimising Lamda1 (h1)*********************
deltah= 0.0001
h1b= 0.001
h1a= 0
DO 10
chisqa= 0
DO 20 i= 1,17
Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
Ctdiff= Ct - Cs
c= Ctdiff**2/eCt**2
chisqa= chisqa + c(i)
20 END DO
!Use initial h1 value of 0 to calculate start-point chisq.

chisqb= 0
DO 30 i= 1,17
h1b= h1b + deltah
Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
Ctdiff(i)= Ct(i) - Cs(i)
c= Ctdiff(i)**2/eCt(i)**2
chisqb= chisqb + c(i)
30 END DO
!First-step h1 used to find competing chisq for comparison.
WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
!Prints the two calculated chisq values to screen.
dchisq= chisqa - chisqb
IF (dchisq.GT.0) THEN
h1a= h1b
ELSE IF (dchisq.LE.0) THEN
deltah= deltah-((deltah*2)/100)

END IF
IF (chisqb.LE.6618.681) EXIT
10 END DO
WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
!*****************************************************************

END PROGRAM assignmentThanks for the help.
Nathan

EDIT: I should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!
 
Physics news on Phys.org
NathanM said:
Hey guys! Having issues with a code I'm writing. It doesn't fit the typical format and I may have already accidentally posted it in the wrong forum, so here's a copy of what was there with the newest version of the code. It'll (eventually) be a four-parameter least squares minimisation, but for the moment I'm just trying to get the single parameter version to work. I've got the first module done but, when I run it it gets caught up in a loop and doesn't actually minimise anything! Good start huh?

Here's my code, bless me with your insight!
Code:
      PROGRAM assignment
      
      ! A program designed to fit experiemental data, using the method
      ! of least squares to minimise the associated chi-squared and
      ! obtain the four control parameters A,B,h1 and h2.
      !*****************************************************************

      IMPLICIT NONE
      INTEGER i,t0
      DOUBLE PRECISION t(17),t2(17),Ct(17),eCt(17),Ctdiff(17),c(17)
      DOUBLE PRECISION Aa,Ab,Ba,Bb,h1a,h1b,h2a,h2b,chisqa,chisqb,dchisq
      DOUBLE PRECISION deltah,Cs(17)

        OPEN(21, FILE='data.txt', FORM='FORMATTED', STATUS='OLD')
           DO i=1,17
              READ(21,*)t(i),Ct(i),eCt(i)
           END DO
        CLOSE(21)

      !Read in data.txt as three one dimensional arrays.
      !*****************************************************************
      !OPEN(21, FILE='outtest.txt', FORM='FORMATTED', STATUS='NEW')
      !DO i=1,17
      !   WRITE(21,*)t(i),Ct(i),eCt(i)
      !END DO
      !CLOSE(21)
      !
      !Just to check input file is being read correctly.
      !*****************************************************************
      !****************************Parameters***************************
      
      Aa= 0
      Ba= 0
      h1a= 0
      h2a= 0

      !*****************************************************************
      !**********************Minimising Lamda1 (h1)*********************
      deltah= 0.0001
      h1b= 0.001
      h1a= 0
           DO 10
      chisqa= 0
              DO 20 i= 1,17   !  *** Is part of the line below supposed to be commented out?
                   Cs(i)= exp(-h1a*t(i))!*Aa !+ Ba*exp(-h2a*t(i))
                     Ctdiff= Ct - Cs  ! *** This is wrong. All three variables are arrays.
                       c= Ctdiff**2/eCt**2
                         chisqa= chisqa + c(i)
 20           END DO
      !Use initial h1 value of 0 to calculate start-point chisq.

      chisqb= 0
              DO 30 i= 1,17
                   h1b= h1b + deltah   ! *** In the line below, are you commenting out part of the line?
                     Cs(i)= exp(-h1b*t(i))!*Ab !+ Bb*exp(-h2b*t(i))
                       Ctdiff(i)= Ct(i) - Cs(i)
                         c= Ctdiff(i)**2/eCt(i)**2
                           chisqb= chisqb + c(i)
 30           END DO
      !First-step h1 used to find competing chisq for comparison.
                 WRITE(6,*) 'Chi-squared a=', chisqa,'for Lamda1=',h1a
                 WRITE(6,*) 'Chi-squared b=', chisqb,'for Lamda1=',h1b
      !Prints the two calculated chisq values to screen.
                 dchisq= chisqa - chisqb
              IF (dchisq.GT.0) THEN
                 h1a= h1b
                 ELSE IF (dchisq.LE.0) THEN
                 deltah= deltah-((deltah*2)/100)

              END IF
           IF (chisqb.LE.6618.681) EXIT
 10        END DO
           WRITE(6,*) 'Chi-squared is', chisqb,'for Lamda1=',h1b
      !*****************************************************************
           
      END PROGRAM assignment

Thanks for the help.
Nathan

EDIT: I should be getting a chi-squared of 6618.681 from this, but it's just stuck between 6921.866 and 6920.031. Help!

I added some comments where I think you're going wrong. My comments look like this: ! ***
In one comment, you have three assignment statements where you are using array variables, but without indexes. This is wrong. Everywhere you have an array variable in a calculation, it should have an index. E.g., it should be Ct(i), not just Ct.

In my other comments, you have a line of code that might be wrong. It looks like you commented it out, but I'm not sure what you did.

Your variable names are terrible! For example, Aa,Ab,Ba,Bb,h1a,h1b. Can't you come up with some names that are self-explanatory? Also, why does your code meander all across the page? This makes it harder to see which statements are inside IF blocks and DO loops and such.
 

Similar threads

  • · Replies 13 ·
Replies
13
Views
2K