# [R] Problem with mle

From: Rainer M KRug <RMK_at_krugs.de>
Date: Sat 03 Jun 2006 - 01:02:05 EST

Hi

I have two problems with mle - probably I am using it the wrong way so please let me know.

I want to fit different distributions to an observed count of seeds and in the next step use AIC or BIC to identify the best distribution.

But when I run the script below (which is part of my original script), I get one error message for the first call of mle:

Error in solve.default(oout\$hessian) : Lapack routine dgesv: system is exactly singular

1: NaNs produced in: dweibull(x, shape, scale, log)
2: NaNs produced in: dweibull(x, shape, scale, log)
3: NaNs produced in: dweibull(x, shape, scale, log)
4: NaNs produced in: dweibull(x, shape, scale, log)

What can I do to avoid this problem?

The second one is following the second call of mle: summary(fit.NE) gives me the summary of the fit, but

profile(fit.NE) gives me the following error message: Error in onestep(step) : profiling has found a better solution, so original fit had not converged

1: NaNs produced in: dexp(x, 1/rate, log)
2: NaNs produced in: dexp(x, 1/rate, log)
3: NaNs produced in: dexp(x, 1/rate, log)
4: NaNs produced in: dexp(x, 1/rate, log)
5: NaNs produced in: dexp(x, 1/rate, log)
6: NaNs produced in: dexp(x, 1/rate, log)

What can I do in this case - which parameters do I have to change that it converges?

Rainer

######################################

library(stats4)

SumSeeds <- c(762, 406, 288, 260, 153, 116, 57, 33, 40, 44, 36, 24, 35, 23, 36, 25, 35, 30)
X <- c(1.74924, 3.49848, 5.24772, 6.99696, 8.74620, 17.49240, 26.23860, 34.98480, 43.73100, 52.47720, 61.22340, 69.96960, 71.71880, 73.46810, 75.21730, 76.96650, 78.71580, 80.46500 )

#Sum of log of values in vector
SumLog <- function (x)
{

sum(log(x))
}

#-ln(L) with Poisson Likelihood estimator NeglogPoisL <- function (obs, est)
{

sel <- est != 0
SumLogObs <- SumLog(obs[sel])
-sum( obs[sel] * log(est[sel]) - est[sel]- SumLogObs )
}

#-ln(L) for Negative Exponential
NENeglogPoisL <- function(a=0.2, rate=0.04)
{

NeglogPoisL( SumSeeds, a * dexp( X, rate) ) }

#-ln(L) for Weibull
WBNeglogPoisL <- function(a=1000000, shape=0.12, scale=1e+37)
{

NeglogPoisL( SumSeeds, a * dweibull(X, shape, scale) ) }

fit.WB <- mle(

WBNeglogPoisL
, start=list(a=1, shape=0.1, scale=1)
, fixed=list()
, control=list(maxit=10000)
)

fit.NE <- mle(
NENeglogPoisL
, start=list(a=1, rate=1)
, fixed=list()
, control=list(maxit=10000)
)
######################################

______________________________________________
R-help@stat.math.ethz.ch mailing list