Использование вещания numpy с scipy truncnorm

Я хочу оценить одностороннее усеченное нормальное распределение для разных значений квантиля и разных значений неусеченного среднего. Для эффективности я хочу использовать numpy трансляцию, а не цикл Python.

В качестве минимального воспроизводимого примера представьте, что три квантили, которые я хочу оценить, равны [3.0, 2.0, 1.0], соответствующие неусеченные средние значения равны [6.0, 5.0, 4.0], нижняя граница равна 1.5, а неусеченное стандартное отклонение равно 3.0.

Оценка их по отдельности работает, как и ожидалось. Если я побегу

import numpy as np
from scipy.stats import truncnorm
print truncnorm.logpdf(3.0, a=(1.5-6.0)/3.0, b=np.inf, loc=6.0, scale=3.0)
print truncnorm.logpdf(2.0, a=(1.5-5.0)/3.0, b=np.inf, loc=5.0, scale=3.0)
print truncnorm.logpdf(1.0, a=(1.5-4.0)/3.0, b=np.inf, loc=4.0, scale=3.0)

я получил

-2.44840736626
-2.3878150686
-inf

(Последнее значение равно -inf, потому что 1.0 меньше порогового значения). Использование широковещательной рассылки numpy для двух значений одновременно также работает должным образом. Если я побегу

print truncnorm.logpdf(
    np.array([3.0, 2.0]),
    a=(1.5-np.array([6.0, 5.0]))/3.0,
    b=np.inf,
    loc=np.array([6.0, 5.0]),
    scale=3.0
)
print truncnorm.logpdf(
    np.array([2.0, 1.0]),
    a=(1.5-np.array([5.0, 4.0]))/3.0,
    b=np.inf,
    loc=np.array([5.0, 4.0]),
    scale=3.0
)

я получил

[-2.44840737 -2.38781507]
[-2.38781507        -inf]

Однако, если я попытаюсь оценить три значения за раз, запустив:

print truncnorm.logpdf(
    np.array([3.0, 2.0, 1.0]),
    a=(1.5-np.array([6.0, 5.0, 4.0]))/3.0,
    b=np.inf,
    loc=np.array([6.0, 5.0, 4.0]),
    scale=3.0
)

Я получаю сообщение об ошибке:

Traceback (most recent call last):
  File "truncnorm_error.py", line 25, in <module>
    scale=3.0
  File "C:\Python27\lib\site-packages\scipy\stats\_distn_infrastructure.py", line 1701, in logpdf
    place(output, cond, self._logpdf(*goodargs) - log(scale))
  File "C:\Python27\lib\site-packages\scipy\stats\_continuous_distns.py", line 4853, in _logpdf
    return _norm_logpdf(x) - self._logdelta
ValueError: operands could not be broadcast together with shapes (2,) (3,)

Что мне не хватает? Я использую Python 2.7, numpy 1.13 и scipy 0.19.


python numpy python-2.7 scipy array-broadcasting
person tcquinn    schedule 25.09.2017    source источник
comment
Похоже на ошибку. Вы можете создать для этого задачу на странице github.com/scipy/scipy/issues (нажмите кнопку большая зеленая кнопка «Новый выпуск»).   -  person Warren Weckesser    schedule 26.09.2017

Ответы (2)


arrow_upward
0
arrow_downward

Причина, по которой это не работает, заключается в том, что logpdf проверяет квантили, чтобы убедиться, что они больше, чем отсечение. Если у вас есть значение меньше, чем усечение, по-видимому, это работает для размеров 1 и 2, но не для 3. Так что это может быть ошибка.

Если вы предоставляете значения, превышающие усечение, все работает нормально. Например, это работает, когда я изменил 1,0 на 1,6 для квантилей:

print truncnorm.logpdf(
    np.array([3.0, 2.0, 1.6]),
    a=(1.5-np.array([6.0, 5.0, 4.0]))/3.0,
    b=np.inf,
    loc=np.array([6.0, 5.0, 4.0]),
    scale=3.0)
person Gerges    schedule 26.09.2017
comment
да. Я нашел то же самое. Поведение срабатывает для длин векторов больше 2, когда один из квантилей падает ниже отсечки. Однако я действительно не хочу добавлять дополнительную логику в свой собственный код для обработки отсечки (потому что это добавит вычислительных затрат), и странно, что такое поведение происходит только для длин векторов, превышающих 2. - person tcquinn; 26.09.2017
comment
Глядя на код scipy, похоже, что функция logpdf заполняет свой выходной вектор -inf и вычисляет только те значения, которые находятся в пределах диапазона усечения. Когда у вас есть скаляр, он возвращается напрямую. Если у вас есть массив со значениями как в пределах диапазона, так и за его пределами, у него есть функция, которая выбирает значения в пределах диапазона. По-видимому, эта функция каким-то образом дает сбой, когда размер массива больше двух. - person Gerges; 26.09.2017

arrow_upward
0
arrow_downward

Спасибо всем. Тем временем я накатил свой:

def left_truncnorm_logpdf(x, untruncated_mean, untruncated_std_dev, left_cutoff):
    f = np.array(np.subtract(stats.norm.logpdf(x, loc=untruncated_mean, scale=untruncated_std_dev),
                             np.log(1 - stats.norm.cdf(left_cutoff, loc=untruncated_mean, scale=untruncated_std_dev))))
    f[x < left_cutoff] = -np.inf
    return f

Это неэлегантно, и я уверен, что у него есть проблемы, но, похоже, он работает для моих целей (например, он правильно передает векторные аргументы для x и untruncated_mean).

person tcquinn    schedule 26.09.2017