NBD with t distributed beta

beta_nbd.R

# NBD with t distributed beta

beta.nbd.density <- function(x, r, alpha, a, b)
{
  # params:
  #   x: count we are evaluating density at
  #   r: nbd r parameter
  #   alpha: nbd alpha parameter
  #   a: beta alpha parameter
  #   b: beta beta parameter
  # output:
  #   probability density at x
  
  k <- gamma(x + r) * alpha^r / (gamma(r) * factorial(x) * beta(a, b))
  f <- function(t) (alpha + t)^-(r + x) * t^(alpha + x - 1) * (1 - t)^(b - 1)
  i <- integrate(f, lower = 0, upper = 1)
  return(k*i$value)
}