# finds the bias of the Method of Moments Estimator of lambda in the
# gamma family

nmlambbias <- function(r,lamb,n,nreps) {
   lambhat <- vector(length=nreps)
   unfudge <- (n-1) / n
   for (i in 1:nreps) {
      x <- rgamma(n,shape=r,rate=lamb)
      xbar <- mean(x)
      s2 <- var(x) * unfudge
      lambhat[i] <- xbar / s2
   }
   mean(lambhat) - lamb
}