Python:如何拟合曲线

Question

我有以下代码来针对三组不同的时间范围绘制三组数据(计数率与时间):

#!/usr/bin/env python

from pylab import rc, array, subplot, zeros, savefig, ylim, xlabel, ylabel, errorbar, FormatStrFormatter, gca, axis
from scipy import optimize, stats
import numpy as np
import pyfits, os, re, glob, sys

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 > tmin-5) & (time < tmax))
time=time[cond]
rate=rate[cond]
error=error[cond]

errorbar(time, rate, error, fmt='r.', capsize=0)
gca().xaxis.set_major_formatter(FormatStrFormatter('%5.1f'))

axis([tmin-10,tmax,-0.00,0.45])
xlabel('Time, MJD')
savefig("sync.eps",orientation='portrait',papertype='a4',format='eps')

因为这样的话剧情太困惑了，我就想拟合一下曲线。我尝试使用 UnivariateSpline，但这完全弄乱了我的数据。请问有什么建议吗？我应该首先定义一个函数来拟合这些数据吗？我还寻找“最小二乘”:这是解决这个问题的最佳方案吗？

Answer 1

这就是我解决的方法:

#!/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
pi = 3.14
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 > tmin-5) & (time < tmax))
time=time[cond]
rate=rate[cond]
error=error[cond]

gauss_fit = lambda p, x: p[0]*(1/(2*pi*(p[2]**2))**(1/2))*exp(-(x-p[1])**2/(2*p[2]**2))+p[3]*(1/sqrt(2*pi*(p[5]**2)))*exp(-(x-p[4])**2/(2*p[5]**2)) #1d Gaussian func
e_gauss_fit = lambda p, x, y: (gauss_fit(p, x) -y) #1d Gaussian fit
v0= [0.20, 56210.0, 1, 0.40, 56234.0, 1] #inital guesses for Gaussian Fit, just do it around the peaks
out = leastsq(e_gauss_fit, v0[:], args=(time, rate), maxfev=100000, full_output=1) #Gauss Fit
v = out[0] #fit parameters out
xxx = arange(min(time), max(time), time[1] - time[0])
ccc = gauss_fit(v, xxx) # this will only work if the units are pixel and not wavelength
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("plotfitting.png")

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

来自here .

如果我想适应曲线的上升和衰减部分的不同函数怎么办？