###########################################
###########################################
########## TP 09 : VAR à densité ##########
###########################################
###########################################

## Importations

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

## Exercice 1 : Lois usuelles

#1.
def unif(a,b) :
    return (b-a)*rd.random() + a

#2.
def expo(mu) :
    return -np.log(rd.random()) / mu

#3.
def normale(mu, sigma) :
    return rd.gauss(mu, sigma)

## Exercice 2 : histogramme de simulations

N = 100000
n = 50
L = [rd.gauss(0,1) for _ in range(N)]
plt.hist(L,n,density = True)
plt.show()

## Exercice 3 : minimum de 2 lois exponentielles indépendantes

#1.
def Z(lam, mu) :
    X = expo(lam)
    Y = expo(mu)
    return min(X,Y)

#2/#3. cf le cours (exercice 15)

#4.
N = 100000
n = 500
lam, mu = 0.1, .3
L = [Z(lam,mu) for _ in range(N)]
plt.hist(L,n,density = True)
limite = 7*np.log(2) / (lam+mu)
X = np.linspace(0,limite,200)
Y = (lam+mu)*np.exp(-(lam+mu)*X)
plt.plot(X,Y,'--', color = 'red')
plt.show()

# on reconnaît la densité d'une loi exponentielle de paramètre lambda + mu

## Exercice 4 : loi donnée par une densité

#1. on détermine une primitive de x |--> x**(-3/2)
# On trouve que f est une densité si et seulement si :
#  lambda == 1/2

#2. si x < 0, alors [Y < x] est quasi-impossible
# si x > 0, alors [Y < x] = [X < exp(x)]
# dont on calcule la probabilité grâce
# à la densité donnée pour X :
# P(Y < x) = intégrale de 1 à exp(x) de f
# P(Y < x) = 1 - exp(-x/2)   après calculs
# On vérifie que FY est continue sur R,
# et de classe C1 sauf éventuellement en 0.
# On en déduit que Y est une VAR à densité.
# On détermine enfin une densité de Y en dérivant
# sa fonction de répartition FY. On trouve :
# g(x) = exp(-x/2) / 2 si x > 0
# et g(x) = 0 sinon

#3. Simulation de Y : pour 0 < u < 1,
# on résout FY(x) = u, soit : 1 - exp(-x/2) = u
# on trouve : x = -2 ln(1-u)
# d'où la simulation suivante :

def Y() :
    u = unif(0,1) # c'est-à-dire u = rd.random()
    return -2*np.log(1-u)
    # rq : u et 1-u suivent la même loi, donc
    # peut aussi renvoyer : -2*np.log(u)

## Exercice 5 : loi de Cauchy

#1. effectuer le changement de variable x = a*t

#2. exprimer pour x réel : P(Z < x) = P(tan(pi.U+pi/2) < x)
# FZ(x) = P(0 < U < 1/pi * Arctan(x) + 1/2) = 1/pi * Arctan(x) + 1/2
# Par opérations, FZ est de classe C1 sur R (tout entier)
# donc Z admet une densité. Une densité de Z s'obtient en dérivant FZ :
# pour tout réel x, fZ(x) = 1/pi * 1/(1+x**2) = f1(x)

#3. simulation de f1
def X() :
    u = unif(0,1) # u = rd.random()
    return np.tan(np.pi * u + np.pi/2)

#4. estimation d'espérance
def espX(N) :
    s = 0
    for _ in range(N) :
        s += X()
    return s/N

L = []
for _ in range(10) :
    m = espX(10000)
    m = int(100*m)/100 # arrondi à 2 décimales
    L.append(m)
print(L)

# exemple de liste obtenue : [-2.01, -0.2, 11.33, 0.0, -2.84, 73.4, 4.16, 2.12, 0.17, -0.42]
# il semble que la fonction espX ne converge pas.
# On montre en effet que : x.fa(x) = ax / (pi*(x²+a²)) ~ a / (pi*x)
# et on sait que l'intégrale sur [1, +infini[ de 1/x diverge
# donc par équivalent pour des fonctions positives,
# l'intégrale sur [1, +infini[ de x.fa(x) diverge.
# X n'admet pas d'espérance !

#5. Dilatation de X

# on montre grâce à sa fonction de répartition que Y = aX admet
# pour densité la fonction fa.

def aX(a) :
    return a*X()


## Exercice 6 : loi double exponentielle

#1. on calcule séparément l'intégrale de f_ab sur [a,+infini[
# et sur ]-infini,a] en simplifiant la valeur absolue.

#2. On détermine la fonction de répartition FX de : X = a + bUV.
# Soit x un réel. FX(x) = P(X <= x) : on ne peut pas a priori
# écrire P(X < x) car on ne sait pas si X est à densité...
# V suit une loi uniforme sur {-1;1} donc [V=-1],[V=1] forment
# un système complet d'événements. On applique la formule des probabilités
# totales associée à ce SCE :
# P(X <= x) = P(V=-1)*P(X <= x | V=-1) + P(V=1)*P(X <= x | V=1)
#           = (1/2)*P(a-bU <= x | V=-1) + (1/2)*P(a+bU <= x | V=1)
#           = (1/2)*P(U >= (a-x)/b | V=-1) + (1/2)*P(U <= (x-a)/b | V=1)
#           = (1/2)*P(U >= (a-x)/b) + (1/2)*P(U <= (x-a)/b)
# par indépendance de U et V.
# 1er cas : x >= a. Alors :
# P(X <= x) = (1/2)*1 + (1/2)*FU((x-a)/b)
#           = 1/2 + (1/2)*(1-exp(-(x-a)/b))
#           = 1 - (1/2)*exp(-(x-a)/b)
# 2eme cas : x < a. Alors :
# P(X <= x) = (1/2)*(1-FU((a-x)/b)) + (1/2)*0
#           = (1/2)*exp(-(a-x)/b)
# On étudie ensuite FX :
# Par opérations, FX est continue et de classe C1
# sur ]-infini,a[ et sur ]a,+infini[.
# De plus, lim FX(x) quand x --> a vaut 1/2 à droite et à gauche
# donc FX est continue sur R tout entier.
# On a donc prouvé que X est une VAR à densité.
# On dérive FX (sauf en x=a) pour trouver une densité :
# x > a : FX'(x) = (1/2b)*exp(-(x-a)/b)
# x < a : FX'(x) = (1/2b)*exp(-(a-x)/b)
# On constate qu'on peut rassembler ces deux formules en :
# FX'(x) = (1/2b)*exp(-|x-a|/b) = f_ab

#3. Simulation de X

def X2(a,b) :
    V = rd.choice([-1,1])
    U = expo(1)
    return a + b*U*V
    
def espX2(a,b,N) :
    s = 0
    for _ in range(N) :
        s += X2(a,b)
    return s/N

def varianceX2(a,b,N) :
    L = [X2(a,b) for _ in range(N)]
    moyenne = sum(L)/N
    L2 = [(L[i]-moyenne)**2 for i in range(N)]
    return sum(L2)/N
# rq : on peut montrer que E(X2(a,b)) = a et V(X2(a,b)) = 2b²