% ======================================================
%  Algoritmo di Hastings-Metropolis: 
%  genera una catena di Markov sugli stati {1,2,...,m}
%  con probabilita' limite p_i=C*exp(-i^2/T);
%  matrice R di transizione uniforme.
% ======================================================
%
%   function X=metropolis(n,m,T,X0)
%   n: numero di passi da simulare
%   m: numero di stati della catena
%   T: temperatura nella distribuzione di Boltzmann
%   X0: stato iniziale
%   X: vettore contenente la lista degli stati visitati
% 
% m=100; T=100;
% n=5000;
% X=metropolis(n,m,T,ceil(m*rand));
% stairs(X);
% m=100; T=100; n=1e6;
% X=metropolis(n,m,T,ceil(m*rand));
% j=[1:m];
% P=exp(-j.^2/T); % distribuzione di Boltzmann
% p=P/sum(P); % distribuzione normalizzata
% theta_ex=(j.^8)*p' % valore esatto
% Y=X(k:end);
% abs(mean(Y.^8)-theta_ex)/theta_ex; % errore relativo
function X=metropolis(n,m,T,X0)

X=zeros(1,n+1);
X(1)=X0;
Xn=X0;
for i=1:n
    Y=ceil(m*rand);  % variabile uniforme su {1,2,...,m}
    U=rand;
    if U<min(exp((Xn^2-Y^2)/T),1) % criterio di rigetto
       Xn=Y;
    end
    X(i+1)=Xn;

end

return