pymc 一次拟合多个观察值

Question

我正在尝试使用 pymc 来适应振荡数据的时间演变。在这里，每个时间步我不仅有一个点，而且还有几个。

我根本找不到在 pymc3 中进行这项工作的有效方法，因为它总是会引发一些输入值错误。所以我想知道是否有一个很好的解决方案。我附上了代码，但也可以在这里找到 ipython notebook here。

# coding: utf-8

# # Two level oscillation tests

# In[2]:

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import numpy.linalg as npl
from IPython.core.pylabtools import figsize
get_ipython().magic('matplotlib inline')


# create an artificial data set with several points per time value

# In[3]:

import pymc3 as pm
sigma =0.2;
Omega =0.5;
Nt = 20;
tmax =2;
Nrep = 5;
tlin = np.linspace(0,tmax,Nt);
t_1 = tlin[:];
t_2 = tlin[:];
n1_simu = np.sin(2*np.pi*Omega*tlin)**2;
n2_simu = 1 - n1_simu;

n1_noise = 0.2*np.random.randn(Nt);
n2_noise = 0.2*np.random.randn(Nt);

n1_exp = n1_simu+n1_noise;
n2_exp = n2_simu+n2_noise;

for jj in np.arange(Nrep):
    n1_noise = 0.2*np.random.randn(Nt); 
    n2_noise = 0.2*np.random.randn(Nt);
    n2_exp = np.r_[n2_exp, n2_simu+n2_noise]
    n1_exp = np.r_[n1_exp, n1_simu+n1_noise]
    t_1 = np.r_[t_1, tlin]
    t_2 = np.r_[t_2, tlin]

nt_exp = np.r_[n1_exp, n2_exp];
t_all = np.r_[t_1, t_2];
plt.figure(1)
plt.clf;
plt.plot(t_1,n1_exp, 'o');
plt.plot(t_1,n2_exp, 'o');
plt.xlabel('t')
plt.ylabel('population')


# Now that we have the simulated datas let us simulate them with pymc. 
# 
# The key is to put the mean value function into the Deterministic     symbol, then pymc unstands that it is supposed to be a variable.

# In[4]:

basic_model = pm.Model()


with basic_model:
    # Priors for unknown model parameters
    sigma = pm.HalfNormal('sigma', sd=1)
    Omega = pm.Normal('omega', mu=0.55, sd=0.1)
    amp = pm.Normal('Amplitude', mu=0.55, sd=0.1)


    # Expected value of outcome
    n1 = amp*pm.sin(2*np.pi*Omega*t_1)**2
    n2 = 1-n1
    Nval = len(nt_exp)
    Nswitch = len(n1_exp)
    idx = np.arange(Nval)
    if n1.shape:
        print(n1.shape)
        rate = pm.switch(Nswitch>= idx, np.r_[n1, n1], np.r_[n2, n2])

    # Likelihood (sampling distribution) of observations
    Y_obs = pm.Normal('Y_obs', mu=rate, sd=sigma, observed=nt_exp)


    # now sample it

    # In[ ]:

    Nsamples  =5000
    with basic_model:
    # obtain starting values via MAP
    start = pm.find_MAP()

    # instantiate sampler
    step = pm.NUTS(scaling=start)

    # draw 500 posterior samples
    trace = pm.sample(Nsamples, step, start=start)

Answer 1

关键是为你拥有的每个 x 值输入一个 y 值。如果每个 x 值有多个点，请放入多个 x,y 对。

我在这里修复的模型还有一些其他问题：https://gist.github.com/3ab2cba82e3cb5291e38

希望我能正确理解您的尝试。本质上这是一个受https://github.com/pymc-devs/pymc3/blob/master/pymc3/examples/gaussian_mixture_model.ipynb启发的混合模型。