将 numpy 广播与 scipy truncnorm 一起使用答案

【问题标题】：Using numpy broadcasting with scipy truncnorm将 numpy 广播与 scipy truncnorm 一起使用
【发布时间】：2018-03-07 01:14:23
【问题描述】：

我想针对分位数的不同值和未截断均值的不同值评估单侧截断正态分布。为了提高效率，我想使用numpy 广播而不是 Python 循环。

对于一个最小可重复的例子，假设我要评估的三个分位数是[3.0, 2.0, 1.0]，相应的未截断平均值为[6.0, 5.0, 4.0]，下限为1.5，未截断的标准差为@ 987654326@.

单独评估这些可以按预期工作。如果我跑

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。

【问题讨论】：

看起来像一个错误。您可以在github.com/scipy/scipy/issues 上为此创建一个问题（点击绿色的“新问题”大按钮）。

标签： python python-2.7 numpy scipy array-broadcasting

【解决方案1】：

这不起作用的原因是，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)

【讨论】：

是的。我发现了同样的事情。当其中一个分位数低于截止值时，将针对大于 2 的向量长度触发该行为。但是，我真的不想在我自己的代码中添加额外的逻辑来处理截断（因为它会增加计算开销），而且奇怪的是，这种行为只发生在大于 2 的向量长度上。
查看scipy 代码，看起来函数logpdf 用-inf 填充其输出向量，并且只计算截断范围内的值。当你有一个标量时，它直接返回。如果您有一个包含范围内和范围外值的数组，它有一个函数可以选择范围内的值。显然，当数组大小大于 2 时，此函数会以某种方式出错。

【解决方案2】：

谢谢大家。与此同时，我推出了自己的：

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 广播矢量参数）。

【讨论】：