# Poisson counts and exponential waiting times
# in case of tire ruptures

W=30000; # width of time window (in Km rather than time)
mu=5000  # mu of time r.v

# start with generation 0 ('g.0') tire (g.0)

pr.end.in.g.0 = exp(-W/mu) # prob of not having exited g.0 by time W
pr.end.in.g.0

dpois(0:12,W/mu) # compare with above

1-pr.end.in.g.0 # prob of having exited g.0 by time W

density.g.1.entries = function(t) (1/mu) * exp(-t/mu)

pr.end.in.g.1.if.enter.g.1 = function(t) exp( -(W - t)/mu )

product = function(t) density.g.1.entries(t) * pr.end.in.g.1.if.enter.g.1(t)

pr.end.in.g.1.if.ENTER.g.1 = integrate(product, 0, W)$value/
                             integrate(density.g.1.entries, 0, W)$value
                             
pr.end.in.g.1.if.ENTER.g.1  # over all entry times from 0 to W

(1 - pr.end.in.g.0 )*pr.end.in.g.1.if.ENTER.g.1 # prob(enter g.1) x prob(stay in g.1 of enter it)

dpois(0:12,W/mu) # compare PoissonProb(1 | mu=6) with above

## Can continue on iteratively, to gen.2 gen.3 etc.....

#  graphically (different coloour for each gen.) 

y = 0:12 ; t = seq(0,W,1000)

m=matrix(NA,nrow=length(y),ncol=length(t))
m[1,]  = rep(0,length(t)) 

for(Y in y[-1]) m[Y+1,]  = ppois(Y-1, t/mu) 

par(mfrow=c(1,1),mar = c(5,5,0.1,0.1))
plot( c(0,1.1*max(t)), c(0,1.15), col="white",xlab="Km",ylab="Proportion in g 0, 1, 2, ... ")
for(Y in y[-1]){
	polygon(c(t,-sort(-t)), c(m[Y,], m[Y+1,length(t):1] ) , col=Y,border=NA)
	text(max(t), (m[Y,length(t)]+m[Y+1,length(t)])/2,toString(Y-1), col=Y, adj=c(-0.5,0.5), cex=0.6)
} 
for(T in seq(2500,30000,2500)){
 text(T, 1.0,paste("E[events]=",toString(T/mu)), adj=c(0,0), cex=0.75,srt=45)
 segments(T, 0,T,1, cex=0.5,col="white",lwd=0.5)
} 
text(1.04*max(t), 0.5, "generation", srt=90)
# event = transition from 1 state or generation to another 
# (failure / tire rupture/ birth / death / marriage / whatever ... )