2013-03-21 1 views
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 

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]) 
+2

수식을'p [0] * exp (-p [2] * (x-p [1]))' – Evert

+0

으로 바꿔 줘서 고마워! 플롯이 corecct를 생성하는 것으로 보입니다. 그러나 "exp에서 오버 플로우가 발생했습니다"라는 메시지가 나타납니다. –

+1

오버 플로우는 분명히 큰 숫자를 지수로 삽입하기 때문에 발생합니다 (약 56200). 최적의 결과를 얻으려면 점을 항상 1 (x 및 y)의 순서로 축척/시프트 한 다음 결과를 다시 축척/이동하십시오. – Evert

답변

2
귀하의 문제는 병이 배열 timesexp(-a*time)에서 사용하는 경우에지도하는 err 기능이 rate 배열 0.에 또한 가까운 작은 값을 포함하고 있기 때문에 어떤 트릭, 0.에 가까운 값을주는 것을 큰 숫자가 포함되어 있기 때문에 조절한다

작은 오류. 즉, 지수 함수로 높은 a이 좋은 해결책을 제공합니다.

를 해결하려면 그 작업을 수행 할 수 있습니다

  • 변화 초기 시간을 포함하도록 붕괴 기능 :
    exp(-a*(time-time0))
  • 변경하여 입력 데이터는 적은 수에서 시작 :
    time -= time.min()

두 옵션 모두 초기 추정을 변경해야합니다. v0v0=[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() 

기부 :

enter image description here

아직도, 최종 결과는 v0에 따라된다, 예를 들면, v0=[1.,1.,time.min()]으로 너무 빠르며 최적을 찾지 못합니다.