#  Code for simulating Gamma data and finding the Maximum Likelihood
#  Estimates by Newton-Raphson.

#  The code uses the "digamma" and "trigamma" functions.  These are
#  the first and second derivative of the logarithm of the gamma function.
#  It also uses "lgamma", which is log(gamma(x)).

gamma.llik <- function(y,alpha,beta){
  n <- length(y)
  return(n*alpha*log(beta)+(alpha-1)*sum(log(y))-beta*sum(y)-n*lgamma(alpha))
}

gamma.nr.step <- function(y,alpha.j,beta.j){
  n <- length(y)
  s1 <- n*log(beta.j)+sum(log(y))-n*digamma(alpha.j)
  s2 <- n*alpha.j/beta.j-sum(y)
  s <- c(s1,s2)
  h11 <- -1*n*trigamma(alpha.j)
  h12 <- n/beta.j
  h22 <- -1*n*alpha.j/beta.j^2
  h <- matrix(c(h11,h12,h12,h22),ncol=2)
  return(c(solve(h,s)))
}


y <- rgamma(40,10,2)

alpha.0 <- mean(y)^2/var(y)
beta.0 <- mean(y)/var(y)

print(alpha.0)
print(beta.0)

gamma.llik(y,alpha.0,beta.0)

step.1 <- gamma.nr.step(y,alpha.0,beta.0)

alpha.1 <- alpha.0-step.1[1]
beta.1 <- beta.0-step.1[2]

print(alpha.1)
print(beta.1)

gamma.llik(y,alpha.1,beta.1)

step.2 <- gamma.nr.step(y,alpha.1,beta.1)

alpha.2 <- alpha.1-step.2[1]
beta.2 <- beta.1-step.2[2]

print(alpha.2)
print(beta.2)

gamma.llik(y,alpha.2,beta.2)

step.3 <- gamma.nr.step(y,alpha.2,beta.2)

alpha.3 <- alpha.2-step.3[1]
beta.3 <- beta.2-step.3[2]

print(alpha.3)
print(beta.3)

gamma.llik(y,alpha.3,beta.3)