2

指数関数的減衰に続く時間に分布するいくつかのデータを当てはめようとしています。ウェブ上のいくつかの適切な例に従おうとしましたが、私のコードはデータに適合しません。直線のみがフィットから得られます。初期パラメータに何か問題があるのでしょうか？これまで、同じ方法を使用してガウスとライン フィットのみを使用してきましたが、この場合には正しくない可能性があります。コードは Web からデータを取得するため、直接実行できます。質問: コードが適合しないのはなぜですか? よろしくお願いします。

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

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

right_exp = lambda p, x: p[0]*exp(-p[1]*x)
err = lambda p, x, y:(right_exp(p, x) -y)
v0= [0.20, 56210.0, 1]
out = leastsq(err, v0[:], args = (time, rate), maxfev=100000, full_output=1)
v = out[0] #fit parameters out
xxx = arange(min(time), max(time), time[1] - time[0])
ccc = right_exp(v, xxx)
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("right exp.png")

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

1 に答える 1

3

で使用すると に近い値を与えるtimes大きな数値が配列に含まれているため、問題は悪条件です。これは、配列に に近い小さな値も含まれているため、関数をだまして小さなエラーにつながります。つまり、指数関数が高いほど、適切な解が得られます。exp(-a*time)0.errrate0.a

これを修正するには、次のことができます。

  • 初期時間を含むように減衰関数を変更します。
    exp(-a*(time-time0))
  • 入力データを小さい数値から開始するように変更します。
    time -= time.min()

v0どちらのオプションでも、最初の推測を変更する必要がありますv0=[0.,0.]。time最初のソリューションはより堅牢に見え、アレイの変更を管理する必要はありません。の適切な初期推定値time0は次のtime.min()とおりです。

right_exp = lambda p, x: p[0]*exp(-p[1]*(x-p[2]))
err = lambda p, x, y:(right_exp(p, x) -y)
v0= [0., 0., time.min() ]
out = leastsq(err, v0, args = (time, rate))
v = out[0] #fit parameters out
xxx = arange(min(time), max(time), time[1] - time[0])
ccc = right_exp(v, xxx)
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
fig.show()

与える:

ここに画像の説明を入力

それでも、最終的な結果は に依存しますv0。たとえば、v0=[1.,1.,time.min()]減衰が速すぎて最適なものが見つかりません。

于 2013-06-02T18:43:53.333 に答える