close all
clear all
clc
% N=number of MH draws
N=20000;

x=zeros(1,N);
like=zeros(1,N);

%initial x
x(1)=0;

% initial likelihood
like(1)=normpdf(x(1),0,1);

ratio=0;

% M-H loop
for j=2:N,
    % the proposals
    xp=x(j-1)+random('norm',0,2.5);
    likep=normpdf(xp,0,1);
        
    % MH step
    if(random('unif',0,1)<likep/like(j-1)),
        x(j)=xp;
        like(j)=likep;
        ratio=(ratio*(j-1)+1)/j;
    else
        x(j)=x(j-1);
        like(j)=like(j-1);
        ratio=ratio*((j-1)/j);
    end
end

%checking the results...mean should be closed to 0 and std closed to 1
disp('Mean and Std. They should be closed to 0 and 1')
mean=mean(x)
std=std(x)
% ratio should be around 40 % if it is not adjust the variance of the proposal normal
disp('Acceptance Ratio, it should be closed to 40%')
ratio*100

% calculating the cdf probabilities

y(1)=0;
y(2)=0.67;
y(3)=1.28;
y(4)=1.65;
y(5)=2.32;
y(6)=2.58;
y(7)=3.09;
y(8)=3.72;

test=zeros(1,8);

for i=1:8;
    for j=1:N;
        if (x(j)<y(i));
            test(i)=test(i)+1;
        end
    end
    test(i)=test(i)/N;
end

disp('result from MH')
test
disp('Exact numbers')
test0=normcdf(y)
disp('Difference between the exact and the MH results')
abs(test-test0)