Mle of Weibull parameters using optim function.

[hasibul212003]

I was trying to obtain mles of the scale and shape parameter for the weibull distribution using the optim function.
The code is as follows:
weibull_data= rweibull(1000,shape=2,scale=3)
log_lik_weibull=function(par,obs_data){
n=length(obs_data)
a=par[1]
b=par[2]
n*(log(a)-log(b))+(a-1)(sum(log(obs_data))-nlog(b))-(a/b)*
sum(obs_data)
}

optim(c(1,1),log_lik_weibull,obs_dat=weibull_data,control = list(fnscale=-1))$par
This gives parameter estimates of 0.86, 2.63 which is significantly different from my shape and scale parameters.
I feel the log likelihood function is correct however I am not able to obtain reliable estimates.
Any help will be appreciated.

[Yarnabrina]

1. You missed a * in between (a - 1) and the next part in parenthesis.
2. I hope that you are aware of MASS::fitdistr.
3. Your log-likelihood function is incorrect. The last term will not be (a / b) * sum(obs_data), as a is in the power of (x / b).

Check this:

simulated_data_Weibull <- rweibull(n = 1000,
                                   shape = 2,
                                   scale = 3)

log_likelihood_Weibull <- function(parameters, observed_data)
{
  a <- parameters[1]
  b <- parameters[2]
  n <- length(x = observed_data)
  return((n * log(x = a)) + ((a - 1) * sum(log(x = observed_data))) - (a * n * log(x = b)) - (sum(observed_data ^ a) / (b ^ a)))
}

optimisation_results <- optim(par = c(1, 1),
                              fn = log_likelihood_Weibull,
                              observed_data = simulated_data_Weibull,
                              control = list(fnscale = -1))

optimisation_results$par
#> [1] 2.006208 2.982978

^{Created on 2019-08-22 by the reprex package (v0.3.0)}

Hope this helps.

[hasibul212003]

Thank you so much. My log likelihood function was wrong. I had the power in the wrong place. I corrected it:
weibull_data= rweibull(1000,shape=2,scale=3)
log_lik_weibull=function(par,obs_data){
n=length(obs_data)
a=par[1]
b=par[2]
sumobs=sum(obs_data)
sumlog_obs= sum(log(obs_data))
n*(log(a)-log(b))+(a-1)(sumlog_obs-nlog(b))-(1/b^a)*sum(obs_data^a)
}

optim(c(1,1),log_lik_weibull,obs_dat=weibull_data,method='BFGS',control = list(fnscale=-1))$par

Got 1.93 & 2.96 now. Didn't know about the fitdistr. I will look into it.
About missing the * it's in my code somehow when I copy the code here it's not showing.