close all;
clear all;

% This code presents a simple example in Monte Carlo integration
% Problem: calculate expectation value of h(x) with distribution f(x)
%          y = integral_a-b {h(x)f(x)}
%          interval: a=0, b=1
%          h(x) = x;  f(x) = exp(x)
%
% Analytical solution: y = 1

% Integral will be calculated with 4 different ways:
% (1) Plain MC
% (2) Inverse Transform MC
% (3) Accept - Reject MC
% (4) Importance Sampling MC

% number of points (random variables) to be created
n = 1000;

% number of iterations to improve statistics 
num_iters = 10;

% define functions
h = @(x) x;
f = @(x) exp(x);


% (1) plain MC: create x_i uniform in [0,1]
for j=1:num_iters
sum_plainMC(j) = 0;
for i=1:n
   x = rand(1);
   sum_plainMC(j) = sum_plainMC(j) + h(x)*f(x);
end
sum_plainMC(j) = sum_plainMC(j)/n;
end
y_plainMC = mean(sum_plainMC);
y_plainMC_var = var(sum_plainMC);


% (2) Inverse transform MC: use inverse function 
%  Note: exp(x) is not a proper PDF. To make it in [0,1] we normalize
%        with the value of integral of exp(x) in [0,1]; i.e. (e-1).
%        Then we calculate the inverse function.
f1 = @(x) exp(x)/(exp(1)-1);        % PDF of the distribution
f1inv = @(x) log(1+(exp(1)-1)*x);   % Inverse CDF of the distribution
const = (exp(1)-1)/n;
for j=1:num_iters
sum_invMC(j) = 0;
for i=1:n
   z = rand(1);     % create z_i uniform in [0,1]
   x = f1inv(z);    % calculate x_i
   sum_invMC(j) = sum_invMC(j) + h(x);
end
sum_invMC(j) = const*sum_invMC(j);
end
y_invMC = mean(sum_invMC);
y_invMC_var = var(sum_invMC);


% (3) Accept - Reject MC: 
%   we use as g(x) distribution the uniform in [0,1]
f1 = @(x) exp(x)/(exp(1)-1);        % PDF of the distribution
g_ar = @(x) 1;       % PDF of the g(x) (=U[0,1])
M = 1.6;           % the maximum of f1 in [0,1] is <1.6
const = (exp(1)-1)/n;
for j=1:num_iters
    sum_acrjMC(j) = 0;
    for i=1:n
        u = 2; x = 1;  % initial values
        % reject x_i for which condition is not valid 
        while (u > f1(x)/(M*g_ar(x)))
            x = rand(1);    % draw x_i from g(x) (=U[0,1])
            u = rand(1);    % draw u uniform in [0,1]
        end
        sum_acrjMC(j) = sum_acrjMC(j) + h(x);
    end
    sum_acrjMC(j) = const*sum_acrjMC(j);
end
y_acrjMC = mean(sum_acrjMC);
y_acrjMC_var = var(sum_acrjMC);


% (4) Importance Sampling MC: 
f1 = @(x) exp(x)/(exp(1)-1);        % PDF of the distribution
mu = 1; sigma = 0.5;   % mean, variance of proposal distribution 
% proposal distribution: semi normal [-inf, 1]
g = @(x) 4/(sqrt(2*pi))*exp(-2*(x-1)^2); 
%const = 1/n;
for j=1:num_iters
    sum_isMC(j) = 0; sum_isMC_w(j) = 0; % initialize variables
    for i=1:n
        x = 2; % initial value
        while (x>1)  % draw x1_i from normal keep only values up to 1
            x = normrnd(mu,sigma); 
        end
        if (x>0)
            w0 = f1(x)/g(x); % compute weights
            sum_isMC(j) = sum_isMC(j) + h(x)*w0;
            sum_isMC_w(j) = sum_isMC_w(j) + w0;
        end
    end
% Estimator  1
%    const = (exp(1)-1)/n;
%    sum_impsamMC(j) = const*sum_impsamMC(j);
% Estimator  2
    const = (exp(1)-1);
    sum_isMC(j) = const*sum_isMC(j)/sum_isMC_w(j);
end
y_isMC = mean(sum_isMC);
y_isMC_var = var(sum_isMC);

% print out all results
fprintf ('\n Total Results (Sampling MC):');
fprintf ('\n Plain:         y=%8.6f, var=%9.3e',y_plainMC,y_plainMC_var);
fprintf ('\n Inverse:       y=%8.6f, var=%9.3e',y_invMC,y_invMC_var);
fprintf ('\n Accept-Reject: y=%8.6f, var=%9.3e',y_acrjMC,y_acrjMC_var);
fprintf ('\n Import. Samp.: y=%8.6f, var=%9.3e\n',y_isMC,y_isMC_var);

%stop;
% plot functions
figure(1);
hold on;
fplot(f1,[0 1],'m');
g_ar1 = @(x) g_ar(x)*M;
fplot(g_ar1,[0 1],'r');
xlabel('x'), ylabel('PDF');
legend ('f*(x)','M*g')
title('Accept - Reject Sampling');
hold off;

figure(2);
hold on;
fplot(f1,[-0.5 1],'m');
fplot(g,[-0.5 1],'r');
xlabel('x'), ylabel('PDF');
legend ('f*(x)','g')
title('Importance Sampling');
hold off;
