Saturday, January 5, 2008

moment generating functions in R

##Binomial
mbinom <- function(t, n, p) (p * exp(t) + (1 - p))^n

##Discrete uniform
mdunif <- function(t, n) apply(exp(t %o% 1:n), 1, sum)

##Geometric
mgeom <- function(t, p) { et <- exp(t); ifelse(t < -log(1 - p), p * et / (1 - (1 - p) * et), NA) }

##Negative binomial
mnbinom <- function(t, p, r) ifelse(t < -log(1 - p), (p / (1 - (1 - p) * exp(t))^r, NA)

##Poisson
mpois <- function(t, r) exp(r * (exp(t) - 1))

##Gamma
mgamma <- function(t, a, b) (1/(1 - b * t))^a

##Normal
mnorm <- function(t, mean, var) exp(mean * t + var * t * t / 2)

##Uniform
munif <- function(t, min, max) (exp(max * t) - exp(min * t)) / ((max - min) * t)

The formulas are copied from Casella and Berger. I haven't checked the code carefully for errors. I'm puzzled why there is an extra exp(t) term in the denominator for the geometric distribution since Geometric(p) = NegBinom(p, r=1). I'll check why later. Now, thank goodness, I'm ready to fall asleep.

No comments: