Re: [R] practical help ... solving a system...

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Mon 20 Jun 2005 - 09:05:29 EST

          "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 Mon Jun 20 09:08:45 2005

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