//////////////////////////
//Calcul de la mesure invariante d'un processus de Markov a sauts
//Q : generateur infinitesimal du semi-groupe
//mu : loi invariante (vecteur ligne)
//////////////////////////

function [mu]=mesureinvariante(Q)
    mu=kernel(Q');
    mu=mu./(sum(mu));
    mu=mu';
endfunction



//////////////////////////
//evolution temporelle de la loi de Xt en fonction du temps t
//Q : generateur infinitesimal du semi-groupe
//mu0 : loi de X0 (vecteur ligne)
//T : temps final
//////////////////////////

function []=mut(Q,mu0,T)
    n=prod(size(mu0));
    h=T/1000;
    Ph=expm(h.*Q);
    mu=mu0;
    aux=mu0;
    for i=1:1000
        aux=aux*Ph;
        mu=[mu;aux];
    end
    plot2d((0:1000).*h,mu);
    hl=legend(string(1:n))
endfunction



//////////////////////////
//Calcul de la matrice de transition de la chaine incluse
//Q : generateur infinitesimal du semi-groupe
//P : matrice de transition de la chaine incluse
//////////////////////////

function [P]=continu2discret(Q)
    d=size(Q);
    d=d(1,1);
    P=zeros(d,d);
    for i=1:d
        if (Q(i,i)==0) then
            //P(i,:)=zeros(1,d);
            P(i,i)=1;
        else
            P(i,:)=-Q(i,:)./Q(i,i);
            P(i,i)=0;
        end
    end
endfunction




//////////////////////////
//Calcul d'une trajectoire d'un processus de Markov a sauts
// Q : generateur infinitesimal du semi-groupe (matrice d*d)
// x0 : position de depart (entier entre 1 et d)
// T : temps final
// tps : vecteur donnant les temps de saut entre 0 et T
// pos : vecteur donnant les positions apres les sauts
//////////////////////////

function [tps,pos]=trajectoire(Q,x0,T)
    d=size(Q);
    d=d(1,1);
    tps=0;
    pos=x0;
    P=continu2discret(Q);
    if (Q(x0,x0)==0) then
        tpsaux=T+1;
        posaux=x0;
    else
        tpsaux=grand(1,1,"exp",-1/Q(x0,x0));
        posaux=grand(1,"markov",P,x0);
    end
    while (tpsaux<T)
        tps=[tps, tpsaux];
        pos=[pos, posaux];
        if Q(posaux,posaux)==0 then
            tpsaux=T+1;
            posaux=posaux;
        else
            tpsaux=tpsaux+grand(1,1,"exp",-1/Q(posaux,posaux));
            posaux=grand(1,"markov",P,posaux);
        end
    end
    tps=[tps, T];
    pos=[pos, pos($)];
endfunction



//////////////////////////
// fonction affichant une trajectoire d'un processus de Markov
// tps,pos : temps de sauts et positions après les sauts (renvoyes par la fonction trajectoire)
// valeurs : vecteur donnant la correspondance entre les entiers entre 1 et d (taille de l'espace des valeurs du processus) et les valeurs associees
//////////////////////////

function tracetrajectoire(tps,pos,valeurs)
    d=size(Q);
    d=d(1,1);
    for i=1:d;
        nb=sum(1.*(pos==i));
        posnew(pos==i)=valeurs(i).*ones(1,nb);
    end
    plot2d2(tps,posnew);
endfunction





//////////////////////////
//fonction illustrant la convergence en temps du theoreme ergodique vers la mesure invariante
// Q : generateur du semi-groupe (matrice d*d)
// x0 : point de depart (entier entre 1 et d)
//T : temps final
//approxmesureinv : approximation d ela mesure invariante renvoyee par le theoreme ergodique au temps T
//////////////////////////

function [approxmesureinv]=thergodique(Q,x0,T)
    [tps,pos]=trajectoire(Q,x0,T);
    diffpos=tabul(pos,"i");
    diffpos=diffpos(:,1);
    nbdiffpos=prod(size(diffpos));
    T=tps($);
    h=T/1000;
    approximationmesure=[];
    for i=1:1000;
        t=i*h;
        selection=(tps<=t);
        tpsaux=tps(selection);
        tpsaux=[tpsaux t];
        posaux=pos(selection);
        posaux=[posaux posaux($)];
        for j=1:nbdiffpos
          posjaux=(posaux(1:($-1))==diffpos(j));
          tpsjaux=tpsaux(2:$);
          tpsjaux=tpsjaux(posjaux)-tpsaux(posjaux);
          tpsjaux=sum(tpsjaux); 
          approximationmesure(j,i)=tpsjaux/t; 
        end
        
    end
    plot2d((1:1000).*h,approximationmesure');
    hl=legend(string(diffpos));
    approxmesureinv=approximationmesure(:,$)';
endfunction




//////////////
// exemple d'utilisation
//////////////

disp('generateur infinitesimal')
Q=[-0.5 0.5 0 0 0 0;1 -3 0 0 1 1;0 1 -2 0 0 1;0 0 1 -2 1 0;1 1 1 1 -5 1;0 1 0 0 0 -1]
disp('mesure invariante')
mesureinvariante(Q)

clf(0);
scf(0);
T=40;
[tps,pos]=trajectoire(Q,1,T);
tracetrajectoire(tps,pos,(-2):3);


disp('mesure invariante approchee par le theoreme ergodique')
clf(1);
scf(1);
T=1000
thergodique(Q,1,1000)

