3

假设我想在两个变量abPyMC 中放置一个自定义先验,例如:

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

(有关这种先验选择背后的动机,请参阅此答案

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

例如,我想在下面的模型中定义这样ab先验。

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. 我应该导入什么才能正确继承Continuous
  2. 在 PyMC2 中,我还需要坚持 Theano 表示法吗?
  3. 最后,我以后如何告诉我的Beta分布alphabeta从这个多元分布中获得先验?

    导入 pymc.Multivariate.Continuous

    类 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)
    
4

2 回答 2

4

是的!这是很有可能的,实际上也很简单。

如果您使用的是 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)))
于 2014-04-22T17:11:23.720 回答
3

在 PyMC2 中,诀窍是将ab参数放在一起:

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

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

于 2014-04-24T01:21:21.597 回答