Plots of Likelihood and log-likelihood function of Bernoulli random variable-CodePudding

data(quilpie, package = "GLMsData"); names(quilpie)
mu = seq(0.01, 0.99, length.out = 100L)
loglik = NULL
lik =NULL
for (i in 1:length(mu)) {
  loglik[i] <- sum(dbinom(quilpie$y, size = 1, prob=mu[i], log=TRUE))
  lik[i] <- sum(dbinom(quilpie$y, size=1, prob=mu[i], log=FALSE))
}
plot(mu, lik, type = "l", lwd = 2,
     xlab = expression(paste(,mu ,)),
     main = "Likelihood function",
     ,ylab = expression(paste("loglik (", mu ,")")))
abline(v = mu[which(lik==max(lik))], lwd = 3, lty =2)
abline(v = optimize(lik, lower = 0, upper = 1, maximum = TRUE)$max)
optimize(lik, lower = 0.01, upper = 0.99, maximum = TRUE)
mu[which(lik==max(lik))]

plot(mu, loglik, type = "l", lwd = 2,
     xlab = expression(paste(,mu ,)),
     main = "Log-likelihood function",
     ,ylab = expression(paste("loglik (", mu ,")")))
abline(v = mu[which(loglik==max(loglik))], lwd = 3, lty =2)
abline(v = optimize(loglik, lower = 0, upper = 1, maximum = TRUE)$max)
optimize(loglik, lower = 0.01, upper = 0.99, maximum = TRUE)
mu[which(loglik==max(loglik))]

I am trying to obtain plots for both the likelihood and the log-likelihood. The log-likelihood plot works well but the likelihood function plot doesn't show anything. Also, the optimize function doesn't work. I need help with this one.

CodePudding user response：

Rewrite the computations of loglik and lik as functions, vectorize them, then optimize. And plot with curve.

Note that the functions are to be optimized on the parameter prob, here as first argument, named x.
Argument y is the data.

floglik <- function(x, y){
  sum(dbinom(y, size = 1, prob = x, log = TRUE))
}
floglik <- Vectorize(floglik, "x")

max_loglik <- optimize(floglik, interval = c(0, 1), y = quilpie$y, maximum = TRUE)

curve(floglik(x, quilpie$y), from = 0, to = 1,
      lwd = 2,
      xlab = expression(paste(mu)),
      main = "Log-likelihood function",
      ylab = expression(paste("loglik (", mu ,")")))
abline(v = max_loglik$maximum)

And the same for the likelihood.

flik <- function(x, y){
  prod(dbinom(y, size = 1, prob = x, log = FALSE))
}
flik <- Vectorize(flik, "x")
max_lik <- optimize(flik, interval = c(0, 1), y = quilpie$y, maximum = TRUE)

curve(flik(x, quilpie$y), from = 0, to = 1,
      lwd = 2,
      xlab = expression(paste(mu)),
      main = "Likelihood function",
      ylab = expression(paste("lik (", mu ,")")))
abline(v = max_lik$maximum)