FloatingPointError от PyMC при выборке из распределения Дирихле

После безуспешно используя декораторы для определения стохастического объекта «логарифма экспоненциальной случайной величины», я решил вручную написать код для этого нового распределения, используя pymc.stochastic_from_dist. Модель, которую я пытаюсь реализовать, доступна здесь (первая модель): введите здесь описание изображения

Теперь, когда я пытаюсь выполнить выборку журнала (альфа) с помощью MCMC Metropolis и с нормальным распределением в качестве предложения (как указано на следующем рисунке в качестве метода выборки), я получаю следующую ошибку:

  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)

Когда я проверяю значения тета-распределения, переданные в гамма-распределение в строке 978 Distributions.py, я вижу, что там не нулевые, а небольшие значения! Итак, я не знаю, как предотвратить эту ошибку с плавающей запятой?


python random floating-point sampling pymc
person Cupitor    schedule 17.10.2013    source источник
comment
Я думаю, что в основном коде Distributions.py должен быть обработчик исключений, чтобы PyMC обрабатывал случаи, когда gammas[0] равен нулю!   -  person Cupitor    schedule 17.10.2013
comment
какой будет результат, если вы добавите np.seterr(divide='ignore') сразу после импорта в код тестовой выборки?   -  person alko    schedule 24.10.2013
comment
@alko, спасибо, но все та же ошибка.   -  person Cupitor    schedule 29.10.2013

Ответы (3)


arrow_upward
1
arrow_downward

Я нашел этот самородок в их документации:

Отсечение стохастической переменной не может быть меньше наибольшего элемента D, иначе плотность D была бы равна нулю. Стандартный пошаговый метод Метрополиса может справиться с этим случаем без проблем; время от времени он будет предлагать недопустимые значения, но они будут отклонены.

Это заставило бы меня поверить, что dtype=np.float (что важно, имеет тот же диапазон, что и float), возможно, не тот метод, который вы хотите использовать. В документации говорится, что это должен быть пустой тип dtype, но он просто должен быть типом dtype, который преобразуется в объект numpy 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
>>>
person Xinthral    schedule 02.07.2020
comment
Эй, Джесси, спасибо за ответ, но в последнее время я почти ничего не помню о PyMC. Хотелось бы лучшего обслуживания. Я пошел дальше, ха-ха! - person Cupitor; 03.07.2020

arrow_upward
0
arrow_downward

Если вы получаете маленькие числа, они могут быть просто слишком малы для числа с плавающей запятой. Как правило, это также то, для чего используются логарифмы. Что делать, если вы используете dtype=np.float64?

person jaap    schedule 28.10.2013
comment
Дело в том, что нужно поиграть с исходным кодом PyMC! Чего я предпочитаю нет! Но спасибо. - person Cupitor; 29.10.2013

arrow_upward
0
arrow_downward

Как вы предположили в конце своего вопроса, проблема связана со слишком маленькими числами, которые приводятся к 0 с плавающей запятой. Одним из решений может быть небольшая настройка исходного кода и замена деления, например, на np.divide и в Условие «где», чтобы добавить явное приведение к малым значениям к заданному порогу.

person EPH    schedule 30.03.2020