# How to pick starting parameters for MASS:fitdist() with the beta distribution?

I have data set of ~700k yes/no events that I want to first aggregate on various features (e.g. group by average), always resulting in a 34 length vector. From there, I want to fit a beta distribution to the resulting vector. Below is an example of one possible vector:

dput(0.0955104277250779, 0.0782381918284555, 0.109683584625186,
0.10115721657354, 0.102377369846524, 0.0691699604743083, 0.0940254652301665,
0.078494747777906, 0.0824474231569216, 0.087598944591029, 0.114703196347032,
0.0966910484151863, 0.0995783979120659, 0.090916800949877, 0.111081696404935,
0.101724137931034, 0.0880689367194321, 0.106699751861042, 0.0772264084644516,
0.109373466383118, 0.0732657833203429, 0.139175257731959, 0.0902068720589541,
0.0900240916465903, 0.103787190855189, 0.0846888376687521, 0.0901782178217822,
0.0943396226415094, 0.0778525161933234, 0.057852491291837, 0.0902365266227454,
0.0940445474171976, 0.0834914611005693, 0.0886513916653852)


From there I want to fit a beta distribution (this is leading up to empiricial bayesian estimation btw). See http://varianceexplained.org/r/empirical_bayes_baseball/ for the general idea.

When I try and fit a beta distribution however, I receive errors as shown below.

library(MASS)

m2 <- MASS::fitdistr(vector, densfun = "beta",  start = list(shape1 = 1, shape2 = 10))

"Error in stats::optim(x = c(0.0955104277250779, 0.0782381918284555, 0.109683584625186,  :
non-finite finite-difference value [1]:"


Although I have been successful in fitting other vectors with the fitdist() function, I have been unable to fit this specific vector, although the data appears rational and well formed. I have encountered a wide variety of errors when aggregating/subsetting the data in different ways as well. In general, I'd like to be able to set the start parameters in such a way as to guarantee a fit, even if the resulting shape is less than ideal. The goal being to remove the needed trial and error tweaking of the start values.

I've read the fitdist documentation thoroughly, and searched for a similiar question/issue with no luck.

Happy to provide additional example vectors as requested.

• Have you tried different starting values, e.g., ones based on methods-of-moments estimates? Otherwise, try changing the initial parameters to be non-integer valued, e.g., $(1.1, 10.1)$ instead of $(1,10)$. Nov 12 '18 at 19:22

Using fitdistrplus works well in this case:

library(fitdistrplus)

x <- c(0.0955104277250779, 0.0782381918284555, 0.109683584625186,
0.10115721657354, 0.102377369846524, 0.0691699604743083,
0.0940254652301665, 0.078494747777906, 0.0824474231569216,
0.087598944591029, 0.114703196347032, 0.0966910484151863,
0.0995783979120659, 0.090916800949877, 0.111081696404935,
0.101724137931034, 0.0880689367194321, 0.106699751861042,
0.0772264084644516, 0.109373466383118, 0.0732657833203429,
0.139175257731959, 0.0902068720589541, 0.0900240916465903,
0.103787190855189, 0.0846888376687521, 0.0901782178217822,
0.0943396226415094, 0.0778525161933234, 0.057852491291837,
0.0902365266227454, 0.0940445474171976, 0.0834914611005693,
0.0886513916653852)

fit <- fitdist(x, "beta")

Fitting of the distribution ' beta ' by maximum likelihood
Parameters:
estimate Std. Error
shape1  35.25746   8.511331
shape2 344.99152  83.815718

plot(fit, las = 1) # Checking fit


The function fitdist is able to chose "reasonable" starting values on its own, whereas fitdistr (MASS) struggles. We can get fitdistr to run without errors by supplying it reasonable starting values (but I'd recommend using the fitdistr package anyway):

library(MASS)

fitdistr(x, "beta", start = list(shape1 = 35, shape2 = 344))

shape1      shape2
35.15658   343.98481
(  8.48699) ( 83.57280)


As @jbowman suggested in the comments, you can use the method-of-moments estimates as starting values. For the beta distribution, the corresponding estimators are:

\begin{align*} \hat{\alpha} &= \bar{x}\left(\frac{\bar{x}(1 - \bar{x})}{s^2} - 1\right) \\ \hat{\beta} &= (1 - \bar{x})\left(\frac{\bar{x}(1 - \bar{x})}{s^2} - 1\right) \end{align*}

Where $$\bar{x}$$ is the sample mean and $$s^2$$ is the sample variance. In R:

beta_mom <- function(x) {

m_x <- mean(x, na.rm = TRUE)
s_x <- sd(x, na.rm = TRUE)

alpha <- m_x*((m_x*(1 - m_x)/s_x^2) - 1)
beta <- (1 - m_x)*((m_x*(1 - m_x)/s_x^2) - 1)

return(list(alpha = alpha, beta = beta))

}

beta_mom(x)

$alpha [1] 33.97962$beta
[1] 332.4864


These are excellent starting values for the optimization algorithms.

• Thank you for the helpful comment. I noted that the fitdist() function still fails on many possible inputs for method = mle, but seems to work better with method = mge. Im doing thorough testing and will update (and close if appropriate) this question shortly. Nov 15 '18 at 22:52
• @LMico And what about the proposed method of estimating the parameters first with the method-of-moments and then supply those to fitdist? Nov 16 '18 at 7:21