Plot the likelihood of weibull

Question

I would like to plot the likelihood function of a size 1000 weibull sample with a sequence of shape parameter theta. I have used standardised weibull so the scale lambda is 1. However the output is a horizontal straight line.

n<-1000
lik <- function(theta, x){
  K<- length(theta)
  n<- length(x)
  out<- rep(0,K)
  for(k in 1:K){
    out[k] <- prod(dweibull(x, shape= theta[k], scale=1))   
  }
  return(out)
}
theta<-seq(0.01, 10, by = 0.01)
x <- rweibull(n, shape= 0.5, scale= 1)
plot(theta, lik(theta, x), type="l", lwd=2)

Try plotting the log0likelihood; sum(dweibull(x, shape= theta[k], scale=1, log=TRUE)) and keep theta in a more plausible range seq(0.01, 1, by = 0.01) — user20650

user20650 user20650 · Accepted Answer · 2020-05-23T17:03:31

There is nothing really wrong about what you have done but computers struggle to calculate the product of many small numbers and so can end up as zero (even 0.99^1000 = 4^-5). And so it is easier to log transform and then sum. (As the log transform is a monotonic increasing function maximising the log-likelihood is the same as maximising the likelihood).Thus change

prod(dweibull(x, shape= theta[k], scale=1))

to

sum(dweibull(x, shape= theta[k], scale=1, log=TRUE))

The other minor change is to plot the likelihood witihin a reasonable range of theta so that you can see the curve.

Working code:

set.seed(1)
n<-1000
lik <- function(theta, x){
  K <- length(theta)
  n <- length(x)
  out <- rep(0,K)
  for(k in 1:K){
    out[k] <- sum(dweibull(x, shape= theta[k], scale=1, log=TRUE))   
  }
  return(out)
}

popTheta = 0.5
theta = seq(0.01, 1.5, by = 0.01)
x = rweibull(n, shape=popTheta, scale= 1)
plot(theta, lik(theta, x), type="l", lwd=2)
abline(v=popTheta)

theta[which.max( lik(theta, x))]

Plot the likelihood of weibull

1 Answers