[R] Warnings in fitdistr() from MASS.

[Rolf Turner]

For a lark, I experimented a bit with the data from Ted Byers' recent
postings.  The result of fitdistr() seemed sensible, but I was bothered
by the warnings about NaNs that arose.  Warnings always make me nervous.

Explicitly this is what I did:

TXT <- "0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0  
0 0 0 0 0
         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1  
1 1 1 1 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1  
1 1 1 1 1
         1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2  
2 2 2 2 2
         2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3  
3 3 3 3 3
         3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5  
5 5 5 5 5
         6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7  
7 7 7 7 7
         7 7 8 8 8 8 9 9 9 9 9 9 9 9 9 10 10 10 10 10 10 10 10 10 11  
11 11 11
         11 11 11 11 11 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12  
12 12 12
         12 18 18 18 18 18 18 18 18 18 18 18 18 18 18"

x <- scan(textConnection(TXT))
closeAllConnections()
try.x <- fitdistr(x,"negative binomial")

Two warnings about NaNs being produced resulted.

Digging into the code with browser() revealed that in the  
optimization process
negative values of "size" were tried on occasion, and this was giving  
the NaNs.

Basically I'm sending this out so that maybe those who are like me  
and are
made nervous by warnings will be able to search the archives and find  
reassurance
that all is actually well.

To keep the warnings from the door, one can set an argument "lower"  
in the
call to fitdistr, e.g.

	eps <- sqrt(.Machine$double.eps)
	fitdistr(x,"negative binomial",lower=c(eps,eps))

Note that setting lower=c(0,0) doesn't work --- you get an *error* to  
the
effect that L-BFGS-B needs finite values of 'fn'!!!

I also tried building my own local version of fitdistr() which had

	if(distname == "negative binomial" & is.null(Call$lower))
         	Call$lower <- rep(sqrt(.Machine$double.eps),2)

just after the assignment ``Call$hessian <- TRUE''.  This *seemed*
to work (i.e. prevent those nervous-making warnings and still give
the right answer).

	cheers,

		Rolf Turner


######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}