% Il modo più semplice per simulare un processo di Poisson consiste nel
% campionare le variabili esponenziali Xi, una alla volta e arrestarsi 
% non appena Sn = somma_i=1,...,n Xi supera l’istante finale T; 
% questo può essere fatto con un semplice ciclo while.
% Il processo è interamente descritto dalla sequenza dei tempi di arrivo Si, 
% e dal numero di eventi verificatesi tra 0 e gli istanti Si (il numero di 
% eventi è una funzione crescente che incrementa di 1 ogni volta che si 
% verifica un evento). E' pertanto opportuno memorizzare un
% vettore S contenete tutti i tempi di arrivo e un vettore N contenente il 
% numero cumulativo di eventi: N(i) = {# eventi tra 0 e S(i)}. 
%% PER LANCIARE LA FUNCTION
% lambda=2; T=30;
% La visualizzazione della funzione N(t) può essere ottenuta con il
% comando (La figura mostra la realizzazione del processo di Poisson ottenuta)
% stairs(S,N)
% title(['Processo di Poisson con rate ' num2str(lambda) ' in [0,' num2str(T) ']'])
% xlabel('t'); ylabel('numero di successi');

%% [S N] = poisson(lambda, T)
%   lambda: parametro del processo (tasso di arrivo)
%   T:   istante finale
%   S:      vettore contenente i tempi di arrivo -- S(1) = 0,  S(end)=T
%   N:      contatore degli eventi -- N(1) = 0, N(end): # eventi in [0,T]

function [S N] = poisson(lambda, T)

% estraggo il primo incremento (quindi il primo punto) da distribuzione
% esponenziale; 1/lambda è il valor medio della distribuzione esponenziale
S(1) = exprnd(1/lambda); % o X=-log(rand)/lambda; -> campionamento di esponenziale
i = 1;
% continuo ad estrarre nuovi incrementi finché non ho superato la lunghezza
% dell'intervallo
N=0;    j=0;
while S(i)<T
    i = i+1;    
    S(i) = S(i-1) + exprnd(1/lambda);   % nuovo tempo di arrivo
    j = j + 1;                          % nuovo evento
    N=[N, j];                           % aggiungi evento alla lista N
end

S(end) = T;
N(end) = N(end-1);

