[R] Problem with fitdistr for beta

[Prof Brian Ripley]

In this example shrinking by (1 - 2e-16) leads to a significant change in 
the distribution: see my probability calculation.  And you can't shrink by 
much less.  A beta(0.1, 0.1) is barely a continuous distribution.

On Fri, 4 Jul 2003, Spencer Graves wrote:

> My standard work-around for the kind of problem you identified is to 
> shrink the numbers just a little towards 0.5.  For example:
> 
>  > library(MASS)
>  > a <- rbeta(100,0.1,0.1)
>  > fitdistr(x=a, "beta", start=list(shape1=0.1,shape2=0.1))
> Error in optim(start, mylogfn, x = x, hessian = TRUE, ...) :
> 	Function cannot be evaluated at initial parameters
> 
>  > r.a <- range(a)
>  > c0 <- 0
>  > c1 <- 1
>  > if(r.a[1]==0)c0 <- min(a[a>0])
>  > if(r.a[2]==1)c1 <- max(a[a<1])
>  > c. <- c(c0, 1-c1)
>  > if(any(c.>0))c. <- min(c.[c.>0])
>  > c.
> [1] 2.509104e-14
>  > fitdistr(x=0.5*c.[1] + (1-c.[1])*a, "beta", 
> start=list(shape1=0.1,shape2=0.1))
>       shape1       shape2
>    0.08728921   0.10403875
>   (0.01044863) (0.01320188)
>  >
> hope this helps.  spencer graves
> 
> Prof Brian Ripley wrote:
> > rbeta(100,0.1,0.1) is generating samples which contain 1, an impossible 
> > value for a beta and hence the sample has an infinite log-likelihood.
> > It is clearly documented on the help page that the range is 0 < x < 1.
> > However, that is not so surprising as P(X > 1-1e-16) is about 1% and hence 
> > values will get rounded to one.
> > 
> > The same would happen for a value of 0.
> > 
> > Your code is syntactically incorrect, at least as received here.
> > 
> > On Fri, 4 Jul 2003, Agustin Lobo wrote:
> > 
> > 
> >>I have the following problem:
> >>
> >>I have a vector x of data (0<x<=1 ) with
> >>a U-shaped histogram and try to fit a beta
> >>distribution using  fitdistr. In fact,
> >>hist(rbeta(100,0.1,0.1)) looks a lot like
> >>my data.
> >>
> >>The equivalent to
> >>the example in the manual
> >>sometimes work:
> >>
> >>
> >>>a <- rbeta(100,0.1,0.1)
> >>>fitdistr(x=a, "beta", start=list(shape1=0.1,shape2=0.1))1)
> >>>     shape1       shape2
> >>
> >>  0.09444627   0.12048753
> >> (0.01120670) (0.01550129)
> >>
> >>but sometimes does not:
> >>
> >>>a <- rbeta(100,0.1,0.1)
> >>>fitdistr(x=a, "beta", start=list(shape1=0.1,shape2=0.1))1)
> >>>Error in optim(start, mylogfn, x = x, hessian = TRUE, ...) :
> >>
> >>        Function cannot be evaluated at initial parameters
> >>
> >>Unfortunately, my data fall in the second case
> >>
> >>I've searched for any weird value that be present in the
> >>cases in which fitdistr exits with the error message, but
> >>could not find any.
> >>
> >>Any help?
> >>(please if anyone answers be sure to answer to my address as well,
> >>I cannot subscribe to the list)
> >>
> >>Thanks
> >>
> >>Agus
> >>
> >>Dr. Agustin Lobo
> >>Instituto de Ciencias de la Tierra (CSIC)
> >>Lluis Sole Sabaris s/n
> >>08028 Barcelona SPAIN
> >>tel 34 93409 5410
> >>fax 34 93411 0012
> >>alobo at ija.csic.es
> >>
> >>______________________________________________
> >>R-help at stat.math.ethz.ch mailing list
> >>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> >>
> > 
> >
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> 

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595