From: Spencer Graves <spencer.graves_at_pdf.com>

Date: Tue 21 Jun 2005 - 02:02:50 EST

There were 15 warnings (use warnings() to see them)

> qchisq(0.95, 1)

[1] 3.841459

> 2*(max(lglk)-lglk)

Carsten Steinhoff wrote:

> Hello Spencer,

>

> thank you for your reply. As I've seen you have optimized only one

value of

* > the binomial distribution by using Fitdistr.
*

* > My problem is that I have to find both size AND prob.
*

> The error message always is: error in using optim. non-finite finite

> difference value [1]

> So I thought to do a MLE "by hand" ...

> Do you have an idea what do modify on the error above?

>

> Greetings, Carsten

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Jun 21 02:17:59 2005

Date: Tue 21 Jun 2005 - 02:02:50 EST

The problem is that "size" in "dbinom" must be an integer:

> dbinom(x=2, size=4.5, prob=.5)

[1] NaN

Warning message:

NaNs produced in: dbinom(x, size, prob, log)

Consider the following:

> library(MASS)

* > set.seed(123)
** > x <- rbinom(100, 10, 0.4)
** > nx <- max(x)
** > fit <- fitdistr(x, dbinom, start=list(prob=0.5),
*

+ size = nx)

Warning message:

one-diml optimization by Nelder-Mead is unreliable: use optimize in:
optim(start, mylogfn, x = x, hessian = TRUE, ...)

> attributes(fit)

$names

[1] "estimate" "sd"

$class

[1] "fitdistr"

#### From this, we learn that fit$est is the MLE for prob given size
#### Try different a range of values for "size"

> lglk <- rep(NA, 12)

#### Start with nx, since it can't be less than that.

> names(lglk) <- nx:(nx+11)

> for(i in 1:12){

+ fit <- fitdistr(x, dbinom, start=list(prob=0.5), + size = nx+i-1) + lglk[i] <- sum(dbinom(x, size=nx+i-1, + prob=fit$estimate, log=TRUE)) + }

There were 15 warnings (use warnings() to see them)

> qchisq(0.95, 1)

[1] 3.841459

> 2*(max(lglk)-lglk)

8 9 10 11 12 13 1.314439853 0.000000000 0.004998298 0.434260975 1.000371029 1.594030787

14 15 16 17 18 19 2.170873826 2.713847684 3.217152046 3.680715184 4.106493564 4.497508846

# 2*log(likelihood ratio) is approximately chi-square, # so a 95% confidence interval for "size" is 8:17.

If I needed this for a scientific research paper, I'd do more research on how to estimate "size" for the binomial -- if I really thought it should be estimated. If I needed an answer today, I'd probably go with what I just did here.

hope this helps, spencer graves ##################

Carsten Steinhoff wrote:

> Hello Spencer,

>

> thank you for your reply. As I've seen you have optimized only one

value of

> The error message always is: error in using optim. non-finite finite

> difference value [1]

>

> So I thought to do a MLE "by hand" ...

>

> Do you have an idea what do modify on the error above?

>

> Greetings, Carsten

Spencer Graves wrote:

> "PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html". If you had provided more > information about what you tried, it might be easier for someone else to > offer effective help. In particular, what did you try with fitdistr, > and why do you think it didn't work (e.g., what error message did you get)? > > The following is a minor modification of examples in the > "fitdistr" help page: > > > library(MASS) > > set.seed(123) > > x <- rbinom(100, 10, 0.4) > > (fit <- fitdistr(x, dbinom, start=list(prob=0.5), size=10)) > prob > 0.39804687 > (0.01547943) > Warning message: > one-diml optimization by Nelder-Mead is unreliable: use optimize in: > optim(start, mylogfn, x = x, hessian = TRUE, ...) > > > The difference between prob=0.4 that I simulated and the > 0.39804687 estimated is shockingly small. (No, I didn't run 100 > simulations with different seeds and pick the best one.) > > If I were concerned about the warning, I could list "fitdist" and > read the code. I might even copy it into a script file and walk through > it line by line. When checked the code for "fitdistr" just now, I > learned that it calls "optim", and "optim" provides other estimation > methods, e.g., BFGR and CG. The following are the results from trying > those two: > > (fit.BFGS <- fitdistr(x, dbinom, start=list(prob=0.5), > + size=10, method="BFGS")) > prob > 0.39800028 > (0.01547882) > Warning messages: > 1: NaNs produced in: dbinom(x, size, prob, log) > 2: NaNs produced in: dbinom(x, size, prob, log) > 3: NaNs produced in: dbinom(x, size, prob, log) > 4: NaNs produced in: dbinom(x, size, prob, log) > 5: NaNs produced in: dbinom(x, size, prob, log) > > (fit.CG <- fitdistr(x, dbinom, start=list(prob=0.5), > + size=10, method="CG")) > prob > 0.39800001 > (0.01547881) > Warning messages: > 1: NaNs produced in: dbinom(x, size, prob, log) > 2: NaNs produced in: dbinom(x, size, prob, log) > 3: NaNs produced in: dbinom(x, size, prob, log) > 4: NaNs produced in: dbinom(x, size, prob, log) > 5: NaNs produced in: dbinom(x, size, prob, log) > > If I'm concerned about these warnings, I can replace "dbinom" by a > local copy that prints the value of "prob" each time it's called: > > dbinom. <- function (x, size, prob, log = FALSE){ > cat(prob, "") > .Internal(dbinom(x, size, prob, log)) > } > (fit.BFGS <- fitdistr(x, dbinom., start=list(prob=0.5), > size=10, method="BFGS")) > 0.5 0.501 0.499 -407.5005 -81.10011 -15.82002 -2.764004 -0.1528009 > 0.3694398 0.3704398 0.3684398 0.3996073 0.4006073 0.3986073 0.3980445 > 0.3990445 0.3970445 0.2133659 0.3611088 0.3906574 0.3965671 0.3977490 > 0.3979854 0.3989854 0.3969854 0.3980003 0.45997 0.4103942 0.4004791 > 0.3984960 0.3980994 0.3980201 0.3980043 0.3980011 0.3980004 0.3980003 > 0.3980003 0.3980003 0.3980003 0.3980003 0.4000003 0.3980003 0.3980003 > 0.3960003 prob > 0.39800028 > (0.01547882) > Warning messages: > 1: NaNs produced in: dbinom(x, size, prob, log) > 2: NaNs produced in: dbinom(x, size, prob, log) > 3: NaNs produced in: dbinom(x, size, prob, log) > 4: NaNs produced in: dbinom(x, size, prob, log) > 5: NaNs produced in: dbinom(x, size, prob, log) > > > It's no wonder dbinom produced NaNs: It does that when prob < 0. > > hope this helps. > spencer graves > > p.s. MASS = "Modern Applied Statistics with S" by Venables and Ripley.. > Have you seen this book? I've learned a lot from it. > > Carsten Steinhoff wrote: > >> Hello, >> >> I want to estimate the parameters of a binomial distributed rv using MLE. >> Other distributions will follow. >> The equation system to solve is not very complex, but I've never done >> such >> work in R and don't have any idea how to start... >> >> The system is: >> >> (1) n*P = X >> >> (2) [sum {from j=0 to J-1} Y{j} /(n-j)] = -n * ln (1-X / n) >> >> >> where * only X is given (empirical mean) >> * J is maximum observed >> * Y is the number of observations above j >> * n and P are the parameters of the binomial distribution to >> find.... >> >> Who could help me with an example-solution for my "first" distribution >> ? I >> also need a hint how to make the sum-element. >> >> Maybe there's another - more simple - way to estimate the parameters... >> first I tried via "FITDISTR" but without success. >> >> Thanx a lot for your help. >> >> Carsten >> >> [[alternative HTML version deleted]] >> >> ______________________________________________ >> 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 > > ______________________________________________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 Tue Jun 21 02:17:59 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:32:54 EST
*