"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