# Vector of sample y-values:

y <- c(9.18,18.33,12.39,5.74,13.51,8.66,18.63,7.82,9.16,10.48,11.58,7.54,8.45,5.43,7.09,2.15,2.01,18.70,5.63,7.71,21.31,4.68,15.82,27.53,16.46,20.40,4.76,5.99,19.51,8.63)


#Getting method of moments estimates as a quick answer:

sum.ysq <- sum(y^2)

# method of moments estimates (see MME example in notes)

alpha.mm <- ( 30*(mean(y))^2 ) / ( sum.ysq - 30*(mean(y))^2 )
beta.mm <- ( sum(y^2) - 30*(mean(y))^2 ) / ( 30*mean(y) )

# Sufficient statistics:
sum.ln.y <- sum(log(y))
sum.y <- sum(y)

# log-likelihood function (to be maximized)
# param[1] = alpha
# param[2] = beta
loglike<-function(param){
alpha <- param[1]
beta <- param[2]
-30*log(gamma(alpha)) - 30*alpha*log(beta) + sum.ln.y*(alpha-1) - sum.y/beta
}

# Use "optim" function to maximize the log-likelihood function
# Use MMEs as starting values:
# The "control=list(fnscale = -1)" part makes it a maximization problem

mle<-optim( par=c(alpha.mm,beta.mm),fn=loglike, control=list(fnscale = -1) )

alpha.mle <- mle$par[1]
beta.mle <- mle$par[2]


# Printing results:

print(paste("MME of alpha:", round(alpha.mm,3), ", MME of beta:", round(beta.mm,3)))

print(paste("MLE of alpha:", round(alpha.mle,3), ", MLE of beta:", round(beta.mle,3)))