27

在使用装饰器定义“指数随机变量的对数”的随机对象失败pymc.stochastic_from_dist后,我决定使用. 我尝试实现的模型可在此处获得(第一个模型): 在此处输入图像描述

现在,当我尝试使用 MCMC Metropolis 并以正态分布作为建议对日志(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 行中传递给 gamma 分布的 theta 值时,我看到值不是零而是小值!所以我不知道如何防止这个浮点错误?

4

3 回答 3

1

我在他们的文档中找到了这个块:

随机变量截止不能小于 D 的最大元素,否则 D 的密度将为零。标准的 Metropolis step 方法可以毫无问题地处理这种情况;它偶尔会提出非法值,但这些会被拒绝。

这会让我相信dtype=np.float(本质上与浮点数具有相同的范围)可能不是您想要的方法。文档说它需要是一个 numpy dtype,但它只需要是一个转换为 numpy dtype 对象的 dtype,并且在 Python2 中(如果我错了,请纠正我)数字 dtypes 是固定大小的类型,这意味着它们是相同的. 也许使用Decimal模块是一种选择。通过这种方式,您可以设置精度级别以封装预期值范围,并将其传递给您的扩展随机方法,在该方法中进行转换。

from decimal import Decimal, getcontext
getcontext().prec = 15
dtype=Decimal

我不知道一旦 numpy 库得到它,它是否仍然不会被截断,或者它是否会尊重继承的精度级别。我没有准确的测试方法,但试一试,让我知道它对你有用。

编辑:我测试了精确继承的概念,它似乎成立:

>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 10
>>> Decimal(1) / Decimal(7)
Decimal('0.1428571429')
>>> np.float(Decimal(1) / Decimal(7))
0.1428571429
>>> getcontext().prec = 15
>>> np.float(Decimal(1) / Decimal(7))
0.142857142857143
>>> 
于 2020-07-02T04:37:43.090 回答
0

如果你确实得到了小数字,它可能对浮点数来说太小了。这通常也是对数要避免的。如果你用dtype=np.float64呢?

于 2013-10-28T17:20:44.810 回答
0

正如您在问题末尾所建议的那样,问题在于浮点数太小,被浮点数转换为 0。一种解决方案可能是稍微调整源代码并将除法替换为例如 np.divide 和“where”条件为给定阈值的小值添加一些显式转换。

于 2020-03-30T22:50:06.640 回答