python - 用于简单多元伯努利推理的多链绘制

Question

我想对具有多个链的多元伯努利(维度 D)进行简单的推理。下面的代码可以工作并正确推断出唯一链的参数值。 我怀疑我错误地定义了我的模型。我没有找到任何简单伯努利推理的简单例子。

返回的错误是: ValueError: Dimension must be 3 but is 2 for 'mcmc_sample_chain/simple_step_size_adaptation___init__/_bootstrap_results/mh_bootstrap_results/hmc_kernel_bootstrap_results/maybe_call_fn_and_grads/value_and_gradients/mcmc_sample_chain_simple_step_size_adaptation___init____bootstrap_results_mh_bootstrap_results_hmc_kernel_bootstrap_results_maybe_call_fn_and_grads_value_and_gradients_Samplemcmc_sample_chain_simple_step_size_adaptation___init____bootstrap_results_mh_bootstrap_results_hmc_kernel_bootstrap_results_maybe_call_fn_and_grads_value_and_gradients_Independentmcmc_sample_chain_simple_step_size_adaptation___init____bootstrap_results_mh_bootstrap_results_hmc_kernel_bootstrap_results_maybe_call_fn_and_grads_value_and_gradients_Bernoulli/log_prob/transpose' (op: 'Transpose') with input shapes: [1,5000 ,2], [2].

这是一个简单的示例，其中 D=2 且 N = 5000(训练集中的样本数)。

import numpy as np 
import tensorflow as tf
import tensorflow_probability as tfp
import functools
tfd = tfp.distributions

# ---------- DATA Generator ------------#

def generate_bernouilli(N,p):
    return np.array([np.random.binomial(size=N, n=1, p = probability) for probability in p ]).T

D = 2
N = 5000
p = np.sort(np.random.random(D))

observations = generate_bernouilli(N,p)

# ---------- Model ------------#

def make_likelihood(theta):
    one_y = tfd.Independent(
        distribution = tfd.Bernoulli(probs=theta),
        reinterpreted_batch_ndims=1)
    y = tfd.Sample(one_y,
          sample_shape=(N,))
    return y

def joint_log_prob(observations, theta):
    return (tf.reduce_sum(make_likelihood(theta).log_prob(observations)))

posterior_log_prob = functools.partial(joint_log_prob, observations)


# ---------- MCMC sampling  ------------#

num_results = int(10e3)
num_burnin_steps = int(1e3)
n_chains = 5

adaptive_hmc = tfp.mcmc.SimpleStepSizeAdaptation(
    tfp.mcmc.HamiltonianMonteCarlo(
        target_log_prob_fn=posterior_log_prob,
        num_leapfrog_steps=3,
        step_size=1.),
    target_accept_prob=tf.constant(.8),
    num_adaptation_steps=int(num_burnin_steps * 0.8))


@tf.function
def run_chain():
# Run the chain (with burn-in).
    samples, is_accepted = tfp.mcmc.sample_chain(
    num_results=num_results,
    num_burnin_steps=num_burnin_steps,
    current_state=tf.ones([n_chains,2])/10,
    kernel=adaptive_hmc,
    trace_fn=lambda _, pkr: pkr.inner_results.is_accepted)

    is_accepted = tf.reduce_mean(tf.cast(is_accepted, dtype=tf.float32))
    return samples, is_accepted


# ---------- Run  ------------#
with tf.device('/CPU:0'):
    samples, is_accepted = run_chain()

如果我们将 current_state 替换为 current_state=tf.ones([2])/10 (从而删除独立链采样)，则代码可以完美运行。

我有几个问题，我将非常感谢您的帮助: + 我的模型是否正确实现？ + 有没有办法在 tf 中调试此类错误？ python 调试器没有太大帮助。

提前致谢!

Answer 1

首先，我显然不是 tensorflow 概率方面的专家，所以这个答案很可能不是最佳实践，我只是利用我对库的有限知识来使其发挥作用，同时尝试了解更多 tensorflow 概率我。

其次，我只是想回答关于模型及其实现的问题部分，关于调试tensorflow的答案，要么谷歌一下，看看是否有一些关于它的教程，要么问另一个问题我觉得这是一个完全不同的问题。

关于模型，它看起来实现得很好，而且我不需要太多改变就能让它工作，但是，出于两个原因，我建议明确使用 theta 上的先验。第一个是，即使您不设置它，也会使用先验(一般来说，统一是一个常数，在这种情况下绝对是统一先验，并且是无界统一先验)，您可能不知道它是哪一个是或假设您正在使用与已实现的模型不同的模型。第二个是，当使用不适合当前问题的先验方法时，您可能会遇到意想不到的问题。例如，这里 theta 是维度 D 的向量，它必须在 0 和 1 之间，但是，在您的实现theta可以采用此范围之外的值；幸运的是，如果 tfd.Bernoulli 的参数在 (0,1) 之外，tensorflow 只会返回 nan，但这可能并不总是在这种情况下，它可能会抛出一个错误(这将在 theta 位于 (0,1) 之外的随机迭代中触发)，或者您可能会得到难以理解的结果，其中概率头数为1.3。

因此，我在代码中添加了一个事先并修改了以下几点:

• 我为观察添加了一个额外的维度，以便可以正确广播
• 我使用了 distribution.log_prob() 而不是 tfd.Sample 的 log_prob。我尝试直接使用 log_prob ，但我无法理解 tfd.Sample 的工作原理以及它如何影响原始发行版的 log_prog 所以我遵循了我更了解的观点。
• 我将轴设置为tf.reduce_sum。这不会给出任何错误，因为 log_prob 之前已执行并失败，但它会出现错误，因为使用多个链时，每个链都是独立的，因此每个链都有其对数后验概率。 posterior_log_prob 必须返回长度为 n_chains 的张量，而不再是标量。

以下是省略未修改部分的结果代码:

observations = generate_bernouilli(N,p)[:, None, :]

# ---------- Model ------------#

def make_prior(D):
    one_theta = tfd.Independent(
        distribution=tfd.Uniform(low=tf.zeros(D)),
        reinterpreted_batch_ndims=1
    )
    return one_theta

def make_likelihood(theta):
    one_y = tfd.Independent(
        distribution = tfd.Bernoulli(probs=theta),
        reinterpreted_batch_ndims=1
    )
    y = tfd.Sample(
          one_y,
          sample_shape=(N,)
    )
    return y

def joint_log_prob(observations, D, theta):
    return (
        make_prior(D).log_prob(theta) + 
        tf.reduce_sum(
            make_likelihood(theta).distribution.log_prob(observations), 
            axis=0
        )
    )

posterior_log_prob = functools.partial(joint_log_prob, observations, D)

# Small comment, for coherence I would also modify the following line
current_state=tf.ones([n_chains,D])/10, 
# otherwise, D != 2 would not work