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

因为，这样的话，情节太混乱了，我想去拟合曲线。 我试过用单变量pline，但这完全弄乱了我的数据。 有什么建议吗？ 我应该先定义一个函数来适应这些数据吗？ 我还寻找“最小平方”：这是解决这个问题的最佳方法吗？在

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。在

如果我想用不同的函数来拟合曲线的上升和下降部分呢？在

Answer 2

我用这个来装衣服。它是从网上某处改编的，但我忘了在哪儿。在

from __future__ import print_function
from __future__ import division
from __future__ import absolute_import

import numpy

from scipy.optimize.minpack import leastsq

### functions ###

def eq_cos(A, t):
    """
    4 parameters
    function: A[0] + A[1] * numpy.cos(2 * numpy.pi * A[2] * t + A[3])
    A[0]: offset
    A[1]: amplitude
    A[2]: frequency
    A[3]: phase
    """
    return A[0] + A[1] * numpy.cos(2 * numpy.pi * A[2] * t + numpy.pi*A[3])

def linear(A, t):
    """
    A[0]: y-offset
    A[1]: slope
    """
    return A[0] + A[1] * t  

### fitting routines ###

def minimize(A, t, y0, function):
    """
    Needed for fit
    """
    return y0 - function(A, t)

def fit(x_array, y_array, function, A_start):
    """
    Fit data

    20101209/RB: started
    20130131/RB: added example to doc-string

    INPUT:
    x_array: the array with time or something
    y-array: the array with the values that have to be fitted
    function: one of the functions, in the format as in the file "Equations"
    A_start: a starting point for the fitting

    OUTPUT:
    A_final: the final parameters of the fitting

    EXAMPLE:
    Fit some data to this function above
    def linear(A, t):
        return A[0] + A[1] * t  

    ### 
    x = x-axis
    y = some data
    A = [0,1] # initial guess
    A_final = fit(x, y, linear, A)
    ###

    WARNING:
    Always check the result, it might sometimes be sensitive to a good starting point.

    """
    param = (x_array, y_array, function)

    A_final, cov_x, infodict, mesg, ier = leastsq(minimize, A_start, args=param, full_output = True)

    return A_final



if __name__ == '__main__':

    # data
    x = numpy.arange(10)
    y = x + numpy.random.rand(10) # values between 0 and 1

    # initial guesss
    A = [0,0.5]

    # fit 
    A_final = fit(x, y, linear, A)

    # result is linear with a little offset
    print(A_final)