#########################################
######## TP informatique bio-spé ########
########        2025-2026        ########
#########################################


#########################################
####### Thème 13 : VAR à densité ########
#########################################

## Importations

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

## Q1a. Loi uniforme sur [a,b]

def unif(a,b) :
    r = rd.random()
    r = r*(b-a)+a
    return r

N = 100000
n = 50
L = [unif(5,15) for _ in range(N)]
plt.hist(L,n,density=True)
plt.show()

## Q1b. Loi exponentielle de paramètre µ

def expo(mu) :
    return (-1/mu)*np.log(rd.random())

N = 100000
n = 500
L = [expo(1/2) for _ in range(N)]
plt.hist(L,n,density=True)
plt.show()

## vérification pour l'espérance

def moyenne(loi,N) :
    s = 0
    for _ in range(N) :
        s += loi()
    return s/N

def expo2() :
    mu = 1/2
    return expo(mu)


## Q1c. Lois normales

## Loi normale standard

def standard() :
    return rd.gauss(0,1)

## Loi normale (mu,sigma)

def normale(mu,sigma) :
    std = standard()
    return sigma*std+mu

N = 100000
n = 500
L1 = [standard() for _ in range(N)]
plt.hist(L1,n,density=True)
L2 = [normale(20,3) for _ in range(N)]
plt.hist(L2,n,density=True)
plt.show()

## Q3. minimum de 2 lois exponentielles indépendantes
## cf : DM9

def Z(mu1,mu2) :
    X = expo(mu1)
    Y = expo(mu2)
    return min(X,Y)

N = 100000
n = 500
mu1, mu2 = 1,3
L1 = [Z(mu1,mu2) for _ in range(N)]
plt.hist(L1,n,density=True)
plt.show()
L2 = [expo(mu1+mu2) for _ in range(N)]
plt.hist(L2,n,density=True)
plt.show()

## Q4. loi puissance

# a) on vérifie que :
#   * f est positive sur R
#   * f est continue sauf peut-être en 1
#   * l'intégrale généralisée de f sur R,
#     donc sur [1,+infini[ converge et
#     vaut 1. On utilise une primitive
#     F : x |---> -2.lambda.x**(-1/2)
#     de limite nulle en +infini, donc
#     l'intégrale de f sur R converge,
#     et lim_{+infini} F - F(1) = 2.lambda
#     donc f est une densité de probabilité
#     si et seulement si : lambda = 1/2.

# b) Loi de Y = ln(X)
#   * Y(Omega) = [0,+infini[ : Y n'est pas discrète
#   * pour x >= 0, FY(x) = P(Y<=x) = P(ln X <= x)
#     donc FY(x) = P(X <= exp(x))
#                = intégrale de 1 à exp(x) de f
#                = 1 - exp(-x/2)
# On reconnaît la fonction de répartition d'une loi
# exponentielle de paramètre 1/2.

def X() :
    Y = expo(1/2)
    return np.exp(Y)

## Q5. Loi de Cauchy

# a) fait en TD (ex.109)

# b) Loi de Z = tan(pi*U+pi/2)
#   * Z(Omega) = R : Z n'est pas discrète
#   * pour x réel,
#   FZ(x) = P(Z <= x) = P(pi*U+pi/2 <= Arctan(x)+pi)
#         = P(U <= Arctan(x)/pi + 1/2) = Arctan(x)/pi + 1/2
#   On vérifie que :
#       1) FZ est  continue sur R
#       2) FZ est de classe C1 sur R
#   Ainsi Z est une VAR à densité, et une densité de Z
#   s'obtient en dérivant FZ : f(x)=FZ'(x) = 1/pi * 1/(1+x**2)
#   Conclusion : f = f1.

def Cauchy_1() :
    U = rd.random()
    return np.tan(np.pi*U+np.pi/2)

N = 100000
print(moyenne(Cauchy_1,N))
# En calculant plusieurs moyennes sur un grand nombre
# de simulations, on a l'impression qu'il n'y a pas
# convergence : une loi de Cauchy n'admet pas d'espérance.
# cf : ex.109 des TD.

def Cauchy(a) :
    return a*Cauchy_1()

## Q6. Variation sur la loi exponentielle

# a) f_{a,b} est une densité de probabilité
# * elle est positive sur R
# * elle est continue sur R
# * son intégrale sur R se décompose en 2 intégrales :
#     - sur ]-infini,a[ on utilise une primitive
#     F1(x) = 1/2*exp((x-a)/b) de limite nulle en -infini
#     - sur [a,+infini[ on utilise une primitive
#     F2(x) = -1/2*exp(-(x-a)/b) de limite nulle en +infini
#     - on trouve alors : F1(a) - F2(a) = 1/2 - (-1/2) = 1

# b) pffff

def exp_R(a,b) :
    U = expo(1)
    V = rd.choice([-1,1])
    return a+b*U*V

N = 100000
a,b = 3,10

def expRab() :
    return exp_R(a,b)

mu = moyenne(expRab,N)
print('Moyenne empirique pour a,b = ',a,',',b)
print(mu)

def carre() :
    return (exp_R(a,b)-mu)**2

v = moyenne(carre,N)
print('Variance empirique pour a,b = ',a,',',b)
print(v)

# on trouve expérimentalement :
# - une espérance égale à a
# - une variance égale à 2b**2