//////////////////////////
//generateur loi exponentielle
//lambda : parametre de la loi
//m,n : taille de la matrice de sortie
//////////////////////////

function X=exponentielle(lambda,m,n)
    X=-log(rand(m,n))./lambda;
endfunction




//////////////////////////
//simulation d'une trajectoire d'un processus de Poisson
//lambda : intensité du processus
//T : temps final pour la trajectoire
//////////////////////////

function poisson(lambda,T)
    Tps=[0];
    N=0;
    S=exponentielle(lambda,1,1);
    while (S<T)
        Tps=[Tps S];
        N=N+1;
        S=S+exponentielle(lambda,1,1);
    end
    plot2d2([Tps T],[0:N N]);
endfunction





//////////////////////////
//illustration de la convergence p.s. de l'estimateur de l'intensite et intervalle de confiance asymptotique
//lambda : intensité du processus
//T : temps final pour la trajectoire
//alpha : risque pour l'intervalle de confiance asymptotique
//////////////////////////

function estimationpoisson(lambda,T,alpha)
    Tps=[0];
    N=0;
    grille=linspace(0,T,1000);
    estimateur=zeros(1,1000);
    S=exponentielle(lambda,1,1);
    while (S<T)
        Tps=[Tps S];
        N=N+1;
        S=S+exponentielle(lambda,1,1);
        estimateur(grille>=S)=estimateur(grille>=S)+1;
    end
    estimateur=estimateur(2:$)./grille(2:$);
    plot2d([0,T],[lambda lambda]);
    plot2d(grille(2:$),estimateur,2);
    quantile=cdfnor("X",0,1,1-alpha,alpha);
    intsup = estimateur+quantile*sqrt(estimateur./grille(2:$));
    intinf = estimateur-quantile*sqrt(estimateur./grille(2:$));
    plot2d(grille(2:$),intsup,3);
    plot2d(grille(2:$),intinf,3);
endfunction



//////////////////////////
//trace la fonction de répartition empirique
//X : échantillon
//////////////////////////

function []=frempirique(X)
    Y=-gsort(-X);
    n=prod(size(Y));
    Y=[Y(1),Y];
    plot2d2(Y,(0:n)/n);
endfunction



//////////////////////////
//Illustration de la convergence en loi de l'estimateur renormalisé
//lambda : intensité du processus
//T : temps final pour la trajectoire
//n : nombre de trajectoires simulées
//////////////////////////

function convloiPoisson(lambda,T,n)
    Tps=zeros(1,n);
    N=zeros(1,n);
    S=exponentielle(lambda,1,n);
    while (~prod(S>=T))
        indices=(S<T);
        N(indices)=N(indices)+1;
        k=sum(indices);
        S(indices)=S(indices)+exponentielle(lambda,1,k);
    end
    R=sqrt(T/lambda)*(N.*(1/T)-lambda);
    frempirique(R);
    X=linspace(-3,3,1000);
    Y=cdfnor("PQ",X,zeros(1,1000),ones(1,1000));
    plot2d(X,Y,2);
endfunction




//////////////////////////
//illustration de la convergence vers un mouvement brownien
//lambda : intensité du processus
//t : temps final pour la trajectoire après renormalisation
//u : taux de compression, i.e. u=T/t avec T le temps final avant renormalisation
//////////////////////////

function convfonct(lambda,t,u)
    Tps=[0];
    N=0;
    grille=linspace(0,t,1000);
    proc=zeros(1,1000);
    S=exponentielle(lambda*u,1,1);
    while (S<t)
        Tps=[Tps S];
        N=N+1;
        S=S+exponentielle(lambda*u,1,1);
        proc(grille>=S)=proc(grille>=S)+1;
    end
    proc(2:$)=sqrt(t)*(proc(2:$)-lambda*u*grille(2:$))/sqrt(lambda*u);
    plot2d(grille,proc);
endfunction




//////////////////////////
//illustration de la prop de décomposition d'un processus de Poisson (cf exercice 6)
//lambda : intensité du processus
//T : temps final pour la trajectoire
//p : proba de choisir le type A
//////////////////////////

function decomp_poisson(lambda,T,p)
    Tps=[0];
    TpsA=[0];
    TpsB=[0];
    N=0;
    NA=0;
    NB=0;
    S=exponentielle(lambda,1,1);
    while (S<T)
        Tps=[Tps S];
        N=N+1;
        if rand(1)<p then
            TpsA=[TpsA S];
            NA=NA+1;
        else
            TpsB=[TpsB S];
            NB=NB+1;
        end
        S=S+exponentielle(lambda,1,1);
    end
    plot2d2([Tps T],[0:N N]);
    plot2d2([TpsA T],[0:NA NA],2);
    plot2d2([TpsB T],[0:NB NB],3);
endfunction




