terça-feira, 12 de dezembro de 2017

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.

sábado, 26 de setembro de 2015

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))



segunda-feira, 21 de setembro de 2015

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$.

sexta-feira, 29 de maio de 2015

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.

sábado, 25 de outubro de 2014

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. 









terça-feira, 16 de setembro de 2014

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)$.