# Graphing the likelihood and the log Likelihood of Gamma(shape=alpha,scale=beta)
#pdf("/home/david/Documents/Spring2009/math_442/march/gammaMLE1.pdf", width=6, height=6)

xx <- c(2.799, 5.857, 4.018, 7.461, 3.518, 3.560, 2.304, 4.114, 3.316, 7.430)

n <- length(xx)

# Define the likelihood function.
# Compute the required constants, sumXvec, sumlnXvec, n, before
# as this function will be executed many times.

# note gamma function is not defined at zero
logL = function(a,b,x)
{
n<-length(x)
-n*log(gamma(a))-n*a*log(b)+(a-1)*sum(log(x))-sum(x)/b
}
L = function(a,b,x)
{
n<-length(x)
exp(-n*log(gamma(a))-n*a*log(b)+(a-1)*sum(log(xx))-sum(x)/b)
}

par(mfrow=c(1,2))
# generate the data
A <- seq(0.01,20,length=100)
B <- seq(0.01,4,length=100)
Z1 <- outer(A,B,logL,x=xx)
Z2<- outer(A,B,L,x=xx)

persp(A,B,Z1, theta=45, phi=45, xlab="alpha", ylab="beta",ticktype="detailed", expand=0.75,
	zlab="Log Likelihood", main="Log Likelood Function",col = "lightblue") 

persp(A,B,Z2, theta=45, phi=50, xlab="alpha", ylab="beta",ticktype="detailed", expand=0.75,
	zlab="Likelihood", main="Likelihood Function",col = "lightblue")