###########################################
###########################################
## TP 05 : Diagonalisation et VAR
###########################################
###########################################

## Importations
import numpy as np
import numpy.linalg as la
import random as rd
import matplotlib.pyplot as plt

###########################################
## I : Diagonalisation
###########################################

A1 = np.array([[-3,2],[-4,3]])
A2 = np.array([[2,5,2],[0,-3,-2],[0,4,3]])
A3 = np.array([[-1,1],[-4,3]])
A4 = np.array([[0,0,1],[-1,1,1],[-2,0,3]])
A5 = np.array([[2,-1],[0,-2]])
A6 = np.array([[1,-2,-2],[-2,1,-2],[-2,-2,1]])
A7 = np.array([[-3,5],[0,2]])
A8 = np.array([[4,3,2],[-3,-2,-2],[-2,-2,-1]])

def diagonalisable(A) :
    (n,p) = A.shape
    if n!= p :
        return 'matrice non carrée'
    (L,P) = la.eig(A)
    try :
        Pm1 = la.inv(P)
        D = np.dot(Pm1, np.dot(A,P))
        for i in range(n) :
            for j in range(n) :
                if i != j and abs(D[i,j]>10**(-12)) :
                    return 'Matrice non diagonalisable'
        print('Matrice diagonalisable, semblable à :')
        print(D)
        print('\nMatrice de passage : P =')
        return P
    except :
        return 'Matrice non diagonalisable'


###########################################
## II : Simulation de VAR
###########################################

## Exercice 1 :

def pileOuFace() :
    return rd.randint(0,1)

def frequenceFace(N) :
    compteur = 0
    for _ in range(N) :
        compteur += pileOuFace()
    return compteur / N

print('\n***  Exercice 1 ***\n')
print('Fréquence de "Face" avec 1000 lancers :')
for _ in range(10) :
    print(frequenceFace(1000))

print('Fréquence de "Face" avec 10 000 lancers :')
for _ in range(10) :
    print(frequenceFace(10000))

print('La probabilité théorique est : 0.5')

## Exercice 2 :

def tirage() :
    sac = [1]*5 + [0]*3
    return rd.choice(sac)

def probaRouge(N) :
    succes = 0
    for _ in range(N) :
        succes += tirage()
    return succes / N

print('\n***  Exercice 2 ***\n')
print("Fréquence d'une boule rouge avec 10 000 tirages :")
for _ in range(10) :
    print(probaRouge(10000))
print("Fréquence théorique = 5/8 :")
print(5/8)

## Exercice 3 :

def bernoulli(p) :
    return int(rd.random() < p)

def frequenceSuccesBernoulli(N,p) :
    succes = 0
    for _ in range(N) :
        succes += bernoulli(p)
    return succes/N

print('\n***  Exercice 3 ***\n')
print("Fréquence des succès lors d'épreuves de Bernoulli de paramètre 0.3 :")
for _ in range(10) :
    print(frequenceSuccesBernoulli(1000,.3))

print("Fréquence théorique = 0.3 (p)")

## Exercice 4 :

def binomiale(n,p) :
    succes = 0
    for _ in range(n) :
        succes += bernoulli(p)
    return succes

def probaKsucces(N,n,p,k) :
    compteur = 0
    for _ in range(N) :
        compteur += binomiale(n,p) == k
    return compteur/N

print('\n***  Exercice 4 ***\n')
print("Estimation de la probabilité d'avoir 3 succès quand n = 10 et p = 0.2 :")
for _ in range(10) :
    print(probaKsucces(1000,10,0.2,3))
print("Probabilité théorique : (k parmi n).p^k.(1-p)^(n-k)")
print(10*9*8/(3*2*1)*0.2**3*0.8**7)

## Exercice 5 :

def frequenceDe(N) :
    Freq = [0]*6
    for _ in range(N) :
        de = rd.randint(1,6)
        Freq[de-1] += 1/N
    return Freq

def probaDeux6consecutifs(n,N) :
    compteur = 0
    for _ in range(N) :
        succes = 0
        for _ in range(n) :
            de = rd.randint(1,6)
            if de == 6 and succes == 1 :
                compteur += 1
                break
            elif de == 6 and succes == 0 :
                succes = 1
            else :
                succes = 0
    return compteur / N

print('\n***  Exercice 5 ***\n')
print('Fréquences des résultats lors de 1000 lancers de dé :')
print(frequenceDe(1000))

print('\nFréquences théoriques :')
print([1/6]*6)

print("\nLa probabilité d'obtenir deux SIX consécutifs en jetant 20 dés est d'environ :")
print(probaDeux6consecutifs(20,10000))
print('Probabilité théorique :')
a,b = (3-5**.5)/30, (3+5**.5)/30
r1,r2 = (5+3*5**.5)/12, (5-3*5**.5)/12
n = 20
p1 = a*r1**2*(1-r1**(n-1))/(1-r1)
p2 = b*r2**2*(1-r2**(n-1))/(1-r2)
print(p1+p2)

## Exercice 6 :

def distributionSommeDeuxDes(N) :
    Freq = [0]*11
    for _ in range(N) :
        de1 = rd.randint(1,6)
        de2 = rd.randint(1,6)
        Freq[de1+de2-2] += 1/N
    return Freq

print('\n***  Exercice 6 ***\n')
print("Distribution des sommes de deux dés lors de 1000 essais :")
print(distributionSommeDeuxDes(1000))
print("\nDistribution théorique :")
print([1/36,2/36,3/36,4/36,5/36,6/36,5/36,4/36,3/36,2/36,1/36])

## Exercice 7 :

def geometrique(p) :
    rang = 1
    while rd.random() > p :
        rang += 1
    return rang

def esperanceGeometrique(p,N) :
    E = 0
    for _ in range(N) :
        E += geometrique(p)
    return E/N

print('\n***  Exercice 7 ***\n')
print("Sur 1000 essais, on observe une moyenne de :")
print(esperanceGeometrique(0.3,1000))
print("\nMoyenne théorique = 1/p :")
print(1/0.3)

## Exercice 8 :

def riemeSucces(p,r) :
    rang = 0
    succes = 0
    echec = True
    while echec :
        rang += 1
        a = rd.random()
        if a < p :
            succes += 1
            if succes == r :
                echec = False
    return rang

def esperanceRiemeSucces(p,r,N) :
    S = 0
    for _ in range(N) :
        S += riemeSucces(p,r)
    return S/N

print('\n***  Exercice 8 ***\n')
print("Sur 1000 essais, rang d'apparition du 5ème succès consécutif quand p = 0.3 :")
print(esperanceRiemeSucces(0.3,5,1000))
print('Espérance théorique : E = r * (1/p) ')
print(5*1/0.3)

## Exercice 9 :

def PPF() :
    rang = 0
    echec = True
    piles = 0
    while echec :
        rang += 1
        piece = rd.randint(0,1) # 0: face, 1: pile
        if piece == 1 :
            piles += 1
        else :
            if piles > 1 :
                echec = False
            else :
                piles = 0
    return rang

def esperancePPF(N) :
    S = 0
    for _ in range(N) :
        S += PPF()
    return S/N

print('\n***  Exercice 9 ***\n')
print("Moyenne du temps d'attente de PilePileFace sur 1000 essais :")
print(esperancePPF(1000))
print('Espérance théorique : 8 ')

## Exercice 10 :

def longueurSerieUnicolore() :
    urne = [0,1]
    couleur = rd.choice(urne)
    longueur = 1
    while rd.choice(urne) == couleur :
        longueur += 1
    return longueur

def longueurMoyenne(N) :
    L = 0
    for _ in range(N) :
        L += longueurSerieUnicolore()
    return L/N
    
print('\n***  Exercice 10 ***\n')
print("Longueur moyenne d'une série unicolore avec 1000 essais :")
print(longueurMoyenne(1000))
print('Espérance théorique : 2 (Longueur suit une loi géométrique de paramètre 1/2) ')

## Exercice 11 :

def repetitionBernoulli(n,p) :
    L = []
    for _ in range(n) :
        L.append( int(rd.random() < p) )
    return L

def attenteDeuxSuccesConsecutifs(p) :
    rang = 0
    succes = 0
    while succes < 2 :
        rang += 1
        r = int(rd.random() < p)
        if r == 1 :
            succes += 1
        else :
            succes = 0
    return rang
    
def attenteMoyenne(p,N) :
    S = 0
    for _ in range(N) :
        S += attenteDeuxSuccesConsecutifs(p) 
    return S/N

print('\n***  Exercice 11 ***\n')
print("Succession de 20 épreuves de Bernoulli de paramètre p = 0.5 :")
print(repetitionBernoulli(20,0.5))
print('\nAttente moyenne de 2 succès consécutifs avec p = 0.5 :')
print(attenteMoyenne(0.5,1000))
print('Espérance théorique : 6 ')

## Exercice 12 :

def marche(n) :
    L = [0]
    for _ in range(n) :
        L.append( L[-1] + rd.choice([-1,1]) )
    return L

def affichage(n) :
    X = [k for k in range(n+1)]
    Y = marche(n)
    plt.plot(X,Y,'--')
    plt.plot(X,Y,'.')
    plt.title('Marches aléatoires sur une droite')
    plt.xlabel('temps')
    plt.ylabel('abscisse du point')
    plt.show()

print('\n***  Exercice 12 ***\n')
print("Positions du point au cours de 20 déplacements :")
print(marche(20))

affichage(100)