

%%

% define two test integrands
f1 = @(x) x + sin(x+pi/2).^2 + cos(4*x+1.6*pi/2) ; 
f2 = @(x) 1 + cos(x).^2 + sin(4*x) ;
% Reference integrals of f1 and f2 in [0,1] using Wolfram 
I1 = 1.13739 ; 
I2 = 2.1407 ; 


% define importance function (e.g. normal distribution)
% we truncate the domain to [0,1]
mu = 0.3;
sig = 0.33; 
gdist = makedist('Normal', mu, sig) ;
tdistr = truncate(gdist, 0,1) ;
g = @(x) tdistr.pdf(x) ;

% Use N samples to approximate integrals
% use iter iterations to estimate variance of N-sample estimator
N = 50; 
iter = 50;

% integrate using random 
xr = rand(N,iter) ;
Ir1 = mean(f1(xr)) ;
Ir2 = mean(f2(xr)) ;

% integrate using importance sampling
% sample random values from tdistr 
xn = random(tdistr, ns, iter); 
In1 = mean(f1(xn)./g(xn)) ;
In2 = mean(f2(xn)./g(xn)) ;


% repeating this for many different numbers of samples 
% so that we can plot convergence as n increases
% this time n is the number of samples
Ns = [5, 100, 1500, 3000] ;
% Ns = [] ;
err1 = zeros(2, length(Ns)) ;
err2 = zeros(2, length(Ns)) ;
var1 = zeros(2, length(Ns)) ;
var2 = zeros(2, length(Ns)) ;

for i=1:length(Ns) 
    ns = Ns(i) ;
    
    % integrate using random 
    xra = rand(ns,iter) ;    
    Ir1a = mean(f1(xra)) ;
    Ir2a = mean(f2(xra)) ;

    % integrate using importance sampling
    % xna are generated from truncated distribution
    xna = random(tdistr, ns, iter); 
    Nfound(i) = length(xna) ;
    In1a = mean(f1(xna)./g(xna)) ;
    In2a = mean(f2(xna)./g(xna)) ;

    err1(1,i) = sqrt((mean(Ir1a)-I1).^2) ; % f1 error - random
    err1(2,i) = sqrt((mean(In1a)-I1).^2) ; % f1 error - IS
    
    err2(1,i) = sqrt((mean(Ir2a)-I2).^2) ; % f2 error - random
    err2(2,i) = sqrt((mean(In2a)-I2).^2) ; % f2 error - IS       

    var1(1,i) = var(Ir1a) ; % f1 variance - random
    var1(2,i) = var(In1a)  ; % f1 variance - IS
    var2(1,i) = var(Ir2a) ; % f2 variance - random
    var2(2,i) = var(In2a)  ; % f2 variance - IS
end


% plot the functions
figure(1);
m = 5;
n = 2; 
Nd = 100; 
x = linspace(0,1,Nd) ;
subplot(m,n,1); 
plot(x,f1(x), 'b') ;
hold on;
plot(x, normpdf(x,mu, sig), ':k') ;
hold off; 
axis([0,1,0,4]) ;
title(['integrand f1'])
subplot(m,n,2); 
plot(x,f2(x), 'r') ;
hold on;
plot(x, normpdf(x,mu, sig), ':k') ;
hold off; 
axis([0,1,0,4]) ;
title(['integrand f2'])

% plot errors ---- random
subplot(m,n,3) ;
errorbar(N, mean(Ir1), sqrt(var(Ir1)), '.b') ;
hold on; 
    plot(N, I1, '*b') ;
hold off;
axis([0,2*N,0,3]) ;
title(['using ' num2str(N) ' random samples'])
subplot(m,n,4); 
errorbar(N, mean(Ir2), sqrt(var(Ir2)), '.r') ;
hold on; 
    plot(N, I2, '*r') ;
hold off;
axis([0,2*N,0,3]) ;
title(['using ' num2str(N) ' random samples'])

% plot errors ---- importance sampled
subplot(m,n,5) ;
errorbar(N, mean(In1), sqrt(var(In1)), '.b') ;
hold on; 
    plot(N, I1, '*b') ;
hold off;
axis([0,2*N,0,3]) ;
title(['using ' num2str(N) ' importance samples'])
subplot(m,n,6); 
errorbar(N, mean(In2), sqrt(var(In2)), '.r') ;
hold on; 
    plot(N, I2, '*r') ;
hold off;
axis([0,2*N,0,3]) ;
title(['using ' num2str(N) ' importance samples'])

% plot error convergence ---- importance sampled
subplot(m,n,7) ;
loglog(Ns, err1(1,:), 'k') ;
hold on;
loglog(Ns, err1(2,:), '--k') ;
hold off;
legend({'rand', 'is'});
ylabel('log-error') ;
xlabel('log-#samples') ;
title(['convergence plot'])

subplot(m,n,8) ;
loglog(Ns, err2(1,:), 'k') ;
hold on;
loglog(Ns, err2(2,:), '--k') ;
hold off;
legend({'rand', 'is'});
ylabel('log-error') ;
xlabel('log-#samples') ;
title(['convergence plot'])

% plot variance convergence ---- importance sampled
subplot(m,n,9) ;
loglog(Ns, var1(1,:), 'k') ;
hold on;
loglog(Ns, var1(2,:), '--k') ;
hold off;
legend({'rand', 'is'});
ylabel('log-var') ;
xlabel('log-#samples') ;
title(['convergence plot'])

subplot(m,n,10) ;
loglog(Ns, var2(1,:), 'k') ;
hold on;
loglog(Ns, var2(2,:), '--k') ;
hold off;
legend({'rand', 'is'});
ylabel('log-var') ;
xlabel('log-#samples') ;
title(['convergence plot'])

% plot sample pdfs
figure(2) ;
subplot(2,1,1); 
nbins = 10; 
hr = zeros(1,nbins) ;
for i = 1:N
    hr = hr + histcounts(xr(i,:), nbins) ;
end
bar(linspace(0,1,nbins),hr); 
axis([0,1,0,1.2*N*iter/nbins])

subplot(2,1,2); 
nbins = 10; 
hn = zeros(1,nbins) ;

for i = 1:N
    hn = hn + histcounts(xn(i,:), nbins) ;
end
bar(linspace(0,1,nbins),hn); 
axis([0,1,0,2*N*iter/nbins])

figure(1) ;

