################################################################## # Poisson process simulations ################################################################## ################################################################## #simulation of one realization of the Poisson process #event time intervals are exponential with rate beta ################################################################## beta=10; EventIntTimeSize=50 EventTimeIntervals = rexp(EventIntTimeSize, beta) #time is cumulative sum of the interarrival times t_cumsum = cumsum(EventTimeIntervals) t = 1 #Nt-number of events until time t=1 Nt = max(which(t_cumsum< t)) ################################################################# ################################################################## #now simulate several Poisson Pr. realizations #at different discrete time points and make plots #For loop goes through two different sample sizes of the #event time intervals #so we repeat the simulation once with sample size 50 and then #with sample size 100 (the length of event interval times) ################################################################## par(mfrow=c(1,2)) for (EventIntTimeSize in c(50,100)){ EventTimeIntervals = rexp(EventIntTimeSize, beta) #time is cumulative sum of the interarrival times t_cumsum = cumsum(EventTimeIntervals) #plot the Poisson process against time #initalize a vector Nt_time that will keep #Poisson Pr. realizations up to times 1:6 Nt=rep(NA,6) for (t in (1:6)){ Nt[t] = max(which(t_cumsum< t)) } #make a plot of Posson Pr. realizations in times (1:6) plot(Nt,xlab='time',main=paste('Poisson Pr. over time, s. size', EventIntTimeSize),col='purple') } ################################################################## ################################################################# #Next we would like to check if #the length t increment of the PP, N(t+s)-N(s) follows Poisson(beta*t) #((t=t+s-s) is the length of the interval) #i.e the number of event in any interval of length t is Poisson(beta*t) #For that we will need many replicates of the Poisson #process at times t+s and t #For simplicity let t=1 and s=1 #so will we plot the N2-N1 and try to match Poisson(beta) ################################################################# ################################################################# #simulPP : more understandable version to simulate PP # using for loop #inputs: #n - is number of times we are simulating Poisson process #t - vector of length 2 corresponding to times of the # increment N(t2=t+s)-N(t1=s) # default is t=c(1,2) correposnding to t=(t1=1,t2=2) #beta - the rate of event interval times #EventIntTimeSize - how many samples of the event interval times # we would like to simulate #output - a nx2 matrix realizations of Poisson process # the first column of the matrix correponds to # n realizations of the PP up to time t1, # the second column of the matrix correponds to # n realizations of the PP up to time t2 ################################################################# simulPP = function(n,t=c(1,2),beta=10,EventIntTimeSize=500) { #initialize a vector N that will keep #n simlated realizations of the Poisson Process N1=rep(NA,n) N2=rep(NA,n) #Repeat n times #simulate event time intervals and #then find the count of the events up to time n #store it in N[i] for ( i in (1:n)){ EventTimeIntervals = rexp(EventIntTimeSize, beta) #time is cumulative sum of the interarrival times t_cum = cumsum(EventTimeIntervals) #print(i) #number of events until time t N1[i] = max(which(t_cum< t[1])) N2[i] = max(which(t_cum< t[2])) } #N1 and N2 are two vectors of the same length so we can form #a matrix N of two columns by cbind #N1 n realizations of the Poisson process up to time t1 #N2 n realizations of the Poisson process up to time t2 N = cbind(N1,N2) return(N) } ################################################################## # Simulations and plots ################################################################## #Simulate nx2 realizations of Poisson process up to time t=(t1=1,t2=2) pp = simulPP(n=5000) #plot PP against iterations plot(pp[,1],main='Realizations of Poisson process up to time 1') plot(pp[,2],main='Realizations of Poisson process up to time 2') #Plot N2-N1 and see how well it matches the Poisson(beta) (true) distribution par(mfrow=c(1,1)) N21=pp[,2]-pp[,1] hist(N21,main='Number of events in the interval [t1,t2] of length t2-t1=1',prob=TRUE ,xlab='N[t2=2]-N[t1=1]') lines(density(rpois(5000, beta)),col='red') ################################################################## #Simulate n realizations of Poisson process all up to time t=(t2=2,t3=3) pp=simulPP(n=5000,t=c(2,3)) par(mfrow=c(1,1)) N32=pp[,2]-pp[,1] hist(N32,main='Number of events in the interval [t2,t3] of length t3-t2=1',prob=TRUE, xlab='N[t3=3]-N[t2=2]') lines(density(rpois(5000, beta)),col='red') ################################################################## #Simulate n realizations of Poisson process all up to time t=(t3=3,t4=4) pp=simulPP(n=5000,t=c(3,4)) par(mfrow=c(1,1)) N43=pp[,2]-pp[,1] hist(N43,main='Number of events in the interval [t3,t4] of length t4-t3=1',prob=TRUE, xlab='N[t4=4]-N[t3=3]') lines(density(rpois(5000, beta)),col='red') ################################################################## ################################################################## #Simulate nx2 realizations of Poisson process all up to time t=(t4=4,t5=5) pp=simulPP(n=5000,t=c(4,5)) par(mfrow=c(1,1)) N54=pp[,2]-pp[,1] hist(N54,main='Number of events in the interval [t4,t5] of length t5-t4=1',prob=TRUE, xlab='N[t5=5]-N[t4=4]') lines(density(rpois(5000, beta)),col='red') ################################################################## ################################################################## #This is more advanced #and efficient version of the simulPP #that uses vectorized for by utilizing apply function #inputs and outputs are the same as simulPP ################################################################## simulPP_adv = function(n,t=c(1,2),beta=10,EventIntTimeSize=500) { #each row of the matrix x corresponds to #a single realization of the Poisson process x = matrix(rexp(n*EventIntTimeSize,beta),n,EventIntTimeSize) #call a wrapper function apply that executes #simulateN on each row (1, if we wanted the apply to #each column we would supply 2) of the matrix x #the result is number of event until time t #per each row of x (there are n rows in x) #i.e. n realizations of Poissson Pr. N1 = apply(x,1,simulateN,t=t[1]) N2 = apply(x,1,simulateN,t=t[2]) #N1 and N2 are two vectors of the same length so we can form #a matrix N of two columns by cbind #N1 n realizations of the Poisson process up to time t1 #N2 n realizations of the Poisson process up to time t2 N = cbind(N1,N2) return(N) } ############################################################# ################################################################# #simulateN: Simulate count of events up to time t # given the event interval times #This function is needed for the apply wrapper function #in the simulPP_adv function #input: #x - vector of event interval times #t - time in N(t) #output : one realization of Possion Pr. ################################################################# simulateN = function(x,t=1) { t_cumsum=cumsum(x) N=max(which(t_cumsum