#------------------------------------------------
# simula un processo di Poisson
# analizzandone poi i diversi "punti di vista"
#
# GdA, PED 27/11/2010
#-------------------------------------------------

# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

r   = 1     # intensita' del processo in s^-1
p.i = 0.01  # (piccola) probabilita' che l'evento accada in dt
dt  = p.i/r # tempo di 'quantizzazione' del fenomeno ai fini della simulazione
dtm = dt/2  # la sua meta'

t.tot = dt #tempo totale da prendere in considerazione


n.evts = 10000 # nr di eventi da generare (non troppi perche' 
               # la simulazione e' inefficiente)

t.i = rep(0, n.evts) # vettore dei tempi in cui accadono gli eventi
t = 0       # istante iniziale
dtm = dt/2  # meta' di dt
n = 0       # n-esimo intervallino
n.gen = 0   # nr di conteggi effettivamente generati 
for (i in 1:n.evts) {
  ok = FALSE
  while (!ok) {
    n = n + 1
    if(runif(1) < p.i) {
      ok = TRUE
      n.gen = n.gen + 1
      t.i[n.gen] = n*dt - dtm 
    }
  }
}

# plot dei tempi dei primi 100 conteggi
plot(t.i[1:100], pch='.', cex=2)

pausa()

# plot dei tempi dei primi 1000 conteggi
plot(t.i[1:1000], pch='.', cex=2)

pausa()

# analizziamo le differenze di tempo fra conteggi successivi
delta.t.i = rep(0, length(t.i) - 1)
for (i in 2:length(t.i) ){
  delta.t.i[i] = t.i[i] - t.i[i-1]
}

hist(delta.t.i, nc=100, xlab="differenze di tempo fra conteggi successivi")

pausa()

# analizziamo in nr. di conteggi in un intervallo di tempo prefissato
Dt = 1  # secondi

nDt = floor (t.i[length(t.i)] / Dt) # nr di intervalli di tempo

n.counts.in.Dt = rep(0, nDt) 
for (i in 1:nDt) {
  t2 = i*Dt
  t1 = t2 - Dt
  n.counts.in.Dt [i] = length(t.i[t.i > t1 & t.i <= t2])
}

conta <- table(n.counts.in.Dt)
barplot(conta)

# provare a cambiare Dt e ripetere soltanto quest'ultima parte