[R] Problem with mle
Rainer M KRug
RMK at krugs.de
Fri Jun 2 17:02:05 CEST 2006
R 2.3.0
Linux, SuSE 10.0
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
In addition: Warning messages:
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
In addition: Warning messages:
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)
)
######################################
More information about the R-help
mailing list