Método dos momentos - Distribuição Uniforme(a, b)

O método dos momentos consiste em solucionar um sistema de equações obtidas ao igualarmos
os momentos populacionais aos momentos amostrais. O primeiro momento momento amostral é dado pela própria média amostral, isto é,
$$m_1=\sum_{i=1}^{n} \frac{X_i}{n}$$
De maneira geral, para $k=1, 2, \dots$, temos o $k$-ésimo momento amostral:
$$m_k=\sum_{i=1}^{n} \frac{X_i^k}{n}$$
Assim, o sistema é obtido ao fazermos:
$$m_k=E(X^k)$$
e resolvendo para encontrar as estimativas dos parâmetros da distribuição desejada.

Exemplo:
Para a distribuição Uniforme(a,b), temos dois parâmetros: $a$ e $b$. Serão necessárias duas equações para encontrarmos $\tilde{a}$ e $\tilde{b}$, os estimadores obtidos por meio do método dos momentos. As equações que utilizaremos são:
$$m_1=E(X) \Rightarrow \sum_{i=1}^{n} \frac{X_i}{n} = \frac{\tilde{a}+\tilde{b}}{2}$$
$$m_2=E(X^2) \Rightarrow \sum_{i=1}^{n} \frac{X_i^2}{n} = \frac{\tilde{a}^2+\tilde{a}\tilde{b}+\tilde{b}^2}{3}$$
Resolvendo a equação $m_1=E(X)$, obtemos $\tilde{a}=2 \bar{X} - \tilde{b}$, em que $\bar{X}$ é a média amostral.
Com a segunda equação, $m_2=E(X^2)$, obtemos $\tilde{b}=\bar{X} + \sqrt{3(m_2-\bar{X}^2)}$.
Assim, temos que os estimadores pelo método dos momentos para a distribuição Uniforme(a,b) são dados por:
$$\tilde{a}=2 \bar{X} - \tilde{b}= \bar{X} -  \sqrt{3(m_2-\bar{X}^2)}$$
$$\tilde{b}=\bar{X} + \sqrt{3(m_2-\bar{X}^2)}$$
O gráfico a seguir mostra o resultado para uma amostra proveniente de uma distribuição Uniforme(5,9), de tamanho $n=1000$. Nas linhas vermelhas estão os valores verdadeiros dos parâmetros ($a=5$ e $b=9$) e nas linhas verdes estão as estimativas obtidas pelos estimadores discutidos nesta postagem ($\tilde{a}=4,97$ e $b=8,90$).

Uma reflexão importante sobre a distribuição uniforme é que seu estimador obtido pelo método dos momentos não é função de estatística suficiente.

Por exemplo, no caso da distribuição Uniforme($0, \theta$) o estimador obtido pelo método dos momentos é igual a $\tilde{\theta}=2\bar{X}$, que não é função da estatística suficiente desta distribuição, que é o máximo, $X_{(n)}$.

Questão de Concurso:
Perceba como essa reflexão foi cobrada no último concurso do Supremo Tribunal Militar (STM), na prova para o cargo de Estatístico, elaborada pela banca CESPE, em 2011:

Item: "Qualquer estimador obtido pelo método dos momentos é uma função de estatística suficiente."
Resposta: O exemplo acima mostra um contraexemplo. Item errado.

Método de aceitação rejeição - Com exemplo feito em R

O método de aceitação rejeição (ou método de rejeição) é usado para simular valores de uma variável aleatória, geralmente quando o método da inversa não pode ser feito. O método de rejeição consiste no seguinte. Considere uma v.a. com densidade $f(x)$ da qual é necessário gerar valores e $g(x)$ a densidade de outra v.a., da qual é possível gerar valores. No método, é gerado um valor $W$ da v.a. com densidade $g(x)$ e ele é aceito com probabilidade proporcional a $f(W)/g(W)$, com
$$\frac{f(w)}{g(w)} \leq c \text{ para todo $w$.}$$
O algoritmo pode ser escrito como:

1. Gere $W$ com densidade g;
2. Gere um número $U \sim \mathcal{U}(0,1)$;
3. Se $U \leq \frac{f(W)}{c g(W)}$, aceite o valor W. Caso contrário, retorne ao passo 1.

Exemplo no R. 

Suponha que queiramos gerar valores de uma v.a. com densidade

$$f(x) = 12 x^2 (1-x), \ 0 < x < 1.$$

Note que $X\sim B(3,2)$. Como esta variável pertence ao intervalo (0,1), podemos considerar o método com $g(x)=1$, para $x \in (0,1)$, i.e, uma distribuição uniforme. Para determinar o menor valor de $c$ tal que $f(x)/g(x) \leq c$, maximizamos 

$$\frac{f(x)}{g(x)} = 12 x^2 (1-x),$$

obtendo $x=2/3$. Assim, $f(x)/g(x) \leq 12 (2/3)^2 (1/3)=16/9=c$. Assim, temos que 

$$\frac{f(x)}{c \ g(x)} = \frac{27}{4} x^2 (1-x).$$

Assim, para este exemplo o algoritmo fica:

1. Gere $u_1$ e $u_2$ com distribuição $\mathcal{U}(0,1)$;
2. Se $u_2 \leq \frac{27}{4} u^2_1 (1-u_1)$, aceite o valor $u_1$. Caso contrário, retorne ao passo 1.

O código em R, adaptado de uma postagem do blog "Playing with R" (para ler, clique aqui),
pode ser escrito como:

x = runif(10000,0,1) #gerando os valores que podem ser aceitos ou não

accept = c()

for(i in 1:length(x)){

  U = runif(1)

  if(U <= (27/4)*(x[i]^2)*(1-x[i])) {

    accept[i] = 'Yes'

  }

  else if(U > (27/4)*(x[i]^2)*(1-x[i])) {

    accept[i] = 'No'

  }

}  

T = data.frame(x, accept = factor(accept, levels= c('Yes','No')))

xseq=seq(0,1,0.01)

hist(T[,1][T$accept=='Yes'], breaks = seq(0,1,0.01), ylim=c(0,2), freq = FALSE, main = 'Histogram of X', xlab = 'X')

lines(xseq, dbeta(xseq,3,2))  

library(ggplot2)

print(qplot(x, data = T, geom = 'density', color = accept))

print(qplot(x, data = T, geom = 'histogram', fill = accept, binwidth=0.01))

Ordens de magnitude estocástica

Considere $(X_n)_{n \geq 1}$ e $(Y_n)_{n \geq 1}$ sequências de variáveis aleatórias definidas em $(\Omega, \mathcal{A}, \mathcal{P})$. Supomos que $Y_n \neq 0$ q.c. $[\mathcal{P}]$, i.e, $P(Y_n \neq 0)=1$, $n=1, 2, \dots$.

Definição: $X_n = O_p(Y_n)$ se, e somente se, para todo $\eta>0$ existirem um número $k=k(\eta)>0$ e $n_0=n_0(\epsilon)$ natural tal que, para todo $n \geq n_0$,

$$P \left( \left| \frac{X_n}{Y_n} \right| \geq k  \right) \leq \eta$$

Em outras palavras, $\left(\frac{X_n}{Y_n}\right)_{n \geq 1}$ é limitada em probabilidade para $n$ suficientemente grande.

Definição: $X_n = o_p(Y_n)$ se, e somente se, para quaisquer números reais $\epsilon$ e  $\eta$ positivos existir um número $n_0=n_0(\epsilon, \eta)$ natural tal que:

$$P \left( \left|\frac{X_n}{Y_n}\right| \geq \epsilon  \right) \leq \eta, \text{ para todo $n \geq n_0$.}$$

Ou, de maneira equivalente, se

$$\lim_{n \rightarrow \infty} P \left( \left|\frac{X_n}{Y_n}\right| \geq \epsilon  \right) =0.$$  

Exemplo: Seja $(X_n)_{n \geq 1}$ uma sequência de variáveis aleatórias definidas no espaço
$(\Omega, \mathcal{A}, \mathcal{P})$. Prove que $O_p(O_p(X_n))=O_p (X_n)$.

Solução:
Escreva $Y_n=O_p(X_n)$ e $W_n=O_p(Y_n)$. Queremos então mostrar que $W_n=O_p(X_n)$.
Como $Y_n=O_p(X_n)$ e $W_n=O_p(Y_n)$, então para todo $\eta>0$ existem $k_1$ e $k_2$, e também $n_0=n_0(k_1)$ e $n_1=n_1(k_2)$ tais que

$$P\left( \left|\frac{Y_n}{X_n}\right| \geq k_1  \right) \leq \eta/2, \text{ para todo $n \geq n_0$, e}$$
$$P\left( \left|\frac{W_n}{Y_n}\right| \geq k_2  \right) \leq \eta/2, \text{ para todo $n \geq n_1$.}$$

Temos assim que, para todo $n \geq$ max$(n_0, n_1)$,

$$\begin{eqnarray}P\left( \left|\frac{W_n}{X_n}\right| \geq k_1 k_2  \right) &=& P\left( \left|\frac{W_n Y_n}{Y_n X_n}\right| \geq k_1 k_2  \right) \nonumber\\ &\leq& P\left( \left|\frac{Y_n}{X_n}\right| \geq k_1  \right) + P\left( \left|\frac{W_n}{Y_n}\right| \geq k_2  \right) \\ &\leq& \eta/2 +\eta/2 = \eta \nonumber \end{eqnarray}$$

Assim, $W_n=O_p(X_n)$ ou $O_p(O_p(X_n))=O_p (X_n)$.

Note que a equação (1) ocorre pois para $X$ e $Y$ variáveis aleatórias e $a$ e $b$ números reais positivos, temos que

  $$\begin{eqnarray}P(XY>cd) &=&P(XY>cd, X>c, Y>d)+ P(XY>cd, X>c, Y<d) \nonumber \\  &+& P(XY>cd, X<c, Y>d) + P(XY>cd, X<c, Y<d) \nonumber\\  &=&P(\{ X>c \text{ ou } Y>d \} \text{ e } XY>cd) \nonumber \\  &\leq& P(X>c \text{ ou } Y>d)\\  &\leq& P(X>c) + P(Y>d) \end{eqnarray}$$

Note que a equação (2) ocorre porque a probabilidade da interseção é menor ou igual que a probabilidade de um dos eventos. Já a equação (3) resulta do fato de que a probabilidade da união é menor ou igual que a soma das probabilidades.

Ordens de magnitude de sequências de números reais

No estudo de teoria assintótica, muitas vezes é interessante comparar sequências de variáveis aleatórias quando $n$, o tamanho da amostra, tende a infinito. Antes de conhecer a notação $O_p$ e $o_p$ usada para comparar sequências de variáveis aleatórias, será apresentada a notação $O$ e $o$, usada para comparar sequências de números reais.

Considere duas sequências de números reais, $(a_n)_{n \geq 1}$  e $(b_n)_{n \geq 1}$, com $(b_n) \neq 0$ para todo $n$ suficientemente grande.

Definição: $a_n=O(b_n)$ quando $n \rightarrow \infty$ se houver um número real $k>0$ e $n_0=n_0 (k)$ natural tal que

$$\left| \frac{a_n}{b_n} \right| \leq k, \text{ para todo $n \geq n_0$.}$$

Em outras palavras, dizemos que a ordem de magnitude de $(a_n)_{n \geq 1}$  é no máximo igual à ordem de magnitude de $(b_n)_{n \geq 1}$  para todo $n$ suficientemente grande.

Note também que $a_n=O(b_n)$ se e somente se $\left(\frac{a_n}{b_n}\right)_{n \geq 1}$ é uma sequência limitada, e $a_n=O(1)$ se e somente se $(a_n)_{n \geq 1}$ é uma sequência limitada.

Definição: $a_n=o(b_n)$ quando $n \rightarrow \infty$ se, e somente se, para todo $\epsilon>0$ existir $n_0=n_0 (\epsilon)$ natural tal que

$$\left| \frac{a_n}{b_n} \right| \leq \epsilon, \text{ para todo $n \geq n_0$.}$$

Em outras palavras, dizemos que a ordem de magnitude de $(a_n)_{n \geq 1}$  é menor do que a ordem de magnitude de $(b_n)_{n \geq 1}$  para todo $n$ suficientemente grande.

Note também que $a_n=o(b_n)$ se e somente se $\lim_{n \rightarrow \infty} \frac{a_n}{b_n}=0$, e $a_n=o(1)$ se e somente se $\lim_{n \rightarrow \infty} a_n=0$.

Melhor estimador não viciado (ENVVUM) de $\theta$ no caso da Distribuição Uniforme($0, \theta$)

Comecemos com uma amostra aleatória de tamanho $n$, $X_1, X_2, \dots, X_n$, de uma população com distribuição Uniforme(0, $\theta$).

Objetivo: Queremos encontrar um estimador não viciado de variância uniformemente mínima (ENVVUM) para $\theta$.

Como faremos isso: Há várias formas de se obter um ENVVUM. Aqui iremos:

1. Buscar uma estatística suficiente e completa;
2. Encontrar um estimador não viciado que seja função da estatística obtida em (1) e 
3. Através do Teorema de Lehmann-Sheffé, obteremos o  ENVVUM.

A densidade para uma observação é dada por: 

$$f(\mathbf{x}|\theta)=\prod_{i=1}^{n} \frac{1}{\theta} I_{(0,\theta)} {(x_i)} =I_{(0,x_n)} {(x_1)} \frac{1}{\theta^n} I_{(0,\theta)} {(x_n)}$$

Como separamos a densidade $f(\mathbf{x}|\theta)$ em uma parte que depende apenas da amostra e outra que depende do parâmetro através de $X_n$, temos que $T=X_n$ é uma estatística suficiente, pelo Critério de Fatoração de Fisher-Neyman.

Para verificar a completude da estatística $T=X_n$, note que se uma estatística é completa, temos que $E(g(T))=0  \iff  g(T)=0$. Lembre de uma postagem anterior que calculamos a densidade de $T=X_n$, obtendo $f_T (x)=\frac{n}{\theta^n} x^{n-1}$. Assim,

$$E(g(T))=\int_{0}^{\theta} g(T) \frac{n}{\theta^n} x^{n-1} dx = 0 \iff g(T)=0$$ ,

Usando o teorema fundamental do cálculo, temos que $g(T)=0$.

Mostramos então que $X_n$ é uma estatística suficiente e completa. Calculamos anteriormente a esperança de $T=X_n$, que é $\frac{n \theta}{(n+1)}$. Assim, considere
$T^*=\left(\frac{n+1}{n}\right)X_n$. Temos naturalmente que $E(T^*)=\theta$. Então $T^*$ é um estimador não viciado para $\theta$ que é função de $X_n$, uma estatística suficiente e completa. Logo, pelo Teorema de Lehmann-Scheffé, $T^*$ é ENVVUM.

Quociente de duas distribuições gama independentes

O quociente de duas distribuições gama independentes é realizado nesta postagem. No primeiro quadro, obtém-se a densidade para $V=X_1/X_2$, em que $X_i \sim \text{ Gama }(\alpha_i, \beta_i, \gamma_i), i=1, 2$. No segundo quadro, a distribução acumulada é encontrada. Tando a função densidade de probabilidade como a função de distribuição acumulada são escritas em termos da função H, cuja definição pode ser encontrada em ambos os quadros. 

Distribuição do produto de duas variáveis exponenciais

Nesta postagem será apresentada a obtenção da função densidade de probabilidade (f.d.p.) da variável $Z=XY$, em que $X$ e $Y$ seguem distribuição exponencial com parâmetros $a$ e $b$, respectivamente, e as duas variáveis são independentes. Para isso, usaremos a transformada de Mellin (para saber mais, clique aqui).

A transformada de Mellin

Quando $f(z)$ se trata de uma f.d.p., a transformada de Mellin de $f$  com argumento $s \in \mathbb{C}$ é dada por:

$$F(s)=E(Z^{s-1})=\int_0^{\infty} z^{s-1} f(z) dz,$$

em que $\alpha \leq s \leq \beta$.

A transformada inversa de Mellin

A partir da expressão para o $r-$ésimo momento de uma distribuição, pode-se conhecer a densidade através da transformada inversa de Mellin, dada por:

$$f(z)=\frac{1}{2 \pi i} \int_L z^{-s} E(Z^{s-1}) ds,$$

em que $i=\sqrt{-1}$.

Distribuição do produto de duas v.a. exponenciais independentes

Problema: $X \sim Exp(a)$ e $Y \sim Exp(b)$. Encontre a f.d.p de $Z=XY$ assumindo que $X$ e $Y$ são independentes.

Demonstração: Através da suposição de independência e da fórmula para o $r$-ésimo momento de uma v.a. exponencial, temos que:

$$E(Z^{s-1})=E(X^{s-1}Y^{s-1})=E(X^{s-1})E(Y^{s-1})=\frac{\Gamma(s)}{a^{s-1}} \frac{\Gamma(s)}{b^{s-1}}$$

Assim, usando a transformada inversa de Mellin, temos que

$$f(z)=\frac{1}{2 \pi i} \int_L z^{-s}E(z^{s-1}) ds = \frac{ab}{2 \pi i} \int_L (z a b)^{-s} \Gamma^2(s) ds$$

Pelo teorema de Cauchy,

$$f(z)=\frac{ab}{2 \pi i} 2 \pi i \sum_{r=1}^{n} Res,$$
em que $Res$ significa resíduos.

Para resolver essa integral complexa, precisamos achar os pólos (i.e., os pontos de singularidade). Para este caso, como tem-se o quadrado de uma função gama, os pólos são de segunda ordem e iguais  a $s = -r$, onde $r=0,1,2,\dots$. Assim, temos que os resíduos são dados por

$$Res = \lim_{s \rightarrow -r} \frac{\partial^{2-1}}{\partial s^{2-1}} \left\{ (s+r)^2 (z a b)^{-s}  \Gamma^2(s) \right\}$$

Note que
$$(s+r) \Gamma(s)=\frac{\Gamma(s+r+1)}{\prod_{j=0}^{r-1} (s+j)} = \frac{1}{r!}$$

Seja $A(s)=(s+r)^2 (z a b)^{-s}  \Gamma^2(s)$. $A(s)$, usando o fato anterior, pode ser escrito como

$$A(s)= \frac{\Gamma^2(s+r+1)}{\left(\prod_{j=0}^{r-1} (s+j) \right)^2} (z a b)^{-s}$$

$$\lim_{s \rightarrow -r}  \log(A(s)) = 2 \log(\Gamma(1)+r \log(z a b) - 2 \log \left(\prod_{j=0}^{r-1} (s+j) \right)$$

Usando o fato de que $\log(A(s))=(1/A(s)) A'(s)$,

$$\log(A(s)) = \left( 2 \Psi(1)-\log(z a b) -2\sum_{j=0}^{r-1} \frac{1}{j-r} \right) \frac{(zab)^r}{(r!)^2} = A'(s)(1/A(s)),$$

Em que $\Psi(u)=\frac{d}{d u} \log(\Gamma(u))=\frac{\Gamma(u)}{u}$.
Como estamos interessados em $Res=A'(s)= \left( 2 \Psi(1)-\log(z a b) -2\sum_{j=0}^{r-1} \frac{1}{j-r} \right)$ voltamos ao cálculo da densidade de $Z$ com

$$f(z)=\frac{ab}{2 \pi i} 2 \pi i \sum_{r=1}^{n} Res = \frac{ab} \sum_{r=1}^{n}  \left( 2 \Psi(1)-\log(z a b) -2\sum_{j=0}^{r-1} \frac{1}{j-r} \right),$$
Com $z \in (0, \infty)$.