指数衰减拟合答案

【问题标题】：exponential decay fitting指数衰减拟合
【发布时间】：2013-06-03 16:35:26
【问题描述】：

我正在尝试拟合一些在指数衰减后的时间内分布的数据。我试图在网上遵循一些合适的例子，但我的代码不适合数据。拟合只产生一条直线。也许初始参数有问题？到目前为止，我只使用了高斯拟合和线拟合，使用相同的方法，这对于这种情况可能不正确。该代码从网络获取数据，因此可以直接执行。问题：为什么代码不适合？非常感谢。

#!/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])

【问题讨论】：

将公式更改为p[0]*exp(-p[2]*(x-p[1]))
非常感谢埃弗特！该情节似乎导致corecct，但我收到消息：“exp中遇到溢出”。
溢出显然是因为您将大量数字插入指数（大约 56200）。为了获得最佳拟合结果，请始终将您的点缩放/移动到 1（x 和 y）的顺序，然后将结果缩放/移动回来。

标签： python numpy scipy

【解决方案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()] 衰减太快，找不到最佳值。

【讨论】：