指数衰减拟合

Question

我正在尝试对一些数据进行拟合，这些数据在时间上呈现指数衰减的分布。我在网上找了一些拟合的例子，但我的代码并没有成功拟合这些数据。结果只得到了一条直线。可能是初始参数设置有问题吧？到目前为止，我只用过高斯拟合和直线拟合，使用的方法可能不适合这种情况。 这段代码是从网上获取的数据，所以可以直接运行。 问题是：为什么代码没有得到任何拟合结果？ 非常感谢！

#!/usr/bin/env python

import pyfits, os, re, glob, sys
from scipy.optimize import leastsq
from numpy import *
from pylab import *
from scipy import *

rc('font',**{'family':'serif','serif':['Helvetica']})
rc('ps',usedistiller='xpdf')
rc('text', usetex=True)
#------------------------------------------------------

tmin = 56200
tmax = 56249

data=pyfits.open('http://heasarc.gsfc.nasa.gov/docs/swift/results/transients/weak/GX304-1.orbit.lc.fits')
time  = data[1].data.field(0)/86400. + data[1].header['MJDREFF'] + data[1].header['MJDREFI']
rate  = data[1].data.field(1)
error = data[1].data.field(2)
data.close()

cond = ((time > 56210) & (time < 56225))
time = time[cond]
rate = rate[cond]
error = error[cond]

right_exp = lambda p, x: p[0]*exp(-p[1]*x)
err = lambda p, x, y:(right_exp(p, x) -y)
v0= [0.20, 56210.0, 1]
out = leastsq(err, v0[:], args = (time, rate), maxfev=100000, full_output=1)
v = out[0] #fit parameters out
xxx = arange(min(time), max(time), time[1] - time[0])
ccc = right_exp(v, xxx)
fig = figure(figsize = (9, 9)) #make a plot
ax1 = fig.add_subplot(111)
ax1.plot(time, rate, 'g.') #spectrum
ax1.plot(xxx, ccc, 'b-') #fitted spectrum
savefig("right exp.png")

axis([tmin-10, tmax, -0.00, 0.45])

Answer 1

你的问题有点棘手，因为你的数组 times 里有一些很大的数字，当这些数字用在 exp(-a*time) 这个公式里时，结果会接近 0.，这就让 err 函数产生了误导，因为你的 rate 数组里也有一些接近 0. 的小值，导致错误很小。换句话说，指数函数里的高 a 值能给出一个好的结果。

要解决这个问题，你可以：

• 修改你的衰减函数，加入一个初始时间：
  exp(-a*(time-time0))
• 把你的输入数据改成从一个较小的数字开始：
  time -= time.min()

对于这两种方法，你都需要改变初始猜测 v0，比如说 v0=[0.,0.]。第一种解决方案看起来更稳妥，而且你不需要去处理 time 数组的变化。一个好的初始猜测对于 time0 是 time.min()：

right_exp = lambda p, x: p[0]*exp(-p[1]*(x-p[2]))
err = lambda p, x, y:(right_exp(p, x) -y)
v0= [0., 0., time.min() ]
out = leastsq(err, v0, args = (time, rate))
v = out[0] #fit parameters out
xxx = arange(min(time), max(time), time[1] - time[0])
ccc = right_exp(v, xxx)
fig = figure(figsize = (9, 9)) #make a plot
ax1 = fig.add_subplot(111)
ax1.plot(time, rate, 'g.') #spectrum
ax1.plot(xxx, ccc, 'b-') #fitted spectrum
fig.show()

这样就能得到：

不过，最终的结果还是依赖于 v0，比如说用 v0=[1.,1.,time.min()] 的话，衰减得太快，找不到最佳结果。