2012-12-16 15 views
13

Czy istnieje funkcja wbudowana w scipy/numpy do uzyskania PMF wielomianu? Nie jestem pewien, czy binom generalizuje w prawidłowy sposób, np.wielomianowy pmf w pythonie scipy/numpy

# Attempt to define multinomial with n = 10, p = [0.1, 0.1, 0.8] 
rv = scipy.stats.binom(10, [0.1, 0.1, 0.8]) 
# Score the outcome 4, 4, 2 
rv.pmf([4, 4, 2]) 

Jaki jest prawidłowy sposób to zrobić? dzięki.

Odpowiedz

9

Nie mam wbudowanej funkcji, którą znam, a prawdopodobieństwa dwumianowe nie generalizują się (musisz znormalizować inny zestaw możliwych wyników, ponieważ suma wszystkich zliczeń musi wynosić n, co nie będzie pod opieką niezależnych dwumianów). Jednakże, jest to dość proste do wdrożenia siebie, na przykład:

import math 

class Multinomial(object): 
    def __init__(self, params): 
    self._params = params 

    def pmf(self, counts): 
    if not(len(counts)==len(self._params)): 
     raise ValueError("Dimensionality of count vector is incorrect") 

    prob = 1. 
    for i,c in enumerate(counts): 
     prob *= self._params[i]**counts[i] 

    return prob * math.exp(self._log_multinomial_coeff(counts)) 

    def log_pmf(self,counts): 
    if not(len(counts)==len(self._params)): 
     raise ValueError("Dimensionality of count vector is incorrect") 

    prob = 0. 
    for i,c in enumerate(counts): 
     prob += counts[i]*math.log(self._params[i]) 

    return prob + self._log_multinomial_coeff(counts) 

    def _log_multinomial_coeff(self, counts): 
    return self._log_factorial(sum(counts)) - sum(self._log_factorial(c) 
                for c in counts) 

    def _log_factorial(self, num): 
    if not round(num)==num and num > 0: 
     raise ValueError("Can only compute the factorial of positive ints") 
    return sum(math.log(n) for n in range(1,num+1)) 

m = Multinomial([0.1, 0.1, 0.8]) 
print m.pmf([4,4,2]) 

>>2.016e-05 

Moja implementacja współczynnika wielomianu jest nieco naiwne, i działa w przestrzeni dziennika, aby zapobiec przepełnieniu. Należy również pamiętać, że n jest zbędny jako parametr, ponieważ jest podany przez sumę zliczeń (i ten sam zestaw parametrów działa dla dowolnego n). Co więcej, ponieważ szybko zniknie w przypadku umiarkowanej n lub dużej wymiarowości, lepiej działaj w przestrzeni dziennika (logPMF również tutaj!)