Python: fit data with gaussian rising and exponential decay

Question

I am trying to fit some data that are distributed in the time following a rising gaussian curve, and then decaying exponentially. I have found this example on the web, that is very similar to my case, but I just started to fit with python, and the example seems quite confusing to me. Nonetheless, I have tryied to adapt the example to my script and data, and in the following is my progress:

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

I am sure to have confused the whole thing, so sorry in advance for that, but at this point I don't know how to go further... The code take the data from the web, so it is directly executable. At the moment the code runs without any error, but it doesn't produce any plot. Again, my goal is to fit the data with those two functions, how can I improve my code to do that? Any suggestion is really appreciated.

The best way to reduce confusion is to split the problem in small parts. First, try to make a simple plot - a straight line, or a sine wave or something. Then, find a plotting function and make sure you understand it by trying it first on something simple - a linear function or a sine wave. Use your experience with plotting to make sure you do it correctly. Then, and only then, try to make a fit of something complicated. You're welcome to ask questions for all these steps. By breaking it down in small bits you make it easier for us to help you. — Robbert, Mar 20 '13 at 11:55
Hello, Robbert, thanks for the reply. I can already do some fits, like the first gaussian part of the data. My problem is that I can not "connect" the functions. Furthermore, I found that example on the web, and I am trying to apply it, but I have no other references to get an idea of what I am doing is done in a good way or whatever. — Py-ser, Mar 20 '13 at 15:28
@Py-ser If I'm reading your code correctly, it doesn't seem to get to any of the plot methods... — silvado, Mar 20 '13 at 20:33

score 1 · Answer 1 · edited May 23 '17 at 10:25

Similarly to your other question, here also I would use a trigonometric function to fit this peaK:

enter image description here

The following code works if pasted after your code:

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()

Python: fit data with gaussian rising and exponential decay

1 Answers1