Method of Moments Estimation Using R

Finding the MOM estimators using R can be boiled down to finding the root of a function.

If we can frame the problem in terms of a single function we can use uniroot(),

uniroot(f, interval,...)

which numerically approximates a solution to f\( (x) = 0 \) with \( x \) ithin the interval specified by interval.

MOM equations (from p. 182): \[ \begin{aligned} \frac{\widehat{\alpha}}{\widehat{\alpha} + \widehat{\beta}} &= \overline{X}_n\\ \dfrac{\widehat{\alpha}\widehat{\beta}}{(\widehat{\alpha} + \widehat{\beta})^2(\widehat{\alpha} + \widehat{\beta} + 1)} &= \dfrac{1}{n} \sum_{i=1}^n (X_i - \overline{X}_n)^2 = \dfrac{n-1}{n}S^2_n = \widehat{\mu}\prime_2 \end{aligned} \]

We can re-express \( \widehat{\mu}\prime_2 \) as

\[ \widehat{\mu}\prime_2 = \dfrac{R \widehat{\alpha}}{ \left( \dfrac{\widehat{\alpha}}{\overline{X}_n} \right)^2 \left( \dfrac{\widehat{\alpha}}{\overline{X}_n} \right)} \]

where \( R = \dfrac{1}{\overline{X}_n} - 1 = \dfrac{\widehat{\beta}}{\overline{X}_n} \)

beta.mom <- function(x, lower = 0.01, upper = 100) {
    x.bar <- mean(x)
    n <- length(x)
    v <- var(x) * (n - 1)/n
    R <- 1/x.bar - 1

    f <- function(a) {
        # note: undefined when a=0
        R * a^2/((a/x.bar)^2 * (a/x.bar + 1)) - v
    }

    u <- uniroot(f, c(lower, upper))

    return(c(shape1 = u$root, shape2 = u$root * R))
}

x <- rbeta(50,2,5) # random draws from a beta
head(x)

[1] 0.2188 0.2225 0.5897 0.1938 0.2558 0.5422

beta.mom(x)

shape1 shape2 
 3.333  8.062 

Of course, it's not always that easy, but this isn't a course in programming so uniroot is useful enough for now.

We will consider a few optimization tools in R when we get to maximum likelihood estimation.

To assess the fit of our model, a good place to start is to compare what we have observed to what we expect.

This can be done using a quantile-quantile (Q-Q) plot.

library(fastR) # load for plotting function

x <- rbeta(50, 2, 5)
(mome <- beta.mom(x))

shape1 shape2 
 2.725  7.585 

# The quantile function
qbeta.plot <- function(x) { qbeta(x, shape1 = mome[1], shape2 = mome[2]) }

Let \( X \) denote the proportion of allotted time that a randomly selected student spends working on a certain aptitude test. Suppose the PDF of \( X \) is

\[ f(x|\theta) = (\theta + 1) x^\theta, \quad 0 \le x \le 1 \]

where \( -1 < \theta \). A random sample of ten students yields the following data:

x <- c(.92, .79, .90, .65, .86, .47, .73, .97, .94, .77)

We wish to use MOM to estimate \( \theta \).

First notice that

\[ E(X_1) = \displaystyle \int_0^1 (\theta + 1) x^{\theta+1} dx = \frac{\theta + 1}{\theta + 2} \]

Now finding the MOM estimator of \( \theta \) is quite easy

\[ \frac{\widehat{\theta} + 1}{\widehat{\theta} + 2} = \overline{X}_n \Longrightarrow \widehat{\theta} = \frac{2 \overline{X}_n - 1}{1-\overline{X}_n} \]

xbar <- mean(x)
theta_hat <- (2 * xbar - 1) / (1 - xbar)
theta_hat

[1] 3

Now, notice that

\[ F(x) = \int_0^x (\theta + 1) t^\theta dt = x^{\theta + 1}, \quad 0 \le x \le 1 \]

So we the quantile function is

\[ x_p = p^{\frac{1}{\theta + 1}} \]

Which is easy to program in R

qaptitude <- function(p, theta = 3) {p^(1/(theta + 1))}