conditional proability as gamma

lbd <- 1.05^((1:50))

n.iter = 10000
gm <- 0
for (i in 1:length(lbd)){
  gm<-gm+rexp(n.iter,lbd[i])
}

hist(gm,prob=TRUE)

l <- sum(1/lbd)/sum(1/lbd^2) # beta
n <- sum(1/lbd)^2/sum(1/lbd^2) # alpha
x = 1:max(gm)
lines(x,dgamma(x,n,l),col="black")

## poisson

spois <- function(rate,increase=0,stop=15){
  n = 0 
  sum = 0
  while(sum < stop){
    inc = rexp(1,rate=rate)
    sum = sum + inc
    rate = rate+increase
    n = n + 1
  }
  return(n)
}

n.iter = 100000
sim = vector(mode='numeric',length = n.iter)
for(i in 1:n.iter){
  sim[i] = spois(rate=0.8)
}


hist(sim,prob=T)
x = 1:max(sim)
lines(x,dpois(x,lambda = 0.8*15))