Python: 用高斯上升和指数衰减拟合数据

Question

我正在尝试处理一些数据，这些数据在时间上呈现出一个先上升的高斯曲线，然后再指数衰减。我在网上找到了一些与我的情况非常相似的例子，这个例子，但我刚开始用Python进行拟合，觉得这个例子有点复杂。尽管如此，我还是尝试把这个例子调整到我的脚本和数据上，下面是我的进展：

#!/usr/bin/env python

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

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 > 56200) & (time < 56220))
time=time[cond]
rate=rate[cond]
error=error[cond]

def expGauss(x, pos, wid, tConst, expMod = 0.5, amp = 1):
    expMod *= 1.0
    gNorm = amp * N.exp(-0.5*((x-pos)/(wid))**2)
    g = expBroaden(gNorm, tConst, expMod)
    return g, gNorm

def expBroaden(y, t, expMod):
    fy = F.fft(y)
    a = N.exp(-1*expMod*time/t)
    fa = F.fft(a)
    fy1 = fy*fa
    yb = (F.ifft(fy1).real)/N.sum(a)
    return yb

if __name__ == '__main__':

# Fit the first set
#p[0] -- amplitude, p[1] -- position, p[2] -- width
    fitfuncG = lambda p, x: p[0]*N.exp(-0.5*(x-p[1])**2/p[2]**2) # Target function
    errfuncG = lambda p, x, y: fitfuncG(p, x) - y # Distance to the target function
    p0 = [0.20, 56210, 2.0] # Initial guess for the parameters
    p1, success = optimize.leastsq(errfuncG, p0[:], args=(time, rate))
    p1G = fitfuncG(p1, time)
   # P.plot(rate,  'ro', alpha = 0.4, label = "Gaussian")
   # P.plot(p1G, label = 'G-Fit')

def expGauss(x, pos, wid, tConst, expMod = 0.5, amp = 1):
  #p[0] -- amplitude, p[1] -- position, p[2] -- width, p[3]--tConst, p[4] -- expMod  
    fitfuncExpG = lambda p, x: expGauss(x, p[1], p[2], p[3], p[4], p[0])[0]
    errfuncExpG = lambda p, x, y: fitfuncExpG(p, x) - y # Distance to the target function  
    p0a = [0.20, 56210, 2.0] # Initial guess for the parameters
    p1a, success = optimize.leastsq(errfuncExpG, p0a[:], args=(time, rate))
    p1aG = fitfuncExpG(p1a, time)
    print type(rate), type(time), len(rate), len(time)  
    P.plot(rate, 'go', alpha = 0.4, label = "ExpGaussian")  
    P.plot(p1aG, label = 'ExpG-Fit')

    P.legend()  
    P.show()

我知道我可能搞混了整个过程，所以提前跟你们说声抱歉，但在这个阶段我不知道该怎么继续……这段代码从网上获取数据，所以可以直接运行。目前这段代码没有报错，但也没有生成任何图表。再次强调，我的目标是用这两个函数来拟合数据，我该如何改进我的代码来实现这个目标呢？任何建议都非常感谢。

Answer 1

和你之前的问题类似，这里我也会用一个三角函数来拟合这个峰值：

下面的代码可以在你的代码后面粘贴使用：

import numpy as np
from scipy.optimize import curve_fit
x = time
den = x.max() - x.min()
x -= x.min()
y_points = rate
def func(x, a1, a2, a3):
    return  a1*sin(1*pi*x/den)+\
            a2*sin(2*pi*x/den)+\
            a3*sin(3*pi*x/den)
popt, pcov = curve_fit(func, x, y_points)
y = func(x, *popt)
plot(time,rate)
plot(x,y, color='r', linewidth=2.)
show()