I've been trying to estimate the parameters of a reliability distribution called Kumaraswamy Inverse Weibull (KumIW), which can be found in the package 'RelDists'. I tried to use the function optim, but it doesn't matter the start that I give, the data generated by 'rKumIW' or the sample size, the mu parameter always converge to one, although the others parameters converge to the correct values. What can be done? Should I try to use another function from some different package?
library(RelDists)
library(gamlss)
ll <- function(par, x){
mu <- par[1]
sigma <- par[2]
nu <- par[3]
ll = n*log(sigma) + n*log(nu) + n*log(mu) -
(mu+1)*sum(log(x)) - sigma*sum(x^-mu) +
(nu-1)*sum(log(1-exp(-sigma*(x^-mu))))
return (ll)
}
n = 10000
set.seed(123)
y = rKumIW(n, mu = 2, sigma = 2.3, 3.5)
start = c(4, 2, 3)
fit1 <- optim(start, ll, x = y, method = "L-BFGS-B", hessian = FALSE, control = list(fnscale = -1),
lower = c(0, 0, 0), upper = c(Inf, Inf, Inf))
fit1
$par
[1] 0.9991335 2.3170352 3.5874946
$value
[1] -12489.38
$counts
function gradient
27 27
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
There's a bug in the quantile function qKumIW(), which is also used inside the random deviate generator rKumIW(), which leads to mu always being set to 1 regardless of the value specified ... see here. If you install my patched version (from a clean R session to make sure the new version gets loaded) everything seems to work consistently ...
## install patched version
remotes::install_github("bbolker/RelDists")
library(RelDists)
## bbmle package isn't necessary, but allows a more
## compact example ...
library(bbmle)
n = 10000
set.seed(123)
y = rKumIW(n, mu = 2, sigma = 2.3, 3.5)
dd <- data.frame(y = y)
m1 <- mle2(y ~ dKumIW(mu, sigma, nu),
data = dd,
start = list(mu = 4, sigma = 2, nu = 3),
method = "L-BFGS-B",
lower = c(0,0,0))
coef(m1)
## mu sigma nu
## 1.998208 2.317107 3.587746