- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
在 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 行检查传递给 gamma 分布的 theta 值时,我看到不是零而是小值!所以我不知道如何防止这个浮点错误?
最佳答案
我找到了这个 nugget在他们的文档中:
The stochastic variable cutoff cannot be smaller than the largest element of D, otherwise D’s density would be zero. The standard Metropolis step method can handle this case without problems; it will propose illegal values occasionally, but these will be rejected.
这会让我相信 dtype=np.float
(本质上与 float 具有相同的范围)可能不是您想要的方法。文档说它需要是一个 numpy dtype,但它只需要是一个转换为 numpy dtype 对象的 dtype,在 Python2 中(如果我错了请纠正我)数字 dtype 是固定大小的类型,这意味着它们是相同的.也许利用 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
>>>
关于python - 来自 PyMC 的 FloatingPointError 从 Dirichlet 分布中抽样,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/19434726/
我想从列表中抽取项目样本,但我想设置每个项目被包含的概率,而不是要抽取的项目总数(所以 random.sample( ) 不起作用)。我用下面的代码得到了我想要的效果(其中 p 是包含的概率,item
我正在使用 Google Analytics Reporting API,但即使指定日期范围内的 session 远少于 500K limit,我也会得到抽样结果。 .我一个月只有约 4K 次 ses
我是一名优秀的程序员,十分优秀!