如何在pymc3中创建容器

Question

我正在尝试为 Langevin 方程（谐波势中的布朗粒子）的特定结果的似然函数建立模型：

这是我在 pymc2 中似乎有效的模型： https://github.com/hstrey/BayesianAnalysis/blob/master/Langevin%20simulation.ipynb

#define the model/function to be fitted.
def model(x): 
    t = pm.Uniform('t', 0.1, 20, value=2.0)
    A = pm.Uniform('A', 0.1, 10, value=1.0)

    @pm.deterministic(plot=False)
    def S(t=t):
       return 1-np.exp(-4*delta_t/t)

    @pm.deterministic(plot=False)
    def s(t=t):
       return np.exp(-2*delta_t/t)

    path = np.empty(N, dtype=object)

    path[0]=pm.Normal('path_0',mu=0, tau=1/A, value=x[0], observed=True)
    for i in range(1,N):
        path[i] = pm.Normal('path_%i' % i,
                        mu=path[i-1]*s,
                        tau=1/A/S,
                        value=x[i],
                        observed=True)
        return locals()

mcmc = pm.MCMC( model(x) )
mcmc.sample( 20000, 2000, 10 )

基本思想是每个点都依赖于链（马尔可夫链）中的前一个点。顺便说一句，x是一个数据数组，N是它的长度，delta_t是时间步长=0.01。知道如何在 pymc3 中实现这个吗？我试过了：

# define the model/function for diffusion in a harmonic potential
DHP_model = pm.Model()
with DHP_model:
    t = pm.Uniform('t', 0.1, 20)
    A = pm.Uniform('A', 0.1, 10)

    S=1-pm.exp(-4*delta_t/t)

    s=pm.exp(-2*delta_t/t)

    path = np.empty(N, dtype=object)

    path[0]=pm.Normal('path_0',mu=0, tau=1/A, observed=x[0])
    for i in range(1,N):
        path[i] = pm.Normal('path_%i' % i,
                        mu=path[i-1]*s,
                        tau=1/A/S,
                        observed=x[i])

不幸的是，我一尝试运行模型就崩溃了。我在我的机器上尝试了一些 pymc3 示例（教程），这很有效。

提前致谢。我真的希望 pymc3 中的新采样器能帮助我处理这个模型。我正在尝试将贝叶斯方法应用于单分子实验。

Answer 1

因为我没有看到我的问题的答案，所以让我自己回答。我提出了以下解决方案：

# now lets model this data using pymc
# define the model/function for diffusion in a harmonic potential
DHP_model = pm.Model()
with DHP_model:
    D = pm.Gamma('D',mu=mu_D,sd=sd_D)
    A = pm.Gamma('A',mu=mu_A,sd=sd_A)

    S=1.0-pm.exp(-2.0*delta_t*D/A)

    ss=pm.exp(-delta_t*D/A)

    path=pm.Normal('path_0',mu=0.0, tau=1/A, observed=x[0])
    for i in range(1,N):
        path = pm.Normal('path_%i' % i,
                        mu=path*ss,
                        tau=1.0/A/S,
                        observed=x[i])

    start = pm.find_MAP()
    print(start)
    trace = pm.sample(100000, start=start)

不幸的是，此代码需要 6 小时到 2 天的时间来编译 N=50。我运行正在一台速度相当快的 PC (24Gb RAM) 运行正在 Ubuntu。我尝试使用 GPU，但运行速度稍慢。我怀疑内存问题，因为它在运行ning 时使用了 99.8% 的内存。我用 Stan 尝试了相同的计算，只需 2 分钟即可达到运行.

Answer 2

无需在循环中创建许多单独的正态分布的一维变量，您可以创建自定义分布（通过扩展 Continuous），该分布知道用于计算整个路径的对数似然的公式。您可以 bootstrap 这个似然公式脱离 pymc3 已经知道的正态似然公式。有关示例，请参阅 the built-in AR1 class。

因为你的粒子遵循马尔可夫属性，你的可能性看起来像

import theano.tensor as T

def logp(path):
    now = path[1:]
    prev = path[:-1]

    loglik_first = pm.Normal.dist(mu=0., tau=1./A).logp(path[0])
    loglik_rest = T.sum(pm.Normal.dist(mu=prev*ss, tau=1./A/S).logp(now))
    loglik_final = loglik_first + loglik_rest

    return loglik_final

我猜您想在每个时间步都为 ss 绘制一个值，在这种情况下您应该确保指定 ss = pm.exp(..., shape=len(x)-1)，以便 prev*ss 在上面的块被解释为逐元素乘法。

然后你可以用

指定你的观察结果

path = MyLangevin('path', ..., observed=x)

这应该运行快很多。

如何在pymc3中创建容器

How to create a container in pymc3

pymc3