Homogener Poisson-Prozess

In [1]:
options(repr.plot.width=5, repr.plot.height=5)
# P(N(t)=n|lambda)
Poissonpf <- function(n, t, lambda) {
    exp(-lambda * t) * (lambda*t)^n/factorial(n)
}
In [2]:
# Zeitintervall und Intensität festlegen
t = 90
lambda = 1/30
In [3]:
(p <- Poissonpf(1:10, t, lambda))
plot(1:10, p, type='o', xlab='n', ylab='P(N) = n')
  1. 0.149361205103592
  2. 0.224041807655388
  3. 0.224041807655388
  4. 0.168031355741541
  5. 0.100818813444924
  6. 0.0504094067224622
  7. 0.0216040314524838
  8. 0.00810151179468143
  9. 0.00270050393156048
  10. 0.000810151179468143
In [4]:
# E(N(t)) und SD(N(t))
c(lambda * t, sqrt(lambda*t))
  1. 3
  2. 1.73205080756888
In [5]:
# Simulation of imax trials and plot example events
# Watch the clustering of events
Poisson <- function(t, lambda) {
    N <- 0
    S <- 0
    T <- c()
    weiter = TRUE
    while (weiter) {
        S <- S-log(runif(1))/lambda
        if (S > t) weiter <- FALSE
        else {
            N <- N+1
            T <- c(T, S)
        }
    }
    return(c(N, T))
}

imax <- 20
plot(1:10, type='n', xlim=c(0, t), ylim=c(0.2, imax+0.8),
    xlab='Duration [sec]', ylab='Trial')
for (i in 1:imax) {
    x <- Poisson(t, lambda)
    N <- x[1]
    T <- x[-1]
    points(T, rep(i, N))
    lines(c(0, t), c(i, i))
}
In [6]:
# Run Nsim simulations and plot relative frequencies
Nsim <- 10000
N <- c()
for (i in 1:Nsim) {
    N <- c(N, Poisson(t, lambda)[1])
}
In [14]:
(tbl <- table(N))
plot(1:10, p, pch=4, col='blue', xlab='n', ylab='P(N)=n')
points(1:10, tbl[as.character(1:10)]/Nsim, type='h', col='red')
points(1:10, tbl[as.character(1:10)]/Nsim, col='red')
N
   0    1    2    3    4    5    6    7    8    9   10   11 
 488 1493 2268 2284 1706  980  455  218   77   22    7    2 
In [ ]: