LAB 3. LIKELIHOOOD

setwd("C:/Inference Course Spring 2006/Labs/Likelihood course")

mu.true<-1
k.true<-0.4
x<-rnbinom(50,mu=mu.true,size=k.true)

#Take a look at the data using some graphics. 

plot (table(factor(x,levels=0:max(x))),ylab="Frequency",xlab="x")

#Now we will build a function that calculates the negative 
# log-likelihood for this distribution given a set a parameters. 
# Arguments for this function are p, the vector of parameters, 
# and dat the vector of data: 

NLLfun1 <- function(p, dat = x) {
mu = p[1]
k = p[2]
 -sum(dnbinom(x, mu = mu, size = k, log = TRUE))
 }

# First we will calculate the negative log-likelihood with the true values of the distribution. We have to combine these values into a vector to be able to pass them to the NLL function:

nll.true <- NLLfun1(c(mu = mu.true, k = k.true))


m = mean(x)
 v = var(x)
 mu.mom = m
 k.mom = m/(v/m - 1)

# And then obtain the MLE using optim

sol1<-optim(fn=NLLfun1,par=c(mu=mu.mom,k=k.mom),hessian=TRUE)

muvec<-seq(0.4,3,by=0.05)
kvec<-seq(0.01,0.7,by=0.01)
resmat<-matrix(nrow=length(muvec),ncol=length(kvec))

##Now we calculate the nll for all these values and store them in the matrix:

for (i in 1:length(muvec)) {
 for (j in 1:length(kvec)) {
 resmat[i, j] = NLLfun1(c(muvec[i], kvec[j]))
 }
 }

# Now plot this.  
contour(muvec,kvec, resmat, xlab="mu",ylab="k")

#And to add more resolution in the middle:

contour(muvec,kvec,resmat,levels=70:80,lty=2,add=TRUE)


# Confidence intervals

NLLfun.mu<-function(p,mu){
    	k<-p[1]
	-sum(dnbinom(x,mu=mu,size=k,log=TRUE))
	}
mu.profile<-matrix(ncol=2,nrow=length(muvec))

# Run optimization for mu
for (i in 1:length(muvec)) {
	Oval = optim(fn = NLLfun.mu, par = sol1$par["k"], method = "L-BFGS-B",
	lower = 0.002, mu = muvec[i])
 	mu.profile[i, ] = c(Oval$par, Oval$value)
 }

colnames(mu.profile) = c("k", "NLL")

# Do the same process for k:
NLLfun.k = function(p, k) {
mu = p[1]
-sum(dnbinom(x, mu = mu, size = k, log = TRUE))
}
k.profile = matrix(ncol = 2, nrow = length(kvec))
 for (i in 1:length(kvec)) {
	Oval = optim(fn = NLLfun.k, par = sol1$par["mu"], method = "L-BFGS-B",
	lower = 0.002, k = kvec[i])
	k.profile[i, ] = c(Oval$par, Oval$value)
 	}

colnames(k.profile) = c("mu", "NLL")

#Redraw the contour plot with profiles added:

contour(muvec, kvec, resmat, xlab = "mu", ylab = "k")
contour(muvec, kvec, resmat, levels = 70:80, lty = 2, add = TRUE)
lines(muvec, mu.profile[, "k"], lwd = 2)
lines(k.profile[, "mu"], kvec, lwd = 2, lty = 2)

#univariate profiles
plot(muvec, mu.profile[, "NLL"], type = "l", xlab = "mu",
ylab = "Negative log-likelihood")

# Support limits
cutoff<-sol1$value + qchisq(0.95,1)/2

relheight = function(mu) {
 	O2 <- optim(fn = NLLfun.mu, par = sol1$par["k"], method = "L-BFGS-B",
 	lower <- 0.002, mu = mu)
  	O2$value -cutoff
 	}

#We know the lower limit is somewhere around 0.7, so going on either side
#should give us values that are negative/positive.

relheight(mu = 0.6)

relheight(mu = 0.8)

#Using R’s uniroot function, which takes a single-parameter function and
#searches for a value that gives zero:

lower = uniroot(relheight, interval = c(0.4, 1))$root
upper = uniroot(relheight, interval = c(1.2, 5))$root
ci.uniroot = c(lower, upper)
plot(muvec, mu.profile[, "NLL"], type = "l", xlab = "mu", ylab = "Negative log-likelihood")
abline(v = ci.uniroot, lty =3)
cutoffs = c(0, qchisq(c(0.95), 1)/2)
nll.levels = sol1$value + cutoffs
abline(h = nll.levels, lty = 1:2) 


##REEF FISH EXERCISE##########################################################

data<-read.table("Reef_fish.txt", header=T, sep="\t")
attach(data)

plot(settlers, recruits)
# so that you can guess starting values for a and b

hist(settlers)

#For the exercise
adult.recruits<-function(a,b,settlers){a*settlers/(1+(a/b)*settlers)}

plot(settlers, recruits)
# so that you can guess starting values for a and b


#likelihood function for a zero inflated binomial:

NLLfun<-function(a,b)
{
recprob=a/(1+(a/b)*settlers)
-sum(dbinom (recruits, prob=recprob, size=settlers,log=TRUE), na.rm=TRUE)
}

m4<-mle(minuslogl=NLLfun,start=list(a=0.5,b=10),method="L-BFGS-B",lower=0.003)

fish.data<-cbind(recruits,settlers)
recruits=recruits[settlers>0]
settlers=settlers[settlers>0]

m4<-mle(minuslogl=NLLfun,start=list(a=0.5,b=10),method="L-BFGS-B",lower=0.003)

# You can now plot the fits that you got from the analyses:

a<-coef(m4)["a"]
b<-coef(m4)["b"]
plot(settlers,recruits)
curve(a*x/(1+(a/b)*x), add=TRUE)

# and calculate the likelihod
NL<-NLLfun(a,b)

## Is there another option?

#############################################################################
## Sapling growth exercise

# Using R functions and the stats4 package
growth<- read.table("Sapling_Growth.txt", header=TRUE)

NLLfun<-function(a,b,scale)
{
pred.growth<-a*light/(1+(a/b)*light)
-sum(dgamma(Growth,shape=pred.growth/scale,scale=scale,log=TRUE),na.rm=TRUE)
}

library(bbmle)
m3<-mle(minuslogl=NLLfun,start=list(a=0.5,b=10,scale=2),method="L-BFGS-B",lower=0.01)

a<-coef(m3)["a"]
b<-coef(m3)["b"]
scale1<-coef(m3)["scale"] # I added the 1 because scale is an R function and it messes up the syntax
plot(light,Growth)
x<-light
curve(a*x/(1+(a/b)*x), add=TRUE)


# Check the fit
 Pred.Growth<-a*x/(1+(a/b)*x)
 plot(Growth,Pred.Growth,xlim=c(0,15),ylim=c(0,15))
 abline(a=0,b=1)


############################################################################################3


## We will use the likelihood package. 

library(likelihood)

# Write the list of parameters you will need. Remember that these must include those of your pdf.

pred.growth<-function(a,b){a*growth$light/(1+(a/b)*growth$light)}

par<-list()
par$a<-0.5
par$b<-5

#How do we provide the parameters of the pdf? 


par$log<-TRUE
par_lo<-list(a=0,b=0,scale=0)
par_hi<-list(a=5,b=10,scale=10)
par_step<-list(a=5,b=10,scale=10)


var<-list()
var$x<-"Growth"
var$mean<-"Predicted"
var$scale<-1
wrap<-function(mean, scale){mean/scale}
var$shape<-wrap



sol2<-anneal(model=pred.growth, par=par, dep_var="Growth", var=var,pdf=dgamma, source_data=growth)