下面的代码使用beta作为优先项来预测二项分布的p。有些时候,我得到了毫无意义的结果(接受率= 0)。当我用pymc3编写相同的逻辑时,我没有问题。我看不出我在这里错过了什么。
import numpy as np
import tensorflow as tf
import tensorflow_probability as tfp
import edward2 as ed
from pymc3.stats import hpd
import numpy as np
import seaborn
import matplotlib.pyplot as plt
p_true = .15
N = [10, 100, 1000]
successN = np.random.binomial(p=p_true, n=N)
print(N)
print(successN)
def beta_binomial(N):
p = ed.Beta(
concentration1=tf.ones( len(N) ),
concentration0=tf.ones( len(N) ),
name='p'
)
return ed.Binomial(total_count=N, probs=p, name='obs')
log_joint = ed.make_log_joint_fn(beta_binomial)
def target_log_prob_fn(p):
return log_joint(N=N, p=p, obs=successN)
#kernel = tfp.mcmc.HamiltonianMonteCarlo(
# target_log_prob_fn=target_log_prob_fn,
# step_size=0.01,
# num_leapfrog_steps=5)
kernel = tfp.mcmc.NoUTurnSampler(
target_log_prob_fn=target_log_prob_fn,
step_size=.01
)
trace, kernel_results = tfp.mcmc.sample_chain(
num_results=1000,
kernel=kernel,
num_burnin_steps=500,
current_state=[
tf.random.uniform(( len(N) ,))
],
trace_fn=(lambda current_state, kernel_results: kernel_results),
return_final_kernel_results=False)
p, = trace
p = p.numpy()
print(p.shape)
print('acceptance rate ', np.mean(kernel_results.is_accepted))
def printSummary(name, v):
print(name, v.shape)
print(np.mean(v, axis=0))
print(hpd(v))
printSummary('p', p)
for data in p.T:
print(data.shape)
seaborn.distplot(data, kde=False)
plt.savefig('p.png')图书馆:
pip install -U pip
pip install -e git+https://github.com/google/edward2.git@4a8ed9f5b1dac0190867c48e816168f9f28b5129#egg=edward2
pip install https://storage.googleapis.com/tensorflow/linux/cpu/tensorflow-2.0.0-cp37-cp37m-manylinux2010_x86_64.whl#egg=tensorflow
pip install tensorflow-probability有时,我看到以下情况(当接受rate=0时):

而且,有时我看到以下内容(当接受rate>.9时):

发布于 2019-12-28 03:31:38
当我在贝叶斯推理中得到不稳定的结果时(我使用mis,但它也使用螺母),这通常是因为先验和可能性都是错误指定的,或者超参数对这个问题没有好处。
第一个图表显示,采样器从未偏离对答案的最初猜测(因此接受率为0)。它也让我担心绿色分布似乎是正确的0。β(1,1)在0处有正概率,但p=0可能是一个不稳定的解?(就像在这里一样,采样器可能无法计算那个点的导数并返回一个NaN,所以不知道下一步在哪里取样??)完全猜到了)。
您是否可以强制初始条件为0,并查看这是否总是造成失败的抽样?
除此之外,我还会尝试调整超参数,如步骤大小、迭代次数等.
此外,您可能希望通过只使用一个N来简化示例,这可能有助于您进行诊断。祝好运!
发布于 2019-12-28 22:07:18
random.uniform's maxval默认值为None。我将其更改为1,结果就稳定了。
random.uniform(( len(N) ,), minval=0, maxval=1)https://stackoverflow.com/questions/59507603
复制相似问题