pyMCMC/pyMC를 사용하여 데이터/관측치에 비선형 함수를 맞 춥니 다.

Question

일부 데이터를 가우시안 (더 복잡한) 함수로 맞추려고합니다. 아래에 작은 예제를 만들었습니다.

내 첫 번째 질문은 일까요?

두 번째 질문은 입니다. x 방향, 즉 관측치/데이터의 x 위치에 오류를 어떻게 추가합니까?

pyMC에서 이러한 종류의 회귀를 수행하는 방법에 대한 훌륭한 가이드를 찾는 것은 매우 어렵습니다. 아마도 최소 자승법이나 비슷한 접근법을 사용하기 쉽기 때문에 결국에는 많은 매개 변수를 갖게되고 pyMC가이를 어떻게 제한하고 비교할 수 있는지 알아야합니다. pyMC는이를위한 좋은 선택 인 것 같습니다.

import pymc 
import numpy as np 
import matplotlib.pyplot as plt; plt.ion() 

x = np.arange(5,400,10)*1e3 

# Parameters for gaussian 
amp_true = 0.2 
size_true = 1.8 
ps_true = 0.1 

# Gaussian function 
gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps 
f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true) 

# add noise to the data points 
noise = np.random.normal(size=len(x)) * .02 
f = f_true + noise 
f_error = np.ones_like(f_true)*0.05*f.max() 

# define the model/function to be fitted. 
def model(x, f): 
    amp = pymc.Uniform('amp', 0.05, 0.4, value= 0.15) 
    size = pymc.Uniform('size', 0.5, 2.5, value= 1.0) 
    ps = pymc.Normal('ps', 0.13, 40, value=0.15) 

    @pymc.deterministic(plot=False) 
    def gauss(x=x, amp=amp, size=size, ps=ps): 
     e = -1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.)) 
     return amp*np.exp(e)+ps 
    y = pymc.Normal('y', mu=gauss, tau=1.0/f_error**2, value=f, observed=True) 
    return locals() 

MDL = pymc.MCMC(model(x,f)) 
MDL.sample(1e4) 

# extract and plot results 
y_min = MDL.stats()['gauss']['quantiles'][2.5] 
y_max = MDL.stats()['gauss']['quantiles'][97.5] 
y_fit = MDL.stats()['gauss']['mean'] 
plt.plot(x,f_true,'b', marker='None', ls='-', lw=1, label='True') 
plt.errorbar(x,f,yerr=f_error, color='r', marker='.', ls='None', label='Observed') 
plt.plot(x,y_fit,'k', marker='+', ls='None', ms=5, mew=2, label='Fit') 
plt.fill_between(x, y_min, y_max, color='0.5', alpha=0.5) 
plt.legend() 

필자는 더 많은 반복을 실행하고, 마지막으로 번인하고 묽게 사용해야 할 수도 있음을 알고 있습니다. 데이터와 적합성을 플로팅 한 그림이 아래에 나와 있습니다.

pymc.Matplot.plot (MDL) 수치 좋게 나타내는 분포 정점과 같다. 이거 좋지, 그렇지?

Answer 1

내 첫 번째 질문은, 내가 바로 그 일을하고 있습니까?

예! 당신은 알고있는 번 - 인 (burn-in) 기간을 포함시켜야합니다. 나는 샘플의 전반부를 버리고 싶다. 씬닝을 할 필요는 없지만 때로는 포스트 -MCMC가 처리하기가 더 빨라지고 저장할 때 작아집니다.

내가 추천하는 유일한 다른 것은 임의의 시드를 설정하여 결과가 "재현 가능"하도록하는 것입니다. np.random.seed(12345)이 트릭을 수행합니다.

아, 정말로 조언을 많이 주면 라고 말하면서 matplotlib 결과가 조금 더 아름답게 보일 것입니다.

두 번째 질문은 x 방향, 즉 관측치/데이터의 x 위치에 오류를 어떻게 추가합니까?

한 가지 방법은 각 오류에 대한 잠재 변수를 포함하는 것입니다. 이것은 당신의 예에서 효과적이지만 더 많은 관찰이 있다면 가능하지 않을 것입니다. 나는이 길을 시작하는 약간의 예를주지 : 여기

입니다 : 당신이 x 및 y에 잡음이있는 경우 좋은 답변을 얻기 어려울 수 있습니다처럼

# add noise to observed x values 
x_obs = pm.rnormal(mu=x, tau=(1e4)**-2) 

# define the model/function to be fitted. 
def model(x_obs, f): 
    amp = pm.Uniform('amp', 0.05, 0.4, value= 0.15) 
    size = pm.Uniform('size', 0.5, 2.5, value= 1.0) 
    ps = pm.Normal('ps', 0.13, 40, value=0.15) 

    x_pred = pm.Normal('x', mu=x_obs, tau=(1e4)**-2) # this allows error in x_obs 

    @pm.deterministic(plot=False) 
    def gauss(x=x_pred, amp=amp, size=size, ps=ps): 
     e = -1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.)) 
     return amp*np.exp(e)+ps 
    y = pm.Normal('y', mu=gauss, tau=1.0/f_error**2, value=f, observed=True) 
    return locals() 

MDL = pm.MCMC(model(x_obs, f)) 
MDL.use_step_method(pm.AdaptiveMetropolis, MDL.x_pred) # use AdaptiveMetropolis to "learn" how to step 
MDL.sample(200000, 100000, 10) # run chain longer since there are more dimensions 

같습니다 a notebook collecting this all up.

Answer 2

누군가가 그것을 필요로 할 때를 대비하여 Abraham Flaxman의 대답을 PyMC3으로 변환했습니다. 약간의 사소한 변화가 있지만 혼란 스러울 수 있습니다.

우선 결정 성 데코레이터 @Deterministic이 분포 형 호출 함수 var=pymc3.Deterministic()으로 대체되었습니다. 정규 분포 확률 변수의 벡터를 생성 할 때 둘째

rvs = pymc2.rnormal(mu=mu, tau=tau) 

가

rvs = pymc3.Normal('var_name', mu=mu, tau=tau,shape=size(var)).random() 

로 치환 전체 코드는 다음과 같다 :

결과

import numpy as np 
from pymc3 import * 
import matplotlib.pyplot as plt 

# set random seed for reproducibility 
np.random.seed(12345) 

x = np.arange(5,400,10)*1e3 

# Parameters for gaussian 
amp_true = 0.2 
size_true = 1.8 
ps_true = 0.1 

#Gaussian function 
gauss = lambda x,amp,size,ps: amp*np.exp(-1*(np.pi**2/(3600.*180.)*size*x)**2/(4.*np.log(2.)))+ps 
f_true = gauss(x=x,amp=amp_true, size=size_true, ps=ps_true) 

# add noise to the data points 
noise = np.random.normal(size=len(x)) * .02 
f = f_true + noise 
f_error = np.ones_like(f_true)*0.05*f.max() 

with Model() as model3: 
    amp = Uniform('amp', 0.05, 0.4, testval= 0.15) 
    size = Uniform('size', 0.5, 2.5, testval= 1.0) 
    ps = Normal('ps', 0.13, 40, testval=0.15) 

    gauss=Deterministic('gauss',amp*np.exp(-1*(np.pi**2*size*x/(3600.*180.))**2/(4.*np.log(2.)))+ps) 

    y =Normal('y', mu=gauss, tau=1.0/f_error**2, observed=f) 

    start=find_MAP() 
    step=NUTS() 
    trace=sample(2000,start=start) 

# extract and plot results 
y_min = np.percentile(trace.gauss,2.5,axis=0) 
y_max = np.percentile(trace.gauss,97.5,axis=0) 
y_fit = np.percentile(trace.gauss,50,axis=0) 
plt.plot(x,f_true,'b', marker='None', ls='-', lw=1, label='True') 
plt.errorbar(x,f,yerr=f_error, color='r', marker='.', ls='None', label='Observed') 
plt.plot(x,y_fit,'k', marker='+', ls='None', ms=5, mew=1, label='Fit') 
plt.fill_between(x, y_min, y_max, color='0.5', alpha=0.5) 
plt.legend() 

y_error

X의 오류에 대한 10

(변수에 'X'접미사 주) : 결과

# define the model/function to be fitted in PyMC3: 
with Model() as modelx: 

    x_obsx = pm3.Normal('x_obsx',mu=x, tau=(1e4)**-2,shape=40).random() 

    ampx = Uniform('ampx', 0.05, 0.4, testval= 0.15) 
    sizex = Uniform('sizex', 0.5, 2.5, testval= 1.0) 
    psx = Normal('psx', 0.13, 40, testval=0.15) 

    x_pred = Normal('x_pred', mu=x_obsx, tau=(1e4)**-2*np.ones_like(x_obsx),testval=5*np.ones_like(x_obsx),shape=40) # this allows error in x_obs 

    gauss=Deterministic('gauss',ampx*np.exp(-1*(np.pi**2*sizex*x_pred/(3600.*180.))**2/(4.*np.log(2.)))+psx) 

    y = Normal('y', mu=gauss, tau=1.0/f_error**2, observed=f) 

    start=find_MAP() 
    step=NUTS() 
    tracex=sample(20000,start=start) 

:

x_error_graph

마지막으로 관찰 할 때이다

traceplot(tracex[100:]) 
plt.tight_layout(); 

(결과가 표시되지 않음), 우리는 sizex이 '감쇠'또는 '회귀 dilu'로 고통 받고있는 것으로 보입니다의 측정 오류로 인해 '오류'가 발생했습니다.