#MLE for pi
################################################################################
# Example to find maximum of likelihood function
w <- 4
n <- 10
L <- function(pi, w, n) {
pi^w * (1 - pi)^(n - w)
}
logL <- function(pi, w, n) {
w * log(pi) + (n - w) * log(1 - pi)
}
curve(expr = L(pi = x, w = w, n = n), xlim = c(0,1), xlab = expression(pi), ylab = "Likelihood", col = "red")
abline(v = w/n, lty = "dotted")
curve(expr = logL(pi = x, w = w, n = n), xlim = c(0,1), xlab = expression(pi), ylab = "log(L)", col = "red")
abline(v = w/n, lty = "dotted")
pi<-seq(from = 0.3, to = 0.5, by = 0.01)
data.frame(pi = pi, L = L(pi = pi, w = w, n = n), logL = logL(pi, w = w, n = n))
################################################################################
#Suppose we could not simply use w and only had y to work with
y<-c(1, 0, 0, 1, 0, 1, 0, 0, 0, 1) #Possible y vector
sum(y)
L2 <- function(pi, y) {
#prod(dbinom(x = y, size = 1, prob = pi)) #alternative
prod(pi^y * (1 - pi)^(1 - y))
}
logL2 <- function(pi, y) {
#sum(dbinom(x = y, size = 1, prob = pi, log = TRUE)) #alternative
sum(y * log(pi) + (1 - y) * log(1 - pi))
}
L2(pi = sum(y)/n, y = y)
logL2(pi = sum(y)/n, y = y)
L(pi = w/n, w = w, n = n)
logL(pi = w/n, w = w, n = n)
################################################################################
#Work with first derivative
partial.logL <- function(pi, w, n) {
w/pi - (n - w)/(1 - pi)
}
pi<-seq(from = 0.3, to = 0.5, by = 0.01)
data.frame(pi = pi, L = partial.logL(pi = pi, w = w, n = n))
curve(expr = partial.logL(pi = x, w = w, n = n), xlim = c(0,1), xlab = expression(pi),
ylab = "Partial derivative of log(L)", col = "red")
abline(v = w/n, lty = "dotted")
abline(h = 0, lty = "dotdash")
################################################################################
# Part of the code is from Appendix B of Bilder and Loughin (2014)
#Initialize some values
epsilon<-0.0001 #Convergence criteria
pi.hat<-0.3 #Start value
save.pi.hat<-pi.hat #Save the results for each iteration here and put pi.hat as first value
change<-1 #Initialize change for first time through loop
#Loop to find the MLE (uses Newton-Raphson algorithm)
while (abs(change) > epsilon) {
num<-(w-n*pi.hat) / (pi.hat*(1-pi.hat))
den<- -w/pi.hat^2 - (n-w)/(1-pi.hat)^2
pi.hat.new<-pi.hat - num/den #pi^(r+1) = pi^(r) - (1st der log(L))/(2nd der log(L))
change<-pi.hat.new-pi.hat #-num/den
pi.hat<-pi.hat.new #Same for next time through loop
save.pi.hat<-c(save.pi.hat, pi.hat) #Keeps iteration history
}
#Print iteration history
data.frame(iteration = 1:length(save.pi.hat), save.pi.hat)
#Log-likelihood function plot
curve(expr = L(pi = x, w = w, n = n), xlim = c(0,1), xlab = expression(pi), ylab = "Likelihood", col = "red")
abline(v = w/n, lty = "dotted")
points(x = save.pi.hat, y = L(pi = save.pi.hat, w = w, n = n), pch = 1) #Each pi^(r)
#Zoomed in Log-likelihood function plot - use when first pi.hat = 0.3
curve(expr = logL(pi = x, w = w, n = n), from = 0.25, to = 0.5,
xlab = expression(pi), ylab = "Log-likelihood function")
points(x = save.pi.hat, y = logL(pi = save.pi.hat, w = w, n = n), pch = 1) #Each pi^(r)
segments(x0 = save.pi.hat, y0 = -10, x1 = save.pi.hat, y1 = logL(pi = save.pi.hat, w = w, n = n),
lty = "dotted") #Could choose a different y0 to make this more general
#first derivation of log-likelihood function
curve(expr = partial.logL(pi = x, w = w, n = n), from = 0.2, to = 0.6,
xlab = expression(pi), ylab = "First derivative of log-likelihood function")
abline(h = 0, col = "blue", lty = "dotted") #We are trying to find the value of pi that makes the first der of L(pi) equal to 0
points(x = save.pi.hat, y = partial.logL(pi = save.pi.hat, w = w, n = n), pch = 1) #Each pi^(i)
#Linear approximation using 1st order Taylor's series expansion: f(x) = f(x0) + (x-x0)*fp(x0) where x0 is what is being expanded
# about and fp() is first derivative of f()
x0<-0.3 #Start value
fp.x0<- -( (1-x0)^2*w + x0^2*(n-w) ) / (x0^2*(1-x0)^2) #Slope of tangent line at x0
f.x0<-(w-n*x0) / (x0*(1-x0))
yint<-f.x0 - x0*fp.x0
abline(a = yint, b = fp.x0, lty = "dashed", col = "red")
legend(x =0.35, y = 14, bty = "n", legend = c(expression(paste(frac(d, paste(d,pi)), " ", paste("log[L(", pi, "|", bold(y), ")]"))),
expression(paste("Tangent line at ", pi==0.3)), expression(paste(pi, " estimates"))), lty = c("solid", "dashed", NA), col = c("black", "red"),
pch = c(NA, NA, 1))
#Where the tangent line intersects x-axis = 0 corresponds to pi^(r+1)
##############################################
#Additional examples
#Note necessarily meant for these simple of problems
optim(par = 0.3, fn = logL, w = w, n = n, method = "BFGS", control = list(fnscale = -1))
#A function that works for one parameter
optimize(f = logL, interval = c(0,1), w = w, n = n, maximum = TRUE)
#