###########################################
####   Poisson Process
###########################################
###  Graficas 11.50 x 6.95 inches // 7.3 6.95

### Exponencial Variables  

xexpo <- 1:600/100
dexpo1 <- dexp(xexpo,0.5)
dexpo2 <- dexp(xexpo,1)
dexpo3 <- dexp(xexpo,2)
dexpo4 <- dexp(xexpo,5)
dfexpo1 <- pexp(xexpo,0.5)
dfexpo2 <- pexp(xexpo,1)
dfexpo3 <- pexp(xexpo,2)
dfexpo4 <- pexp(xexpo,5)
par(mfrow=c(1,2))
clrs <- rainbow(50)
plot(xexpo,dexpo1,type='l', main='Exponential Density',ylim=c(0,5),
     lwd=3,col=clrs[49],ylab = 'density',xlab = '')
lines(xexpo,dexpo2,col=clrs[1],lwd=3,lty=2)
lines(xexpo,dexpo3,col=clrs[3],lwd=3,lty=3)
lines(xexpo,dexpo4,col=clrs[5],lwd=3,lty=4)
legend('topright',c('0.5','1','2','5'),lty= 1:4,
       lwd=rep(3,4),col=c(clrs[49],clrs[1],clrs[3],clrs[5]))
plot(xexpo,dfexpo1,type='l', main='Exponential Distribution Function',ylim=c(0,1),
     lwd=3,col=clrs[49],ylab='distribution function',xlab='')
lines(xexpo,dfexpo2,col=clrs[1],lwd=3,lty=2)
lines(xexpo,dfexpo3,col=clrs[3],lwd=3,lty=3)
lines(xexpo,dfexpo4,col=clrs[5],lwd=3,lty=4)
legend('bottomright',c('0.5','1','2','5'),lty=1:4,
       lwd=rep(2,4),col=c(clrs[49],clrs[1],clrs[3],clrs[5]))
par(mfrow=c(1,1))

## Variables with Gamma distribution, rate=1
xgam <- 1:1000/100
dgam1 <- dgamma(xgam,0.5);dfgam1 <- pgamma(xgam,0.5)
dgam2 <- dgamma(xgam,1); dfgam2 <- pgamma(xgam,1)
dgam3 <- dgamma(xgam,2); dfgam3 <- pgamma(xgam,2)
dgam4 <- dgamma(xgam,4); dfgam4 <- pgamma(xgam,4)
dgam5 <- dgamma(xgam,6); dfgam5 <- pgamma(xgam,6)
dgam6 <- dgamma(xgam,8); dfgam6 <- pgamma(xgam,8)
par(mfrow=c(1,2))
plot(xgam,dgam1,type='l', main='Gamma Density',ylim=c(0,1),
     lwd=3,col=clrs[29],ylab = 'density',xlab = '')
lines(xgam,dgam2,col=clrs[31],lwd=3,lty=2)
lines(xgam,dgam3,col=clrs[33],lwd=3,lty=3)
lines(xgam,dgam4,col=clrs[35],lwd=3,lty=4)
lines(xgam,dgam5,col=clrs[37],lwd=3,lty=5)
lines(xgam,dgam6,col=clrs[39],lwd=3,lty=6)
legend('topright',c('k=0.5','k=1','k=2','k=4','k=6','k=8'),
       lty= 1:6,lwd=rep(3,6),
       col=c(clrs[29],clrs[31],clrs[33],clrs[35],clrs[37],clrs[39]))

library(latex2exp)
plot(xgam,dfgam1,type='l', main='Gamma Distribution Function',ylim=c(0,1),
     lwd=3,col=clrs[29],ylab = 'density',xlab = '')
lines(xgam,dfgam2,col=clrs[31],lwd=3,lty=2)
lines(xgam,dfgam3,col=clrs[33],lwd=3,lty=3)
lines(xgam,dfgam4,col=clrs[35],lwd=3,lty=4)
lines(xgam,dfgam5,col=clrs[37],lwd=3,lty=5)
lines(xgam,dfgam6,col=clrs[39],lwd=3,lty=6)
legend('bottomright',c('k=0.5','k=1','k=2','k=4','k=6','k=8'),
       lty= 1:6,lwd=rep(3,6),
       col=c(clrs[29],clrs[31],clrs[33],clrs[35],clrs[37],clrs[39]))

###  Gamma distribution for different values of k and lambda
dgam11 <- dgamma(xgam,2,0.5);dfgam11 <- pgamma(xgam,2,0.5)
dgam12 <- dgamma(xgam,2,1);  dfgam12 <- pgamma(xgam,2,1)
dgam13 <- dgamma(xgam,2,2);  dfgam13 <- pgamma(xgam,2,2)
dgam14 <- dgamma(xgam,2,4);  dfgam14 <- pgamma(xgam,2,4)
dgam15 <- dgamma(xgam,5,0.5); dfgam15 <- pgamma(xgam,5,0.5)
dgam16 <- dgamma(xgam,5,1);  dfgam16 <- pgamma(xgam,5,1)
dgam17 <- dgamma(xgam,5,2);  dfgam17 <- pgamma(xgam,5,2)
dgam18 <- dgamma(xgam,5,4);  dfgam18 <- pgamma(xgam,5,4)
plot(xgam,dgam11,type='l', main='Gamma Density',ylim=c(0,1.4),
     lwd=3,col=clrs[29],ylab = 'density',xlab = '')
lines(xgam,dgam12,col=clrs[31],lwd=3,lty=2)
lines(xgam,dgam13,col=clrs[33],lwd=3,lty=3)
lines(xgam,dgam14,col=clrs[35],lwd=3,lty=4)
lines(xgam,dgam15,col=clrs[47],lwd=3,lty=1)
lines(xgam,dgam16,col=clrs[49],lwd=3,lty=2)
lines(xgam,dgam17,col=clrs[1],lwd=3,lty=3)
lines(xgam,dgam18,col=clrs[3],lwd=3,lty=4)
legend('topright',c(TeX('k=2 $\\lambda$=0.5$'),TeX('k=2 $\\lambda$=1$'),
                    TeX('k=2 $\\lambda$=2$'),TeX('k=2 $\\lambda$=4$'),
                    TeX('k=5 $\\lambda$=0.5$'),TeX('k=5 $\\lambda$=1$'),
                    TeX('k=5 $\\lambda$=2$'),TeX('k=5 $\\lambda$=4$')),
       lty= rep(1:4,2),lwd=rep(3,8),
       col=c(clrs[29],clrs[31],clrs[33],clrs[35],clrs[47],clrs[49],clrs[1],clrs[3]))

plot(xgam,dfgam11,type='l', main='Gamma Distribution Function',ylim=c(0,1),
     lwd=3,col=clrs[29],ylab = 'distribution function',xlab = '')
lines(xgam,dfgam12,col=clrs[31],lwd=3,lty=2)
lines(xgam,dfgam13,col=clrs[33],lwd=3,lty=3)
lines(xgam,dfgam14,col=clrs[35],lwd=3,lty=4)
lines(xgam,dfgam15,col=clrs[47],lwd=3,lty=1)
lines(xgam,dfgam16,col=clrs[49],lwd=3,lty=2)
lines(xgam,dfgam17,col=clrs[1],lwd=3,lty=3)
lines(xgam,dfgam18,col=clrs[3],lwd=3,lty=4)
legend('bottomright',c(TeX('k=2 $\\lambda$=0.5$'),TeX('k=2 $\\lambda$=1$'),
                       TeX('k=2 $\\lambda$=2$'),TeX('k=2 $\\lambda$=4$'),
                       TeX('k=5 $\\lambda$=0.5$'),TeX('k=5 $\\lambda$=1$'),
                       TeX('k=5 $\\lambda$=2$'),TeX('k=5 $\\lambda$=4$')),
       lty= rep(1:4,2),lwd=rep(3,8),
       col=c(clrs[29],clrs[31],clrs[33],clrs[35],clrs[47],clrs[49],clrs[1],clrs[3]))



library(RColorBrewer) # For setting color scales
blues <-  brewer.pal(9,"Blues")
reds <-  brewer.pal(9,"Reds")
greens <-  brewer.pal(9,"Greens")
greys <-  brewer.pal(9,"Greys")
oranges <- brewer.pal(9,"Oranges")
###########################################


### Sum of exponencial variables
### Start by fixing the number of events,  
### which is the number of exponential variables

set.seed(999)
par(mfrow=c(3,1))
nev <- 20  # Number of events
tTimes <- c(0,cumsum(rexp(nev,1)))
plot(tTimes,0:nev, type='s',xlim=c(0,25),cex.lab=1.5,
     xlab='Time',ylab='Events',lwd=2, cex.main=1.8,
     main='Poisson process with intensity 1')
abline(h=0,col='red')

tTimes <- c(0,cumsum(rexp(nev,0.5)))
plot(tTimes,0:nev, type='s',xlim=c(0,25),cex.lab=1.5,
     xlab='Time',ylab='Events',lwd=2,cex.main=1.8,
     main='Poisson process with intensity 0.5')
abline(h=0,col='red')

tTimes <- c(0,cumsum(rexp(nev,0.1)))
plot(tTimes,0:nev, type='s',xlim=c(0,25),cex.lab=1.5,
     xlab='Time',ylab='Events',lwd=2,cex.main=1.8,
     main='Poisson process with intensity 0.1')
abline(h=0,col='red')

set.seed(999)
par(mfrow=c(3,1))
tTimes <- c(0,cumsum(rexp(nev,1)))
plot(tTimes,0:nev, type='s',ylim=c(0,nev),cex.lab=1.5,
     xlab='Time',ylab='Events',lwd=2,cex.main=1.8,
     main='Poisson process with intensity 1')
abline(h=0,col='red')

tTimes <- c(0,cumsum(rexp(nev,0.5)))
plot(tTimes,0:nev, type='s',ylim=c(0,nev),cex.lab=1.5,
     xlab='Time',ylab='Events',lwd=2,cex.main=1.8,
     main='Poisson process with intensity 0.5')
abline(h=0,col='red')

tTimes <- c(0,cumsum(rexp(nev,0.1)))
plot(tTimes,0:nev, type='s',ylim=c(0,nev),cex.lab=1.5,
     xlab='Time',ylab='Events',lwd=2,cex.main=1.8,
     main='Poisson process with intensity 0.1')
abline(h=0,col='red')
par(mfrow=c(1,1))

### Graphs on the same panel
set.seed(111)
nev <- 40
tTimes <- c(0,cumsum(rexp(nev,1)))
plot(tTimes,0:nev, type='s',ylim=c(0,nev),cex.lab=1.4,
     xlab='Time',ylab='Events',lwd=2,cex.main=1.7,
     main='Poisson process with intensity 1',col=blues[9])
abline(h=0,col='red')
for (i in 2:5)
  lines(c(0,cumsum(rexp(nev,1))),0:nev,type='s',
        lwd=2,col=blues[10-i])

tTimes <- c(0,cumsum(rexp(nev,0.5)))
plot(tTimes,0:nev, type='s',ylim=c(0,30),cex.lab=1.4,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,60),cex.main=1.7,
     main='Poisson process with intensity 0.5',col=reds[9])
abline(h=0,col='red')
for (i in 2:5)
  lines(c(0,cumsum(rexp(nev,0.5))),0:nev,type='s',
        lwd=2,col=reds[10-i])

tTimes <- c(0,cumsum(rexp(nev,0.1)))
plot(tTimes,0:nev, type='s',ylim=c(0,30),cex.lab=1.4,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,250),cex.main=1.7,
     main='Poisson process with intensity 0.1',col=oranges[9])
abline(h=0,col='red')
for (i in 2:5)
  lines(c(0,cumsum(rexp(nev,0.1))),0:nev,type='s',
        lwd=2,col=oranges[10-i])

### Poisson distributon + Uniform
### Another method is to generate a Poisson variable
### with parameter lambda*t
### to get the total number of events 
### and then place these events using  
### a uniform  distribution.
### In this scheme, the time we want to simulate
### is fixed.

set.seed(1235)
par(mfrow=c(3,1))
tt <- 30  # Time of observation
lbd <- 1  # Intensity
nev <- rpois(1,tt*lbd) # Number of events 
tTimes <-c(sort(runif(nev,max=tt)),tt) # Position of the events
plot(c(0,tTimes),c(0:nev,nev), type='s',cex.lab=1.5,cex.main=1.7,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),ylim=c(0,30),
     main='Poisson process with intensity 0.5')
abline(h=0,col='red')

lbd <- 0.5
nev <- rpois(1,tt*lbd)
tTimes <-c(sort(runif(nev,max=tt)),tt)
plot(c(0,tTimes),c(0:nev,nev), type='s',cex.lab=1.5,cex.main=1.7,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),ylim=c(0,30),
     main='Poisson process with intensity 0.5')
abline(h=0,col='red')

lbd <- 0.1
nev <- rpois(1,tt*lbd)
tTimes <-c(sort(runif(nev,max=tt)),tt)
plot(c(0,tTimes),c(0:nev,nev), type='s',cex.lab=1.5,cex.main=1.7,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),ylim=c(0,30),
     main='Poisson process with intensity 0.1')
abline(h=0,col='red')
par(mfrow=c(1,1))

### Graphs on the same panel
set.seed(1357)
tt <- 30
lbd <- 1
nev <- rpois(1,tt*lbd)
tTimes <-c(sort(runif(nev,max=tt)),tt)
plot(c(0,tTimes),c(0:nev,nev), type='s',cex.lab=1.5, cex.main=1.7,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),
     main='Poisson process with intensity 1',col=blues[9])
abline(h=0,col='red')
for (i in 2:5)
{
  nev <- rpois(1,tt*lbd)
  lines (c(0,c(sort(runif(nev,max=tt)),tt)),c(0:nev,nev),type='s',col=blues[10-i],lwd=2)
}
  
lbd <- 0.5
nev <- rpois(1,tt*lbd)
tTimes <-c(sort(runif(nev,max=tt)),tt)
plot(c(0,tTimes),c(0:nev,nev), type='s',cex.lab=1.5, cex.main=1.7,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),
     main='Poisson process with intensity 0.5',col=reds[9])
abline(h=0,col='red')
for (i in 2:5)
{
  nev <- rpois(1,tt*lbd)
  lines (c(0,c(sort(runif(nev,max=tt)),tt)),c(0:nev,nev),type='s',col=reds[10-i],lwd=2)
}

lbd <- .1
nev <- rpois(1,tt*lbd)
tTimes <-c(sort(runif(nev,max=tt)),tt)
plot(c(0,tTimes),c(0:nev,nev), type='s',cex.lab=1.5, cex.main=1.7,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),
     main='Poisson process with intensity 0.1',col=oranges[9])
abline(h=0,col='red')
for (i in 2:5)
{
  nev <- rpois(1,tt*lbd)
  lines (c(0,c(sort(runif(nev,max=tt)),tt)),c(0:nev,nev),type='s',col=oranges[10-i],lwd=2)
}


### Compound Poisson Process
### To each event of a Poisson process 
### we associate a N(0,3) normal variable
nev <- 25
lbd <- 0.5
set.seed(999111)
tTimes <- c(0,cumsum(rexp(nev,lbd)))
pcomp <- c(0,cumsum(rnorm(nev,sd=sqrt(3))))
plot(tTimes,pcomp,type='s',lwd=2,xlab='Times',xlim=c(0,35),
     main='Compound Poisson Process',ylab='',
     ylim=c(-15,15),cex.lab=1.5, cex.main=1.7)
abline(h=0,col='red')
for (i in 2:5)
{
  tTimes <- c(0,cumsum(rexp(nev,lbd)))
  pcomp <- c(0,cumsum(rnorm(nev,sd=sqrt(3))))
  lines(tTimes,pcomp,type='s',col=i,lwd=2)
}


### Decomposition of a Poisson process
nev <- 25
lbd <- 0.5
set.seed(9119)
tTimes <- cumsum(rexp(nev,lbd))
plot(c(0,tTimes),0:nev,type='s',lwd=2,xlab='Times',
     main='Decomposition of a Poisson process',
     ylab='',cex.main=1.7,cex.lab=1.5)
abline(h=0,col='red')
mark <- rbinom(nev,1,0.25)
lines(c(0,tTimes[mark==1],tTimes[nev]),
      c(0:length(mark[mark==1]),length(mark[mark==1])),
        type='s',lwd=2,col='darkblue')
lines(c(0,tTimes[mark==0],tTimes[nev]),
      c(0:length(mark[mark==0]),length(mark[mark==0])),
      type='s',lwd=2,col='red')

### Superposition of Poisson processes.
set.seed(8264)
tt <- 30  # Simulation Time 
lbd1 <- 0.6  # Intensity for the first process
lbd2 <- 0.3  # Intensity for the second process
nev1 <- rpois(1,tt*lbd1) #Number of Events process 1
nev2 <- rpois(1,tt*lbd2) #Number of Events process 2
tt1 <-sort(runif(nev1,max=tt)) # Location of events P1
tt2 <-sort(runif(nev2,max=tt)) # Location of events  P2
tTimes <- sort(c(tt1,tt2))
plot(c(0,tTimes,tt),c(0:(nev1+nev2),(nev1+nev2)), 
     type='s',xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),
     main='Superposition of Poisson processes',
     cex.main=1.7,cex.lab=1.5)
abline(h=0,col='red')
lines(c(0,tt1,tt),c(0:nev1,nev1), type='s',col='blue',lwd=2)
lines(c(0,tt2,tt),c(0:nev2,nev2), type='s',col='red',lwd=2)


### Non-homogeneous Poisson processes
### Start with a simple non-homogeneous Poisson
### process with intensity 
### lambda(t) = 0.5 + 0.5*t/30
### Simulate 30 time units of this process
### by thinning a homogeneous process with
### parameter 1

set.seed(3838)
tt <- 30  # Time of observation
lbd <- 1  # Intensity
nev <- rpois(1,tt*lbd) # Number of events 
tTimes <-c(sort(runif(nev,max=tt)),tt) # Position of the events
plot(c(0,tTimes),c(0:nev,nev), type='s',cex.lab=1.5,cex.main=1.7,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),ylim=c(0,30),
     main='Non-homogeneous Poisson process',col=blues[9])
abline(h=0,col='red')
mark1 <- rbinom(nev,1,0.5+(0.5*tTimes/30))
lines(c(0,tTimes[mark1==1],30),
      c(0:length(mark1[mark1==1]),length(mark1[mark1==1])),
      type='s',lwd=2,col=blues[7])
mark2 <- rbinom(nev,1,(tTimes/30))
lines(c(0,tTimes[mark2==1],30),
      c(0:length(mark2[mark2==1]),length(mark2[mark2==1])),
      type='s',lwd=2,col=blues[5])
segments(3,27,8,27,col = blues[9])
segments(3,26,8,27,col = blues[7])
segments(3,25,8,27,col = blues[5])
segments(3,25,8,25,col='red')
segments(3,25,3,28,col='red')
text(3,24.2,0);text(8,24.2,30);text(2.7,27,1)
text(5.5,28.5,'Intensities')


###  Second example
###  Mixture of homogeneous and non-homogeneous
###  Poisson process

### Intensity function
xx <- (0:800)/100
plot(xx,rep(20,801), type='n', main='Intensity function',
     xlab='Time (h)',ylab = 'customers/hour',ylim=c(0,25),
     cex.main=1.7, cex.lab=1.5)
lines((0:200)/100, 20*sqrt((0:200)/200),type='l',
      lwd=2,col=blues[7])
lines((200:600)/100, rep(20,401),type='l',lwd=2,col=blues[7])
lines((600:800)/100, 20*sqrt((800-(600:800))/200),
      type='l',lwd=2,col=blues[7])
abline(h=0,col=reds[7])


### Homogeneous process
set.seed(3838)
tt <- 8  # Time of observation
lbd <- 20  # Intensity
nev <- rpois(1,tt*lbd) # Number of events 
tTimes <-sort(runif(nev,max=tt)) # Position of the events
###  Graph of the homogeneous process
plot(c(0,tTimes,tt),c(0:nev,nev), type='s',cex.lab=1.5,cex.main=1.7,
     xlab='Time',ylab='Events',lwd=2,xlim=c(0,tt),ylim=c(0,150),
     main='Non-homogeneous Poisson process',col=blues[9])
abline(h=0,col='red')
### Non-homogeneous process
tt1 <- tTimes[tTimes<=2] # events during the first 2 hrs  
ll1 <- length(tt1) # number of events in the first 2 hrs
tt2 <- tTimes[tTimes>=6] # events during the last 2 hrs
ll2 <- length(tt2) # number of events in the last 2 hrs
mark1 <- rbinom(ll1,1,sqrt(tt1/2)) # thinning at the beginning
lmk1 <- sum(mark1) # Number of events in the first 2 hrs
lines(c(0,tt1[mark1==1]), # graph of first 2 hrs
      c(0:lmk1),type='s',lwd=2,col=blues[7])
lines(c(tTimes[ll1:(nev-ll2)],6), # graphs of homogeneous interval
      c(lmk1:(nev-ll2-ll1+lmk1),(nev-ll2-ll1+lmk1)), 
      type='s',lwd=2,col=blues[7])
mark2 <- rbinom(ll2,1,sqrt((8-tt2)/2)) # thinning at the end
lmk2 <- sum(mark2) # number of events in the last 2 hrs
lines(c(6,tt2[mark2==1],8), # graph of last 2 hrs
      c((nev-ll2-ll1+lmk1),
        (nev-ll2-ll1+lmk1):(nev-ll2-ll1+lmk1+lmk2-1),
        (nev-ll2-ll1+lmk1+lmk2-1)),
      type='s',lwd=2,col=blues[7])
abline(v=c(2,6), col=reds[7],lwd=2)