MCMC

Doing MCMC is fun! Thanks to KnuthClass

aa

% probability step
if x(n) < 0
    px = 3/5;
else
    px = 2/5;
end


% Metropolis Hastings MCMC

p = min([1, p(y)/p(x(n))]);
r = rand;

n = n+1;
if r < p 
    x(n) = y;
else
    x(n) = x(n-1);
end


%%

% do this 100 times

x(1) = 1;   % start in the low prob region
px(1) = 2/5;

for n = 1:10000
    
    % take a step
    y = (-1)*x(n);
    
    % probability step
    if y < 0
        py = 3/5;
    else
     py = 2/5;
    end
    
    % Metropolis Hastings MCMC
    p = min([1, py/px(n)]);
    r = rand;

    %n = n+1;
    if r < p 
        x(n+1) = y;
        px(n+1) = py;
    else
        x(n+1) = x(n);
        px(n+1) = px(n);
    end
    
    x(n)
    
end
    

%%
clear

% Gaussian probability density with zero mean and unit variance
prob = @(z)((sqrt(2*pi))^(-1)*exp(-z^2/2));

x(1) = 5*(2*rand-1);
px(1) = prob(x);

for n = 1:100000
    step = 0.01*randn;
    y = x(n)+step;
    
    py = prob(y);
    
    p = min([1,py/px(n)]);
    r = rand;

    if r < p 
        x(n+1) = y;
        px(n+1) = py;
    else
        x(n+1) = x(n);
        px(n+1) = px(n);
    end
end

 aa