PyMC 中的自定义先验 [英] Custom priors in PyMC

问题描述

假设我想在 PyMC 中的两个变量 a 和 b 上放置一个自定义先验，例如:

Say I want to put a custom prior on two variables a and b in PyMC, e.g.:

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

(关于这种先验选择背后的动机，请参阅这个答案)

(for the motivation behind this choice of prior, see this answer)

这可以在 PyMC 中完成吗?如果是这样怎么办?

Can this be done in PyMC? If so how?

举个例子，我想在下面的模型中在 a 和 b 上定义这样的先验.

As an example, I would like to define such prior on a and b in the model below.

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)，但我的问题是:

Update

Following John Salvatier's advice, I try the following (note that I am in PyMC2, although I would be happy to switch to PyMC3), but my questions are:

1. 我应该导入什么才能正确继承Continuous?
2. 在 PyMC2 中，我还需要坚持使用 Theano 表示法吗?

3. 最后，我以后如何告诉我的 Beta 分布 alpha 和 beta 有一个来自这个多变量分布的先验?

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?

导入 pymc.Multivariate.Continuous

import pymc.Multivariate.Continuous

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

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)

推荐答案

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

In PyMC2, the trick is to put the a and b parameters together:

# 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]

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

这篇关于PyMC 中的自定义先验的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！