# -*- coding: utf-8 -*-
"""
Created on Sat Mar 22 10:44:31 2025

@author: arnau
"""

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

#%% Intro




#%% LGN

def piece(p):
    if rd.rand()<p:
        return 1
    else :
        return 0

def valeur_approchee(p,n=1000):
    s = 0
    for k in range(n):
        s+=piece(p)
    return s/n
      
    
    
def test_moyenne(e,p,n=100,N=1000):
    s=0
    for i in range(N):
        moy=0
        for k in range(n):
            if piece(p)==1:
                moy+=1/n
        if np.abs(moy-p)>e:
            s+=1
    return s/N

#%% LGN fonction de répartition empirique

a = 2

def repart_exp(a,x):
    if x<0:
        return 0
    else:
        return 1-np.exp(-a*x)
    
def repart_exp_emp(a,x):
    expo = rd.exponential(1/a,1000)
    s = 0
    for simu in expo:
        if simu<=x:
            s+=1
    return s/1000


X = np.linspace(-5,10,1000)
plt.plot(X, [repart_exp(a,x) for x in X], label='fct repartion E('+str(a)+')')
plt.plot(X, [repart_exp_emp(a,x) for x in X],'--', label='fct repartion empirique E('+str(a)+')')
plt.legend()
#%% Approximation des lois de Poisson 


lambda_ = 2
n = 30

def poisson(l):
    return [l**k/math.factorial(k)*np.exp(-l) for k in range(n+1)]

def binom(n,p):
    return [ math.factorial(n)/(math.factorial(k)*math.factorial(n-k))*p**k*(1-p)**(n-k) for k in range(n+1)]

plt.bar(range(n+1),poisson(lambda_),width = 0.4, align="edge", label = 'Poisson')
plt.bar(range(n+1),binom(n,lambda_/n),width = 0.4, label = 'Binomiale')
plt.legend()

'''
Diagramme en barres ( != histogramme ). Abscisses : entiers k ; ordonnées : P(X=k)
'''

#%% TCL

def Mstar_poisson(n,lambda_):
    '''
    Renvoie Mn* pour une suite de variables indépendantes de loi P(lambda_)

    '''
    temp = rd.poisson(lambda_,n)
    return (sum(temp)-n*lambda_)/np.sqrt(n*lambda_)


def Mstar_geo(n,p):
    '''
    Renvoie Mn* pour une suite de variables indépendantes de loi G(p)

    '''
    temp = rd.geometric(p,n)
    return (sum(temp)-n/p)/np.sqrt(n*(1-p)/p**2)

def Mstar_exp(n,a):
    '''
    Renvoie Mn* pour une suite de variables indépendantes de loi E(a)

    '''
    temp = rd.exponential(1/a,n)
    return (sum(temp)-n/a)/np.sqrt(n/a**2)


def densite_normal(x):
    return 1/(np.sqrt(2*np.pi))*np.exp(-x**2/2)


t = np.linspace(-5,5,1000)
densite = densite_normal(t)
plt.plot(t,densite)


N = 10000
n = 1000
L = []
for _ in range(N):
    '''
    simule N = 10 000 réalisation de Mn*
    '''
    L.append(Mstar_exp(n,2))

plt . hist (L,50, density = True, align ='mid')

#%% Moivre Laplace

def moivre(n,p):
    x = rd.binomial(n,p)-n*p
    return x/(np.sqrt(n*p*(1-p)))

L = [moivre(300,1/3) for k in range(10000)]

plt.hist(L,20, density=True , align='mid')
densite = densite_normal(t)
plt.plot(t,densite,color='red')