PyMC中的自定义先验分布

Question

假设我想在PyMC中对两个变量a和b设置自定义先验分布，例如：

p(a,b)∝(a+b)^(−5/2)

(有关此先验选择的动机，请参见此答案)

是否可以在PyMC中实现这个功能？如果可以，应该如何操作？

例如，在下面的模型中，我希望定义这样的先验分布来控制a和b。

import pymc as pm

# ...
# Code that defines the prior: p(a,b)∝(a+b)^(−5/2)
# ...

theta   = pm.Beta("prior", alpha=a, beta=b)

# Binomials that share a common prior
bins = dict()
for i in xrange(N_cities):
    bins[i] = pm.Binomial('bin_{}'.format(i), p=theta,n=N_trials[i],  value=N_yes[i], observed=True)

mcmc = pm.MCMC([bins, ps])

更新

根据John Salvatier的建议，我尝试了以下方法（请注意，我使用的是PyMC2，但我很乐意切换到PyMC3），但我的问题是：

1. What should I import so that I can properly inherit from Continuous?
2. In PyMC2, do I still need to stick to Theano notation?

3. Finally, how can I later tell my Beta distribution that alpha and beta have a prior from this multivariate distribution?

   import pymc.Multivariate.Continuous

   class CustomPrior(Continuous): """ p(a,b)∝(a+b)^(−5/2)

   :Parameters:
       None
   
   :Support:
       2 positive floats (parameters to a Beta distribution)
   """
   def __init__(self, mu, tau, *args, **kwargs):
       super(CustomPrior, self).__init__(*args, **kwargs)
   
   def logp(self, a,b):
   
   
       return np.log(math.power(a+b),-5./2)
   

- Amelio Vazquez-Reina

---

哦，抱歉，我以为你在使用PyMC3。虽然你可以在pymc2中完成这个任务，但是方法不同。我会更新我的回答。 - John Salvatier

Answer 1

是的！这是完全可能的，事实上也很简单。

如果您使用的是PyMC 2，请查看随机变量创建的文档。

@pymc.stochastic(dtype=int)
def switchpoint(value=1900, t_l=1851, t_h=1962):
    """The switchpoint for the rate of disaster occurrence."""
    if value > t_h or value < t_l:
        # Invalid values
        return -np.inf
    else:
        # Uniform log-likelihood
        return -np.log(t_h - t_l + 1)

如果你正在使用PyMC 3，请查看multivariate.py。请记住传递给init和logp的值都是theano变量，而不是numpy数组。这足以让你开始吗？

例如，这是多元正态分布。

class MvNormal(Continuous):
    """
    Multivariate normal

    :Parameters:
        mu : vector of means
        tau : precision matrix

    .. math::
        f(x \mid \pi, T) = \frac{|T|^{1/2}}{(2\pi)^{1/2}} \exp\left\{ -\frac{1}{2} (x-\mu)^{\prime}T(x-\mu) \right\}

    :Support:
        2 array of floats
    """
    def __init__(self, mu, tau, *args, **kwargs):
        super(MvNormal, self).__init__(*args, **kwargs)
        self.mean = self.median = self.mode = self.mu = mu
        self.tau = tau

    def logp(self, value):
        mu = self.mu
        tau = self.tau

        delta = value - mu
        k = tau.shape[0]

        return 1/2. * (-k * log(2*pi) + log(det(tau)) - dot(delta.T, dot(tau, delta)))

- John Salvatier

---

非常感谢@John。我认为我理解了上面的定义，但是我尝试按照步骤操作却没有成功。我已经在原帖中更新了更具体的问题。 - Amelio Vazquez-Reina

谢谢@John。这非常有帮助。唯一剩下的问题是我如何使用多元先验来获得单个变量（例如二项分布中的alpha和beta，如OP中的代码）的先验。我如何从多元先验中获得单个变量的先验？ - Amelio Vazquez-Reina

哦，你想在具有不同名称的变量上使用多元先验分布吗？我想没有办法这样做。您可以创建一个多元变量并对其进行索引以单独获取组件，然后对它们进行命名。这是你想要的吗？ - John Salvatier

感谢@John。在我发帖的代码示例中，我定义了theta = pm.Beta("prior", alpha=a, beta=b)。我想要做的是将我的先验分布定义为p(a,b)∝(a+b)^(−5/2)。看起来我可以使用你在回答中提到的@pymc.stochastic在PyMC2中定义p(a,b)。但我之后需要传递a和b给上面的Beta调用。这是我不知道如何做的部分。 - Amelio Vazquez-Reina

Answer 2

在 PyMC2 中，诀窍是将 a 和 b 参数放在一起：

# Code that defines the prior: p(a,b)∝(a+b)^(−5/2)
@pm.stochastic
def ab(power=-2.5, value=[1,1]):
    if np.any(value <= 0):
        return -np.Inf
    return power * np.log(value[0]+value[1])

a = ab[0]
b = ab[1]

这个笔记本有一个完整的例子。

- Abraham D Flaxman

---

谢谢！这正是我一直在寻找的。我已经寻找了很长时间这个问题的答案。顺便问一下，你知道如何在PyMC3中实现吗？ - Amelio Vazquez-Reina

不客气。我已经更新了一个错误修复（返回对数概率）和 PyMC3 版本的开端。 - Abraham D Flaxman