%% Problem Sheet 10
% Stochastic Simulation

close all; clear all; clc;

%% Exercise 3. a) 1.

lambda = @(t) 3*(cos(t)+1);
lambda_star = 6;
t_max = 10;
interarrivals = [];

% homogeneous Poisson process
mysum = 0;
while mysum < t_max
    X = -(1/lambda_star)*log(rand);
    mysum = mysum + X;
    interarrivals = [interarrivals; X];
end
% remove last jump (after t_max)
interarrivals = interarrivals(1:end-1);
arrivals = cumsum(interarrivals);

% thinning
keep = rand(length(arrivals), 1) <= lambda(arrivals)/lambda_star;
arrivals = arrivals(keep==1);

t = [0; arrivals; t_max];
N_t = 0:length(arrivals);

figure(1);
clf;
hold on;
handles = zeros(1,3);
for i = 1:(length(t)-1);
    handles(1) = line([t(i), t(i+1)], [N_t(i), N_t(i)]);
end;
x = linspace(0, 10, 100);
handles(2) = plot(x, lambda(x), 'green');
handles(3) = plot(arrivals, zeros(length(arrivals),1), 'ro'); % (not asked)
hold off;
title({'realization of an inhomogeneous Poisson process', ...
       'simulation using thinning'});
legend(handles, {'Poisson process', 'rate function', 'arrival times'}, 'Location', 'NorthWest');


%% Exercise 3: a) 2.
% (inhomogeneous Poisson process, grid)

h = 0.01;
t = 0:h:t_max; % from 0 to t_max by h

arrivals = (rand(length(t)-1, 1) <= lambda(t(2:end))'*h);
N_t = [0; cumsum(arrivals)];

figure(2);
clf;
handles = zeros(1,3);
for i=0:max(N_t)
    handles(1) = line([min(t(N_t==i)), max(t(N_t==i))], [i, i]);
end
hold on;
x = linspace(0, 10, 100);
handles(2) = plot(x, lambda(x), 'green');
handles(3) = plot(t(arrivals==1)+h, zeros(sum(arrivals),1), 'ro');
hold off;
title({'realization of an inhomogeneous Poisson process', ...
       'simulation using a grid'});
legend(handles, {'Poisson process', 'rate function', 'arrival times'}, 'Location', 'NorthWest');


%% Exercise 3: b)
% (inhomogeneous Poisson process, estimation)

lambda = @(t) 3*(cos(t)+1);
lambda_star = 6;
t_max = 3*pi;

N = 10^4;

N_3pi = zeros(N,1);
for i=1:N
    
    % simulate a homogeneous Poisson process
    interarrivals = [];
    mysum = 0;
    while mysum < t_max
        X = -(1/lambda_star)*log(rand);
        mysum = mysum + X;
        interarrivals = [interarrivals; X];
    end
    arrivals = cumsum(interarrivals);
    % remove last arrival (after t_max)
    arrivals = arrivals(1:end-1);
    
    % thin it
    keep = rand(length(arrivals), 1) <= lambda(arrivals)/lambda_star;
    N_3pi(i) = sum(keep==1);
end

exp_est = mean(N_3pi);
exp_true = 9*pi;
RE_est = std(N_3pi)/(sqrt(N)*exp_true);

fprintf('The estimated expecation is %f,\n', exp_est);
fprintf('the true expectation is %f\n', exp_true);
fprintf('and the estimated relative error is %f\n\n', RE_est);


%% Exercise 4 a)

close all; clear all;

lambda = 1.5;
df = 4;
t_max = 5;

% generate the homogeneous PP (when a fish bites)
interarrivals = [];
mysum = 0;
while mysum < t_max
    X = -(1/lambda)*log(rand);
    mysum = mysum + X;
    interarrivals = [interarrivals; X];
end
arrivals = cumsum(interarrivals(1:end-1));

% generate the jump sizes (weight of the fish)
jumpsizes = (1/4)*chi2rnd(df, length(arrivals), 1);

t = [0; arrivals; t_max];
X_t = [0; cumsum(jumpsizes)];

% plot the realization
figure(1);
clf;
for i=1:length(t)-1
    line([t(i), t(i+1)], [X_t(i), X_t(i)]);
end


%% Exercise 4 b)

N = 10^4;
morethan6 = zeros(N,1);

for i=1:N
    % generate the homogeneous PP
    interarrivals = [];
    mysum = 0;
    while mysum < t_max
        X = -(1/lambda)*log(rand);
        mysum = mysum + X;
        interarrivals = [interarrivals; X];
    end
    arrivals = cumsum(interarrivals(1:end-1));
    % generate the jump sizes
    jumpsizes = (1/4)*chi2rnd(df, length(arrivals), 1);
    morethan6(i) = (sum(jumpsizes)>6);
end

fprintf('The estimated probability is %.4f\n', mean(morethan6));