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

[Spencer Graves]

	  "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 at 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