!******************************************************* !* Program to demonstrate the Parafit subroutine * !* --------------------------------------------------- * !* Reference: BASIC Scientific Subroutines, Vol. II * !* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * !* * !* F90 version by J-P Moreau, Paris * !* (www.jpmoreau.fr) * !* --------------------------------------------------- * !* SAMPLE RUN: * !* * !* Parametric least squares fit * !* * !* The input data are: * !* * !* X( 1) = 1.000000 Y( 1) = 0.033702 * !* X( 2) = 2.000000 Y( 2) = 0.249029 * !* X( 3) = 3.000000 Y( 3) = 0.944733 * !* X( 4) = 4.000000 Y( 4) = 1.840089 * !* X( 5) = 5.000000 Y( 5) = 1.840089 * !* X( 6) = 6.000000 Y( 6) = 0.944733 * !* X( 7) = 7.000000 Y( 7) = 0.249029 * !* X( 8) = 8.000000 Y( 8) = 0.033702 * !* X( 9) = 9.000000 Y( 9) = 0.002342 * !* X(10) = 10.000000 Y(10) = 0.000084 * !* * !* The coefficients are: * !* 2.000000000000000 * !* 4.500000000000000 * !* 3.000000000000000 * !* * !* The standard deviation of the fit is 0.000000 * !* * !* The number of iterations was 49 * !******************************************************* PROGRAM Parafit parameter(LMAX=10,SIZE=25) real*8 A(LMAX),X(SIZE),Y(SIZE); real*8 d,e,ee1 integer i,l,m,n n = 10; l = 3; m = 0; d = 0.d0 print *, 'PARAMETRIC LEAST SQUARES FIT' print *, ' ' print *, 'The input data are:' print *, ' ' do i = 1, n X(i) = dfloat(i) Y(i) = 2.d0 * dexp(-(X(i) - 4.5d0) * (X(i) - 4.5d0) / 3.d0); write(*,50) i,X(i),i,Y(i) end do e = 0.1d0; ee1 = 0.5d0; A(1) = 10.d0; A(2) = 10.d0; A(3) = 10.d0 call Param_LS(l,m,n,d,e,ee1,A,X,Y) print *, ' ' print *,'The coefficients are:' print *, A(1) print *, A(2) print *, A(3) print *, ' ' write (*,60) d write (*,70) m print *, ' ' stop 50 format(' X(',I2,') = ',F9.6,' Y(',I2,') = ',F9.6) 60 format(' The standard deviation of the fit is ',F10.6) 70 format(' The number of iterations was ',I2) End !main program !Function subroutine Subroutine S500(l,x,y,A) real*8 x,y,A(l) y = A(1) * dexp(-(x - A(2)) * (x - A(2)) / A(3)) return end !Residual generation subroutine Subroutine S200(l,l2,n,d,A,X,Y) parameter(SIZE=25) real*8 d,l2,A(l),X(SIZE),Y(SIZE) integer l,n,j real*8 xx,yy l2 = 0.d0 do j = 1, n xx = X(j) ! Obtain function call S500(l,xx,yy,A) l2 = l2 + (Y(j) - yy) * (Y(j) - yy) end do d = dsqrt(l2 / dfloat(n - l)) return end !*************************************************************** !* Parametric least squares curve fit subroutine. This program * !* least squares fits a function to a set of data values by * !* successively reducing the variance. Convergence depends on * !* the initial values and is not assured. * !* n pairs of data values, X(i), Y(i), are given. There are l * !* parameters, A(j), to be optimized across. * !* Required are initial values for the A(l) and e. Another * !* important parameter which affects stability is e1, which is * !* initially converted to e1(l)for the first intervals. * !* The parameters are multiplied by (1 - e1(i)) on each pass. * !*************************************************************** Subroutine Param_LS(l,m,n,d,e,ee1,A,X,Y) !Labels: 50,100 parameter(SIZE=25) real*8 d,e,ee1,A(l),X(SIZE),Y(SIZE) real*8 E1(l) integer i,l,m,n real*8 a0,l1,l2,m0,m1 do i = 1, l E1(i) = ee1 end do ; !Set up test residual l1 = 1.d6 !Make sweep through all parameters 50 do i = 1, l a0 = A(i) !Get value of residual A(i) = a0 100 call S200(l,l2,n,d,A,X,Y) !Store result in m0 m0 = l2 !Repeat for m1 A(i) = a0 * (1.d0 - E1(i)) call S200(l,l2,n,d,A,X,Y) m1 = l2 !Change interval size if called for !If variance was increased, halve E1(i) if (m1 > m0) then E1(i) = -E1(i) / 2.d0 end if !If variance was reduced, increase step size by increasing E1(i) if (m1 < m0) then E1(i) = 1.2d0 * E1(i) end if !If variance was increased, try to reduce it if (m1 > m0) then A(i) = a0 end if if (m1 > m0) then goto 100 end if end do !i loop !End of a complete pass !Test for convergence m = m + 1 if (l2.eq.0.d0) then return end if if (dabs((l1 - l2) / l2) < e) then return end if !If this point is reached, another pass is called for l1 = l2 goto 50 end ! Param_LS() ! End file parafit.f90