Simulations de lois de probabilités

Lois de probabilités

import random as rd
import matplotlib.pyplot as plt
import numpy as np

## Affichage loi

def affiche_loi(X):
    univ,X = X
    plt.clf()
    plt.bar(univ,X)
    plt.draw(); plt.show()

def compare_lois(X,Y):
    univ,X = X
    leg=['Loi théorique', 'Loi simulée']
    plt.clf()
    plt.bar(univ,X,width=.3,align='edge')
    plt.bar(univ,Y,width=-.3,align='edge')
    plt.legend(leg,loc=2,bbox_to_anchor=(0,1.15))
    plt.draw(); plt.show()

## Espérance et variance

def esperance(X):
    univ,X = X
    s=0
    for i in range(len(univ)):
        s += univ[i]*X[i]
    return s

## Loi uniforme 

def loi_unif(a,b):
    return (range(a,b+1),[1/(b-a+1) for i in range(a,b+1)])

def simul_loi_unif(a,b,N):
    res = [0 for i in range(a,b+1)]
    for _ in range(N):
        k = rd.randint(0,b-a)
        res[k] += 1/N
    return res

## Loi binomiale

def arrangement(n,p):
    k=1
    for j in range(p):
        k *= (n-j)
    return k

def fact(n):
    return arrangement(n,n)

def combinaison(n,p):
    return arrangement(n,p)//fact(p)

def loi_binomiale(n,p):
    univ = range(n+1)
    X=[combinaison(n,k) * p**k *(1-p)**(n-k) for k in range(n+1)]
    return (univ,X)

def simul_bernoulli(p):
    if rd.random() < p:
        return 1
    else:
        return 0

def simul_binomiale(n,p,N):
    res = [0 for i in range(n+1)]
    for _ in range(N):
        exp = [simul_bernoulli(p) for i in range(n)]
        s = sum(exp)
        res[s] += 1/N
    return res

## Loi hypergéométrique

def loi_hypergeometrique(N,n,R):
    univ = range(n+1)
    X = [combinaison(R,k) * combinaison(N-R,n-k) / combinaison(N,n) for k in univ]
    return (univ,X)

def simul_hypergeometrique(N,n,R,K):
    res = [0 for i in range(n+1)]
    for _ in range(K):
        population = [1 for _ in range(R)] + [0 for _ in range(N-R)]
        cmpt = 0
        for i in range(n):
            r = rd.randint(0,N-i-1)
            cmpt += population.pop(r)
        res[cmpt] += 1/K
    return res

from math import *
def compare_h_b(N,n,R):
    univ,X = loi_binomiale(n,R/N)
    univ,Y = loi_hypergeometrique(N,n,R)
    leg=['Loi binomiale', 'Loi hypergéometrique']
    plt.clf()
    plt.bar(univ,X,width=.3,align='edge')
    plt.bar(univ,Y,width=-.3,align='edge')
    plt.legend(leg,loc=2,bbox_to_anchor=(0,1.15))
    plt.draw(); plt.show()

def compare_hb_anim(N,n,R,K=20):
    for i in range(K):
        compare_h_b(N,n,R)
        plt.pause(.1)
        N*=2; R*=2

## Test

compare_lois(loi_hypergeometrique(30,10,7),simul_hypergeometrique(30,10,7,1000))
compare_lois(loi_binomiale(10,.5),simul_binomiale(10,.5,1000))
compare_hb_anim(30,20,7)