如何决定在 PyMC3 中将哪些先验分布用于参数？
How to decide on what priors distributions to use for parameters in PyMC3?
我正在研究 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')
附带问题
我也一直在研究在使用除高斯分布以外的不同分布时对适应性测试和卡尔曼滤波器实施优度的想法，所以，如果你有时间，我将不胜感激，如果你能看看它们？。这两个问题都可以在这里找到：
拟合优度检验link：Goodness-to-fit test
卡尔曼滤波器link：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 组合在一个图中。
第一个错误：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')
我正在研究 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')
附带问题
我也一直在研究在使用除高斯分布以外的不同分布时对适应性测试和卡尔曼滤波器实施优度的想法，所以，如果你有时间，我将不胜感激，如果你能看看它们？。这两个问题都可以在这里找到：
拟合优度检验link：Goodness-to-fit test
卡尔曼滤波器link：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 组合在一个图中。
第一个错误：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')