从Dirichlet分布中进行采样的PyMC中的FloatingPointError

Question

being unsuccessful in using decorators定义了“指数随机变量的对数”的随机对象之后，我决定使用pymc.stochastic_from_dist手动编写此新分布的代码。我尝试实现的模型在这里可用（第一个模型）：

现在，当我尝试使用MCMC Metropolis以正态分布作为建议来采样log（alpha）（如以下图片中作为采样方法所述时，出现以下错误：

  File "/Library/Python/2.7/site-packages/pymc/distributions.py", line 980, in rdirichlet
    return (gammas[0]/gammas[0].sum())[:-1]

FloatingPointError: invalid value encountered in divide

尽管采样没有出错，但采样直方图与本文中的相匹配。我的分层模型是：

"""
A Hierarchical Bayesian Model for Bags of Marbles

logalpha ~ logarithm of an exponential distribution with parameter lambd
beta ~ Dirichlet([black and white ball proportions]:vector of 1's)
theta ~ Dirichlet(alpha*beta(vector))

"""

import numpy as np
import pymc
from scipy.stats import expon
lambd=1.
__all__=['alpha','beta','theta','logalpha']
#------------------------------------------------------------
# Set up pyMC model: logExponential
# 1 parameter: (alpha)

def logExp_like(x,explambda):
    """log-likelihood for logExponential"""
    return -lambd*np.exp(x)+x
def rlogexp(explambda, size=None):
    """random variable from logExponential"""
    sample=np.random.exponential(explambda,size)
    logSample=np.log(sample)
    return logSample
logExponential=pymc.stochastic_from_dist('logExponential',logp=logExp_like,
                                          random=rlogexp,
                                          dtype=np.float,
                                          mv=False)
#------------------------------------------------------------
#Defining model parameteres alpha and beta.
beta=pymc.Dirichlet('beta',theta=[1,1])
logalpha=logExponential('logalpha',lambd)

@pymc.deterministic(plot=False)
def multipar(a=logalpha,b=beta):
    out=np.empty(2)
    out[0]=(np.exp(a)*b)
    out[1]=(np.exp(a)*(1-b))
    return out
theta=pymc.Dirichlet('theta',theta=multipar)

我的测试抽样代码是：

from pymc import Metropolis
from pymc import MCMC
from matplotlib import pyplot as plt
import HBM
import numpy as np
import pymc
import scipy
M=MCMC(HBM)
M.use_step_method(Metropolis,HBM.logalpha, proposal_sd=1.,proposal_distribution='Normal')
M.sample(iter=1000,burn=200)

[当我检查在distributions.py的第978行中传递给伽玛分布的theta的值时，我发现不存在零而是很小的值！所以我不知道如何防止这种浮点错误？

Answer 1

如果您得到的数字很小，则可能太小而无法浮动。通常，这也是避免对数的原因。如果您使用dtype=np.float64怎么办？

Answer 2

正如您在问题末尾所建议的那样，问题在于数字太小而无法转换为0。一种解决方案是调整一些源代码，并用np.divide和在“ where”条件下，将一些显式强制转换添加到给定阈值的较小值。