test binomiale in Python per numeri molto grandi

Question

ho bisogno di fare un test binomiale in Python che permette di calcolo per 'n' numeri dell'ordine di 10000.

Ho implementato una funzione binomial_test rapido utilizzando scipy.misc. il pettine, tuttavia, è praticamente limitato attorno a n = 1000, immagino perché raggiunge il massimo numero rappresentabile mentre calcola i fattoriali o il combinatorio stesso. Qui è la mia funzione:

from scipy.misc import comb 
def binomial_test(n, k): 
    """Calculate binomial probability 
    """ 
    p = comb(n, k) * 0.5**k * 0.5**(n-k) 
    return p 

Come potrei utilizzare una funzione (o NumPy, SciPy ...) pitone nativo al fine di calcolare che binomiale di probabilità? Se possibile, ho bisogno del codice compatibile di Scipy 0.7.2.

Grazie mille!

Answer 1

Modificato per aggiungere questo commento: si noti che, come menziona Daniel Stutzbach, il "test binomiale" non è probabilmente quello che il poster originale stava chiedendo (sebbene abbia usato questa espressione). Sembra che stia chiedendo la funzione di densità di probabilità di una distribuzione binomiale, che non è quello che sto suggerendo di seguito.

Hai provato scipy.stats.binom_test?

rbp@apfelstrudel ~$ python 
Python 2.6.2 (r262:71600, Apr 16 2009, 09:17:39) 
[GCC 4.0.1 (Apple Computer, Inc. build 5250)] on darwin 
Type "help", "copyright", "credits" or "license" for more information. 
>>> from scipy import stats 
>>> print stats.binom_test.__doc__ 

    Perform a test that the probability of success is p. 

    This is an exact, two-sided test of the null hypothesis 
    that the probability of success in a Bernoulli experiment 
    is `p`. 

    Parameters 
    ---------- 
    x : integer or array_like 
     the number of successes, or if x has length 2, it is the 
     number of successes and the number of failures. 
    n : integer 
     the number of trials. This is ignored if x gives both the 
     number of successes and failures 
    p : float, optional 
     The hypothesized probability of success. 0 <= p <= 1. The 
     default value is p = 0.5 

    Returns 
    ------- 
    p-value : float 
     The p-value of the hypothesis test 

    References 
    ---------- 
    .. [1] http://en.wikipedia.org/wiki/Binomial_test 


>>> stats.binom_test(500, 10000) 
4.9406564584124654e-324 

Piccolo Modifica per aggiungere link alla documentazione: http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom_test.html#scipy.stats.binom_test

BTW: funziona su SciPy 0.7.2, così come sulle attuali 0,8 dev.

Answer 2

Non specificamente una soluzione Python, ma se si può fare con piccoli errori frazionali, si potrebbe provare a utilizzare approssimazione di Stirling per n !:

pettine (n, k) = n!/(K! * (Nk)!), dove n! è approssimativamente sqrt (2 * Pi n) (n/e)^n per grande n.

Per n> 1000 gli errori frazionari dovrebbero essere molto piccoli.

Per il calcolo delle probabilità con n grande, utilizzare logaritmi per i risultati intermedi:

log p = log (pettine (n, k)) - n * log (2)

p = exp (log (p))

Answer 3

Vorrei guardare in GNU Multi-Precision package (gmpy), che permette di eseguire calcoli a precisione arbitraria: probabilmente si potrebbe fare:

comb(n, k, exact=1)/2**k/2**(n-k) 

ma con le lunghe interi di gmpy.

Infatti, se si utilizzano calcoli interi esatti, è possibile raggiungere facilmente n = 10000 per la parte di combinazioni; per questo, è necessario utilizzare:

comb(n, k, exact=1) 

invece della virgola mobile approssimazione comb(n, k), che trabocca.

Tuttavia, come notato nel Poster originale, il numero intero (lungo) restituito potrebbe essere troppo lungo per essere moltiplicato per un valore float!

Inoltre, si verifica rapidamente un altro problema: 0.5**1000 = 9.3 ... e-302 è già molto vicino al float del flusso ...

In sintesi: se hai davvero bisogno di risultati precisi per tutti k per n~10,000, devi utilizzare un approccio diverso rispetto alla formula del post originale, che soffre dai limiti dell'aritmetica in virgola mobile a doppia precisione. Usare gmpy come indicato sopra potrebbe essere una soluzione (non testata!).

Answer 4

Qualsiasi soluzione che assomigli a comb(n, k) * 0.5**k * 0.5**(n-k) non funzioni correttamente per n. Sulla maggior parte delle piattaforme (tutte?), Il valore minimo che un float float può memorizzare è compreso tra 2 ** e 1022. Per il grande n-k o grande k, il lato destro verrà arrotondato a 0. Analogamente, il pettine (n, k) può crescere in modo così grande da non adattarsi a un oggetto mobile.

Un approccio più robusto consiste nel calcolare lo probability density function come differenza tra due punti consecutivi nello cumulative distribution function, che può essere calcolato utilizzando la funzione beta incompleta regolarizzata (consultare il pacchetto "funzioni speciali" di SciPy). Matematicamente:

pdf(p, n, k) = cdf(p, n, k) - cdf(p, n, k-1) 

Un'altra opzione è quella di utilizzare il Normal approximation, che è abbastanza preciso per la grande n. Se la velocità è una preoccupazione, questa è probabilmente la strada da percorrere:

from math import * 

def normal_pdf(x, m, v): 
    return 1.0/sqrt(2*pi*v) * exp(-(x-m)**2/(2*v)) 

def binomial_pdf(p, n, k): 
    if n < 100: 
     return comb(n, k) * p**k * p**(n-k) # Fall back to your current method 
    return normal_pdf(k, n*p, n*p*(1.0-p)) 

Non ho ancora testato il codice, ma che dovrebbe darvi l'idea generale.

Answer 5

GMPY supporta anche calcoli in virgola mobile di precisione estesa. Ad esempio:

>>> from gmpy import * 
>>> 
>>> def f(n,k,p,prec=256): 
...  return mpf(comb(n,k),prec) * mpf(p,prec)**k * mpf(1-p,prec)**(n-k) 
... 
>>> print(f(1000,500,0.5)) 
0.0252250181783608019068416887621024545529410193921696384762532089115753731615931 
>>> 

Ho specificato una precisione in virgola mobile di 256 bit. A proposito, la versione di Source Forge è obsoleta. La versione corrente è mantenuta su code.google.com e supporta Python 3.x. (Disclaimer: io sono l'attuale responsabile della gmpy.)

casevh

Answer 6

# This imports the array function form numpy 

from numpy import array 

    # the following defines the factorial function to be used in the binomial commands/ 
# n+1 is used in the range to include the nth term 

def factorial (n): 
    f=1 
    for x in range(1,n+1): 
     f=f*(x) 
    return f 

# The follwong calculates the binomial coefficients for given values of n & k 
def binomial (n,k): 
    b=1 
    b=(factorial(n)/(factorial(k)*factorial(n-k))) 
    return int(b) 

# the following lines define the pascal triangle , and print it out for 20 rows./ 
# in order to include nth term, the n +1 term needs to be in the range. The commands/ 
# append the next binomial coeficiant to a raw first and then append rows to the triangle/ 
# and prints a 20 row size pascal triangle 
def pascal(T): 
    triangle=[] 
    for n in range(T): 
     r=[] 
     for k in range(n+1): 
      r.append(binomial(n,k)) 
     triangle.append(r) 
    return triangle 

for r in pascal(20): 
    print((r))