#Corrigé du TP de probabilités

import numpy as np
import numpy.linalg as alg
import numpy.random as rd
from math import factorial
import matplotlib.pyplot as plt

np.set_printoptions(precision=3)

##Exercice 1

L=rd.poisson(7,1000)
M=np.zeros(11)
for x in L:
    if x<11:
        M[x]+=1
M=1/1000*M
N=np.array([np.exp(-7)*7**n/factorial(n) for n in range(11)])
print(M)
print(N)

##Exercice 2

def rademacher(n):
    T=rd.binomial(1,0.5,n)
    return(2*T-1)

def parcours(n):
    T=rademacher(n)
    s=0
    P=[0]
    for k in T:
        s+=k
        P.append(s)
    return(P)

n=20
L=[k for k in range(n+1)]
for j in range(3):
    P=parcours(n)
    plt.plot(L,P)
plt.grid()
plt.show()

def simulT():
    c=0
    s=0
    while c==0 or s!=0:
        x=rademacher(1)
        s+=x
        c+=1
    return(c)

S=0
for k in range(100):
    S+=simulT()
print(S/100)


##Exercice 3

def matalea(n,p):
    cpt=0
    for k in range(n):
        A=rd.geometric(p,(2,2))
        if alg.det(A)==0:
            cpt+=1
    return(cpt/n)

print(matalea(1000,2/3))

## Exercice 4

def geomCr(q):
    p=1-q
    C=0
    a=0
    b=rd.geometric(p)
    while a<=b:
        a=b
        b=rd.geometric(p)
        C+=1
    return(C)

def moyenneGeomCr(q):
    S=0
    for k in range(1000):
        S+=geomCr(q)
    return(S/1000)

# Si C a une variance, la loi faible des grands nombres (conséquence de Bienaymé-Tchebychev)
# assure que la moyenne de n répétitions indépendantes de la variable tend vers l'espérance lorsque n tend vers +infini,
# au sens où pour tout epsilon>0 fixé, la probabilité de s'écarter de E(X) de plus de epsilon tend vers 0.

Lq=np.linspace(0.01,0.99,101)
Lc=[moyenneGeomCr(q) for q in Lq]
plt.plot(Lq,Lc)
plt.show()

#On conjecture que E(C) tend vers +infini lorsque q tend vers 0.

print(moyenneGeomCr(0.999))
#Pour la limite lorsque q tend vers 1, on trouve une limite aux alentours de 1.7. On peut montrer que cette limite vaut e-1

## Exercice 5

def simulX():
    U=rd.random()
    return(int(1/U))

def simulY(p):
    X=simulX()
    return(rd. binomial(X,p))

##Exercice 6

def gw(p,lam):
    T=rd.poisson(lam)
    Z=0
    for k in range(T):
        Z+=rd.binomial(1,p)
    return(Z)

def verifgw(p,lam):
    S=0
    for k in range(1000):
        S+=gw(p,lam)
    EZ=S/1000
    return([EZ,p*lam]) #On a E(X)=p et E(T)=lam

print(verifgw(0.6,18))

##Exercice 7

#Le point (x,y) se trouve dans le quart du cercle si et seulement si x^2+y^2<=1.
#La probabilité que cela se produse vaut pi/4 (quotient de l'aire du quart de cercle par celle du carré)
#En calculant 4 fois cette probabilité de façon empirique, on aura donc une valeur approchée de pi.

def montecarlo(n):
    cpt=0
    for k in range(n):
        x,y=rd.random(2)
        if x**2+y**2<=1:
            cpt+=1
    return(4*cpt/n)

# On adapte la fonction pour obtenir un bôdessin

def montecarlographe(n):
    Labsbleu,Lordbleu,Labsrouge,Lordrouge=[],[],[],[]
    for k in range(n):
        x,y=rd.random(2)
        if x**2+y**2<=1:
            Labsbleu.append(x)
            Lordbleu.append(y)
        else:
            Labsrouge.append(x)
            Lordrouge.append(y)
    return(Labsbleu,Lordbleu,Labsrouge,Lordrouge)

Labsbleu,Lordbleu,Labsrouge,Lordrouge=montecarlographe(10000)
plt.plot(Labsbleu,Lordbleu,'.b',Labsrouge,Lordrouge,'.r')
plt.axis('equal')
plt.show()