python - scipy 多项式 pmf 返回 nan

Question

我正在尝试使用 scipy.stats (python) 中的 multinominal.pmf 函数。

当我在输入中所有概率都大于零的情况下使用此函数时，它工作正常。问题是当我想在其中一个概率为零时使用该函数。

下面的例子说明了我的意思:

In [18]: multinomial.pmf([3, 3, 0], 6, [1/3.0, 1/3.0, 1/3.0])
Out[18]: 0.027434842249657095

In [19]: multinomial.pmf([3, 3, 0], 6, [2/3.0, 1/3.0, 0])
Out[19]: nan

可以看出，在所有概率 > 0 的第一次使用该函数没有问题。然而，当我将其中一个概率更改为零时，函数返回 nan，即使函数应该返回 0.21948。

有没有一种方法(在 python 中)可以在其中一个概率为零时计算 pmf？另一种可以处理它的方法，或者解决此功能的方法。

附加信息

示例中的函数应该返回的值是我在 matlab 中使用 mnpdf 函数计算的。然而，由于我的其余代码是用 python 编写的，所以我更喜欢找到一种在 python 中计算它的方法。

Answer 1

好地方!这是 scipy 中的错误。源代码可以找到here .

3031 到 3051 行:

def pmf(self, x, n, p):
    return np.exp(self.logpmf(x, n, p))

第 2997 到 3017 行:

def logpmf(self, x, n, p):
    n, p, npcond = self._process_parameters(n, p)

第 2939 行到 2958 行:

def _process_parameters(self, n, p):

    p = np.array(p, dtype=np.float64, copy=True)
    p[...,-1] = 1. - p[...,:-1].sum(axis=-1)

    # true for bad p
    pcond = np.any(p <= 0, axis=-1)  # <- Here is why!!!
    pcond |= np.any(p > 1, axis=-1)

    n = np.array(n, dtype=np.int, copy=True)

    # true for bad n
    ncond = n <= 0

    return n, p, ncond | pcond

行pcond = np.any(p <= 0, axis=-1)结果 pcond正在true如果 p 有任何值是 <= 0。

然后在logpmf第 3029 行:

return self._checkresult(result, npcond_, np.NAN)

结果为 logpmf和 pmf返回 nan !

请注意，实际结果计算正确(第 3020 行，2994-2995 行):

result = self._logpmf(x, n, p)

def _logpmf(self, x, n, p):
    return gammaln(n+1) + np.sum(xlogy(x, p) - gammaln(x+1), axis=-1)

你的值(value)观:

import numpy as np
from scipy.special import xlogy, gammaln

x = np.array([3, 3, 0])
n = 6
p = np.array([2/3.0, 1/3.0, 0])

result = np.exp(gammaln(n+1) + np.sum(xlogy(x, p) - gammaln(x+1), axis=-1))
print(result)

>>>0.219478737997