// Chaîne de Markov

// A partir d'une matrice de transitions.

function transition1=Q(N) // on implémente une matrice de transition
    q=zeros(N+1,N+1) // attention aux tailles des matrices et aux numéros des états
    for j=1:N+1
    q(j,:)=binomial((j-1)/N,N)
    end
    transition1=q
endfunction

function trajectoire=traj(k,N,n)
    q=Q(N)
    trajectoire=[k,grand(n,"markov",q,k+1)-1]; //on simule une trajectoire d'état initial k, de taille n issue de la CM de matrice de transition Q. Dans les programmes, on doit, en entré, respecter l'espace d'état (ie 0,1,...,N) mais faire attention aux indices dans les codes (ie que les états précédents correspondent aux indices 1,...,N+1)).
    // plot2d(trajectoire)
endfunction

function sim2(N,n1,n2) // On simule n1 trajectoires de taille n2 et on regarde la proportion de fois où on arrive effectivement en N
    res=zeros(1,N+1)
    for j=0:N
        compt=0 // le compteur nous assure que les trajectoires considérées sont bien absorbées en 0 ou en N.
        for k=1:n1
            sim=traj(j,N,n2)
            if sim(n2)==N then res(j+1)=res(j+1)+1 ; compt=compt+1
            end
            if sim(n2)==0 then compt=compt+1
            end
        end
        res(j+1)=res(j+1)/compt
    end
    plot2d([0:N],res,5)
    aux=[0:N]/N // valeurs théoriques
    plot([0:N],aux,'B')
    title('Probabilité que X atteigne N','fontsize',5)
    xlabel('état i de départ','fontsize',5)
    ylabel('proportion en N','fontsize',5)
endfunction

function sim3(N,n1,n2) // On calcule par la méthode de Monte-Carlo le temps d'absorption moyen
    res=zeros(1,N+1)
    for j=1:(N-1)
        compt=0
        for k=1:n1
            sim=traj(j,N,n2)
            l=1
            bool=0
            while and([l<n2,bool==0]), // on s'arrête au bout de la chaine ou quand un état absorbant est atteint.
                if or([sim(l)==N,sim(l)==0]) then // on cherche à connaitre le premier indice à partir duquel la chaîne est absorbée en N ou 0.
                    bool=1;res(j+1)=res(j+1)+l;compt=compt+1
                else l=l+1
                end
            end
        end
        res(j+1)=res(j+1)/compt
    end
    plot([1:N-1],res(2:N),'B')
    aux=[1:N-1]/N
    aux=-2*N*(aux.*log(aux)+(1-aux).*log(1-aux)) // valeurs théoriques
    plot([1:N-1],aux,'R')
    title('Méthode de Monte-Carlo','fontsize',5)
    xlabel('état i de départ','fontsize',5)
    ylabel('estimation de mi','fontsize',5)
endfunction

function res=abso(N) // On calcule la valeur théorique de sim3
    P=zeros(N+1,N+1)
    for j=0:N
        P(j+1,:)=binomial(j/N,N)
    end
    P=eye(N+1,N+1)-P
    X=ones(N+1,1)
    res=P\X
    res=res'
endfunction


// A partir d'une chaîne quelconque vérifiant des propriétés de récurrence et de transience.

function fq(n,k)
    H=genmarkov([4,1,5],2) // on génère une matrice de transitions d'une CM à 12 états, qui possède 3 classes récurrentes de cardinaux 4,1 et 5 et 2 états transients.
    t=grand(n,"markov",H,k)
    m=tabul(t) // tabul retoune un tableau à 2 colonnes: la première liste les valeurs prises par la trajectoire et la seconde compte les occurrences de chaque état.
    plot2d3(m(:,1),m(:,2)/n) // plot2d3 permet de tracer les fréquences empiriques de chaque état.
endfunction

Q1=genmarkov([5],0)
[a,b]=linsolve(Q1'-eye(5,5),zeros(5,1)) // linsolve permet de calculer une mesure invariante. Attention a la transposée ! 
pi=b'/sum(b) // Penser à renormaliser pour avoir une mesure de proba
// disp(pi)

function dist(n) // un calcul de distance pour conforter la convergence.
    P=ones(5,1)*pi
    res=zeros(1,n)
    for i=1:n
        res(i)=log(max(abs(Q1^i-P)))/i
    end
    plot2d(res)
endfunction

function w=f(x)
    w=x^(2)
endfunction

function ergo(n,ini)
    t=grand(n,"markov",Q1,ini)
    plot2d([1:n],cumsum(f(t))./ [1:n])
    p=0
    for j=1:5
        p=p+pi(j)*f(j)
    end
    plot(p*ones(1,n)) // Pour tracer une constante, il en fait tracer un tableau ne contenant que la constante.
endfunction

function cvloi(n,ini) // On veut montrer une convergence en loi, on trace donc des histogrammes ou des fréquences empiriques.
    absc=linspace(0.5,5.5,6) // on demande 5 batons répartis entre 0.5 et 5.5 (on centre un baton sur chaque état).
    t=grand(n,"markov",Q1,ini)
    //m=tabul(t)
    //plot2d3(m(:,1),m(:,2)/n,2)
    plot2d3([1:5],pi(:)) // on compare avec la proba invariante.
    histplot(absc,t,5)
endfunction
