2013-03-20 139 views
1

我試圖擬合一些數據,這些數據在上升高斯曲線之後的時間內分佈,然後以指數形式衰減。 我發現this example on the web,這和我的情況非常相似,但我剛開始適應python,這個例子似乎讓我感到困惑。 但是,我剛纔tryied的例子適應我的腳本和數據,並在下面是我的進步: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() 

我肯定混淆了整個事情,很抱歉提前爲這一點,但在這一點我不知道怎麼走得更遠... 代碼從網絡上獲取數據,所以它是直接可執行的。 目前代碼運行時沒有任何錯誤,但不會產生任何情節。 同樣,我的目標是將這些數據與這兩個函數進行匹配,我該如何改進我的代碼才能做到這一點? 任何建議真的很感激。

+4

的最好方式減少混淆是將問題分解成小部分。 首先,嘗試做一個簡單的情節 - 一條直線,或一個正弦波或其他東西。然後,找到一個繪圖函數,並確保通過首先嚐試一些簡單的線性函數或正弦波來理解它。使用你的經驗與繪圖,以確保你做得對。然後,只有這樣,才能嘗試適應一些複雜的東西。 歡迎您就所有這些步驟提出問題。通過將它分成小塊,您可以更輕鬆地幫助您。 – Robbert 2013-03-20 11:55:32

+0

你好,羅伯特,謝謝你的回覆。 我已經可以做一些擬合,如數據的第一個高斯部分。 我的問題是我不能「連接」功能。 此外,我在網上找到了這個例子,我正在嘗試應用它,但是我沒有其他參考資料可以讓我知道我正在做的事情是以一種好的方式或其他方式完成的。 – 2013-03-20 15:28:42

+0

爲什麼你會期望代碼產生一個情節? – silvado 2013-03-20 15:43:38

回答

1

類似於your other question,這裏也是我會用三角函數,以適應這個峯:

enter image description here

如果你的代碼後粘貼下面的代碼工作:

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