#########################################
######## TP informatique bio-spé ########
########        2025-2026        ########
#########################################


#########################################
####### Thème 14 : Statistiques #########
#########################################

## Importations

import random as rd
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

## Q1. Loi normale standard [centrée réduite]

def phi(t) :
    return np.exp(-t**2/2)/np.sqrt(2*np.pi)

def PHI(x,N=1000) :
    # on néglige l'intégrale de -infini à -3 de phi
    if x < -3 :
        return 0
    S = 0
    pas = (x+3)/N
    x0 = -3
    for _ in range(N) :
        S += phi(x0)
        x0 += pas
    return pas*S

def trace_f(f,a=-3,b=3,N=200,style='-') :
    X = np.linspace(a,b,N)
    Y = [f(x) for x in X]
    plt.plot(X,Y,style)
    plt.show()

trace_f(phi)
trace_f(PHI)
trace_f(norm.cdf,-3,3,200,'--')

## Q2. Théorème de Moivre-Laplace

def B(p=1/2) :
    return int(rd.random()<p)

def Mn(n=30,p=1/2) :
    S = 0
    for _ in range(n) :
        S += B(p)
    return S/n

def run(N=100,n=30,p=1/2) :
    return [Mn(n,p) for _ in range(N)]

def run_star(N=1000,n=30,p=1/2) :
    m = p
    sigma = np.sqrt(p*(1-p)/n)
    return [(Mn(n,p)-m)/sigma for _ in range(N)]

def histo(nb = 20) :
    L = run_star(1000*nb)
    plt.hist(L,nb,density=True,ec='w',alpha=0.8)
    X = np.linspace(-3,3,200)
    Y = [phi(x) for x in X]
    plt.plot(X,Y)
    plt.show()

histo()

## Q3. Convergence en loi

def unif() :
    return rd.randint(1,6)

def expo(mu=1) :
    return -np.log(rd.random())/mu

def run_star2(loi,N,n) :
    R = []
    for _ in range(N) :
        s = 0
        for _ in range(n) :
            s += loi()
        R.append(s/n)
    m = sum(R)/N
    s2 = sum([(R[i]-m)**2 for i in range(N)])/N
    s = np.sqrt(s2)
    return [(R[i]-m)/s for i in range(N)]

def repartition(loi) :
    X = np.linspace(-3,3,200)
    Y = [norm.cdf(x) for x in X]
    plt.plot(X,Y)
    L = run_star2(loi,1000,30)
    plt.hist(L,50,cumulative=True,density=True,histtype='step')
    plt.show()

# choisir la comparaison à effectuer

repartition(B)
# repartition(unif)
# repartition(expo)
# repartition(rd.gauss)