function [X,acc] = MHsimple()

% Computer code to generate numbers drawn from a target distribution.
% Uses the MH algorithm.  Inspiration from:
%
% https://stephens999.github.io/fiveMinuteStats/MH_intro.html
% http://physics.ujep.cz/~mmaly/vyuka/poc_fyz_1/zdroje/MeropolisAlg_clanek.pdf

% *************************************************************************
% Numerical parameters:

burnin = 100; % number of burn-in iterations
nsamp = 10000; % number of samples to draw
sig = 5; % standard deviation of Gaussian proposal
x = 1; % start point

% *************************************************************************

% storage
X = zeros(nsamp,1); % samples drawn from the Markov chain
acc = [0 0]; % vector to track the acceptance rate

% *************************************************************************
% Burn-in:

for i = 1:burnin
    [x,a] = MHstep(x,sig); % iterate chain one time step
    acc = acc + [a 1]; % track accept-reject status
end

% *************************************************************************
% MH routine

for i = 1:nsamp
    [x,a] = MHstep(x,sig); % iterate chain one time step
    acc = acc + [a 1]; % track accept-reject status
    X(i) = x; % store the i-th sample
end

end

% *************************************************************************

function [x1,a] = MHstep(x0,sig)
    xp = normrnd(x0,sig,1); % generate candidate from Gaussian
    accprob = targetdist(xp) / targetdist(x0); % acceptance probability
    u = rand; % uniform random number
    if u <= accprob % if accepted
        x1 = xp; % new point is the candidate
        a = 1; % note the acceptance
    else % if rejected
        x1 = x0; % new point is the same as the old one
        a = 0; % note the rejection
    end
end

% *************************************************************************

function probX = targetdist(x)

    lambda=1;
    
    if(x<0)
        probX=0;
    else
        probX = lambda*exp(-lambda*x);
    end
end

