scipy Multinomial pmf return nan

Question

I'm trying to use the multinominal.pmf function from scipy.stats (python).

When I use this function where all probabilities in the input bigger than zero, it work fine. The problem is when I want to use the function when one of the probabilities is zero.

The following example shows what I mean:

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

As can be seen, in the first time where all probabilities > 0 there is no problem to use the function. However when I change one of the probabilities to zero, the function return nan, even through the function should return 0.21948.

Is there a way (in python) to calculate the pmf when one of the probabilities is zero? either another way that can handle it, or a work around for this function.

additional info

The value that the function in the example should have returned I calculated using mnpdf function in matlab. However since the rest of my code is in python I prefer to find a way to calculate it in python.

Answer 1

Good spot! This is a bug in scipy. The source code can be found here.

Line 3031 to 3051:

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

Line 2997 to 3017:

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

Line 2939 to 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

The line pcond = np.any(p <= 0, axis=-1) results in pcond being true if any value of p is <= 0.

Then in logpmf line 3029:

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

results in logpmf and pmf returning nan!

Note that the actual result is calculated properly (line 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)

With your values:

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