33 Estimação pelo método de momentos (EMM)

Seja \(\boldsymbol{X}_n\) uma a.a. de \(X\) tal que \[ E_\theta (\lvert X \rvert^k) < \infty, k \in \{1,\dots,p\}, \Theta \subseteq \mathbb{R}^P \]

O estimador obtido pelo método de momentos é aquele que satisfaz \[ E_\theta(X^k) = \frac{1}{n} \sum^n_{i=1} X_i^k, k = 1, \dots, p. \]

Note que a estimação pelo método de momentos não utiliza toda a informação contida na função de verossimilhança. Precisamos conhecer a forma dos primeiros \(p\) momentos.

33.1 Exemplos

33.1.1 Exemplo 1

Seja \(\boldsymbol{X}_n\) a.a. de \(X\sim U(0,\theta), \theta > 0\). Encontre o estimador pelo método de momentos.

Calcule o EQM do EMV e do estimador obtido pelo método de momentos.

33.1.1.1 Resposta

Note que \(\Theta = (0,\infty) \subseteq \mathbb{R}\). Portanto, \(p=1\). Precisamos encontrar o valor de “\(\theta\)” que satisfaz \[ E_\theta(X) = \frac{1}{n} \sum X_i \]

Sabemos que \[ E_\theta(X) = \int^\theta_0 x \frac{1}{\theta} dx = \frac{\theta}{2} = \bar{X} \]

Logo, \(\hat{\theta}_{\mathrm{MM}}(\boldsymbol{X}_n) = 2\bar{X}\) é o estimador obtido pelo método de momentos.

Sabemos que \(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n) = \max\{X_1,\dots,X_n\}\), logo, \[ \begin{aligned} \mathrm{EQM}_\theta(\hat{\theta}_{\mathrm{MM}}(\pmb{X}_n), \theta) &= E_\theta((\hat{\theta}_{\mathrm{MM}}(\pmb{X}_n) - \theta)^2) \\ &=\mathrm{Var}_\theta(\hat{\theta}_{\mathrm{MM}}(\pmb{X}_n)) + \mathrm{Viés}(\hat{\theta}_{\mathrm{MM}}(\pmb{X}_n), \theta)^2 \\ \\ E_\theta(\hat{\theta}_{\mathrm{MM}}(\pmb{X}_n)) &= 2 E_\theta(X) = 2 \frac{\theta}{2} = \theta \\ \Rightarrow \mathrm{Viés}_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n),\theta) &= 0, \forall \theta \in \Theta. \\ \\ \mathrm{Var}_\theta(\hat{\theta}_{\mathrm{MM}}(\pmb{X}_n)) = \mathrm{Var}_\theta(2\bar{X}) &\stackrel{\mathrm{i.i.d.}}{=} 4 \frac{\mathrm{Var}_\theta(X)}{n} \stackrel{\mathrm{Unif}}{=} \frac{\theta^2}{3n} \\ \Rightarrow \mathrm{EQM}_\theta(\hat{\theta}_{\mathrm{MM}}(\pmb{X}_n),\theta) = \frac{\theta^2}{3n} \end{aligned} \]

Para o EMV, \[ \mathrm{EQM}_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n), \theta) = \mathrm{Var}_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n)) +\mathrm{Viés}(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n), \theta)^2 \\ \\ \]

Precisaremos encontrar a distribuição do estimador, \[ P_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n) \leq t) = P_\theta(\max\pmb{X}_n \leq t) \stackrel{\mathrm{iid}}{=} \prod P_\theta(X \leq t) = P_\theta(X \leq t)^n \] note que \[ P_\theta(X \leq t) = \left\{\begin{array}{ll} 0, & t \leq 0 \\ \frac{t}{\theta}, & 0 < t \leq \theta \\ 1, & t \geq \theta \\ \end{array}\right. \Rightarrow P_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n)\leq t)^n = \left\{\begin{array}{ll} 0, & t \leq 0 \\ \left(\frac{t}{\theta}\right)^n, & 0 < t \leq \theta \\ 1, & t \geq \theta \\ \end{array}\right. \]

Como \(P_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n))\) é absoulatemente contínua para todo \(\theta > 0\), temos que \(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n)\) é uma variável aleatória contínua cuja f.d.p. é dada por \[ f_\theta^{\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n)}(x) =\left. \frac{d}{dt} P_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n) \leq t) \right\rvert_{t=x} = \left\{\begin{array}{ll} \frac{n x^{n-1}}{\theta^n}, & 0 < x \leq \theta \\ 0, & \mathrm{c.c.} \end{array}\right. . \]

Logo, \[ \begin{aligned} E_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n)) &= \int^{\infty}_{-\infty} w f_\theta^{\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n)}(w) dw \\ &= \int^\theta_0 w \frac{n w^{n-1}}{\theta^n} dw \\ &= \frac{n}{\theta^n} \int^\theta_0 w^n dw = \left.\frac{n}{\theta^n}\frac{w^{n+1}}{n+1} \right\rvert_{0}^\theta = \frac{n}{n+1}\theta, \forall \theta \in \Theta, \\ \Rightarrow E_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n)^2) &= \frac{n\theta^2}{n+2}, \forall \theta \in \Theta, \\ \Rightarrow \mathrm{Var}_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n)) &= \frac{n\theta^2}{n+2} - \left(\frac{n}{n+1}\right)^2\theta^2 \\ \Rightarrow \mathrm{EQM}_\theta(\hat{\theta}_{\mathrm{MV}}(\pmb{X}_n),\theta) &= \frac{n}{n+2}\theta^2 - \frac{n^2}{(n+1)^2}\theta^2 + \left(\frac{n}{n+1}\theta - \theta\right)^2 \\ &= \theta^2\left(\frac{n}{n+2} - \frac{2}{n+1} + 1\right) \end{aligned} \]

33.1.1.2 Simulação

using Random, Distributions, StatsBase, Plots, LaTeXStrings

Random.seed!(76)

M = 10_000
n = 10
theta0 = rand(Poisson(4)) + 1

d = Uniform(0, theta0)

thetaMM = zeros(M)
thetaMV = zeros(M)
for i in 1:M
  x = rand(d, n)
  thetaMM[i] = 2*mean(x)
  thetaMV[i] = maximum(x)
end

# Comparar
println("Média simulada do EMM = $(mean(thetaMM)),
        calculada = $theta0")
println("Média simulada do EMV = $(mean(thetaMV)),
        calculada = $(10/11 * theta0)")
println("Variância simulada do EMM = $(var(thetaMM)),
        calculada = $(theta0^2 /(3*n))")
println("Variância simuladado EMV = $(var(thetaMV)),
        calculada = $(n * theta0^2 / (n+2) - (n/(n+1))^2 * theta0^2)")

densidade(x) = 0 < x < theta0 ? n * x^(n-1) / theta0^n : 0

pmm = histogram(thetaMM, label = "",
                title = "Histograma do EMM", normalize=:pdf)
pmv = histogram(thetaMV, label = "",
                title = "Histograma do EMV", normalize=:pdf)
plot!(minimum(thetaMV):0.001:maximum(thetaMV)-0.01, densidade,
      label="Densidade Teórica")
display(plot(pmm, pmv))

Média simulada do EMM = 4.012626175744703,
        calculada = 4
Média simulada do EMV = 3.6375235806062722,
        calculada = 3.6363636363636362
Variância simulada do EMM = 0.5384935090654293,
        calculada = 0.5333333333333333
Variância simuladado EMV = 0.11001021153100962,
        calculada = 0.11019283746556674

33.1.2 Exemplos adicionais

Seja \(\boldsymbol{X}_n\) a.a. de \(X \sim f_\theta, \theta \in \Theta\). Encontre o estimador pelo método de momentos para os casos abaixo:

\(X\sim \mathrm{Ber}(\theta), \theta \in (0,1)\)
\(X\sim \mathrm{Bin}(m, \theta), \theta \in (0,1)\)
\(X\sim \mathrm{Poiss}(\theta), \theta \in (0,\infty)\)
\(X\sim \mathrm{N}(\mu, \sigma^2), \theta = (\mu, \sigma^2) \in \mathbb{R}\times\mathbb{R}_+\)
\(X\sim \mathrm{Exp}(\theta), \theta > 0\)
\(X\sim \mathrm{Gamma}(\alpha, \beta), \theta = (\alpha, \beta) \in \mathbb{R}_+\times\mathbb{R}_+\)

33.1.2.1 Respostas

\(E_\theta(X) = \theta\). Portanto, \(\hat{\theta}_{\mathrm{MM}}(\boldsymbol{X}_n) = \bar{X}\).
\(E_\theta(X) = m \theta\). Portanto, \(\hat{\theta}_{\mathrm{MM}}(\boldsymbol{X}_n) = \frac{\bar{X}}{m}\).
\(E_\theta(X) = \theta\). Portanto, \(\hat{\theta}_{\mathrm{MM}}(\boldsymbol{X}_n) = \bar{X}\).
\(E_\theta(X) = \theta, E_\theta(X^2) = \sigma^2 + \mu^2\). Note que \[ \begin{aligned} &\left\{\begin{array}{ll} \mu &= \bar{X} \\ \sigma^2 + \mu^2 &= \frac{1}{n} \sum X_i^2 \end{array}\right. \\ \Rightarrow&\left\{\begin{array}{ll} \hat{\mu}_{\mathrm{MM}}(\boldsymbol{X}_n) &= \bar{X} \\ \hat{\sigma^2}_{\mathrm{MM}}(\boldsymbol{X}_n) + \hat{\mu}_{\mathrm{MM}}(\boldsymbol{X}_n)^2 &= \frac{1}{n} \sum X_i^2 \end{array}\right. \\ \Rightarrow&\left\{\begin{array}{ll} \hat{\mu}_{\mathrm{MM}}(\boldsymbol{X}_n) &= \bar{X} \\ \hat{\sigma^2}_{\mathrm{MM}}(\boldsymbol{X}_n) = \frac{1}{n} \sum X_i^2 - \bar{X}^2 \end{array}\right. \\ \Rightarrow& \hat{\theta}_{\mathrm{MM}}(\boldsymbol{X}_n) = \left(\bar{X}, \frac{1}{n} \sum X_i^2 - \bar{X}^2\right) \end{aligned} \]
\(E_\theta(X) = \frac{1}{\theta}\). Portanto, \(\hat{\theta}_{\mathrm{MM}}(\boldsymbol{X}_n) = \frac{1}{\bar{X}}\)
Distribuição Gama

Se \(X \sim \mathrm{Gama}(\alpha,\beta)\), então \[ \left\{\begin{array}{ll} f_\theta^{X}(x) = \frac{\beta^\alpha x^{\alpha-1}\mathrm{e}^{-\beta x}}{\Gamma(\alpha)}, & x > 0\\ 0, & \mathrm{c.c.} \end{array}\right. \] Logo, \[ \begin{aligned} E_\theta(X^k) &= \int^\infty_0 x^k\frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} \mathrm{e}^{-\beta x} dx \\ &= \int^\infty_0 \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha + k - 1} \mathrm{e}^{-\beta x} dx \\ \\ y = \beta x, dy &= \beta dx, x = \frac{y}{\beta} \\ \\ \Rightarrow E_\theta(X^k) &= \frac{\beta^\alpha}{\Gamma(\alpha)} \int^\infty_0 \left(\frac{y}{\beta}\right)^{\alpha + k - 1} \mathrm{e}^{-y} \frac{dy}{\beta} \\ &= \frac{\beta^\alpha}{\Gamma(\alpha)} \frac{1}{\beta^{\alpha+k}} \int^\infty_0 y^{\alpha + k - 1} \mathrm{e}^{-y} dy \\ &= \frac{1}{\Gamma(\alpha) \beta^k} \Gamma(\alpha+k) \\ \Rightarrow E_\theta(X) &= \frac{1}{\Gamma(\alpha) \beta} \Gamma(\alpha +1 ) = \frac{\alpha \Gamma(\alpha)}{\Gamma(\alpha) \beta} = \frac{\alpha}{\beta} \\ \Rightarrow E_\theta(X^2) &= \frac{1}{\Gamma(\alpha) \beta^2} \Gamma(\alpha +2)\\ &= \frac{(\alpha+1) \Gamma(\alpha+1)}{\Gamma(\alpha) \beta} = \frac{(\alpha+1)\alpha\Gamma(\alpha)}{\Gamma(\alpha)\beta^2} = \frac{\alpha(\alpha+1)}{\beta^2} \end{aligned} \]

Seguindo para as equações de momentos, \[ \begin{aligned} &\left\{\begin{array}{ll} \frac{\alpha}{\beta} = \bar{X} \\ \frac{\alpha (\alpha+1)}{\beta^2} = \frac{1}{n} \sum X_i^2 \end{array}\right. \\ \Rightarrow &\left\{\begin{array}{ll} \alpha = \bar{X}\beta \\ \frac{\bar{X}\beta(\bar{X}\beta+1)}{\beta^2} = \frac{1}{n} \sum X_i^2 \end{array}\right. \\ \Rightarrow &\left\{\begin{array}{ll} \alpha = \bar{X}\beta \\ (\bar{X}\beta)^2 + \bar{X}\beta = \beta^2 \frac{1}{n} \sum X_i^2 \end{array}\right. \\ \Rightarrow &\left\{\begin{array}{ll} \alpha = \bar{X}\beta \\ \beta^2\left(\frac{1}{n}\sum X_i^2\right) - (\beta\bar{X})^2 - \beta \bar{X} = 0 \end{array}\right. \\ \Rightarrow &\left\{\begin{array}{ll} \alpha = \bar{X}\beta \\ \beta^2\underbracket{\left(\frac{1}{n}\sum X_i^2 - \bar{X}^2\right)}_{S^2} -\beta\bar{X} = 0 \end{array}\right. \\ \\ \frac{\bar{X} \pm \sqrt{\bar{X}^2 - 4\cdot S^2 \cdot 0}}{2S^2} &= \frac{\bar{X} \pm \sqrt{\bar{X}^2}}{2S^2}\\ &= \left\{\begin{array}{l} 0 \\ \frac{\bar{X}}{S^2} \end{array}\right. \\ \\ \Rightarrow &\left\{\begin{array}{ll} \alpha = \bar{X}\beta \\ \beta^2 S^2 = \beta \bar{X} \end{array}\right. \\ \Rightarrow &\left\{\begin{array}{ll} \hat{\alpha}_{\mathrm{MM}}(\boldsymbol{X}_n) = \frac{\bar{X}^2}{S^2} \\ \hat{\beta}_{\mathrm{MM}}(\boldsymbol{X}_n) = \frac{\bar{X}}{S^2} \end{array}\right. \\ \Rightarrow& \hat{\theta}_{\mathrm{MM}}(\boldsymbol{X}_n) = \left(\frac{\bar{X}^2}{S^2}, \frac{\bar{X}}{S^2}\right). \end{aligned} \]

using Random, Distributions, StatsBase, Plots, LaTeXStrings

Random.seed!(31)

M = 10_000
n = 10 # Tamanho da amostra
alpha0 = 4
beta0 = 6

d = Gamma(alpha0, 1/beta0)

alphaMM = zeros(M)
betaMM = zeros(M)
for i in 1:M
  x = rand(d, n)
  s = var(x, corrected = false)
  alphaMM[i] = mean(x)^2/s
  betaMM[i] = mean(x)/s
end

# Comparar
println("Média simulada do EMM para alpha = $(mean(alphaMM)),
        real = $alpha0")
println("Média simulada do EMM para beta = $(mean(betaMM)),
        real = $beta0")

# Plotar
pmmalpha = histogram(alphaMM, label = "",
                title = "Histograma do EMM para " * L"\alpha", normalize=:pdf,
                     bins = 40)
vline!([alpha0], label="Valor real")
pmmbeta = histogram(betaMM, label = "",
                title = "Histograma do EMM para " * L"\beta", normalize=:pdf,
                    bins = 40)
vline!([beta0], label="Valor real")
display(plot(pmmalpha, pmmbeta))

Média simulada do EMM para alpha = 5.838811157394312,
        real = 4
Média simulada do EMM para beta = 8.994281492594121,
        real = 6