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)
}
# Zeitintervall und Intensität festlegen
t = 90
lambda = 1/30
(p <- Poissonpf(1:10, t, lambda))
plot(1:10, p, type='o', xlab='n', ylab='P(N) = n')
# E(N(t)) und SD(N(t))
c(lambda * t, sqrt(lambda*t))
# 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))
}
# Run Nsim simulations and plot relative frequencies
Nsim <- 10000
N <- c()
for (i in 1:Nsim) {
N <- c(N, Poisson(t, lambda)[1])
}
(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')