Chapter 16 Likelihood

Let’s study the log-likelihood profiles as a function of \(p\) for samples sizes of 10,100 and 1000.

#set a range of recombination probabilities (ranging from 0 to 0.5) for which we want to calculate the likelihood
r <- seq(0, 0.5, by=0.01)

#define a matrix in which we will keep track of likelihood values as a function of probability
l <- matrix(data=NA, nrow=length(r), ncol=2, byrow=TRUE, dimnames=list(c(),c("p","log-likelihood")))

#make a loop to step through the likelihood calculation for each value of the parameter p
for(i in 1:length(r)) { 

print("...", quote=FALSE)
print(paste("step ",i," of loop. Things look good.", sep=""))

#define probability for each step in the loop
p <- r[i]

#print some output each time we go thorugh the loop to tell us the current computation
print(paste("computing likelihood for p =", p, sep=" "))

#write down the likelihood function and assign it to a variable
#to caluculate the log.likelihood for 10x the data change: 10=>100, 2=>20, 8=>80
#to caluculate the log.likelihood for 100x the data change: 10=>1000, 2=>200, 8=>800
log.likelihood <- log(choose(10,2)) + 2*log(p) + 8*log(1-p)

#print the likelihood to the terminal so we can see it in real time
print(paste("log-likelihood", log.likelihood, sep=" "))

#record probability and log-likelihood value
l[i,] <- c(p,log.likelihood)

if(i==51) {print(paste("...","exiting loop","...", sep=" "))}

}

## [1] ...
## [1] "step 1 of loop. Things look good."
## [1] "computing likelihood for p = 0"
## [1] "log-likelihood -Inf"
## [1] ...
## [1] "step 2 of loop. Things look good."
## [1] "computing likelihood for p = 0.01"
## [1] "log-likelihood -5.48408056903387"
## [1] ...
## [1] "step 3 of loop. Things look good."
## [1] "computing likelihood for p = 0.02"
## [1] "log-likelihood -4.17900517962613"
## [1] ...
## [1] "step 4 of loop. Things look good."
## [1] "computing likelihood for p = 0.03"
## [1] "log-likelihood -3.45012696474731"
## [1] ...
## [1] "step 5 of loop. Things look good."
## [1] "computing likelihood for p = 0.04"
## [1] "log-likelihood -2.95766511612812"
## [1] ...
## [1] "step 6 of loop. Things look good."
## [1] "computing likelihood for p = 0.05"
## [1] "log-likelihood -2.59514841243807"
## [1] ...
## [1] "step 7 of loop. Things look good."
## [1] "computing likelihood for p = 0.06"
## [1] "log-likelihood -2.31516217349445"
## [1] ...
## [1] "step 8 of loop. Things look good."
## [1] "computing likelihood for p = 0.07"
## [1] "log-likelihood -2.09242312677392"
## [1] ...
## [1] "step 9 of loop. Things look good."
## [1] "computing likelihood for p = 0.08"
## [1] "log-likelihood -1.9118476703586"
## [1] ...
## [1] "step 10 of loop. Things look good."
## [1] "computing likelihood for p = 0.09"
## [1] "log-likelihood -1.76371416330336"
## [1] ...
## [1] "step 11 of loop. Things look good."
## [1] "computing likelihood for p = 0.1"
## [1] "log-likelihood -1.64139182148038"
## [1] ...
## [1] "step 12 of loop. Things look good."
## [1] "computing likelihood for p = 0.11"
## [1] "log-likelihood -1.54015786665673"
## [1] ...
## [1] "step 13 of loop. Things look good."
## [1] "computing likelihood for p = 0.12"
## [1] "log-likelihood -1.45653155470894"
## [1] ...
## [1] "step 14 of loop. Things look good."
## [1] "computing likelihood for p = 0.13"
## [1] "log-likelihood -1.38787570595085"
## [1] ...
## [1] "step 15 of loop. Things look good."
## [1] "computing likelihood for p = 0.14"
## [1] "log-likelihood -1.33214634085202"
## [1] ...
## [1] "step 16 of loop. Things look good."
## [1] "computing likelihood for p = 0.15"
## [1] "log-likelihood -1.28772891598364"
## [1] ...
## [1] "step 17 of loop. Things look good."
## [1] "computing likelihood for p = 0.16"
## [1] "log-likelihood -1.25332753488452"
## [1] ...
## [1] "step 18 of loop. Things look good."
## [1] "computing likelihood for p = 0.17"
## [1] "log-likelihood -1.22788781962538"
## [1] ...
## [1] "step 19 of loop. Things look good."
## [1] "computing likelihood for p = 0.18"
## [1] "log-likelihood -1.21054187620424"
## [1] ...
## [1] "step 20 of loop. Things look good."
## [1] "computing likelihood for p = 0.19"
## [1] "log-likelihood -1.2005681743982"
## [1] ...
## [1] "step 21 of loop. Things look good."
## [1] "computing likelihood for p = 0.2"
## [1] "log-likelihood -1.19736174561156"
## [1] ...
## [1] "step 22 of loop. Things look good."
## [1] "computing likelihood for p = 0.21"
## [1] "log-likelihood -1.20041167492758"
## [1] ...
## [1] "step 23 of loop. Things look good."
## [1] "computing likelihood for p = 0.22"
## [1] "log-likelihood -1.20928384987723"
## [1] ...
## [1] "step 24 of loop. Things look good."
## [1] "computing likelihood for p = 0.23"
## [1] "log-likelihood -1.22360756342282"
## [1] ...
## [1] "step 25 of loop. Things look good."
## [1] "computing likelihood for p = 0.24"
## [1] "log-likelihood -1.24306498712405"
## [1] ...
## [1] "step 26 of loop. Things look good."
## [1] "computing likelihood for p = 0.25"
## [1] "log-likelihood -1.26738281208371"
## [1] ...
## [1] "step 27 of loop. Things look good."
## [1] "computing likelihood for p = 0.26"
## [1] "log-likelihood -1.29632554843427"
## [1] ...
## [1] "step 28 of loop. Things look good."
## [1] "computing likelihood for p = 0.27"
## [1] "log-likelihood -1.32969010891481"
## [1] ...
## [1] "step 29 of loop. Things look good."
## [1] "computing likelihood for p = 0.28"
## [1] "log-likelihood -1.36730139763174"
## [1] ...
## [1] "step 30 of loop. Things look good."
## [1] "computing likelihood for p = 0.29"
## [1] "log-likelihood -1.40900869380712"
## [1] ...
## [1] "step 31 of loop. Things look good."
## [1] "computing likelihood for p = 0.3"
## [1] "log-likelihood -1.45468267039141"
## [1] ...
## [1] "step 32 of loop. Things look good."
## [1] "computing likelihood for p = 0.31"
## [1] "log-likelihood -1.50421292436223"
## [1] ...
## [1] "step 33 of loop. Things look good."
## [1] "computing likelihood for p = 0.32"
## [1] "log-likelihood -1.55750592310229"
## [1] ...
## [1] "step 34 of loop. Things look good."
## [1] "computing likelihood for p = 0.33"
## [1] "log-likelihood -1.61448329204991"
## [1] ...
## [1] "step 35 of loop. Things look good."
## [1] "computing likelihood for p = 0.34"
## [1] "log-likelihood -1.67508038466687"
## [1] ...
## [1] "step 36 of loop. Things look good."
## [1] "computing likelihood for p = 0.35"
## [1] "log-likelihood -1.73924508796667"
## [1] ...
## [1] "step 37 of loop. Things look good."
## [1] "computing likelihood for p = 0.36"
## [1] "log-likelihood -1.806936826321"
## [1] ...
## [1] "step 38 of loop. Things look good."
## [1] "computing likelihood for p = 0.37"
## [1] "log-likelihood -1.87812573368988"
## [1] ...
## [1] "step 39 of loop. Things look good."
## [1] "computing likelihood for p = 0.38"
## [1] "log-likelihood -1.95279197029709"
## [1] ...
## [1] "step 40 of loop. Things look good."
## [1] "computing likelihood for p = 0.39"
## [1] "log-likelihood -2.03092516446481"
## [1] ...
## [1] "step 41 of loop. Things look good."
## [1] "computing likelihood for p = 0.4"
## [1] "log-likelihood -2.11252396410592"
## [1] ...
## [1] "step 42 of loop. Things look good."
## [1] "computing likelihood for p = 0.41"
## [1] "log-likelihood -2.19759568545622"
## [1] ...
## [1] "step 43 of loop. Things look good."
## [1] "computing likelihood for p = 0.42"
## [1] "log-likelihood -2.2861560491725"
## [1] ...
## [1] "step 44 of loop. Things look good."
## [1] "computing likelihood for p = 0.43"
## [1] "log-likelihood -2.37822899604707"
## [1] ...
## [1] "step 45 of loop. Things look good."
## [1] "computing likelihood for p = 0.44"
## [1] "log-likelihood -2.47384657639288"
## [1] ...
## [1] "step 46 of loop. Things look good."
## [1] "computing likelihood for p = 0.45"
## [1] "log-likelihood -2.57304890871019"
## [1] ...
## [1] "step 47 of loop. Things look good."
## [1] "computing likelihood for p = 0.46"
## [1] "log-likelihood -2.67588420461821"
## [1] ...
## [1] "step 48 of loop. Things look good."
## [1] "computing likelihood for p = 0.47"
## [1] "log-likelihood -2.7824088582735"
## [1] ...
## [1] "step 49 of loop. Things look good."
## [1] "computing likelihood for p = 0.48"
## [1] "log-likelihood -2.89268759964339"
## [1] ...
## [1] "step 50 of loop. Things look good."
## [1] "computing likelihood for p = 0.49"
## [1] "log-likelihood -3.00679371209474"
## [1] ...
## [1] "step 51 of loop. Things look good."
## [1] "computing likelihood for p = 0.5"
## [1] "log-likelihood -3.12480931582913"
## [1] "... exiting loop ..."

#identify the maximum log-likelihood value
l.max <- max(l[,2])

#determine the value of p that results in the maximum likelihood.  This is the maximum likelihood estimation of p.
p.max <- l[l[,2]==l.max,1]

#plot the log-likelihood for all p
plot(l[,1],l[,2],lwd=3,col=2, xlab="recombination fraction", ylab="log-likelihood", type="l")

#indicate where the maximum likelihood lies
abline(v=p.max)