################################################################################ #Bernoulli distribution # ################################################################################ ################################################################################ # Initial work w <- 4 n1 <- 10 #Using n1 rather than n to avoid confusion with n argument in fderiv() L <- function(pi, w, n1) { pi^w * (1 - pi)^(n1 - w) } logL <- function(pi, w, n1) { w * log(pi) + (n1 - w) * log(1 - pi) } curve(expr = L(pi = x, w = w, n1 = n1), xlim = c(0,1), xlab = expression(pi), ylab = "Likelihood", col = "red") abline(v = w/n1, lty = "dotted") curve(expr = logL(pi = x, w = w, n1 = n1), xlim = c(0,1), xlab = expression(pi), ylab = "log(L)", col = "red") abline(v = w/n1, lty = "dotted") ################################################################################ # Derivatives library(pracma) partial.logL <- function(pi, w, n1) { w/pi - (n1 - w)/(1 - pi) } curve(expr = partial.logL(pi = x, w = w, n1 = n1), xlim = c(0,1), xlab = expression(pi), ylab = "Partial derivative of log(L)", col = "red", lwd = 2) abline(v = w/n1, lty = "dotted") abline(h = 0, lty = "dotdash") partial.logL2 <- function(pi, w, n1) { fderiv(f = logL, x = pi, n = 1, w = 4, n1 = 10) } partial.logL2(pi = 0.4, w = 4, n1 = 10) #At MLE partial.logL2(pi = c(0.4,0.5), w = 4, n1 = 10) # Confirms numerical derivative is same curve(expr = partial.logL2(pi = x, w = w, n1 = n1), col = "blue", add = TRUE, xlim = c(0.05, 0.95), lty = "dashed", lwd = 2) uniroot(f = partial.logL2, interval = c(0.01,0.99), w = 4, n1 = 10) #