PyMC (2)를 사용한 강력한 비선형 회귀

Question

이 질문은 PyMC를 사용하여 비선형 회귀를 시도한다는 점에서 Fit a non-linear function to data/observations with pyMCMC/pyMC과 비슷합니다.

그러나 관찰 된 변수를 PyMC를 사용하여 비정규 분포 (즉, T 분포)를 따르는 방법을 아는 사람이 있는지 궁금합니다. 나는 그들이 T 분포를 포함하고 있다는 것을 알고 있지만, 나는 이것을 관찰 된 변수로 포함시키는 방법을 확신하지 못했습니다.

다음은 내가 문제에 직면 해있는 곳의 가짜 데이터를 사용하여 신속하게 데모 한 것입니다. 분명히 이상한 데이터 요소 중 일부를 보호하는 출력 분포를 사용하고 싶습니다.

import numpy as np 
import pymc as pm 
import matplotlib.pyplot as plt 

# For reproducibility 
np.random.seed(1234) 

x = np.linspace(0, 10*np.pi, num=150) 

# Set real parameters for the sinusoid 
true_freq = 0.9 
true_logamp = 1.2 
true_decay = 0.12 
true_phase = np.pi/4 


# Simulate the true trajectory 
y_real = (np.exp(true_logamp - true_decay*x) * 
      np.cos(true_freq*x + true_phase)) 

# Add some noise 
y_err = y_real + 0.05*y_real.max()*np.random.randn(len(x)) 


# Add some outliers 
num_outliers = 10 
outlier_locs = np.random.randint(0, len(x), num_outliers) 
y_err[outlier_locs] += (10 * y_real.max() * 
         (np.random.rand(num_outliers))) 



# Bayesian Regression 

def model(x, y, p0): 

    log_amp = pm.Normal('log_amp', mu=np.log(p0['amplitude']), 
         tau=1/(np.log(p0['amplitude']))) 

    decay = pm.Normal('decay', mu=p0['decay'], 
         tau=1/(p0['decay'])) 

    period = pm.TruncatedNormal('period', mu=p0['period'], 
           tau=1/(p0['period']), 
           a=1/(0.5/(np.median(np.diff(x)))), 
           b=x.max() - x.min()) 

    phase = pm.VonMises('phase', mu=p0['phase'], kappa=1.) 

    obs_tau = pm.Gamma('obs_tau', 0.1, 0.1) 

    @pm.deterministic(plot=False) 
    def decaying_sinusoid(x=x, log_amp=log_amp, decay=decay, 
          period=period, phase=phase): 

     return (np.exp(log_amp - decay*x) * 
       np.cos((2*np.pi/period)*x + phase)) 

    obs = pm.Normal('obs', mu=decaying_sinusoid, tau=obs_tau, value=y, 
        observed=True) 

    return locals() 

p0 = { 
    'amplitude' : 2.30185, 
    'decay'  : 0.06697, 
    'period' : 7.11672, 
    'phase'  : 0.93055, 
} 


MDL = pm.MCMC(model(x, y_err, p0)) 
MDL.sample(20000, 10000, 1) 

# Plot fit 
y_min = MDL.stats()['decaying_sinusoid']['quantiles'][2.5] 
y_max = MDL.stats()['decaying_sinusoid']['quantiles'][97.5] 
y_fit = MDL.stats()['decaying_sinusoid']['mean'] 
plt.plot(x, y_err, '.', label='Data') 
plt.plot(x, y_fit, label='Fit') 
plt.plot(x, y_real, label='True') 

plt.fill_between(x, y_min, y_max, color='0.5', alpha=0.5) 
plt.legend() 

감사합니다!

Answer 1

PyMC2는 pm.T의 t 분포를 가지고 있지만 0 인 중앙에 있으므로이 응용 프로그램에서 직접 사용할 수 없습니다. 대신, x 및 dof 매개 변수 nu에서 로그 확률을 계산하여 사용자 정의 관측 확률을 정의하는 pm.t_like(x, nu) 함수를 사용할 수 있습니다. 관찰 된 변수에 대한 이러한 맞춤 분포를 만드는 것은 간단합니다. 59-60 줄을 다음과 같이 변경하십시오.

@pm.observed 
def obs(mu=decaying_sinusoid, tau=obs_tau, value=y): 
    return pm.t_like(value-mu, tau)