T = 500; % Numero massimo di iterazioni
sigma = 100; % Deviazione standard della Q(x)= Normal(x | mu, sigma)
xmin = -30; xmax = 30; % range per i valori iniziali di x
x = zeros( 1 , T ); % Buffer per i campioni
seed=1; rand( 'state' , seed ); randn('state',seed ); % fissa il random seed
x(1) = unifrnd( xmin , xmax ); % genera i valori di partenza
%% Inizia campionamento
t = 1;
while t < T % Itera fino a T campioni
    t = t + 1;
    % Propone un nuovo valore per x usando Q(x)=Normale
    x_star = normrnd( x(t-1) , sigma ); x_star = x(t-1) + sigma*randn(1);
    % Calcola l'acceptance ratio
    alpha = min( [ 1 cauchy( x_star ) / cauchy( x(t-1) ) ] );
%     alpha = min( [ 1 normpdf ( x_star ) / normpdf( x(t-1) ) ] );
    % Campiona da Unif[ 0 1 ]
    u = rand;
    % Accettare la proposta?
    if u < alpha
        x(t) = x_star; % SI: setta il nuovo stato
    else
        x(t) = x(t-1); % NO: usa il vecchio stato
    end
end
%% Mostra l'istogramma dei campioni
% Campionamento da una distribuzione di Cauchy usando il metodo di 
% Metropolis. La figura in alto mostra l'istogramma dei campioni (500)
% con sovrapposta la densita teorica; la figura in basso, mostra la sequenza di
% campioni in una catena
figure( 1 ); clf;
subplot( 3,1,1 );
nbins = 200;
xbins = linspace( xmin , xmax , nbins );
counts = hist( x , xbins );
bar( xbins , counts/sum(counts) , 'k' );
xlim( [ xmin xmax ] );
xlabel( 'x' ); ylabel( 'p(x)' );
%% Sovrappone la densita' teorica
y = cauchy( xbins );
hold on;
plot( xbins , y/sum(y) , 'r--' , 'LineWidth' , 3 );
set( gca , 'YTick' , [] );
%% Mostra la storia dei campioni (sequenza di campioni in una catena)
subplot( 3,1,2:3 );
stairs( x , 1:T , 'k-' );
ylabel( 't' ); xlabel( 'x' );
set( gca , 'YDir' , 'reverse' );
xlim( [ xmin xmax ] );