# Code to model in 3d the MLE and find estimators of alpha and lambda for 8.46 
# David Read, California State University, East Bay, 10APR09

#NOTE:The rgl package must be installed on this R version for this code to work!!!

library(rgl)

partitions=800

whale = c(0.2252, 0.7407, 0.2825, 1.1494, 0.4367, 0.346, 1.7241, 0.5882, 0.3448, 0.3984, 1.087, 2.0408, 0.2551, 0.3623, 0.2273, 0.4219, 0.2415, 1.1905, 0.3817, 0.1667, 0.2283, 0.7194, 0.3571, 0.4587, 1.0101, 0.5556, 0.5076, 0.2273, 0.5208, 0.4566, 1.7857, 1.0204, 0.4505, 0.3003, 0.1776, 0.3367, 4, 1.2987, 0.2053, 0.3497, 0.7576, 0.3802, 0.1642, 0.3968, 2.8571, 0.5405, 0.7463, 0.3584, 1.13158, 0.7519, 1.4286, 1.1494, 0.3436, 3.7037, 0.2252, 0.2203, 0.2188, 0.289, 1.2658, 0.2146, 0.3021, 0.2525, 0.4587, 0.2398, 0.6803, 1.3699, 0.4237, 0.4032, 0.6369, 0.4525, 0.3448, 0.3145, 0.2525, 0.4717, 0.3268, 0.8696, 0.2222, 0.1912, 0.1984, 0.2421, 0.3636, 0.3268, 0.1675, 0.4785, 0.188, 0.4049, 0.3731, 0.3367, 1.18182, 0.3559, 2.2222, 1.7857, 0.2079, 0.3497, 0.3003, 0.4292, 0.1718, 0.3106, 0.5405, 0.3311, 0.2353, 0.232, 0.2151, 0.2538, 0.346, 0.2123, 0.2169, 1.4925, 2.2727, 0.3165, 0.3021, 0.5181, 1.5152, 0.2342, 0.641, 0.3289, 0.3623, 0.2959, 0.1475, 0.3636, 0.2571, 1.1494, 0.2959, 0.369, 0.1812, 0.4065, 0.1337, 1.7241, 0.5988, 0.5128, 0.4049, 0.3311, 0.2141, 0.4386, 0.2016, 0.4132, 0.2915, 0.2016, 0.2347, 0.2899, 0.2494, 0.3106, 0.1724, 0.463, 2.5, 0.1828, 0.5618, 0.3401, 2.6316, 0.4167, 0.3584, 0.2786, 0.3086, 0.9346, 0.3484, 0.1451, 0.2179, 0.3731, 1.1364, 0.3268, 0.2096, 0.3448, 0.2179, 0.2778, 0.3521, 1.5625, 0.3058, 0.4274, 0.4762, 0.271, 3.2258, 0.4464, 0.667, 1.0638, 0.9434, 0.2865, 0.3086, 0.2538, 0.7407, 0.2119, 0.3497, 0.1466, 0.2128, 0.3344, 0.1686, 0.4065, 4.3478, 0.3185, 0.3279, 0.2028, 0.3571, 0.7246, 0.5882, 0.6369, 0.3155, 0.3802, 2.5641, 0.3636, 0.3012, 0.8696, 1.4925, 0.2364, 0.2404, 0.2907, 0.5917, 0.202, 0.2841, 0.1337, 0.6135, 0.237)

v1<-rep(seq(1.0001, 3, by=(2/partitions)), times=partitions)   #v1 is ALPHA

v2<-seq(1.0001, 3, by=(2/partitions))                          #v3 is LAMBDA
v3<-rep(v2, each=partitions)

#NOTE:The bloody MLE function is  -(n*alp*log(lam)+(alp-1)*sum(log(whale))-lam*sum(whale)-n*log(gamma(alp)))

# These sums taken once and turned into constants save hundreds of millions of calculations!

n=length(whale)
totalwhale=sum(whale)
logtotalwhale=sum(log(whale))

plot3d(v1[1:length(v1)], v3[1:length(v3)], -(n*v1[1:length(v1)]*log(v3[1:length(v3)])+(v1[1:length(v1)]-1)*logtotalwhale-v3[1:length(v3)]*totalwhale-n*log(gamma(v1[1:length(v1)]))), xlab="ALPHA", ylab="LAMBDA", zlab="f(x,y)", main="Plot of MLE")

#Code below determines estimates for alpha and lambda, more partitions means a better estimate, but slower processing!

v4<-numeric(length(v1))

for (i in 1:length(v1))
	{
	Next = -(n*v1[i]*log(v3[i])+(v1[i]-1)*logtotalwhale-v3[i]*totalwhale-n*log(gamma(v1[i])))
	v4[i] = Next
	}

element=which(v4==min(v4))
v1[element]     #ALPHA ESTIMATE
v3[element]     #LAMBDA ESTIMATE