//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////    Gradient stochastique avec intéraction : Equation de Langevin avec intéraction champ moyen                        ////
////////    Cas : V(x)=x^2/2 et alpha>0  Simulation de la mesure invariante par un échantiullonnage de Gibbs                             ////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////



clf();
stacksize("max");


N=5;        //Nombre de particules (N=2 => tracé de la densité de pi et histogramme 3d) Nb : N<10
sigma=1;    //Coefficient de diffusion
alpha=1;    //taux d'interaction autour de la moyenne des particules (alpha>0)
beta=1;     //force du gradient sur chaque particule  
T=100;       //Horizon temporel
n=10000;     //Pas de temps du schéma d'Euler T/n


// Echantillonneur de Gibbs de la mesure invariante

m=zeros(N,N+1); 
s=zeros(N,N+1);
X=zeros(1,N);
Gibbs=[X];

T_gibbs=N*1000;  //Iteration de l'algo de Gibbs

for t_gibbs=2:T_gibbs

P=Gibbs(t_gibbs-1,:);

 for i=1:N
    
    for j=1:N
    m(i,j)= P(i)+sum(P(j)-P); 
    s(i,j)=N*sigma/sqrt(2*alpha);
    end
    s(i,N+1)=sigma/(sqrt(2*beta));
    s(i,i)=(N/(N-1))*(sigma/sqrt(2*alpha));
     m(i,i)=(sum(P)-P(i))/(N-1);
     m(i,N+1)=0;
    v=sum((s(i,:)^{-2}));
    mm=sum(s(i,:)^{-2}.*m(i,:))/v;
    Gibbs(t_gibbs,i)=grand(1,1,"nor",mm,sqrt(1/v));
   end
end

// Affichage de l'évolution de l'échantillonnage de Gibbs des deux premières  
// coordonnées.

clf(0);xset("window",0);show_window();

subplot(121);
 plot2d(1:T_gibbs,Gibbs(:,1),5)
subplot(122);
 plot2d(1:T_gibbs,Gibbs(:,2),5)


clf(1);xset("window",1);show_window(); //Histogramme des valeurs 
subplot(121);
histplot(linspace(min(Gibbs(:,1)),max(Gibbs(:,1)),20),Gibbs(:,1),style=2);
subplot(122);
histplot(linspace(min(Gibbs(:,2)),max(Gibbs(:,2)),20),Gibbs(:,2),style=2);







function f=freq2d(echant,bornex,borney)
  kx=length(bornex)-1;
  ky=length(borney)-1;  

for i=1:kx,
    for j=1:ky,
        f(i,j)=length(find((echant(1,:)>bornex(i))&(echant(1,:)<=bornex(i+1))&(echant(2,:)>bornex(j))&(echant(2,:)<=bornex(j+1))))        
    end;
end;
f=f/sum(f);
endfunction


//Histogramme des états visités (x1,x2) par l'échantillonneur de Gibbs

clf(2);xset("window",2);show_window();


bx=[-1.5:.15:1.5];
by=bx;
U=Gibbs(:,1:2);
f=freq2d(U',bx,by);
hist3d(list(f,bx,by), alpha=89.7, alpha=50);



//Affichage de la densité de pi(x1,x2) (à constante près) dans le cas N=2

clf(3);xset("window",3);show_window(); 
x=linspace(-1.5,1.5,31);
y=linspace(-1.5,1.5,31);
z=(exp(-(((2*beta/sigma^2)+(alpha/4))*x^2))'* exp(-(((2*beta/sigma^2)+(alpha/4))*y^2))).*exp((alpha/sigma^2)*x'*y);
plot3d(x,y,z,flag=[2 4 4], alpha=60, theta=50);



// Simulation du modèle de Langevin (simple puit en 0) à champ moyen (cas V(x)=x^2/2) 

X=zeros(1,N); //Initialisation des processus et de sa moyenne empirique


for t=1:n
  
X=[X;X(t,:)-beta*X(t,:)*(T/n)+epsilon*alpha*(mean(X(t,:))-X(t,:))*(T/n)+sigma*grand(1,N,"nor",0,sqrt(T/n))];

end


//AFFICHAGE DES RESULTATS 

clf(4);xset("window",4);show_window();

for NN=1:N
 plot2d((T/n)*[0:n]',X(:,NN),NN)
end



//Histogramme des états (x1,x2) visités par le processus en champ moyen


clf(5);xset("window",5);show_window(); 


bx=[-1.5:.15:1.5];
by=bx;
U=X(:,1:2);
f=freq2d(U',bx,by);
hist3d(list(f,bx,by), alpha=89.7, alpha=50);