ab0.R

# (a, b, 0) class
gen.dist <- function(a=0, b=10, MAX=100) {
  # generate the distribution for x = 0, 1, 2, ..., MAX
  pdf <- numeric(MAX + 1)
  ratios = as.numeric(lapply(1:MAX, function(k) {force(k); a + b/k  }))
  
  total = 1
  for(k in 1:MAX) total = total + prod(ratios[1:k])
  pdf[1] = 1/total # p0
  
  for(k in 1:MAX) pdf[k+1] = pdf[k]*(a + b/k)
  pdf
}

plot.dist <- function(dist, MAX=100, title="") {
  plot(1:(MAX+1) - 1, dist, type='h', ylab="p(x)", xlab="x", main=title)
}

plot.dist(gen.dist(a=0, b=10), title="Poisson(lambda = 10)")

plot.dist(gen.dist(a=0.9, b=0), title="Geometric(B/(1+B) = 0.9)")

plot.dist(gen.dist(a=0.5, b=10), title="NB(B/(1+B) = 0.5, r = 19)")

plot.dist(gen.dist(a=-0.25, b=100), title="Bin(q = 1/3, m = 399)")