Simulation
==========

Construction de variables aléatoires de loi exponentielle
---------------------------------------------------------

.. math:: X=-\frac{1}{\lambda}\ln U

où :math:`U` est une uniforme sur :math:`[0,1]`

-  la fonction map applique la fonction :math:`\exp(-x)` à tous les
   points de la liste /bins/, ce qui permet de tracer la densité
   théorique sur l’histogramme
-  la variable /bins/ représente le nombre d’intervalles sur lequel est
   construit l’histogramme

.. code:: ipython3

    %matplotlib inline
    from numpy import *
    import numpy as np
    import matplotlib.pyplot as plt
    
    
    n=3
    plt.clf()
    tableau=-log(np.random.rand(10**n))
    counts, bins = np.histogram(tableau,bins=30)
    y=map(lambda x: exp(-x),bins+bins[1]/2)
    plt.hist(bins[:-1], bins, weights=counts,density=True)
    plt.plot(bins,list(y))
    plt.show()




.. image:: output_2_0.png


Loi de Laplace
==============

-  Selon le texte 3 du rapport de la préparation agreg marocaine, la loi
   de Laplace qui a pour densité

   .. math:: \frac12 \exp(-|x|)

se construit comme le produit d’une variable aléatoire de Bernoulli de
paramètre :math:`1/2` et d’une autre variable aléatoire indépendante de
loi exponentielle de paramètre :math:`1`.

-  On utilise ce schéma pour simuler de telles variables et on compare
   l’histogramme empirique ainsi obtenu à la densité théorique.

.. code:: ipython3

    from numpy import *
    import numpy as np
    import matplotlib.pyplot as plt
    
    
    n=5
    plt.clf()
    
    def bernoulli(p):
        if np.random.rand()<=p:
            return 1
        else:
            return -1
        
    v=[bernoulli(0.5) for k in range(10**n)]
    y=-log(np.random.rand(10**n))
    tableau=y*v
    counts, bins = np.histogram(tableau,bins=50)
    y=map(lambda x: 0.5*exp(-abs(x)),bins)
    plt.hist(bins[:-1], bins, weights=counts,density=True)
    plt.plot(bins,list(y))
    plt.show()
    




.. image:: output_4_0.png


Loi gaussienne par rejet
========================

-  On voit assez facilement que

   .. math:: \sup_{x\in\mathbf R} \frac{\exp(-x^2/2)/\sqrt{2\pi}}{\exp(-|x|)/2}=\sqrt{e}

   ce qui revient à dire que

   .. math::  \frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}{2})\le \sqrt{e}\  \frac{1}{2}\exp(-|x|)

   donc on peut appliquer la méthode du rejet pour simuler une
   gaussienne à partir de la simulation d’une loi de Laplace.

-  Comme :math:`\sqrt{e}\simeq 1.65`, il faudra en moyenne :math:`1.65`
   tirages de la loi de Laplace, soit plus de 3 appels au générateur
   aléatoire, pour un tirage de gaussienne. C’est moins performant que
   la méthode de Rayleigh.

.. code:: ipython3

    def laplace():
        return -bernoulli(0.5)*log(np.random.rand())
    
    def f(x):
        return exp(-x**2/2)/sqrt(2*pi)
    
    def g(x):
        return 0.5*exp(-abs(x))
    
    def gaussienne():
        y=laplace()
        u=np.random.rand()
        while(f(y)<=exp(0.5)*u*g(y)):
            y=laplace()
            u=np.random.rand()
        return y
    
    tableau=[gaussienne() for k in range(10**5)]
    counts, bins = np.histogram(tableau,bins=50)
    y=map(lambda x: f(x),bins)
    plt.hist(bins[:-1], bins, weights=counts,density=True)
    plt.plot(bins,list(y))
    plt.show()



.. image:: output_6_0.png


File d’attente Geo/Geo/1
------------------------

-  On considère la file Geo/Geo/1 décrite dans le texte 4. Dans cette
   file, il y a un seul serveur et une file d’attente infinie. À chaque
   instant, il y a une arrivée avec probabilité :math:`p` et un départ
   avec probabilité :math:`q`.

-  Cela signifie que le temps entre deux arrivées suit une loi
   géométrique de paramètre :math:`p` et les temps de service une loi
   géométrique de paramètre :math:`q`.

-  La simulation consiste à partir d’un état :math:`x_0` choisi par
   l’utilisateur et on utilise l’algorithme du poly.

-  Ici, les transitions sont simples

   -  si :math:`x=0`, alors on passe à l’état :math:`1` avec probabilité
      :math:`p` et on reste en :math:`0` sinon
   -  si :math:`x>0`, on passe en :math:`x+1` avec probabilité
      :math:`p(1-q)`, en :math:`x-1` avec probabilité :math:`q(1-p)`,
      sinon on reste en :math:`x`.
   -  Du point de vue de la simulation, on pose

      .. math::  r_1=p(1-q), \ r_2=r1+q(1-p)

      si le nombre aléatoire :math:`u` est inférieur à :math:`r_1`, on
      augmente :math:`x` de :math:`1`. Sinon et si :math:`u` est
      inférieur à :math:`r_2`, on diminue :math:`x` de :math:`1`. Sinon
      on ne fait rien.

Remarques
~~~~~~~~~

-  en faisant varier les valeurs relatives de :math:`p` et :math:`q`, on
   constate aisément deux comportements

   -  si :math:`p>q` alors très rapidement, on quitte l’état :math:`0`
      et la valeur du nombre de clients dans le système semble tendre
      vers l’infini.
   -  si :math:`q<p`, on repasse très souvent par :math:`0`, ce qui
      semble dire qu’il est récurrent.

-  le cas limite où :math:`p=q` est un peu plus compliqué à observer. On
   se rend compte que l’on revient à :math:`0` mais que le temps entre
   deux passages à :math:`0` est assez long. En fait, on montre que dans
   ce cas, la chaîne est récurrente nulle : le temps de retour à
   :math:`0` est fini mais d’espérance infinie.

Explication
^^^^^^^^^^^

On peut réécrire

.. math::  p<q \Longleftrightarrow \frac{1}{q}<\frac{1}{p}

ce qui revient à dire que

/Le terme de service moyen est strictement inférieur au temps moyen
entre deux arrivées/

Il est donc logique après un petit moment de réflexion, que le système
ne sature pas. A contrario, si :math:`q>p`, il arrive plusieurs clients
pendant la durée d’un service

.. code:: ipython3

    from numpy import *
    import numpy as np
    import matplotlib.pyplot as plt
    
    
    def transition(x,p,q):
        if (x==0):
            if (np.random.rand()<=p):
                return 1
            else:
                return 0
        else:
            r1=p*(1-q)
            r2=r1+q*(1-p)
            u=np.random.rand()
            if (u<=r1):
                return x+1
            else:
                if (u<=r2):
                    return x-1
                else:
                    return x
            
            
    x=0
    etats=[]
    N=1000
    for j in range(N):
        etats+=[x]
        x=transition(x,0.5,0.25)    
    plt.scatter(np.asarray(range(N)),np.asarray(etats),s=1)
    plt.show()



.. image:: output_8_0.png


Estimation par les moyennes de Césaro
=====================================

.. code:: ipython3

    def transition(x,p,q):
        if (x==0):
            if (np.random.rand()<=p):
                return 1
            else:
                return 0
        else:    
            r1=p*(1-q)
            r2=r1+q*(1-p)
            u=np.random.rand()
            if (u<=r1):
                return x+1
            else:
                if (u<=r2):
                    return x-1
                else:
                    return x
            
    
    N=100000
    cumul=0
    x=0
    for j in range(N):
        x=transition(x,0.5,0.6)
        if (x==2):
            cumul +=1
    cumul/N




.. parsed-literal::

    0.18486



Estimation directe
==================

pas ce qui était demandé

.. code:: ipython3

    N=100000
    x=0
    K=1000
    cumul=0
    for k in range(K):
        for j in range(N):
            x=transition(x,0.5,0.6)
        if(x==2):
            cumul+=1
    cumul/K





.. parsed-literal::

    0.18



Calcul algébrique de la probabilité stationnaire
================================================

.. code:: ipython3

    
    p=0.5
    pbar=0.5
    q=0.6
    qbar=1-q
    A=np.array([[pbar-1,q*pbar,0,0],[p,p*q+pbar*qbar-1,q*pbar,0],[0,p*qbar,p*q+pbar*qbar-1,q],[1,1,1,1]])
    np.linalg.solve(A,np.array([0,0,0,1]))




.. parsed-literal::

    array([0.24107143, 0.40178571, 0.26785714, 0.08928571])



.. code:: ipython3

    def transition(x,p,q):
        if (x==0):
            if (np.random.rand()<=p):
                return 1
            else:
                return 0
        elif (x==3):
            if (np.random.rand()<q):
                return 2
            else:
                return 3
        else:    
            r1=p*(1-q)
            r2=r1+q*(1-p)
            u=np.random.rand()
            if (u<=r1):
                return x+1
            else:
                if (u<=r2):
                    return x-1
                else:
                    return x
                
    N=1000000
    pi_=np.zeros(4)
    x=0
    for j in range(N):
        x=transition(x,0.5,0.6)
        pi_[x]+=1
    
    pi_/N




.. parsed-literal::

    array([0.241419, 0.401801, 0.267611, 0.089169])



Simulation régénérative
=======================

.. math::  \pi(y) = \frac{\sum_{j=1}^K N_{xy}^j}{\sum_{j=1}^k L^j}

où - :math:`N_{xy}^j` =nb de visites à :math:`y` dans la :math:`j`-ième
excursion - :math:`L^j` = longueur de la :math:`j`-ième excursion

VOir Asmussen, Applied probabilities and queues

https://nextcloud.r2.enst.fr/nextcloud/index.php/s/gYSfF4AFjYX429T

page 174

.. code:: ipython3

    def transition(x,p,q):
        if (x==0):
            if (np.random.rand()<=p):
                return 1
            else:
                return 0
        elif (x==3):
            if (np.random.rand()<q):
                return 2
            else:
                return 3
        else:    
            r1=p*(1-q)
            r2=r1+q*(1-p)
            u=np.random.rand()
            if (u<=r1):
                return x+1
            else:
                if (u<=r2):
                    return x-1
                else:
                    return x
                
                
    
    K=1000
    visites_=np.zeros((K,4))
    longueurs_=np.ones(K)
    k=0
    while(k<K):
        x=0
        y=transition(x,0.5,0.6)
        while(y !=x):
            visites_[k,y]+=1
            longueurs_[k]+=1
            y=transition(y,0.5,0.6)
        k+=1
    
        

.. code:: ipython3

    for i in range(1,4):
            L=np.sum(longueurs_)
            print(i,np.sum(visites_[:,i])/L)



.. parsed-literal::

    1 0.41118944735046625
    2 0.26381623834432566
    3 0.09756652262906527







Estimation de pi
----------------

On tire des points :math:`(U_i,V_i)` uniformément dans le carré
:math:`[0,1]\times [0,1]`. D’après la loi forte des grands nombres

.. math::  \frac{1}{N} \sum_{i=1}^N \mathbf 1_{\{U_i^2+V_i^2\le 1\}}\longrightarrow \frac{\pi}{4}


.. code:: ipython3

    %matplotlib inline
    from numpy import *
    import numpy as np
    import scipy
    from pylab import* 
    import pylab as pylab
    
    
    N=10**2
    u=2*np.random.rand(N)-1
    v=2*np.random.rand(N)-1
    frequence=sum(u**2+v**2<=1)*1.0/N
    variance=frequence*(1-frequence)
    print "estimation de pi ", 4*frequence
    print "\nintervalle de confiance à 95% ",4*(frequence-1.96*sqrt(variance)/sqrt(N)),4*(frequence+1.96*sqrt(variance)/sqrt(N))


.. parsed-literal::

    estimation de pi  2.96
    
    intervalle de confiance à 95%  2.6161107527124465 3.3038892472875534