[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