% 1D Monte carlo integration
% without and with sampling
% Example $\int exp(-x^2) dx from 0 to 1
% the  answer for the integral is 0.74682413

clear  % Clears command window
clc    % Clears command history 
tic
Date = today('datetime');
fprintf(' Date:   %s  \n', Date)
fprintf(' Monte Carlo integration   \n')
fprintf(' Using importance sampling  \n')

rng('default')      % reset RNG
rng('shuffle')      % use seed as current time

Nmax = 100000;
a = 0.0;
b = 1.0;
fprintf('\n N points = %9i  \n', Nmax)

% =====================================
% Method 0: Using Matlab
% =====================================
fun0 = @fint;
I0 = integral(fun0,a,b);
fprintf(' \n MatLab integration \n')
fprintf(' integral = %9.5f \n \n',I0)

% =====================================
% Method 1: the mean value
% =====================================
Sint = 0.0;
fav1 = 0.0;
fav2 = 0.0;
for it=1:Nmax
    x = a+(b-a)*rand;
    Sint = Sint + fint(x);
    fav1 = fav1 + fint(x);
    fav2 = fav2 + (fint(x))^2;
end
Sint = Sint*(b-a)/double(Nmax);
fav1 = fav1/double(Nmax);
fav2 = fav2/double(Nmax);
Serr = (b-a)*sqrt((fav2 - fav1^2)/Nmax);
fprintf(' The mean value \n')
fprintf(' integral = %9.5f \n',Sint)
fprintf(' error    = %9.5f \n',Serr)

% ========================================
% Method 2: Smpling using inverse function
% Sampling function F = A*exp(-x)
% ========================================

% calculating the normalization
global A
fun = @fnorm;
I1 = integral(fun,a,b);
A = 1.0/I1;
%A = 1.5820; % normalization for exp(-x) on [0,1]

Sint = 0.0;
fav1 = 0.0;
fav2 = 0.0;

% getting correct limits of integration for the inverse function
bi = exp(-b);
ai = exp(-a);

for it=1:Nmax
    x = (-1.0)*log(bi+(ai-bi)*rand);
    Sint = Sint + fint2(x);
    fav1 = fav1 + fint2(x);
    fav2 = fav2 + (fint2(x))^2;
end
Sint = Sint/double(Nmax);
fav1 = fav1/double(Nmax);
fav2 = fav2/double(Nmax);
Serr = (b-a)*sqrt((fav2 - fav1^2)/Nmax);
fprintf(' \n Sampling using transform \n')
fprintf(' integral = %9.5f \n',Sint)
fprintf(' error    = %9.5f \n',Serr)

% ==============================================
% Method 3: Sampling using the Metropolis method
% ==============================================

Sint = 0.0;
h = (b-a)/10.;
km = 1;
xm(1) = 0.0;

for k = 1:Nmax
    xtest = xm(km) + h*(2*rand - 1.0);
    xtest = abs(xtest);
    
    if xtest < a || xtest > b
        continue
    end
    r = fint3(xtest)/fint3(xm(km));
    if r < rand
        continue
    else
        km = km + 1;
        xm(km) = xtest;
        Sint = Sint + fint2(xtest);
        fav1 = fav1 + fint2(xtest);
        fav2 = fav2 + (fint2(xtest))^2;
    end
end
%histogram(xm)
Sint = Sint/double(km);
fav1 = fav1/double(km);
fav2 = fav2/double(km);
Serr = (b-a)*sqrt((fav2 - fav1^2)/km);

fprintf(' \n Sampling using Metropolis \n')
fprintf(' integral = %9.5f \n',Sint)
fprintf(' error    = %9.5f \n',Serr)
fprintf(' success  =  %5.2f \n', km/Nmax)

timeE = ceil(toc);
fprintf('\n Elapsed time = %4i sec \n',timeE)

function F = fint(x)
%    F = exp(-x*x);
    F = exp(-x.^2);
end
function F = fint2(x)
global A
    F = exp(-x*x)/(A*exp(-x));
end
function F = fint3(x)
    F = exp(-x);
end

function f = fnorm(x)
f = exp(0.0-x);
end