使用MINPACK进行曲线拟合:实现?

计算科学 优化 最小二乘
2021-12-01 08:31:38

我需要在 Fortran 中实现一个非线性拟合算法,并选择使用 MinPACK 的 Levenberg-Marquardt 算法作为最小二乘法的基础。但是,我似乎误解了有关如何使用这些例程的一些非常基本的东西。我的理解是我必须调用例程 lmdif1 并定义一个函数 fcn 由MINPACK例程调用,该函数计算要拟合的数据与模型函数之间的差异的平方;然后 lmdif1 将返回数组 x 中自由拟合参数的最佳解。

我编写了一个小测试程序(如下所示)来尝试一下。该程序通过使用给定参数a(在这种情况下=-0.5)计算模型函数(在这种情况下为exp(ax))的值并通过一个小的随机“错误”来改变它们来生成一组“数据”。然后用相同的模型函数对该数据集进行拟合,其中 a 是自由参数。但它并没有真正起作用,拟合通常非常糟糕,除非我几乎给出“正确”的解决方案作为开始猜测。

考虑到据说 Levenberg-Marquardt 非常健壮并且看到它在 Gnuplot 中的工作情况如何,我想知道我在这里做错了什么。有任何想法吗?当链接到MINPACK 例程时,下面的代码应该可以工作。

汤姆

! 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
0个回答
没有发现任何回复~