% Presents two lagorithms that simulate a pPoisson random variable

%clear memory and screen
clear all; clc;

%----------
%Parameters
%----------
n_traj = 100000; %number of paths
lambda = 5;     %Expectation of the Poisson distribution

%%
%-----------------------
%First algorithm
%-----------------------
fprintf('First algorithm');
fprintf('\n');

start_time=clock; start_cpu=cputime; %initialize timing variables
N = zeros(n_traj,1);
for i=1:n_traj
    X = -log(rand(1,1))/lambda;
    while X < 1
        N(i,1) = N(i,1) + 1;
        X = X + -log(rand(1,1))/lambda;
    end 
end
% End program and print out the results
moy_N = mean(N)
std_Nbarre = std(N)/sqrt(n_traj)
end_time=clock; end_cpu=cputime;
fprintf('Program executed in %10.5f seconds',etime(end_time,start_time));
fprintf('\n');
fprintf('Using %24.5f CPU time', end_cpu-start_cpu);
fprintf('\n');

%%
start_time=clock; start_cpu=cputime; %initialize timing variables

p = exp(-lambda);
F = p;
N = zeros(n_traj,1);
for i=1:n_traj
    p = exp(-lambda);
    F = p;
    U = rand(1,1);
    while U > F
        N(i,1) = N(i,1) + 1;
        p = p*lambda/N(i,1);
        F = F + p;
    end 
end
% End program and print out the results
moy_N = mean(N)
std_Nbarre = std(N)/sqrt(n_traj)
end_time=clock; end_cpu=cputime;
fprintf('Program executed in %10.5f seconds',etime(end_time,start_time));
fprintf('\n');
fprintf('Using %24.5f CPU time', end_cpu-start_cpu);
fprintf('\n');