In [1]:
import numpy as np
import random
import numpy.random as npr
import matplotlib.pyplot as plt

Processus de Poisson

In [2]:
# Simulation d'un processus de Poisson de paramètre "lambd"

def Simu_Poisson(t,lambd):
    T = 0
    Z = 0 
    Time_list = [T]
    Ech_list = [Z]
    
    while T < t:
        T += -np.log(npr.rand())/lambd # ~Exp(lambd)
        if T < t:
            Time_list.append(T)
            Z += 1
            Ech_list.append(Z)
        else: 
            Time_list.append(t)
            Ech_list.append(Z)  
            
    return(Time_list,Ech_list)
In [3]:
T,Z = Simu_Poisson(10,1/2) 
plt.step(T,Z,'--P', where='post')
plt.show()
In [4]:
lamba_caché = 14
t = 100

T,Z = Simu_Poisson(t,lamba_caché)
#Z[-2]/T[-2] # estimateur de lambda : X_t/t
# autre estimateur (on supprime T[0]=0 de la liste)
# np.mean((np.array(Z[1:])/np.array(T[1:])))

Simulation d'une chaîne de Markov en temps continu (CMTC)

In [5]:
# Simulation d'une CMTC de générateur infinitesimal Q et de mesure initiale mu sur un espace d'états E 

def Simu_init(mu):
    i = random.choices(range(len(mu)), mu)[0]
    return(E[i])


def Simu_CMTC(t,Q,mu):
    N = len(Q[0]) # taille de E
    T = 0
    Time_list = [T]
    Z = Simu_init(mu)
    Ech_list = [Z]
    lambd = -np.diag(Q) # paramètres des exponentielles entre deux sauts
    P = (Q.T/lambd).T + np.identity(N)  # matrice de transition de la chaîne de Markov incluse
    
    while T < t:
        i = E.index(Z)
        T += -np.log(npr.rand())/lambd[i]
        if T < t:
            Time_list.append(T)
            j = random.choices(range(N),P[i])[0]
            Z = E[j]
            Ech_list.append(Z)
        else:
            Time_list.append(t)
            Ech_list.append(Z)
            
    return(Time_list,Ech_list)
In [6]:
# EXEMPLE 

E = [0,1.5,3,7]

mu0 = np.ones(4)*0.25

Q = np.array(
    [[-5, 3, 1, 1],
     [1, -1, 0, 0],
     [2, 1, -4, 1],
     [0, 2, 2, -4]])

# Tests : 
# np.diag(Q)<0
# np.sum(Q,axis=1)==0

# Q2 = np.array(
#    [[-1, 0.3, 0.7, 0],
#     [1, -2, 0.5, 0.5],
#     [1, 1, -3, 1],
#     [0.1, 0.1, 0.3, -0.5]])
In [7]:
T,Z = Simu_CMTC(100,Q,mu0)
plt.step(T,Z,'--P', where='post')
plt.show()
In [8]:
# Mesure d'équilibre

def Mesure_invariante(Q):
    N = len(Q[0])
    A = np.append(Q.T,[np.ones(N)],axis=0)
    b = np.zeros(N+1) 
    b[-1] = 1
    pi = np.linalg.solve(np.dot(A.T,A),np.dot(A.T,b))
    return(pi)
#Test : 
#np.dot(Mesure_invariante(Q),Q)
In [9]:
Mesure_invariante(Q)  # NB : calcul exact = [14/83, 58/83, 6/83, 5/83]
Out[9]:
array([0.1686747 , 0.69879518, 0.07228916, 0.06024096])

Ainsi, pour toute fonction $f : E \to \mathbb R$, on a $$ \lim_{t\to\infty}\frac1t\int_0^t f(X_s) d s = \sum_{x\in E}f(x)\pi_x $$ Par exemple, pour $f(x)=x^2$ on obtient que la moyenne temporelle de la variante de $X_t$ c'est la variance d'une variable de loi $\pi$, donnée ici par :

In [10]:
Var_Pi = np.dot(np.array(E)**2,Mesure_invariante(Q))-np.dot(np.array(E),Mesure_invariante(Q))**2
Var_Pi
Out[10]:
2.329583393816229

Retrouvons à l'aide de ce résultat la mesure invariante (en prenant $f(x)=\boldsymbol 1_{x=y}$ pour $y\in E$)

In [11]:
# Temps moyen passé par état

def frequences_moyennes(t,Q,mu):
    T,Z = Simu_CMTC(t,Q,mu)
    N = len(Q[0])
    Njumps = len(T)
    frequences = np.zeros(N)
    
    for i in range(Njumps):
        j = E.index(Z[i])
        if i==0:
            frequences[j] += T[1]  # on rappelle que T[0]=0
        elif i<Njumps-1:
            frequences[j] += T[i+1]-T[i]
        else:
            frequences[j] += t-T[-1]
    
    return(frequences/t)
In [12]:
frequences_moyennes(int(1e4),Q,mu0)
Out[12]:
array([0.16611697, 0.69873821, 0.07452608, 0.06061874])

Files d'attente

Le generateur $Q$ d'une file d'attente $M/M/K$ de loi d'arrivées $\mathcal(\lambda)$ et de loi de durée de traitement $\mathcal E(\mu)$ est donné par $Q_{0,0}=-Q_{0,1}=-\lambda$ puis, pour $n\geq 1$, $$ Q_{x,x-1}=\min(n,K)\mu,\quad Q_{x,x}=-(\lambda+\min(n,K)\mu),\quad Q_{x,x+1}=\lambda $$ et $Q_{x,y}=0$ autrement. Autrement dit, les taux de saut sont $\lambda(n)=\lambda+\min(n,K)\mu$ pour tout $n\in\mathbb N$ et, la matrice de transition de la chaîne incluse est $P_{0,1}=1$ puis, pour tout $n\geq 1$, $$ P_{n,n-1}=\mu\min(n,K)/(\lambda+\min(n,K)\mu),\qquad P_{n,n+1}=\lambda/(\lambda+\min(n,K)\mu) $$

In [13]:
def MM(K,lambd,mu,t): # inter-arrivées ~ Exp(lambd), temps de traitement ~ Exp(mu), K guichets.
    T = 0
    Time_list = [T]
    Z = 0
    Ech_list = [Z]
    
    while T < t:
        T += -np.log(npr.rand())/(lambd+min(Z,K)*mu)
        if T < t:
            Time_list.append(T)
            if npr.rand()<lambd/(lambd+min(Z,K)*mu):
                Z += 1
            else:
                Z -= 1
            Ech_list.append(Z)
        else:
            Time_list.append(t)
            Ech_list.append(Z)
            
    return(Time_list,Ech_list) 


# Taille moyenne empirique de la file 
def moyEst(T,Z):
    Tm = T[:-1] # on supprime T[-1]=t
    Tp = T[1:] # on supprime T[0]=0
    Tdiff = np.array(Tp)-np.array(Tm)
    t = T[-1]
    Zp = Z[:-1] # on supprime Z[-1] (qui est égal à Z[-2])
    return(Tdiff.dot(Zp)/t)
In [14]:
# Guichet M/M/1, lambd=1, mu=20 
t = int(1e4)
T,Z = MM(1,1,20,t)
plt.step(T,Z,'--P', where='post')
plt.show()
moyEst(T,Z)
Out[14]:
0.0524081289937855
In [15]:
# Guichet M/M/2, lambd=1, mu=10
t = int(1e5)
T,Z = MM(2,1,10,t)
plt.step(T,Z,'--P', where='post')
plt.show()
moyEst(T,Z)
Out[15]:
0.10131975917839088
In [ ]: