Combinaisons & Loi hypergéométrique [Résolu]

[Myshl]

Aux fanas du calcul probabiliste.
Python 3.

Une fonction "combinaisons" fondée sur une des nombreuses propriétés du triangle de Pascal, elle ne récurre que sur p.

def cnp(n,p):
if p>0:
return int(cnp(n,p-1)/p*(n-p+1))
return 1

Et donc la fonction "Loi hypergéométrique"
Np=Nombre population
Sp=Succès population
Ne=Nombre échantillon
Se=Succès échantillon

def lhg(Se,Ne,Sp,Np):
return cnp(Np-Sp,Ne-Se)/cnp(Np,Ne)*cnp(Sp,Se)

Le contrôle de cohérence sera fait avant appel

[Matthieu Brucher]

Sauf que si p et n sont entiers, il peut y avoir des problèmes d'arrondi ! Je mettrai la division en dernier, au moins.

[tyrtamos]

Bonjour,

Le code proposé pour le nombre de combinaisons est très concis et très élégant, mais pose 2 petits problèmes:

- il est effectivement nécessaire d'utiliser la division entière comme le dit Matthieu Brucher (sinon, avec de fortes valeurs de n et p: OverflowError: integer division result too large for a float). Mais, dans la formule, il faut absolument que la multiplication soit faite avant la division, sinon, le résultat peut être faux. La formule corrigée est donc:

Code :

1
2
3
4
def cnp(n,p):
    if p>0:
        return (cnp(n,p-1)*(n-p+1))//p
    return 1

- comme la fonction est récursive, elle est limitée par la pile de Python: env. 1000 pour Python 2.7 (je ne sais pas pour Python 3.x). Donc, un appel comme cnp(10000, 1000) déclenche une exception.

Voilà un calcul non récursif qui n'a pas cet inconvénient:

Code :

1
2
3
4
5
6
7
8
9
10
def combin(n, p):
    """Nombre de combinaisons de n objets pris p a p"""
    if p > n//2:
        p = n-p
    x, y, i = 1, 1, n-p+1
    while i <= n:
        x = (x*i)//y
        y += 1
        i += 1
    return x

Ce code non-récursif est un peu plus rapide que la version récursive (env. -25%).

Avec un tel code, on peut faire un calcul comme combin(1000000, 100000) qui donne un nombre de 141179 chiffres en un peu moins de 2 minutes.

[Myshl]

Merci, j'ai pris à mon compte les bons conseils.

A la fonction "Combinaison" j'ajoute la vérification de la relation d'ordre "0 <= p <= n".
(Par convention Cnp(n,p) = 0 si p > n)

Code :

1
2
3
4
5
6
7
8
9
10
11
12
def cnp(n,p):
    """Nombre de combinaisons"""
    if 0 <= p <= n:
        if p > n//2:
            p = n-p
        r, i, j = 1, 1, n-p+1
        while j <= n:
            r = (r*j)//i
            i += 1
            j += 1
        return r
    return 0

A la fonction "Loi hypergéométrique" j'ajoute aussi la vérification d'ordre "0 <= Ne <= Np" sur la combinaison en dénominateur, pour éviter une division par zéro, et puisque la probabilité est nulle si elle n'est pas calculable, et donc l'évènement impossible.

Code :

1
2
3
4
5
6
7
8
def lhg(Se,Ne,Sp,Np):
    """Loi hypergéométrique"""
    # Np = Nombre population  # Sp = Succès population
    # Ne = Nombre échantillon # Se = Succès échantillon
    r = 0
    if 0 <= Ne <= Np:
        r = cnp(Np-Sp,Ne-Se)/cnp(Np,Ne)*cnp(Sp,Se)
    return r

[Myshl]

Mea culpa!! Je corrige mon erreur dans le contrôle de la relation d'ordre pour la fonction lhg().
L'ensemble des contrôles dans lhg() et cnp() est désormais exaustif et sans redondance des contraintes propres à la loi hypergéométrique. La fonction lhg() retourne 0.0 dans tous les cas où elle n'est pas calculable (Division par 0 notamment), donc quand la probabilité qu'elle calcule est celle d'un évènement impossible. Je fais ici le choix de retourner une probabilité nulle qui me dispense de traiter des erreurs. (Le 0.0 suffit pour le calcul de probabilités composées ET et OU tout en allègeant le programme du traitement d'erreurs.)

Je conserve le contrôle de la relation d'ordre pour la fonction cnp() qui vient en complément indispensable au contrôle de la relation d'ordre dans lhg() pour cnp(Np-Sp,Ne-Se) et cnp(Sp,Se). Modifier le contrôle dans cnp() rend fausse la fonction lhg().

Pour ce qui concerne mon propre besoin (1), mes deux fonctions sont abouties. Elles calculent ce qui est calculable, ou renvoient 0 ou 0.0 (et non pas une erreur à traiter) pour ce qui est incalculable qui traduit une impossibilité et donc un évènement de probabilité nulle. Mes deux fonctions ont été vérifiées et confirmées en faisant varier les paramètres Se, Ne, Sp et Np chacun de 0 à 4, sous Python3 et sous Excel.

Merci à Matthieu Brucher et à Tyrtamos pour le coup de main.
Tyrtamos, merci du while vs récursion.
Si vous passez par ici, giflez moi si j'ai fait une erreur.
Mais bon, j'ai bossé.

Cordialement
Myshl

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def cnp(n,p):
    "Nombre de combinaisons de p dans n"
    if 0 <= p <= n:
        if p > n//2:
            p = n-p
        r, i, j = 1, 1, n-p+1
        while j <= n:
            r = (r*j)//i
            i += 1
            j += 1
        return r
    return 0
 
def lhg(Se,Ne,Sp,Np):
    "Loi hypergéométrique: probabilité de Se dans Ne quand Sp dans Np"
    if 0 <= Se <= Ne <= Np:
        return cnp(Np-Sp,Ne-Se)*cnp(Sp,Se)/cnp(Np,Ne)
    return 0.0

(1) Un cauchemar masochiste: Optimiser (gain >= mise) pour un tirage du kéno (fdj) en variant: La mise totale, le nombre de grilles, le nombre de numéros par grille, la mise par grille, le Multiple Oui/Non, le nombre de numéros communs entre grilles, plus en prenant en compte le montant du JackPot pour le tirage visé. Tout ça uniquement en terme de probabilités, donc sans toucher à de pseudo statistiques.