"""

Thème 12

"""

import random as rd   # indispensable pour les probas
import numpy as np
import matplotlib.pyplot as plt
## Les deux fonctions essentielles

rd.random()        # renvoie un flottant de ]0;1[ choisi au hasard
rd.randint(1,50)   # renvoie un ENTIER tiré au hasard dans [1,50]

# ces deux fonctions suffisent à simuler toutes les
# expériences aléatoires qu'on étuidie. Notamment :

# 1. Jouer à pile ou face
# 2. Lancer un dé
# 3. Tirer une boule dans une urne
# 4. etc.

# Il y a une troisième fonction utile mais qui peut se construire
# avec rd.randint : la fonction rd.choice (prend au hasard un objet
# dans une liste). Elle est dans le formulaire de l'agro.


## DEFINITION DE SIMUMLER UNE VARIABLE ALEATOIRE DISCRETE

# Soit X une variable aléatoire discrète
# Mettons : X(Omega)={x1,x2,x3,.....}
#           et P(X=x1)=p1,P(X=x2)=p2,P(X=x3)=p3,....
# Simuler informatiquement X, c'est programmer une fonction f
# telle que sur un grand nombre d'appels de f, f renvoie la valeur
# x1 à la fréquence p1,  x2 à la fréquence p2, ...

# Par exemple si X suit une loi de Bernoulli de de paramètre 1/4,
# simuler X, c'est programmer une fonction f telle que
# si on appelle f 1000 fois, f a renvoyé environ 250 fois la valeur 1, et le reste des appels a renvoyé 0



## [Q1.] a) Simulation d'un pile ou face


def pileFace():
    return rd.random()<1/2  # voir le cours
    #ou : rd.randint(0,1) choisit équiprobablement entre 0 et 1

# Je teste que la fréqence du pile sur N lancers est voisine de 1/2.
N=10000

# Méthode 1
#----------
compteur = 0        # initialisation du compteur de pile
for k in range(N):  # boucle de N tirages
    compteur+= pileFace()
Freq = compteur/N   # calcul de la fréquence des "pile"
print(Freq)

# Méthode 2 : calcul à l'aide listes
#-----------------------------------
run =[pileFace() for _ in range(N)] # Série de N lancers de pièces

# le nombre de piles dans run est égal à la somme de ses termes :
S = sum(run)
# D'où la fréquence du pile :
freq = S/N
print(freq) # Python affiche 0,5037 : acceptablement proche de 0,5000

##[Q1.] b) Simulation du lancer d'un dé

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

# La proba d'obtenir 2 est approchée par la fréquence d'observation
# du 2 sur un grand nombre N de lancers du dé (d'après la loi faible
# des grands nombres, qu'on verra plus tard dans l'année).

# je fais encore un run de N lancers de dé et je ne retiens que les 2 :
N = 10000
run =[ de()  for _ in range(N)]
Liste =[a for a in run if a==2] # Cette liste contient les 2 de la liste run
# Pour avoir la fréquence du 2 :
freq = len(Liste)/N
print(freq) # devrait valoir environ 1/6 =0.1666...

##[Q1.] c) Simulation du lancer de deux dés et de la somme
N = 10000
run =[de()+de() for _ in range(N)]   # Surtout pas 2*de()
Liste= [a for a in run if a==8]      # je ne retiens que les 8
freq= len(Liste)/N
print(freq) # devrait valoir environ 5/36
print(5/36)

##[Q1.] d) Même chose que le lancer d'un dé à 10 faces ! je copie-colle Q1b)

def boule():
    return rd.randint(1,10)
N = 10000
run = [ boule()  for _ in range(N)]
Liste = [a for a in run if a==1]
freq = len(Liste)/N
print(freq) # devrait valoir environ 1/10

##[Q1.] e) Tirages successifs avec remise dans une urne
# J'applique diviser pour régner : découpage de la tâche en sous-fonctions :
#1. Fonction qui génère une suite de 5 tirages avec remise
#2. Fonction qui pour une liste  renvoie  :
#                        - True si la liste contient un 1
#                        - False sinon
#3. Calcul de la fréquence sur un grand nombre de mains


# 1. Fonction qui me génère une suite de 5 tirages avec remise
def main():
    return [boule() for _ in range(5)]
#2. Fonction qui pour une liste  renvoie  :
#                        - True si la liste contient un 1
#                        - False sinon
def contient(x,L):
    """
    Fonction qui renvoie True ssi x est dans la liste L
    """
    return x in L
#3. Calcul de la fréquence
N = 10000
compteur=0
for k in range(N):
    L = main()
    compteur+=contient(1,L)
freq = compteur/N
print(freq) # s'appelle fréquence empirique


# Remarque: combien vaut la vraie proba (on dit la probabilité théorique) ?
# réponse : Soit X la VAR égale au nombre
# de 1 obtenus dans une série de 5 tirages
# avec remise. On veut P(X>=1).
# Or X suit une loi B(5,1/10)
# Comme P(X>=1) = 1 - P(X=0),
# on déduit que P(X>=1)=1-q**5 où q =9/10
q=0.9
print(f"Proba théorique = {1-q**5}")

##[Q1.] f) Variante de e)
def nbDe3dans(L):
    """ compte le nombre de 3 dans la main L"""
    L2 = [a for a in L if a==3] # sous-liste de L constituée de ses "3"
    return len(L2) # nous donne le nombre de 3 dans L

N = 10000
run=[nbDe3dans(main())  for _ in range(N)] # run = liste de N mains

# Je retiens dans run les mains où le nb de 3 vaut 2 :
liste = [ nbde3 for nbde3 in run if nbde3==2]
freq=len(liste)/N
print(freq)

# Remarque 1. La valeur théorique est P(X=2) avec X qui
# suit  la loi  binomiale de 1e) : p = (2 parmi 5)p²q^3
p = 10*(1/10)**2*(9/10)**3
print(f"Proba théorique={p}")

# Remarque 2. Calcul de fréquence plus classique pour ceux
# qui préfèrent avec un compteur:
N = 10000
compteur = 0
for _ in range(N):
    tirage=main() # je simule 5 tirages avec remise
    compteur += (nbDe3dans(tirage)==2) # ajoute 0 ou 1 au compteur
freq=compteur/N
print(freq)



##[Q1.] g)  Temps d'attente : loi géométrique, vu en cours

def attend10():       # simule la loi géométrique
    essai = 1         # compteur du numéro du 1er succes
    while boule()<10: # tant qu'on échoue ...
        essai+=1      # ... essaie encore !
    return essai      # renvoie le numéro du 1er succès

N = 10000
run = [attend10() for _ in range(N)]
liste = [a for a in run if a>=5]
freq=(len(liste)/N)
# Remarque la valeur théorique est P(X>4) avec X qui suit
# une loi GN*(1/10), donc elle vaut q**4
p = (9/10)**4
print(f"fréquence empirique observée={freq}")
print(f"Proba théorique={p}")



## [Q1.] h) Tirages sucessifs sans remise


# Solution probablement plus proche de ce qui est attendu
# au concours au vu du formulaire Oral info Agro-Veto

def tirage_sans_rem0(L):
    #Contrairement a pop, remove ne renvoie rien
    boule = rd.choice(L) # choisit un elem au hasard dans L
    L.remove(boule)      # le retire de L
    return boule


def exph():
    Urne = [a for a in range(1,11)]
    main= []
    for _ in range(5):
        main.append(tirage_sans_rem0(Urne))
    return main

# Estimation de la proba de A : "la main contient un 10"
#-------------------------------------------------------
N = 10000     # Nombre de répétitions de l'exp.
compteur = 0  # compte le nombre de fois que A est réalisé
for _ in range(N):
    compteur+= 10 in exph()
freq=compteur/N
print(f"Fréquence observée : {freq}")
print(f"Proba théorique : 1/2")

##  Autre solution possible avec pop pour 1.h
def tirage_sans_rem2(L):
    """
    AVEC EFFET DE BORD
    utilise pop
    ne renvoie que l'élément pioché, L est modifié par
    effet
    """
    n = len(L)  # me dit combien d'élem dans L
    i = rd.randrange(n)
    x = L[i] # l'élément pioché dans L
    L.pop(i) # cette action modifie L : elle se vide.
    return x

## Q1.i) Tirages sans remise croissants

# on crée une fonction prenant en entrée une liste
# de flottants L et qui revoie True ssi la liste
# est strictement croissante

def is_croissante(L):
    """
    Teste si la suite des termes de L est strictement
    croissante
    """
    n = len(L)
    for k  in range(n-1):
        if not L[k+1]>L[k]:   # pas rangés dans le bon ordre
            return False
    return True

# Estimation de la probabilité
#-------------------------------------------------------
N = 10000    # Nombre de répétitions de l'exp.
compteur = 0  # compte le nombre de fois que A est réalisé
for _ in range(N):
    compteur+= is_croissante(exph())
freq=compteur/N
print(f"Fréquence observée : {freq}")
# Rem : la proba théorique vaut 1/(5!)
p = 1/120
print(f"Proba théorique : {p}")

## Solution cousue main (passer en 1ere lecture) pour 1.h
def tirage_sans_remise(L):
    """
    SANS EFFET DE BORD
    renvoie un élément de L tiré au hasard
    et la liste L sans cet élément
    """
    n = len(L)  # me dit combien d'élem dans L
    i = rd.randrange(n) # je choisis une place
    x =  L[i]           # je mémorise l'élem de L à cette place
    # je garde tout  ce qui n'est pas en place i dans L
    return x, L[:i]+L[i+1:]     # je renvoie l'élem sélectionné, et la liste vidée de x.

## [Q2] tirages dans une urne de contenu qui évolue

def suite_tirages(c,s):
    B=1
    N=0
    Urne=2*[N]+2*[B]
   # s tirages au max for avec interruption
    for _ in range(s):
        boule = rd.choice(Urne)
        if boule==B:            # si j'ai tiré du blanc
            return 1            # je renvoie 1
        else:
            Urne=Urne+c*[N]     # sinon je rajoute c boules noires
    return 0 # si j'arrive ici, c'est qu je n'ai jamais croisé return
             # dans le for, donc le if n'a jamais été vrai :
             # on n'a jamais eu de blanc

# Estimation de la proba (facile car j'ai simulé la VAR de Bernoulli
# qui vaut 1 si je suis tombé sur une boule blanche et 0 sinon)


def estimation_proba(c,s,N = 10000): # si l'utilisateur ne spécifie
                                     # pas la valeur de N, python
                                     # choisit 10000 par défaut
    compteur = 0
    for _ in range(N):
        compteur+=suite_tirages(c,s)
    freq=compteur/N
    return freq

## [Q3] Brebis galeuses



# je simule le lancer de n pieces. Je m'interesse au cas où le
# nombre X de piles observés vaut soit 1, soit n-1
# La loi de X est une loi B(n,1/2)

def is_galeuse(n):
    """
    cette fonction simule le lancer de n pieces et renvoie
    True si toutes les pièces sont tombées du même côté
    sauf une.
    """
    npile=0
    for _ in range(n):    # Cette boucle simule une loi B(n,1/2).
        npile+=pileFace() # pileFace() a été programmé en [Q1.]a)
    return npile==1 or npile==n-1 # 1 seul pile ou 1 seul face

def brebis(n):
    essai = 1
    while not is_galeuse(n): # Tant qu'il n'y a pas de brebis galeuse
        essai+=1
    return essai            # on a simulé une variale gémoétrique
                            # son paramètre est p=P(X=1 U X=n-1)
                            # où X est une VAR de loi B(n,1/2)

# Calcul du nombre moyen de lancers nécessaires
#----------------------------------------------
# je simule brebis(n) un grand nombre de fois.
# je calcule la moyenne des résultats observés :

N = 10000
n = 15 # Nombre de brebis dans l'assemblée
run =[ brebis(n) for _ in range(N)]
Nmoy=sum(run)/N
print(f"Nombre moyen de lancers des n pièces : {Nmoy}")


# Valeur théorique : la loi du nombre de lancers G nécessaires
# pour observer une Brebis galeuse est une loi géométrique
# de param p = P(X=1 U X=n-1) où X est une  VAR de loi B(n,1/2).
# Le nombre moyen de lancers est donc E(G)=1/p = 2^(n-1)/n
EG= 2**(n-1)/n
print(f"Moyenne  théorique {EG}")

##[Q4.] Fait en classe

#a)
N = 1000  # je tire 10000 points au hasard dans S = [0,1[
X = [rd.random() for _ in range(N)]
# Je trace ces points
Y = [0]*N    # liste de N zéros

plt.plot(X,Y,'.')
plt.show()

# on obtient un nuage de points distribué uniformément dans S

def Bernoulli(p):
    return rd.random()<p

def binomial(n,p):
    """
    on répète Bernoulli(p) et on compte le nombre de succès
    """
    compteur = 0
    for _ in range(n):
        compteur+=Bernoulli(p)
    return compteur
##[Q5.]


# Le plus efficace : stocker les tirages réalisés dans une liste

def simul_X(p):
    main=[ rd.random()<p for _ in range(2)] # 2 tirages
    n=2  # ma main est de longueur 2 : j'initialise n à 2 tirages
    while main[-1]==main[-2]:     # on compare les deux derniers tirages
        boule = Bernoulli(p)      # je tire de nouveau une boule
        main.append(boule)        # je la note dans ma liste
        n+=1
    return n

# REM : la loi du temps d'apparition de la première boule qui détone
# des précédentes N'EST pas géométrique


# Exercice de maths  : calculer la loi de la VAR X

# Je note Bk (Nk) : la kème boule tirée est blanche (noire)

#1. Espace image :  X(Omega)={2,3,4,...}

#2. Calcul des probas : Soit k >=2. J'écris  :
# [X=k]=(B1 Ո B2 Ո...Ո Bk-1 Ո Nk)U(N1 Ո N2 Ո...Ո Nk-1 Ո Bk)
# D'où par additivité finie et mutuelle indépendance des tirages :
# P(X = k)= pq**(k-1) + qp**(k-1)

#3. Il faudrait calculer E(X) pour avoir le nb moyen de
# tirages attendus ou espérés (reconnaître les t.g. de
# séries géométriques dérivées)
# Calcul : E(X)= 1/p + 1/q -1 = 1/(pq)-1

# Calcul empirique de l'espérance :
p = 1/2
q = 1-p
N = 10000
run =[simul_X(p) for _ in range(N)]
Nmoy=sum(run)/N
print(f"Nombre moyen de tirages : {Nmoy}")
EX = 1/(p*q) -1
print(f"Moyenne théorique {EX}")



## Q6. Simulation de la loi binomiale (cf. fiches révisions sup.)

def binomial(N,p):
    X=0
    for k in range(N):
        X+=Bernoulli(p)
    return X

## Q7. La loi du nombre de boules blanches est binomiale
def nblanches():
    N,p=30,5/8
    return binomial(N,p)

## Q8. La loi du no du 1er lancer donnant pile est géométrique

def geom():
    k = 1
    is_pile = Bernoulli(1/2)
    while  not is_pile:
        is_pile = Bernoulli(1/2)
        k+=1
    return k

## Q9, Q10 : vus en cours  (réponses dans le poly. d'info toutes manières !)







## Q11

def sumcum(L):
    s = 0    # initialisation des sommes cumulés
    S =[]    # initialisation de la liste des sommes cumulées
    for elem in L:
        s+=elem
        S.append(s)
    return S

## Q12 Classique et dans les fiches révisions (tombé Agro B 2025)

def simul(X,P):
    # Je calcule les sommes cumulées de la liste P, ce qui me donne
    # les beta_k, que je range dans une liste B.
    B = sumcum(P)
    r  = rd.random()  # je tire un flottant au hasard dans ]0,1[
    # je cherche le premier k tel que r>beta k
    n = len(B)
    for k in range(n): # je peux faire un for avec arrêt car je sais
        if r<B[k]:     # que B[-1]=1
            return X[k]

## Q13 Exemple d'utilisation de Q12

X = [0,1,2]
P = [1/2,1/3,1/6]

# Je vérifie qu'avec 6000 appels de simul(X,P) j'ai
# environ 3000 fois le 0, 2000 fois le 1, et 1000 fois le 2 :

run = [ simul(X,P) for _ in range(6000)]
occur = [ run.count(a) for a in X]
print(occur)

## Plus compliqué : le support de Poisson est infini
# on ne peut pas utiliser la fonction simul car la liste X
# est infinie ! pas le choix : while et calcul de la fonction
# de masse de proche en proche

def  Poisson(mu):
    """
    simule une VAR de loi de Poisson de paramètre mu
    à noter : pk = P(X=k)= exp(_mu)* mu**k/k!
    donc les pk se calculent par récurrence : pk+1 = pk * mu/(k+1)
    """
    pk,k = np.exp(-mu),0  # initialise la loi de X
    s = pk                # initialise les sommes cumulées
    r = rd.random()
    while r >= s:
        pk*=mu/(k+1)
        s+=pk
        k+=1
    return k

def cybele(p,mu):
    # je tire le nombre de visiteurs :
    N = Poisson(mu)
    # je calcule le nombre de mordus :
    return binomial(N,p)

# Essai :
p,mu = 1/3,15
print(cybele(p,mu))