cap4
33 pág.

cap4


DisciplinaInferência Bayesiana27 materiais237 seguidores
Pré-visualização7 páginas
Seja \u3b8 = (\u3b81, . . . , \u3b8k)
\u2032
um vetor aleato´rio de dimensa\u2dco k com func¸a\u2dco de densidade p(\u3b8). Neste caso os
valores gerados sera\u2dco tambe´m vetores \u3b81, . . . ,\u3b8n e o estimador de Monte Carlo
fica
I\u2c6 =
1
n
n\u2211
i=1
g(\u3b8i)
Exemplo 4.2 : Suponha que queremos calcular Pr(X < 1, Y < 1) onde o ve-
tor aleato´rio (X, Y ) tem distribuic¸a\u2dco Normal padra\u2dco bivariada com correlac¸a\u2dco
igual a 0,5. Note que esta probabilidade e´ a integral de p(x, y) definida no inter-
valo acima, portanto simulando valores desta distribuic¸a\u2dco poderemos estimar esta
probabilidade como a proporc¸a\u2dco de pontos que caem neste intervalo. A Figura 4.3
apresenta um diagrama de dispersa\u2dco dos valores simulados e foi obtida usando os
camandos do R abaixo.
\u22123 \u22122 \u22121 0 1 2 3
\u2212
3
\u2212
2
\u2212
1
0
1
2
3
x
y
Figura 4.3: Diagrama de dispersa\u2dco de 1000 valores simulados da distribuic¸a\u2dco N(0,1)
bivariada.
54 CAPI´TULO 4. ME´TODOS APROXIMADOS
\u22123 \u22121 0 1 2 3
0.
0
0.
1
0.
2
0.
3
x
p(x
)
\u22123 \u22121 0 1 2 3
0.
0
0.
1
0.
2
0.
3
y
p(y
)
\u22123 \u22121 0 1 2 3
0.
0
0.
1
0.
2
0.
3
p(x
 | y
<1
)
\u22123 \u22121 0 1 2 3
0.
0
0.
1
0.
2
0.
3
0.
4
p(y
 | x
<1
)
Figura 4.4: Estimativas das densidades marginais e condicionais no Exemplo 4.2.
Uma grande vantagem dos me´todos de simulac¸a\u2dco e´ que apo´s uma amostra
de vetores aleato´rios ser gerada podemos facilmente calcular caracter´\u131sticas das
distribuic¸o\u2dces marginais e condicionais. No Exemplo 4.2, para calcular Pr(X < 1)
basta calcular a freque\u2c6ncia relativa de pontos (xi, yi) tais que xi < 1. Para
calcular a probabilidade condicional Pr(X < 1|Y < 1) basta selecionar somente
aqueles pontos cuja segunda coordenada e´ menor do que 1. Depois calcula-se a
freque\u2c6ncia relativa dos pontos restantes cuja primeira coordenada e´ menor do que
1.
4.4.1 Monte Carlo via Func¸a\u2dco de Importa\u2c6ncia
Em muitas situac¸o\u2dces pode ser muito custoso ou mesmo imposs´\u131vel simular valores
da distribuic¸a\u2dco a posteriori. Neste caso, pode-se recorrer a` uma func¸a\u2dco q(\u3b8) que
seja de fa´cil amostragem, usualmente chamada de func¸a\u2dco de importa\u2c6ncia. O
procedimento e´ comumente chamado de amostragem por importa\u2c6ncia.
Se q(\u3b8) for uma func¸a\u2dco de densidade definida no mesmo espac¸o de variac¸a\u2dco
de \u3b8 enta\u2dco a integral (4.1) pode ser reescrita como
I =
\u222b
g(\u3b8)p(\u3b8)
q(\u3b8)
q(\u3b8)dx = E
[
g(\u3b8)p(\u3b8)
q(\u3b8)
]
4.4. ME´TODO DE MONTE CARLO SIMPLES 55
onde a esperanc¸a agora e´ com respeito a distribuic¸a\u2dco q. Assim, se dispomos de
uma amostra aleato´ria \u3b81, . . . , \u3b8n tomada da distribuic¸a\u2dco q o estimador de Monte
Carlo da integral acima fica
I\u2c6 =
1
n
n\u2211
i=1
g(\u3b8i)p(\u3b8i)
q(\u3b8i)
.
e tem as mesmas propriedades do estimador de Monte Carlo simples.
Em princ´\u131pio na\u2dco ha´ restric¸o\u2dces quanto a` escolha da densidade de importa\u2c6ncia
q, pore´m na pra´tica alguns cuidados devem ser tomados. Pode-se mostrar que
a escolha o´tima no sentido de minimizar a varia\u2c6ncia do estimador consiste em
tomar q(\u3b8) \u221d g(\u3b8)p(\u3b8).
Exemplo 4.3 : Para uma u´nica observac¸a\u2dco X suponha que
X|\u3b8 \u223c N(\u3b8, 1) e \u3b8 \u223c Cauchy(0, 1).
Enta\u2dco,
p(x|\u3b8) \u221d exp[\u2212(x\u2212 \u3b8)2/2] e p(\u3b8) = 1
pi(1 + \u3b82)
.
Portanto, a densidade a posteriori de \u3b8 e´ dada por
p(\u3b8|x) =
1
1 + \u3b82
exp[\u2212(x\u2212 \u3b8)2/2]\u222b
1
1 + \u3b82
exp[\u2212(x\u2212 \u3b8)2/2]d\u3b8
.
Suponha agora que queremos estimar \u3b8 usando func¸a\u2dco de perda quadra´tica. Como
vimos no Cap´\u131tulo 3 isto implica em tomar a me´dia a posteriori de \u3b8 como esti-
mativa. Mas
E[\u3b8|x] =
\u222b
\u3b8p(\u3b8|x)d\u3b8 =
\u222b
\u3b8
1 + \u3b82
exp[\u2212(x\u2212 \u3b8)2/2]d\u3b8\u222b
1
1 + \u3b82
exp[\u2212(x\u2212 \u3b8)2/2]d\u3b8
e as integrais no numerador e denominador na\u2dco te\u2c6m soluc¸a\u2dco anal´\u131tica exata.
Uma soluc¸a\u2dco aproximada via simulac¸a\u2dco de Monte Carlo pode ser obtida usando
o seguinte algoritmo,
1. gerar \u3b81, . . . , \u3b8n independentes da distribuic¸a\u2dco N(x, 1);
2. calcular gi =
\u3b8i
1 + \u3b82i
e g\u2217i =
1
1 + \u3b82i
;
56 CAPI´TULO 4. ME´TODOS APROXIMADOS
3. calcular E\u2c6(\u3b8|x) =
\u2211n
i=1 gi\u2211n
i=1 g
\u2217
i
.
Este exemplo ilustrou um problema que geralmente ocorre em aplicac¸o\u2dces
Bayesianas. Como a posteriori so´ e´ conhecida a menos de uma constante de
proporcionalidade as esperanc¸as a posteriori sa\u2dco na verdade uma raza\u2dco de inte-
grais. Neste caso, a aproximac¸a\u2dco e´ baseada na raza\u2dco dos dois estimadores de
Monte Carlo para o numerador e denominador.
Exercicios
1. Para cada uma das distribuic¸o\u2dces N(0, 1), Gama(2,5) e Beta(2,5) gere 100,
1000 e 5000 valores independentes. Fac¸a um gra´fico com o histograma e
a func¸a\u2dco de densidade superimposta. Estime a me´dia e a varia\u2c6ncia da
distribuic¸a\u2dco. Estime a varia\u2c6ncia do estimador da me´dia.
2. Para uma u´nica observac¸a\u2dco X com distribuic¸a\u2dco N(\u3b8, 1), \u3b8 desconhecido,
queremos fazer infere\u2c6ncia sobre \u3b8 usando uma priori Cauchy(0,1). Gere um
valor de X para \u3b8 = 2, i.e. x \u223c N(2, 1).
(a) Estime \u3b8 atrave´s da sua me´dia a posteriori usando o algoritmo do
Exemplo 4.3.
(b) Estime a varia\u2c6ncia da posteriori.
(c) Generalize o algoritmo para k observac¸o\u2dces X1, . . . , Xk da distribuic¸a\u2dco
N(\u3b8, 1).
4.5 Me´todos de Reamostragem
Existem distribuic¸o\u2dces para as quais e´ muito dif´\u131cil ou mesmo imposs´\u131vel simular
valores. A ide´ia dos me´todos de reamostragem e´ gerar valores em duas etapas.
Na primeira etapa gera-se valores de uma distribuic¸a\u2dco auxiliar conhecida. Na
segunda etapa utiliza-se um mecanismo de correc¸a\u2dco para que os valores sejam
representativos (ao menos aproximadamente) da distribuic¸a\u2dco a posteriori. Na
pra´tica costuma-se tomar a priori como distribuic¸a\u2dco auxiliar conforme proposto
em Smith & Gelfand (1992).
4.5.1 Me´todo de Rejeic¸a\u2dco
Considere uma func¸a\u2dco de densidade auxiliar q(\u3b8) da qual sabemos gerar valores.
A u´nica restric¸a\u2dco e´ que exista uma constante A finita tal que p(\u3b8|x) < Aq(\u3b8).
O me´todo de rejeic¸a\u2dco consiste em gerar um valor \u3b8\u2217 da distribuic¸a\u2dco auxiliar q
e aceitar este valor como sendo da distribuic¸a\u2dco a posteriori com probabilidade
4.5. ME´TODOS DE REAMOSTRAGEM 57
p(\u3b8\u2217|x)/Aq(\u3b8\u2217). Caso contra´rio, \u3b8\u2217 na\u2dco e´ aceito como um valor gerado da pos-
teriori e o processo e´ repetido ate´ que um valor seja aceito. O me´todo tambe´m
funciona se ao inve´s da posteriori, que em geral e´ desconhecida, usarmos a sua
versa\u2dco na\u2dco normalizada, i.e p(x|\u3b8)p(\u3b8).
Podemos enta\u2dco usar o seguinte algoritmo,
1. gerar um valor \u3b8\u2217 da distribuic¸a\u2dco auxiliar;
2. gerar u \u223c U(0, 1);
3. se u < p(\u3b8\u2217|x)/Aq(\u3b8\u2217) fac¸a \u3b8(j) = \u3b8\u2217, fac¸a j = j + 1 e retorne ao passo 1.
caso contra´rio retorne ao passo 1.
Tomando a priori p(\u3b8) como densidade auxiliar a constante A deve ser tal que
p(x|\u3b8) < A. Esta desigualdade e´ satisfeita se tomarmos A como sendo o valor
ma´ximo da func¸a\u2dco de verossimilhanc¸a, i.e. A = p(x|\u3b8\u2c6) onde \u3b8\u2c6 e´ o estimador
de ma´xima verossimilhanc¸a de \u3b8. Neste caso, a probabilidade de aceitac¸a\u2dco se
simplifica para p(x|\u3b8)/p(x|\u3b8\u2c6).
Podemos enta\u2dco usar o seguinte algoritmo para gerar valores da posteriori,
1. gerar um valor \u3b8\u2217 da distribuic¸a\u2dco a priori;
2. gerar u \u223c U(0, 1);
3. aceitar \u3b8\u2217 como um valor da posteriori se u < p(x|\u3b8\u2217)/p(x|\u3b8\u2c6), caso contra´rio
rejeitar \u3b8\u2217 e retornar ao passo 1.
Exemplo 4.4 : Suponha que X1, . . . , Xn \u223c N(\u3b8, 1) e assume-se uma distribuic¸a\u2dco
a priori Cauchy(0,1) para \u3b8. A func¸a\u2dco de verossimilhanc¸a e´,
p(x|\u3b8) =
n\u220f
i=1
(2pi)\u22121/2 exp
{
\u2212(xi \u2212 \u3b8)
2
2
}
\u221d exp
{
\u22121
2
n\u2211
i=1
(xi \u2212 \u3b8)2
}
\u221d exp
{
\u2212n
2
(x\u2212 \u3b8)2
}
e o estimador de ma´xima verossimilhanc¸a e´ \u3b8\u2c6 = x¯. Usando o algoritmo acima,
gera-se um valor da distribuic¸a\u2dco Cauchy(0,1) e a probabilidade de aceitac¸a\u2dco neste
caso fica simplesmente exp[\u2212n(x¯ \u2212 \u3b8)2/2]. A func¸a\u2dco do R a seguir obte´m uma
amostra de tamanho m de \u3b8 e como ilustrac¸a\u2dco vamos gerar 50 observac¸o\u2dces