2017-12-19 9 views
2

scipy.stats (python)의 multinominal.pmf 함수를 사용하려고합니다.scipy Multinomial pmf return nan

입력이 0보다 큰 모든 확률이있는 곳에서이 함수를 사용할 때 제대로 작동합니다. 문제는 확률 중 하나가 0 일 때 함수를 사용하려고 할 때입니다.

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이 기능을 사용하려면 문제가 없습니다 처음에 볼 수 있습니다 :

다음 예

무슨 뜻인지 보여줍니다. 그러나 하나의 확률을 0으로 변경하면 함수는 nan을 반환하고 함수를 통해서도 0.21948을 반환해야합니다.

확률 중 하나가 0 일 때 pmf를 계산하는 방법이 있습니까 (파이썬에서)? 이 기능을 처리 할 수있는 또 다른 방법이거나이 기능에 대한 해결 방법입니다.

추가 정보

예에서 함수 I는 매트랩 mnpdf 함수를 이용하여 계산 한 복귀해야하는 값. 그러나 나머지는 내 코드가 파이썬이므로 파이썬으로 계산하는 방법을 선호합니다.

+1

. 이 문제는 https://github.com/scipy/scipy/issues에서 만들 수 있습니까? (녹색의 "New issue"버튼을 클릭하십시오.) –

+0

예, 지금하겠습니다. – rfire

답변

4

좋은 점! 이것은 scipy의 버그입니다. 소스 코드는 here입니다.

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

라인 2997 3017로 : : 3031 3051에

라인

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에서 pcond = np.any(p <= 0, axis=-1) 결과가있는 경우 truep의 값은입니다. 3029 logpmf 라인이어서 543,210 = 0

: logpmfpmfnan로 되돌아

return self._checkresult(result, npcond_, np.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) 

:

버그처럼 보이는
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 
+0

훌륭한 분석. 잠시라도이 정보를 @rfire가 만든 scipy 문제 (https://github.com/scipy/scipy/issues/8235)에 추가 할 수 있습니까? –

+0

@WarrenWeckesser 이것은 이미 완료되었습니다! –

+0

빨리, 고마워! –