function MCint()
n = 1000000; % the # of samples
figure(1);
clf; hold on;
t = linspace(0,1,200);
plot(t,f1(t),'Linewidth',2,'color','b');
grid;

Iexact = (125 - 125*exp(-5)*cos(1/5) + 5*exp(-5)*sin(1/5))/626;

method = 2;
ntrial = 100;

%% no importance sampling, uniform random variable
if method == 1 
    for k = 1 : ntrial
        u = rand(n,1);
        g = f1(u);
        I(k) = sum(g)/n;
 %       fprintf('n = %d\t I = %d\t Iexact = %d\n',n,I(k),Iexact);
    end
    Imean = sum(I)/ntrial;
    var = sum((I - Imean).^2)/(ntrial - 1);
    sigma = sqrt(var);
    fprintf('ntrial = %d\t Imean = %d\t Imean - Iexact = %d\t var = %d\t sigma = %d\n',...
        ntrial, Imean, Imean - Iexact, var, sigma);
    figure(3);
    hold on; 
    plot([1:ntrial],I - Iexact,'Linewidth',2,'color','b');
end


%% importance sampling, random variable has pdf 5*exp(-5x)/(1 - exp(-5)) on [0,1]
if method == 2
    for k = 1 : ntrial
        u = rand(n,1);
        eta = -log(1 - u*(1 - exp(-5)))/5;
        g = f2(eta);
        fac = (1 - exp(-5))/5;
        I(k) = fac*sum(g)/n;
%        fprintf('n = %d\t I = %d\t Iexact = %d\n',n,I(k),Iexact);
    end
    Imean = sum(I)/ntrial;
    var = sum((I - Imean).^2)/(ntrial - 1);
    sigma = sqrt(var);
    fprintf('ntrial = %d\t Imean = %d\t Imean - Iexact = %d\t var = %d\t sigma = %d\n',...
        ntrial, Imean, Imean - Iexact, var, sigma);
    figure(3);
    hold on; 
    plot([1:ntrial],I - Iexact,'Linewidth',2,'color','m');

end
end

%%
function f = f1(x)
f = cos(x/5).*exp(-5*x);
end

%%
function f = f2(x)
f = cos(x/5);
end
        