python - pymc3 のベイズ因子

Question

PyMC 3 で 2 つのモデルを比較するためにベイズ係数を計算することに興味があります。このウェブサイトによると、PyMC 2 では、手順は比較的簡単に見えます。ベルヌーイ確率変数と、最初のモデルの尤度を返すカスタム尤度関数を含めます。ベルヌーイ変数の値が 0 の場合は 2 番目のモデルの尤度、値が 1 の場合は 2 番目のモデルの尤度です。ただし、PyMC 3 では確率ノードが Theano 変数である必要があるため、事態はさらに複雑になります。

私の 2 つの尤度関数は二項式なので、このクラスを書き直す必要があると思います。

class Binomial(Discrete):
    R"""
    Binomial log-likelihood.
    The discrete probability distribution of the number of successes
    in a sequence of n independent yes/no experiments, each of which
    yields success with probability p.
    .. math:: f(x \mid n, p) = \binom{n}{x} p^x (1-p)^{n-x}
    ========  ==========================================
    Support   :math:`x \in \{0, 1, \ldots, n\}`
    Mean      :math:`n p`
    Variance  :math:`n p (1 - p)`
    ========  ==========================================
    Parameters
    ----------
    n : int
        Number of Bernoulli trials (n >= 0).
    p : float
        Probability of success in each trial (0 < p < 1).
    """
    def __init__(self, n, p, *args, **kwargs):
        super(Binomial, self).__init__(*args, **kwargs)
        self.n = n
        self.p = p
        self.mode = T.cast(T.round(n * p), self.dtype)

    def random(self, point=None, size=None, repeat=None):
        n, p = draw_values([self.n, self.p], point=point)
        return generate_samples(stats.binom.rvs, n=n, p=p,
                                dist_shape=self.shape,
                                size=size)

    def logp(self, value):
        n = self.n
        p = self.p

        return bound(
            binomln(n, value) + logpow(p, value) + logpow(1 - p, n - value),
            0 <= value, value <= n,
            0 <= p, p <= 1)

どこから始めるべきかについて何か提案はありますか?