[R] Problem with mle

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


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)
			 )
######################################

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sat Jun 03 01:18:38 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sat 03 Jun 2006 - 06:10:30 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.