I need to implement a non-linear fitting algorithm in Fortran and chose to use MINPACK's flavor of the Levenberg-Marquardt algorithm as a basis for the least-squares stuff. However, I seem to misunderstand something very basic about how these routines are to be used. My understanding is that I have to call the routine lmdif1 and define a function fcn to be called by the MINPACK routines that calculates the square of the difference between the data to be fitted and the model function; lmdif1 would then return the best solution for the free fitting parameters in an array x.
I have written a small test program (given below) to try this out. The program produces a set of "data" by calculating values from the model function (exp(ax) in this case) with a given parameter a (=-0.5 in this case) and changing them by a small random "error". This dataset is then fitted with the same model function, with a being the free parameter. But it doesn't really work, the fits are usually very bad, unless I almost give the "correct" solution as a starting guess.
Considering that Levenberg-Marquardt is said to be quite robust and seeing how well it works e.g. in Gnuplot, I wonder what I'm doing wrong here. Any ideas? The code below should work when linking to the MINPACK routines.
Tom
! test lmdif
program testanf
implicit none
integer, parameter :: m=10,n=1,lwa=m*n+5*n+m
integer :: info,iwa(n)
double precision :: x(n),fvec(m),wa(lwa),tol=1d-3
external :: fcn
x(1)=-0.5d0
call lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa)
call translate_info(info,n,tol)
write(*,*) 'a=',x(1)
end program testanf
subroutine fcn(m,n,x,fvec,iflag)
implicit none
! m: number of data points
! n: number of fitting parameters
integer, intent(in) :: m,n,iflag
integer :: i
real, allocatable, save :: t(:),fdata(:)
real :: tmin,tmax,dt,rand
double precision, intent(inout) :: x(n)
double precision, intent(out) :: fvec(m)
character(len=1) :: c
logical, save :: init=.true.
if (init) then
! create "measured data"
tmin=0.
tmax=10.
dt=(tmax-tmin)/m
call random_seed(put=(/1,2,3,4,5,6,7,8/))
allocate(fdata(m),t(m))
do i=1,m
t(i)=i*dt
fdata(i)=fmodel(-0.5d0,t(i))
end do
tmax=0.1*maxval(fdata)
do i=1,m
call random_number(rand)
fdata(i)=fdata(i)+(rand-0.5)*tmax ! add measurement errors
write(10,*) t(i),fdata(i)
end do
init=.false.
end if
! minimize the residuals
do i=1,m
fvec(i)=(fdata(i)-fmodel(x(1),t(i)))**2
end do
if (iflag < 0) deallocate(t,fdata)
return
contains
! model function
real function fmodel(a,t)
implicit none
real, intent(in) :: t
double precision, intent(in) :: a
fmodel=exp(a*t)
return
end function fmodel
end subroutine fcn
fvec
, as opposed to just the difference. However, in terms of the quality of the result, that didn't help much - the crucial improvement in this respect seems to be to use double precision throughout, as in your example. If further tests I'll do confirm that, I'd find it quite stunning that it would make such a difference. $\endgroup$ – TomR Aug 5 '13 at 12:09