如何决定 PyMC3 中的参数使用哪些先验分布？

Question

我正在研究 PyMC3 包，我有兴趣在我有多个不同信号且每个信号具有不同幅度的场景中实现该包。

但是，我不知道需要使用哪种类型的先验才能在其中实现 PyMC3 以及实现可能性分布。场景示例如下图所示：

我尝试在这里实现它，但是，每次我都不断收到错误：

pymc3.exceptions.SamplingError: Bad initial energy

我的代码

## Signal 1:
    with pm.Model() as model:
        # Parameters:
        # Prior Distributions:
        # BoundedNormal = pm.Bound(pm.Exponential, lower=0.0, upper=np.inf)
        # c = BoundedNormal('c', lam=10)
        # c = pm.Uniform('c', lower=0, upper=300)
        alpha = pm.Normal('alpha', mu = 0, sd = 10)
        beta = pm.Normal('beta', mu = 0, sd = 1)
        sigma = pm.HalfNormal('sigma', sd = 1)
        mu = pm.Normal('mu', mu=0, sigma=1)
        sd = pm.HalfNormal('sd', sigma=1)

        # Observed data is from a Multinomial distribution:
        # Likelihood distributions:
        # bradford = pm.DensityDist('observed_data', logp=bradford_logp, observed=dict(value=S1, loc=mu, scale=sd, c=c))
        # observed_data = pm.Beta('observed_data', mu=mu, sd=sd, observed=S1)
        observed_data = pm.Beta('observed_data', alpha=alpha, beta=beta, mu=mu, sd=sd, observed=S1)

    with model:
        # obtain starting values via MAP
        startvals = pm.find_MAP(model=model)

        # instantiate sampler
        # step = pm.Metropolis()
        step = pm.HamiltonianMC()
        # step = pm.NUTS()

        # draw 5000 posterior samples
        trace = pm.sample(start=startvals, draws=1000, step=step, tune=500, chains=4, cores=1, discard_tuned_samples=True)

        # Obtaining Posterior Predictive Sampling:
        post_pred = pm.sample_posterior_predictive(trace, samples=500)
        print(post_pred['observed_data'].shape)

    plt.title('Trace Plot of Signal 1')
    pm.traceplot(trace, var_names=['mu', 'sd'], divergences=None, combined=True)
    plt.show(block=False)
    plt.pause(5)  # Pauses the program for 5 seconds
    plt.close('all')

    pm.plot_posterior(trace, var_names=['mu', 'sd'])
    plt.title('Posterior Plot of Signal 1')
    plt.show(block=False)
    plt.pause(5)  # Pauses the program for 5 seconds
    plt.close('all')

附加问题

我也一直在研究在使用除高斯分布以外的其他分布时实现适应性测试和卡尔曼滤波器的优点，所以，如果你有时间，如果你能看看它们，我将不胜感激。这两个问题都可以在这里找到：

拟合优度测试链接：Goodness-to-fit test

卡尔曼滤波器链接：Kalman Filter

---

编辑 1

假设我有大约 5 个信号并且想要实现贝叶斯接口以查看信号 PDF 的差异。我该如何解决这个问题？我是否需要创建多个模型并获得它们的后验分布？就像这张图片：

如果我需要得到后验分布，是否使用下面的代码？

# Obtaining Posterior Predictive Sampling:
post_pred = pm.sample_posterior_predictive(trace, samples=500)

---

编辑 2

如果我有多个信号，我可以通过这种方式实现它以查看所有信号中alpha 和beta 的变化吗？

        observed_data_S1 = pm.Beta('observed_data_S1', alpha=alpha[0], beta=beta[0], observed=S1[0])
        observed_data_S2 = pm.Beta('observed_data_S2', alpha=alpha[1], beta=beta[1], observed=S2[0])
        observed_data_S3 = pm.Beta('observed_data_S3', alpha=alpha[2], beta=beta[2], observed=S3[0])
        observed_data_S4 = pm.Beta('observed_data_S4', alpha=alpha[3], beta=beta[3], observed=S4[0])
        observed_data_S5 = pm.Beta('observed_data_S5', alpha=alpha[4], beta=beta[4], observed=S5[0])
        observed_data_S6 = pm.Beta('observed_data_S6', alpha=alpha[5], beta=beta[5], observed=S6[0])

编辑 3：

如何在一张图中绘制多条轨迹？因为我正在查看多个信号并考虑将所有 alpha 和 beta 组合在一个图中。

Answer 1

第一个错误：Beta 分布的参数alpha 和beta 必须为正数。您在它们上使用了 Normal 先验，它允许 RV 取负值和 0 值。您可以通过在 pm.Normal 分发版上使用 pm.Bound 或使用 pm.HalfNormal 分发版来轻松解决此问题。

第二个错误：另一个不一致的地方是指定mu 和sigma 以及alpha 和beta 参数。 Beta 接受 mu 和 sigma 或 alpha 和 beta 但不能同时接受。默认行为是使用alpha 和beta 参数而不是mu 和sigma 参数。您在推断 mu 和 sigma 时浪费了大量计算能力。

其他评论：您不应该在 3.8 版以后的任何发行版中使用 sd 参数，因为它已被弃用，并将在 3.9 中删除。请改用sigma。

修正版：

import numpy as np
import theano
import theano.tensor as tt
import pymc3 as pm
import matplotlib.pyplot as plt

S1 = np.random.rand(10)

## Signal 1:
with pm.Model() as model:
    # Parameters:
    # Prior Distributions:
    # BoundedNormal = pm.Bound(pm.Exponential, lower=0.0, upper=np.inf)
    # c = BoundedNormal('c', lam=10)
    # c = pm.Uniform('c', lower=0, upper=300)
    alpha = pm.HalfNormal('alpha', sigma=10)
    beta = pm.HalfNormal('beta', sigma=1)

    # Observed data is from a Multinomial distribution:
    # Likelihood distributions:
    # bradford = pm.DensityDist('observed_data', logp=bradford_logp, observed=dict(value=S1, loc=mu, scale=sd, c=c))
    # observed_data = pm.Beta('observed_data', mu=mu, sd=sd, observed=S1)
    observed_data = pm.Beta('observed_data', alpha=alpha, beta=beta, observed=S1)

with model:
    # obtain starting values via MAP
    startvals = pm.find_MAP(model=model)

    # instantiate sampler
    # step = pm.Metropolis()
    step = pm.HamiltonianMC()
    # step = pm.NUTS()

    # draw 5000 posterior samples
    trace = pm.sample(start=startvals, draws=1000, step=step, tune=500, chains=4, cores=1, discard_tuned_samples=True)

    # Obtaining Posterior Predictive Sampling:
    post_pred = pm.sample_posterior_predictive(trace, samples=500)
    print(post_pred['observed_data'].shape)

plt.title('Trace Plot of Signal 1')
pm.traceplot(trace, var_names=['alpha', 'beta'], divergences=None, combined=True)
plt.show(block=False)
plt.pause(5)  # Pauses the program for 5 seconds
plt.close('all')

pm.plot_posterior(trace, var_names=['alpha', 'beta'])
plt.title('Posterior Plot of Signal 1')
plt.show(block=False)
plt.pause(5)  # Pauses the program for 5 seconds
plt.close('all')