Capítulo XI. Monte Carlo Cadeias de Markov (MCMC)



1. Introdução


Lindley (1983) explicou o paradigma Bayesiano da seguinte forma. O interesse aqui é uma quantidade \(\theta\), que é desconhecida. Mas podemos ter alguma ideia pessoal sobre sua distribuição que deve expressar nossa opinião relativa quanto à verossimilhança de vários valores possíveis de \(\theta\) serem o valor verdadeiro. Esta será a distribuição a priori de \(\theta\) e será denotada \(\pi(\theta)\).

Isso representa o estado de nosso conhecimento, de alguma forma, antes de realizar um experimento ou observar os dados. Além disso, existe uma distribuição de probabilidade \(f(x|\theta)\), que descreve a verossimilhança relativa dos valores \(x\), dado que \(\theta\) é o verdadeiro valor do parâmetro, para uma variável aleatória \(X\). Mas não há nada de novo aqui. Esta é a distribuição paramétrica de uma variável aleatória \(X\), que algumas vezes foi denotada por \(f_\theta(x)\). A principal diferença pode ser que \(\theta\) é uma variável aleatória aqui e, portanto, \(f_\theta(x)\) agora é uma distribuição condicional, de \(X\), dada a informação que \(\theta\) é o valor verdadeiro.

Com base nessas duas funções, usamos o Teorema de Bayes para calcular \[ \pi^*(\theta |x) = \dfrac{f(x|\theta)\pi(\theta)}{\displaystyle \int f(x|\theta)\pi(\theta)\mbox{d}\theta}, \] que é a distribuição a posteriori de \(\theta\).

Pode ser visto como uma opinião revisada, uma vez que os resultados do experimento são conhecidos. Uma vez que tenhamos essa distribuição a posteriori, que contém todo o conhecimento que temos, nossa subjetiva a priori e a amostra, que pode ser vista como mais objetiva; podemos calcular qualquer quantidade de interesse. Isso pode ser feito tanto usando fórmulas analíticas, quando as últimas são tratáveis, quanto usando simulações de Monte Carlo para gerar valores possíveis de \(\theta\), dada a amostra observada.

Considere, por exemplo, o caso em que\(\theta\) é uma medida de posição que pode ser a moda, a mediana, a média, etc. Na configuração clássica e padrão, geralmente calculamos intervalos de confiança, para que possamos afirmam que com probabilidade \(1-\alpha\), \(\theta\) pertence a algum intervalo. Na estrutura Bayesiana, calculamos um intervalo de credibilidade de \(1-\alpha\), de modo que \[ P\big(\theta\in [\ell,u] \, | \, X\big) = \int_\ell^u \pi^*(\theta \, | \, X)\mbox{d}\theta = 1-\alpha\cdot \]

Também podemos usar uma versão Bayesiana do Teorema do Limite Central para derivar um intervalo de credibilidade: Sob condições adequadas, a distribuição a posteriori pode ser aproximada com uma distribuição Gaussiana e um intervalo de credibilidade aproximado é então \[ \mbox{E}\big( \theta \, | \, X\big ) \pm \Phi^{-1}(1-\alpha/2)\sqrt{\mbox{Var}\big( \theta \, | \, X\big )}\cdot \]


1.1. Dois tipos de probabilidade


Existem dois tipos de probabilidade: objetiva e subjetiva. Grosso modo, as probabilidades objetivas são definidas pela realidade por meio de leis científicas e processos físicos. As probabilidades subjetivas são criadas pelas pessoas para ajudá-las a quantificar suas crenças e analisar as consequências de suas crenças.

Por exemplo, assume as sentenças da Tabela 1. As frases da primeira linha são objetivas. Se a moeda é justa ou não (50% de probabilidade de cara) ou viciada depende não do que alguém acredita sobre a moeda, mas das propriedades físicas da moeda (sua forma e distribuição de massa).

\[ \begin{array}{l|c|c} & \mbox{Exemplo comum} & \mbox{Exemplo de seguro} \\\hline \mbox{Objetivo} & \mbox{A moeda tem 25% de probabilidade} & \mbox{Um segurado tem 1% de chance} \\ & \mbox{de sair cara duas vezes seguidas.} & \mbox{de perda de propriedade.} \\ \hline \mbox{Subjetivo} & \mbox{Há 90% de chance de que exista vida} & \mbox{Há uma chance de 50% de que os custos} \\ & \mbox{em outros planetas.} & \mbox{ de responsabilidade nacional } \\ & & \mbox{por exposição diminuam no próximo ano.} \end{array} \]

TABELA 1. Probabilidades objetivas e subjetivas.

Por outro lado, as frases da segunda linha parecem menos objetivas. A vida existe em outros planetas ou não; não está claro o que tornaria objetivamente verdade que a probabilidade disso é de 90%. Uma probabilidade como essa varia entre as pessoas e reflete a força da crença.

Embora o limite exato entre probabilidade subjetiva e objetiva seja controverso, a distinção básica em si é um tanto óbvia. É importante para este artigo porque essa distinção é o que distingue a estatística Bayesiana da clássica. Bayesianos acham útil lidar com probabilidades subjetivas. Por outro lado, os estatísticos clássicos pensam que apenas as probabilidades objetivas valem a pena, seja porque as probabilidades subjetivas não fazem sentido ou porque são muito subjetivas.


1.2. Trabalhando com probabilidades subjetivas na vida real


Como probabilidades subjetivas que são apenas sobre crenças teriam algum valor prático? Considere o seguinte argumento:
1. Há 50% de chance de que os custos de responsabilidade por exposição diminuam.
2. Independentemente dos custos de responsabilidade, há 50% de chance de que nosso CEO decida reduzir despesas.
3. Nossa empresa só será lucrativa se, e somente se, os custos de responsabilidade diminuirem e o CEO reduzir as despesas.
4. Portanto, há menos de 25% de chance de fazermos um ajuste profissional.

Embora as afirmações (1) e (2) possam ser probabilidades puramente subjetivas, um argumento como o acima ainda pode ser poderoso. Mesmo que as duas primeiras afirmações possam não ser objetivamente verdadeiras ou falsas, alguém que endossa (1)-(3) e ainda assim rejeita (4) parece estar errado.

As estatísticas Bayesianas realmente funcionam quando você adiciona o princípio da condicionalização:

Se suas crenças estão associadas à função de probabilidade \(P\) e você aprende a evidência \(e\), então suas novas crenças devem agora ser associadas à função de probabilidade \(P_e\), onde para todo \(x\), \[ P_e(x)=P(x \, | \, e)\cdot \] Por definição ou axioma, \[ P(x \, | \, e) =\dfrac{P(x \wedge e)}{P(e)}\cdot \]

Chamamos \(P(x)\) de probabilidade a priori de \(x\) porque é a probabilidade de \(x\) antes de conhecer a evidência \(e\). Da mesma forma, \(P_e(x) = P(x \, | \, e)\) é chamada de probabilidade a posteriori de \(x\).

Com o princípio da condicionalização, a estatística Bayesiana fornece uma sugestão simples, porém poderosa, sobre a dinâmica das crenças: como as crenças de uma pessoa devem mudar quando ela aprende coisas novas. Por exemplo, suponha que adicionamos às sentenças (1)-(4) acima:
5. O CEO decidiu reduzir as despesas.
6. A nova probabilidade de obtermos um lucro é de 50%.
Quando (5) é aprendido, então (6) descreve a nova probabilidade de um lucro através da regra de condicionalização.


2. Conjugados Bayesianos


Como veremos em breve, embora o Bayesianismo tenha uma declaração simples, sua implementação geralmente requer cálculos complexos. No entanto, essa complexidade muitas vezes pode ser evitada usando conjugados Bayesianos, também chamados de distribuições conjugadas.


2.1. Perspectiva histórica


Pierre-Simon Laplace observou, de 1745 a 1784, 770.941 nascimentos: 393.386 meninos e 377.555 meninas. Assumindo que o sexo de um bebê é uma variável aleatória Bernoulli, de modo que \(p\) denota a probabilidade de ter um menino, Pierre-Simon Laplace queria quantificar a credibilidade da hipótese \(p\geq 0.5\). Sua idéia era assumir que a priori, \(p\) poderia ser uniformemente distribuído no intervalo unitário e então, dadas as estatísticas observadas, calcular a probabilidade a posteriori de que \(p\geq 0.5\).

Se \(N\) é o número de meninos de \(n_0\) nascimentos, então \[ f(p \, | \, N=n) = \dfrac{\pi(p)P \big( N=n \, | \, p}{P(N=n)} \approx {n \choose n_0} p^n (1-p)^{n_0-n}\cdot \]

Observe que \[ f(p \, | \, N=n) \approx \underbrace{\pi(p)}_{\mbox{priori}} \times \underbrace{P\big( N=n \, | \, p\big)}_\mbox{verossimilhança}, \] como mecionado na introdução. Um modelo mais geral é assumir que a priori \(p\) tem uma distribuição Beta, \(p\sim B(\alpha,\beta)\).

Então \[ f(p \, | \, N=n) \approx p^{\alpha+n}(1-p)^{\beta+n_0-n}, \]

que é a densidade de uma distribuição Beta, \(B(\alpha+n,\beta+n_0-n)\).

Como consequência, \[ \mbox{E}\big(p \, | \, N=n \big) = \dfrac{\alpha+n}{\alpha+\beta+n_0} \approx \dfrac{n}{n_0} \]

quando \(n\) é grande suficiente e \[ \mbox{Var}\big(p \, | \, N=n\big) = \dfrac{(\alpha+\beta)(\beta+n_0-n)}{(\alpha+\beta+n_0)^2(\alpha+\beta++1)} \approx \dfrac{n(n_0-n)}{n_0^3} = \dfrac{1}{n_0}\dfrac{n}{n_0}\left( 1-\dfrac{n}{n_0}\right)\cdot \]

Neste exemplo, a distribuição da variável de interesse está na família exponencial, aqui uma distribuição Binomial e a distribuição a priori é a distribuição conjugada, aqui uma distribuição Beta.

Considere outro exemplo, onde após discutir com 81 agentes, 51 responderam positivamente a alguma proposição. Gostaríamos de testar se esta proposição excede 2/3, ou não.

Seja \(X = (X_1,\cdots,X_n)\) as respostas dos \(n\) agentes, \(X_i\) é 1 quando eles responderam positivamente, caso contrário 0. Aqui, \(X_i\sim B(\theta)\) e a interpretação Bayesiana do teste é baseada em o cálculo de \(P(\theta > 2/3 \, | \, X)\). Considere a distribuição Beta a priori para \(\theta\), de modo que \(\theta\sim B(\alpha,\beta)\).

Então \[ \theta \, | \, X=x \sim B\left( \alpha+\sum x_i \, , \, \beta+n-\sum x_i\right)\cdot \]

Para visualizar a distribuição a priori e a posteriori, considere o seguinte código: 

n <- 81
sumyes <- 51
priorposterior <- function(a,b){
 u <- seq(0,1,length=251)
 prior <- dbeta(u,a,b)
 posterior <- dbeta(u,sumyes+a,n-sumyes+b)
 plot(u,posterior,type="l",lwd=2, col="red")
 lines(u,prior,lty=2)
 grid()
 abline(v=(sumyes+a)/(n+a+b),col="grey")}

Assumindo como a priori \(\pi(\theta) \approx 1\), corresponde a uma distribuição uniforme, obtida quando \(\alpha=\beta= 1\). Uma alternativa pode ser considerar uma a priori mais informativa, por exemplo, \(\alpha= 5\) e \(\beta= 3\). 

par(mfrow = c(1,2))
priorposterior(1,1)
priorposterior(5,3)

FIGURA 1. Distribuições a priori (linha pontilhada) e a posteriori (linha contínua) para o parâmetro de uma distribuição Bernoulli.


Essas distribuições podem ser visualizadas na figura acima, sendo a linha pontilhada a distribuição a priori e a sólida a distribuição a posteriori. A linha reta vertical representa a média a posteriori.

Se voltarmos à questão inicial, \(P(\theta > 2/3 \, | \, X)\) pode ser obtido facilmente como:

1-pbeta(2/3,1+51,1+81-51)
## [1] 0.2273979
1-pbeta(2/3,5+51,3+81-51)
## [1] 0.2351575
com as duas a priori consideradas.


2.2. Motivação em pequenas amostras


Outro exemplo simples pode motivar a abordagem Bayesiana. Imagine que uma empresa implantou um novo produto há 5 anos e gostaria de uma estimativa prospectiva da sinistralidade para o próximo ano.

Aqui estão os dados disponíveis:

Qual poderia ser a distribuição da sinistralidade do próximo ano?

Podemos modelar esse problema assumindo que o processo aleatório atuarial produz uma taxa de perda lognormalmente distribuída \(L\): \[ \log(L)\sim N(\mu,\sigma)\cdot \]

Por \(\sigma\) representamos nossas crenças usando a constante 0.25, de modo que não importa o que seja \(\mu\), \(l\) geralmente está entre 0.75\(\mu\) e 1.25\(\mu\). Se \(\mu\) estiver próximo de 70% da média, isso implica um intervalo típico de 0.75\(\times\) 70% = 52.5% a 1.25\(\times\) 70% = 87.5%, que é semelhante a 70%\(\pm\) 20% do intervalo de nossas crenças.

Para \(mu\), precisamos reconhecer que a verdadeira sinistralidade é desconhecida, mas acreditamos que seja de 70%\(\pm\) 10%. Uma maneira simples de modelar essa crença é assumindo que a verdadeira sinistralidade mediana, \(e^\mu\) é lognormalmente distribuída com média de 70% e desvio padrão de 10%. Então \[ \log(e^\mu) = \mu \sim N(\mu_0,\sigma_0), \] onde \(\mu_0 = -0.367\) e \(\sigma_0 = 0.142\). Como \(\mu_0\) e \(\sigma_0\) determinam uma distribuição de \(\mu\), e ele próprio parametriza a distribuição, \(\mu_0\) e \(\sigma_0\) são hiperparâmetros. O uso de hiperparâmetros é comum nas estatísticas Bayesianas porque frequentemente temos crenças sobre quais podem ser os parâmetros reais de algum processo.

Usando a terminologia da seção anterior, as equações acima determinam nossa função de probabilidade a priori \(P(x)\). Agora precisamos condicionar em nossa evidência \(e\), as taxas de perda históricas, para determinar nossas crenças posteriores \(P_e(x) = P(x \, | \, e)\).

Para calcular a distribuição a posteriori de \(\mu\), podemos usar o Teorema de Bayes. Tem uma variedade de declarações, mas esta versão é comum para proposições e variáveis discretas, onde \(T\) pode ser pensado como representando a teoria e \(e\) para representar a evidência:

\[ P(T \, | \, e) = \dfrac{P(e \, | \, T)P(T)}{P(e)} = \dfrac{P(e \, | \, T)P(T)}{\displaystyle \sum_{i=1}^n P(e \, | \, T_i)P(T_i)}, \]

para variáveis contínuas, a versão a seguir é mais comum. \(\theta\) pode ser um escalar ou um vetor e geralmente representa parâmetros para uma teoria ou distribuição \[ p(\theta \, | \, e) = \dfrac{p(e \, | \, \theta)P(\theta)}{P(e)} = \dfrac{p(e \, | \, \theta)P(\theta)}{\displaystyle \int p(e \, | \, \phi)P(\phi)\mbox{d}\phi}\cdot \]

Todas as versões do Teorema de Bayes seguem da probabilidade básica e da definição de probabilidade condicional.

Neste caso, o Teorema de Bayes implica que nossa distribuição a posteriori para \(\mu\) é dada por \[ p(\mu \, | \, e) = \dfrac{p(e \, | \, \mu)P(\mu)}{\displaystyle \int p(e \, | \, \phi)P(\phi)\mbox{d}\phi}, \] onde \(e\) é agora os cinco índices de perda históricos observados. Agora podemos resolver teoricamente essa equação porque \(p(e \, | \, \mu)\) e \(p(e \, | \, \phi)\) são dados pela distribuição lognormal e \(p(\phi )\) é normalmente distribuído.

Felizmente, a equação tem uma solução simples neste caso. Como nossa probabilidade a priori para \(\mu\) é normal e a verossimilhança de cada índice de sinistralidade \(e\), em uma escala logarítmica, dada \(\mu\) também é normal, a distribuição a posteriori para \(\mu\) também é normal, não importa qual seja \(e\). Quando é garantido que a priori e a a posteriori estão na mesma família, eles são conhecidos como distribuições conjugadas ou conjugados Bayesianos.

Há duas outras razões pelas quais isso é conveniente:

  1. A a posteriori permanecerá na mesma família que a a priori quando um ponto de dados for observado. Assim, esse raciocínio pode ser repetido para mostrar que a a posteriori permanece na mesma família após qualquer número de observações.
  2. Os parâmetros da distribuição a posteriori são tipicamente funções simples dos parâmetros da distribuição a priori e dos dados observados.

Nosso problema é um exemplo do par conjugado normal/normal. Consulte a Tabela 2 para obter uma tabela de pares conjugados e suas equações. No caso normal/normal, sob a distribuição a posteriori, \[ \mu_1 \sim N\left( \Big(\dfrac{1}{\sigma_0^2} + \dfrac{5}{\sigma^2}\Big)^{-1}\Big( \dfrac{\mu_0}{\sigma_0^2}+\dfrac{1}{\sigma^2\sum_{i=1}^5 \log(r_i)}\Big) \, , \, \Big(\dfrac{1}{\sigma_0^2} +\dfrac{5}{\sigma^2}\Big)^{-1} \right), \]

onde \(r_i\) representa os cinco índices de sinistralidade observados. A média posterior de \(\mu\) é \(\mu_1 = -0.113\) e o desvio padrão a posteriori é \(\sigma_1 = 0.0879\). A situação pode ser representada graficamente.

Tabela 2. A a priori conjugadas para distribuições na família exponencial.

As taxas de perda observadas são marcadas com linhas verticais. Aqui está o conjunto de parâmetros a priori:

sigma <- 0.25
prior.mean <- 0.7
prior.sd <- 0.1
sigma0 <- sqrt(log((prior.sd / prior.mean)^2 + 1))
mu0 <- log(prior.mean) - prior.sd^2/2


Observações de registro estão aqui:

r <- c(0.958, 0.614, 0.977, 0.921, 0.756)
n <- length(r)


Usando a discussão anterior, podemos calcular probabilidades a posteriori

sigma1 <- (sigma0^-2 + n / sigma^2)^-.5
mu1 <- (mu0 / sigma0^-2 + sum(log(r)) / sigma^2) * sigma1^2


Essas distribuições podem ser visualizadas usando

u <- seq (from = 0.3, to = 1.3, by = 0.01)
prior <- dlnorm(u,mu0,sigma0)
posterior <- dlnorm(u,mu1,sigma1)
plot(u,posterior,type="l",lwd=2)
lines(u,prior,lty=2)
grid()


Concluímos este exemplo com algumas observações:

  1. Nossa distribuição em torno do índice de sinistralidade real \(\mu\) aumentou após observar os 5 anos de dados. Isso era esperado porque os índices de sinistralidade observados foram, em geral, maiores do que nossa média anterior de 70%.
  2. A distribuição a posteriori da verdadeira sinistralidade é mais estreita, tem um desvio padrão menor. Isso também era esperado depois de observar os dados, sabemos mais e, portanto, nossa incerteza sobre a verdadeira sinistralidade diminuiu.
  3. A distinção conceitual entre o parâmetro e as distribuições preditivas é bastante importante. As distribuições dos parâmetros são distribuições sobre os parâmetros para outras distribuições, \(\mu\) neste caso. Se quisermos prever novas observações, devemos usar a distribuição preditiva. A distribuição dos parâmetros apenas quantifica o risco do parâmetro.
  4. As distribuições preditivas são significativamente mais incertas do que as distribuições dos parâmetros. Isso também é bastante intuitivo, grande parte da variabilidade nas perdas de seguros é contribuída pelo risco do processo.
  5. Nesse caso, nossa distribuição preditiva a posteriori é realmente mais incerta do que nossa distribuição preditiva anterior (32.9% versus 30.5%). Embora haja menos risco dos parâmetros, há mais risco do processo devido à maior taxa de perda esperada.


2.3. Cisnes Negros e Metodologia Bayesiana


O uso da abordagem de Bayes é interessante quando você não tem nenhuma experiência. Considere o seguinte caso simples. Considere uma amostra Bernoulli \(\{0,0,0,0,0\}\). O que você poderia dizer sobre o parâmetro de probabilidade? Como frequentista, nada. Esta foi realmente a pergunta feita a Longley-Cook na década de 1950: é possível prever a probabilidade de dois aviões colidirem no ar? Nunca houve nenhuma colisão séria de aviões comerciais naquela época. Sem qualquer experiência passada, os estatísticos não poderiam responder a essa pergunta.

Suponha que os \(X_i\)’s sejam variáveis i.i.d. Variáveis \(B(p)\). Com um beta a priori para \(p\), com parâmetros \(\alpha\) e \(\beta\), então \(p\) dado \(X\) tem uma distribuição beta com parâmetros \(\alpha\) e \(\beta+5\). Assim, com a a priori \(\alpha = \beta= 1\),

qbeta(.95,1,6)
## [1] 0.3930378

o intervalo de confiança de 95% para \(p\) é [0%; 40%]. Com a priori não informativa de Jeffrey, consulte a Seção 2.6.,

qbeta(.95,.5,.5+5)
## [1] 0.3057456

o intervalo de confiança de 95% para p é [0%; 30%]. É possível visualizar a transformação posterior, com as informações chegando, na Figura 3.3. Para o uniforme anterior,

par(mfrow=c(1,2))
pmat =persp(0:5,0:1,matrix(0,6,2),zlim=c(0,5), ticktype= "detailed",
            box = FALSE, theta =-30)
title("A priori Uniforme",cex.main=.9)
y=seq(0,1,by=.01)
for(k in 0:5){
 z=dbeta(y,1,1+k)
 indx=which(y<=qbeta(.95,1,1+k))
 xy3d=trans3d(rep(k,2*length(indx)),c(rev(y[indx]),y[indx]),
 c(rep(0,length(y[indx])),z[indx]),pmat)
 polygon(xy3d,col="grey",density=50,border=NA,angle=-50)
 xy3d=trans3d(rep(k,length(y)),y,z,pmat)
 lines(xy3d,lwd=.5)
 xy3d=trans3d(rep(k,length(indx)),y[indx],z[indx],pmat)
 lines(xy3d,lwd=3)}
grid()

pmat =persp(0:5,0:1,matrix(0,6,2),zlim=c(0,7), ticktype= "detailed",
            box = FALSE, theta =-30)
title("A priori não-informativa de Jeffrey",cex.main=.9)
y=seq(0,1,by=.01)
for(k in 0:5){
 z=dbeta(y,.5,.5+k)
 indx=which(y<=qbeta(.95,.5,.5+k))
 xy3d=trans3d(rep(k,2*length(indx)),c(rev(y[indx]),y[indx]),
 c(rep(0,length(y[indx])),z[indx]),pmat)
 polygon(xy3d,col="grey",density=50,border=NA,angle=-50)
 xy3d=trans3d(rep(k,length(y)),y,z,pmat)
 lines(xy3d,lwd=.5)
 xy3d=trans3d(rep(k,length(indx)),y[indx],z[indx],pmat)
 lines(xy3d,lwd=3)}
grid()

FIGURA 3. Distorção da distribuição a posteriori, com a chegada da informação.

Voltando à colisão do avião, em 1955, Longley-Cook previu “qualquer coisa de 0 a 4 […] colisões nos próximos dez anos”, sem nenhuma experiência, portanto, comprar um resseguro pode ser uma boa ideia. Dois anos depois, 128 pessoas morreram no Grand Canyon e 4 anos depois, mais de 133 pessoas morreram na cidade de Nova York.


2.4. Modelos Bayesianos em gestão de portfólio e finanças


Seguindo Frost & Savarino (1986) é possível usar uma metodologia Bayesiana para discutir a seleção de portfólio. Considere um ativo e uma sequência mensal de retornos para um determinado ativo, uma série temporal, \(r = (r_1,\cdots,r_T)\).

Seja \(\mu\) e \(\sigma\) a média e a volatilidade desse retorno. Suponha que os retornos sejam sorteios independentes de uma distribuição gaussiana, \(f(r \, | \, \mu )\). Considere o conjugado a priori, \(\mu\sim N(\mu_0,\tau_0)\).

Então, a distribuição a posteri é gaussiana \(N(\mu_T,\tau_T)\), onde \[ \mu_T = \left( \dfrac{\mu_0}{\tau_T^2}+\dfrac{T}{s^2}\dfrac{1}{T}\sum_{t=1}^T r_i\right)\left( \dfrac{1}{\tau_T^2}+\dfrac{T}{s^2}\right)^{-1} \] e

\[ \tau_T=\left( \dfrac{1}{\tau_T^2}+\dfrac{T}{s^2}\right)^{-1}\cdot \]

A média a posteri \(\mu_T\) é aqui uma média ponderada da média a priori \(\mu_0\) e da média histórica. Observe que quanto menos volátil a informação, maior a sua precisão e à medida que mais dados se tornam disponíveis, as informações a priori se tornam cada vez mais irrelevantes, a média a posteriori se aproxima da média histórica.

Claro, é possível estender esse modelo a dimensões maiores, com vários ativos. Considere uma série de retornos mensais, \((r_1,\cdots,r_T)\). Seja \(\mu\) e \(\Sigma\) o vetor de médias e a matriz de volatilidade desses retornos. Suponha que os retornos sejam retiradas independentes da distribuição Gaussiana, \(f(r \, | \, \mu)\). Considere o conjugado a priori para as médias \(\mu\sim N(\mu_0,\tau_0)\).

Aqui, \(\tau_0\) reflete nossos antecedentes na covariância das médias, não a covariância dos retornos dos ativos. Na verdade, também podemos considerar a priori na covariância dos retornos. A distribuição a priori da matriz de covariâncias possui uma distribuição Wishart inversa, com \(\nu\) graus de liberdade, \(\Sigma\sim IW(\nu,\Sigma_0)\).

Então, a distribuição a posteriori é gaussiana \(N(\mu_T,\Sigma_T)\), onde \[ \mu_T = \dfrac{\left(\dfrac{\mu_0}{\tau_0}+ \widehat{\Sigma}^{-1}T\widehat{\mu}\right)}{\dfrac{1}{\tau_0}+T\widehat{\Sigma}^{-1}} \] e \[ \Sigma_T = \dfrac{T+1+\nu}{T+\nu-d-2}\left( \dfrac{\nu}{\nu+T}\Sigma_0+\dfrac{T}{\nu+T}\widehat{\Sigma}\right), \]

onde \(d\) é o número de ativos considerados. Médias e covariâncias a posteriori são misturas lineares entre as estimativas a priori e empíricas, ou seja, históricas.

Se tivermos a priori não informativa \(\nu= 0\) e \(\tau_0^{-1}=0\), então \[ \mu_T = \widehat{\mu} \qquad \mbox{e} \qquad \Sigma_T = \dfrac{T+1}{T-d-2}\widehat{\Sigma}\cdot \]

Uma extensão é o chamado modelo Black & Litterman, introduzido em Black & Litterman (1992) e discutido em Satchell & Scowcroft (2000). Este modelo é popular entre os atuários e mais geralmente os investidores; porque permite distinguir as previsões da convicção ou crenças, sobre a incerteza associada à previsão.


2.5. Relação com a credibilidade da Bühlmann


Na estrutura de credibilidade, assumimos que as carteiras de seguros são heterogêneas. Existe um fator de risco subjacente, não observável, para todos os segurados e gostaríamos de modelar a frequência de sinistro individual para questões de taxação. Seja \(Y_{i,t}\) o número de sinistros para o ano \(t\) e seguradora \(i\). Aqui gostaríamos de usar observações históricas passadas, ou seja, usaremos \(\mbox{E}(Y_{i,t+1} \, | \, \pmb{Y}_{i,t})\) como uma aproximação para \(\mbox{E}(Y_{i,t+1} \, | \, \theta_i)\), onde \(\pmb{Y}_{i,t} = (Y_{i}, Y_{i,t-1}, Y_{i,t-2},\cdots)\).

Considere um contrato, observado durante \(T\) anos. Seja \(Y_t\) o número de reivindicações e suponha que \(Y_t\in \{0,1\}\), de modo que \(Y_t\sim B(p)\). Da seção anterior, se assumirmos distribuição a priori Beta \(B(\alpha,\beta)\) para \(p\), temos \[ \mbox{E}(p |, | \, \pmb{Y}) = \dfrac{\alpha+T\widehat{Y}_T}{\alpha+\beta+T}=\dfrac{T}{\underbrace{\alpha+\beta+T}_Z} \widehat{Y}_T+\underbrace{\left(1-\dfrac{T}{\alpha+\beta+T} \right)}_{1-Z}\dfrac{\alpha}{\underbrace{\alpha+\beta}_{\mbox{E}(p)}} \] onde \(Z=\in [0,1]\). Porque \[ \mbox{E}(Y_t \, | \, p) =p \qquad \mbox{e} \qquad \mbox{Var}(Y_t \, | \, p) = p(1-p)\cdot \]

Então \[ \dfrac{\mbox{E}\big(\mbox{Var}(Y_t \, | \, p) \big)}{\mbox{Var}\big(\mbox{E}(Y_t \, | \, p) \big)}=\dfrac{\mbox{E}\big(p(1-p) \big)}{\mbox{Var}(p)}=\dfrac{\mbox{E}(p)-\big(\mbox{Var}(p)+\mbox{E}^2(p)\big)}{\mbox{Var}(p)} = \alpha+\beta\cdot \]

Desta forma, \[ Z=\dfrac{T}{\alpha+\beta+T} = \dfrac{T}{k+T}, \qquad \mbox{onde} \qquad k = \dfrac{\mbox{E}\big(\mbox{Var}(Y_t \, | \, p) \big)}{\mbox{Var}\big(\mbox{E}(Y_t \, | \, p) \big)}\cdot \] Isso será conhecido como abordagem de credibilidade de Bühlmann: Se o segurado \(i\) foi observado \(T\) anos, com experiência anterior \(\pmb{Y}_T = (Y_1,\cdots,Y_T)\), então \(\mbox{E}(Y_{T+1} \, | \, \pmb{Y}_T )\) pode ser aproximado por \[ Z\overline{Y}_{i,T}+(1-Z)\mu \qquad \mbox{onde} \qquad \mu = \overline{Y}_T=\dfrac{1}{T}\sum_{t=1}^T Y_{i,t}, \]

e \(\mu\) é a média de toda a carteira. \(Z\) é então a parcela do prêmio que se baseia em informações passadas para um determinado segurado e está relacionada à quantidade de credibilidade que damos a informações passadas.

De Bühlmann (1969), inspirado na abordagem Bayesiana descrita anteriormente, \(Z\) satisfaz \[ Z=\dfrac{T}{k+T} \qquad \mbox{onde} \qquad k = \dfrac{\mbox{E}\big(\mbox{Var}(Y_t \, | \, p) \big)}{\mbox{Var}\big(\mbox{E}(Y_t \, | \, p) \big)}\cdot \]

Como mencionado em Bühlmann & Gisler (2005), em um modelo Bayesiano, desejamos calcular \(\mbox{E}(Y_{T+1} \, | \, \pmb{Y}_T)\), mas aqui nos restringimos a projeções em combinações lineares de experiências passadas. O estimador de credibilidade de Bühlmann é a melhor aproximação linear de mínimos quadrados para o estimador de credibilidade Bayesiano. Quando \(X_{i,t}\)’s são variáveis aleatórias i.i.d. na família exponencial e se a distribuição a priori de \(\theta\) for conjugada para esta família, então o estimador de credibilidade de Bühlmann é igual ao estimador de credibilidade Bayesiano. Este foi o caso anteriormente, onde \(X_{i,t}\)’s eram Bernoulli e \(\theta\) tem um Beta a priori. Outro exemplo popular é obtido quando os \(X_{i,t}\)’s são Poisson e \(\theta\) têm uma Gama a priori.

Para ilustrar a credibilidade de Bühlmann, considere as contagens de sinistros de Norberg (1979), \(Y_{i,t}\) para contratos \(i = 1,\cdots,m\) e tempo \(t = 1,\cdots,T\). 

norberg = matrix( c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    1, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    0, 1, 1, 0, 1, 1, 1, 0, 0, 1,
                    0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                    1, 1, 0, 0, 1, 0, 0, 0, 1, 0,
                    0, 0, 0, 0, 0, 1, 0, 1, 0, 1,
                    0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
                    0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    1, 1, 0, 1, 0, 0, 1, 0, 0, 1,
                    1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0), byrow = TRUE, ncol = 10)
colnames(norberg) = c("year00", "year01", "year02", "year03", "year04",
                       "year05", "year06", "year07", "year08", "year09")

Para cada contrato, desejamos prever o valor esperado \(X_{i,T+1}\) dadas as observações anteriores. A estimativa de Bühlman é \[ Z\times \overline(Y)_i+(1-Z)\times \overline{Y}, \]

onde

\[ Z = \dfrac{T}{k+T}, \qquad \overline{Y}_i=\dfrac{1}{T}\sum_{t=1}^T Y_{i,t}, \qquad \overline{Y} = \dfrac{1}{m}\sum_{i=1}^m \overline{Y}_i \qquad \mbox{e} \qquad k=\dfrac{s^2}{a}, \] sendo \[ s^2= \dfrac{1}{m(T-1)}\sum_{i=1}^m\sum_{t=1}^T (X_{i,t}-\overline{X}_i)^2 \qquad \mbox{e} \qquad a=\dfrac{1}{m-1}\sum_{i=1}^m (\overline{X}_i-\overline{X})^2 -\dfrac{s^2}{n}; \]

ver Herzog (1996) para uma discussão desses estimadores. O código para calcular a estimativa de Bühlman para o ano 11 será:

T<- ncol(norberg)
(m <- mean(norberg))
## [1] 0.145
(s2 <- mean(apply(norberg,1,var)))
## [1] 0.1038889
(a <- var(apply(norberg,1,mean))-s2/T)
## [1] 0.02169006


O fator de credibilidade de Bühlman é baseado em

(Z <- T/(T+s2/a))
## [1] 0.6761462


e as estimativas de Bühlman são

Z*apply(norberg,1,mean)+(1-Z)*m
##  [1] 0.0469588 0.0469588 0.1821880 0.0469588 0.1821880 0.0469588 0.1821880
##  [8] 0.0469588 0.4526465 0.1145734 0.3174173 0.2498027 0.1145734 0.1145734
## [15] 0.0469588 0.0469588 0.3850319 0.1145734 0.1145734 0.0469588


2.6. A priori não informativa


Na primeira parte desta seção, usamos a priori conjugadas devido à facilidade computacional. Como veremos na próxima seção, não ter uma expressão analítica simples não é grande coisa. Portanto, pode ser interessante ter a priori mais críveis, mesmo que não seja um problema simples, como discutido formalmente em Berger et al. (2009).

Um exemplo simples para começar pode ser o caso em que \(X\) dado \(\theta\) tem uma distribuição de Bernoulli. Como vimos anteriormente, sobre nascimentos, a ideia de Pierre Simon Laplace era assumir uma distribuição a priori uniforme, \(\pi(\theta) = 1\). Assumindo que \(\pi(\theta) \propto 1\), para distribuições mais gerais será chamado plano a priori. Se Pierre-Simon Laplace pensou que seria o mais neutro possível ou não informativo, na verdade não é o caso. Uma ideia mais moderna de formalizar a ideia de a priori não informativa (ver Yang & Berger (1998) para uma discussão) é que devemos obter um resultado equivalente ao considerar um parâmetro transformado. Dada alguma transformação um-para-um \(h(\cdot)\), de modo que o parâmetro não se afaste demais de \(\theta\), mas \(\tau= h(\theta)\), queremos aqui ter distribuições invarintes.

Mais formalmente, lembre-se de que a densidade de \(\tau\) é aqui \[ p(\tau)=\pi(\theta) \times \left|\dfrac{\partial \theta}{\partial \tau} \right|\cdot \]

Lembre-se que este Jacobiano tem uma interpretação em termos de informação, pois a informação de Fisher de \(\tau\) é \[ \mbox{I}(\tau)=-\mbox{E}\left( \dfrac{\partial^2 \log\big( f(x \, | \, \tau)\big)}{\partial\tau^2}\right)= \mbox{I}(\theta)\times \left|\dfrac{\partial \theta}{\partial \tau} \right|^2\cdot \]

A a densidade de \(\tau\) pode então ser escrita como \(p(\tau)\sqrt{\mbox{I}(\tau)}=\pi(\theta)\sqrt{\mbox{I}(\theta)}\). Assim, é natural ter uma densidade a priori proporcional à raiz quadrada da informação de Fisher, \[ \pi(\theta)\propto \sqrt{\mbox{I}(\theta)}\cdot \] Isso é chamado de princípio de Jeffrey.

Para a distribuição Poisson, significa que \(\pi(\theta) \propto \theta^{-1/2}\) e se voltarmos ao problema inicial, nos nascimentos, a priori não informativa será uma distribuição Beta com parâmetros 1/2, pois \[ \pi(\theta) \propto \sqrt{\theta\Big(\frac{1}{\theta}-\dfrac{\theta}{1-\theta} \Big)^2+(1-\theta)\Big(\frac{0}{\theta}-\dfrac{\theta}{1-\theta} \Big)^2} = \dfrac{1}{\sqrt{\theta(1-\theta))}}\cdot \]


3. Considerações computacionais


Em muitas aplicações, desejamos calcular \[ \mbox{E}\big( g(\theta) \, | \, x\big) = \dfrac{\displaystyle \int g(\theta)\pi^*(\theta \, | \, x) \mbox{d}\theta}{\displaystyle \int \pi^*(\theta \, | \, x) \mbox{d}\theta}, \]

onde \(\pi^*\) é proporcional à densidade a posteriori. Um método possível é usar a aproximação de Gauss-Hermite, discutida em Klugman (1992) em detalhes. Uma alternativa é baseada na integração de Monte Carlo.


3.1. Maldição da dimensionalidade


Infelizmente, além dos conjugados Bayesianos, a maioria dos problemas Bayesianos são difíceis de calcular. O problema é que a distribuição a posteriori é muitas vezes contínua e multidimensional e tem o formato \[ \dfrac{p(e \, | \, \theta)p(\theta)}{\displaystyle \int p(e \, | \, \phi)p(\phi)\mbox{d}\phi} \] para o vetor \(n\)-dimensional \(\theta\).

Embora \(p(e \, | \, \theta)p(\theta)\) seja frequentemente fácil de calcular, a constante de normalização \[ \int p(e \, | \, \phi)p(\phi)\mbox{d}\phi \] é uma integral \(n\)-dimensional. Normalmente, a menos que sejam usadas distribuições conjugadas, essa integral não é analiticamente solúvel e deve ser aproximada numericamente.

Aqui, a estatística Bayesiana sofre da maldição da dimensionalidade: o rápido aumento na dificuldade computacional de integração à medida que o número de dimensões aumenta.

Em geral, se precisarmos integrar uma função sobre o hipercubo unidade \(d\)-dimensional \([0,1]^d\), então esse pensamento sugere que precisaremos de cerca de \(100^d\) avaliações. Assim, o número de funções ou chamadas aumentará rapidamente com a dimensionalidade. Os computadores são muito mais rápidos do que costumavam ser, mas mesmo problemas moderadamente pequenos tornam-se computacionalmente inviáveis.

Existem algoritmos de integração numérica, como a quadratura gaussiana, que podem reduzir um pouco o número de pontos necessários, mas seu erro ainda diminui na taxa lenta de \(O(d^{-1\frac{1}{n}})\). Para piorar as coisas, as distribuições de probabilidade a posteriori Bayesianas são muitas vezes pontiagudas, ou seja, elas têm grandes áreas com baixa probabilidade e áreas relativamente pequenas com muita probabilidade, como as distribuições mostradas na Figura 1.

Do ponto de vista computacional, nosso interesse é aproximar o valor de uma integral, que é uma área para um função \(\mathbb{R}\to \mathbb{R}\) ou um volume para uma uma função \(\mathbb{R}^d\to\mathbb{R}\). Uma ideia simples é usar uma caixa, ao redor do volume, para particionar a caixa e contar a proporção de pontos no volume considerado. Por exemplo, para calcular a área de um disco com raio 1, considere uma caixa \([-1,+1]\times [-1, +1]\) e calcule a proporção de pontos em uma grade na caixa que pertence ao disco

diskincube <- function(n){
  gridn <- seq(-1,+1,length=n)
  gridcube <- expand.grid(gridn,gridn)
  inthedisk <- apply(gridcube^2,1,sum)<=1
  mean(inthedisk)
}
diskincube(200)
## [1] 0.7766

Aqui, 77.66% dos pontos do quadrado pertencem ao disco e como a área do cubo é \(2^2\), significa que a área do disco deve ser próxima de

diskincube(200)*2^2
## [1] 3.1064


diskincube(2000)*2^2
## [1] 3.138388
o valor verdadeiro é \(\pi\). Observe que sempre subestimamos o valor real com essa técnica.

Este método pode ser estendido para dimensões mais altas. Por exemplo, na dimensão 4,

diskincube.dim4 <- function(n){
  gridn <- seq(-1,+1,length=n)
  gridcube <- expand.grid(gridn,gridn,gridn,gridn)
  inthedisk <- apply(gridcube^2,1,sum)<=1
  mean(inthedisk)
}

O volume da esfera unitária na dimensão 4 é maior do que

diskincube.dim4(40)*2^4
## [1] 4.4488
sendo o valor verdadeiro \(\pi^2/2 \approx 4.9348\).

Mais formalmente, desejamos estimar um volume \(\nu\subset \mathbb{R}^d\), utilizando \(m=k^d\) pontos, para algum inteiro \(k\). Se \(\nu\in [a,b]\), considere para qualquer \(i\) a partição uniforme \([a_i,b_i]\) em \(k\) segmentos. Então, definimos \(x_{i,j}=a_i+(j-1)(b_i-a_i)/(k-1)\), para todo \(j=1,\cdots,k\).

Uma aproximação do volume \(V\) é então \[ V(\nu) \approx \left(\prod_{i=1}^d (b_i-a_i) \right)\times \dfrac{1}{n}\sum_{j_1,\cdots,j_d \, \in \, \{1,\cdots,k\}}\pmb{1}\big( (x_{1,j_1},\cdots,x_{d,j_k})\in \nu \big) \cdot \]

Seja \(V_n(\nu)\) denota esta aproximação. Se \(L(\nu)\) o comprimento do contorno do volume, na dimensão 2, então \[ |V_n(\nu)-V(\nu)| \leq \dfrac{L(\nu)}{k}=\dfrac{L(\nu)}{n^{1/d}}\cdot \]

O ponto aqui é que o erro de aproximação é \(O(1/n^{1/d})\), que aumentará com a dimensão.

Isso é mais ou menos o que fazemos quando calculamos integrais, \[ \int_{[0,1]} h(u)\mbox{d}u \approx \sum_{j=1}^n \dfrac{1}{n}\times h\Big(\frac{j}{n}\Big) \] e em altas dimensões, \[ \int_{[0,1]^d} h(u)\mbox{d}u \approx \sum_{j_1,\cdots,j_d \, \in \, \{1,\cdots,k\}} \Big(\dfrac{1}{k}\Big)^d\times h\Big(\frac{j_1}{k},\cdots,\frac{j_d}{k}\Big)\cdot \]

Claro, esse algoritmo numérico pode ser melhorado, usando trapézios em vez de retângulos ou utilizando a interpolação polinomial de Simpson. Mas a ordem permanecerá inalterada e o erro é da ordem \(O(n^{-1/d})\). Isso é|mais formalmente|o que é chamado de maldição da dimensionalidade.


3.2. Integração Monte Carlo


Agora, se voltarmos ao exemplo inicial, é possível ver o que são os métodos de Monte Carlo: Ao invés de considerar uma grade homogênea, uma ideia pode ser simplesmente desenhar aleatoriamente \(n\) pontos, uniformemente, na caixa. Uma aproximação do volume será \[ V(\nu) \approx \left(\prod_{i=1}^d (b_i-a_i) \right)\times \dfrac{1}{n}\sum_{i=1}^n \pmb{1}\big(X_i\in\nu \big) \cdot \]

Seja \(\widetilde{V}_n(\nu)\) esta aproximação. Para discutir e quantificar erros, observe o fato de que um ponto estar ou não dentro do volume \(\nu\) pode ser modelado usando variáveis aleatórias Bernoulli. É então possível provar que \[ \mbox{E}\big( \widetilde{V}_n(\nu)\big) = V(\nu) \] e \[ \mbox{Var}\big(\widetilde{V}_n(\nu) \big) = \dfrac{1}{n}\times \dfrac{V(\nu)}{\displaystyle \prod_{i=1}^d (b_i-a_i)}\times \left( 1-\dfrac{V(\nu)}{\displaystyle \prod_{i=1}^d (b_i-a_i)}\right)\cdot \]

Agora, essa ideia pode ser usada para aproximar integrais, como \[ \int_{[0,1]^d} h(u)\mbox{d}u = \mbox{E}\big(h(U)\big) \] onde \(U\) tem componentes independentes com \(U_i\sim U([0, 1])\) e do Teorema dos Grandes Números \[ \frac{1}{n}\sum_{i=1}^n h(U_i) \xrightarrow{\mbox{q.c.}} \mbox{E}\big(h(U)\big)\cdot \]

Para quantificar a incerteza, de ne \[ \widetilde{V}_n(\nu) = \dfrac{n}{n-1}\left( \dfrac{1}{n}\sum_{i=1}^n h(U_i)^2- \Big(\dfrac{1}{n}\sum_{i=1}^n h(U_i) \Big)^2\right), \]

e do Teorema do LimiteCcentral, o erro de forma mais probabilística é, então, de ordem \(O(1/\sqrt{n})\). Observe que a ordem não depende da dimensão \(d\) aqui, na verdade, depende, no termo constante. Assim, simulações de Monte Carlo podem ser mais interessantes em alta dimensão do que métodos determinísticos. 


# diskincube.dim4(40^4)*2^4 # 4.943763

O valor verdadeiro é \(\pi^2/2 \approx 4.9348\). Lembre-se de que, ao usar métodos determinísticos, ainda estávamos longe do verdadeiro valor.


Para comparar o algoritmo determinístico e a abordagem Monte Carlo, podemos traçar as aproximações obtidas usando os dois algoritmos, ver Figura 4 abaixo: 

v <- 1:40
plot(v,2^4*Vectorize(diskincube.dim4)(v),type="b",ylim=c(0,5))
#lines(v,2^4*Vectorize(diskincube.dim4)(v^4),type="b",pch=19,cex=.6,lty=2)
abline(h=pi^2/2,col="grey")
grid()

FIGURA 3. Maldição da dimensionalidade, calculando o volume da esfera unitária em \(R^4\), com o método determinístico e usando simulações de Monte Carlo, em função do número de pontos, \(n = 4^k\).


3.3. MCMC (Markov Chain Monte Carlo)


Atualmente, a maneira mais popular na estatística Bayesiana de contornar a maldição da dimensionalidade é usar técnicas de Markov Chain Monte Carlo (MCMC). Curiosamente, a primeira técnica MCMC, o algoritmo Metropolis, foi desenvolvida pelo físico homônimo em 1953. Foi executada em alguns dos primeiros computadores.

Apesar desse início precoce, o algoritmo Metropolis e as técnicas MCMC relacionadas só ganharam popularidade na comunidade estatística na década de 1990. O ressurgimento da estatística Bayesiana desde então deve muito à eficácia e praticidade desses algoritmos.

As técnicas MCMC nem sequer tentam calcular a integral \(\int p(e \, | \, \phi)p(\phi)\mbox{d}\phi\). Em vez disso, elas fornecem amostras de \(\theta\) da distribuição a posteriori de \(p(\theta \, | \, e)\). As técnicas MCMC são Monte Carlo, o que significa que dependem de amostragem aleatória, mas também produzem Cadeias de Markov, o que significa que as amostras que produzem não são independentes. Cada amostra depende da anterior.

O algoritmo Metropolis usa a percepção de que não precisamos calcular \(\int p(e \, | \, \phi)p(\phi)\mbox{d}\phi\) para obter amostras; em vez disso, podemos trabalhar com a distribuição de probabilidade não normalizada \(p(e \, | \, \theta)p(\theta)\), que geralmente é mais fácil de calcular.

O algoritmo pode ser resumido informalmente assim:

  1. Começar com um valor arbitrário na distribuição a posteriori \(\theta_0\).
  2. Dada uma amostra existente \(\theta\), escolha aleatoriamente uma amostra candidata \(\theta'\) próxima. Isso é feito por amostragem de outra distribuição de probabilidade \(Q(\cdot \, | \, \theta_i)\). Freqüentemente, \(Q\) é escolhida para ser multinormal, de modo que \(\theta'\sim N( \theta_i,\Sigma)\), para algum \(\Sigma\) fixo.
  3. Comparar \(p(e \, | \, \theta_i)p(\theta_i)\) e \(p(e \, | \, \theta')p(\theta')\) calculando a taxa de aceitação \[ R = \dfrac{p(e \, | \, \theta')p(\theta')}{p(e \, | \, \theta_i)p(\theta_i)}\cdot \]
  4. Moverse para \(\theta'\) se \(R\) for grande o suficiente: com probabilidade \(\min(1,R)\) defina \(\theta_{i+1} = \theta'\), isso é chamado de aceitação de \(\theta'\); com probabilidade \(1-\min(1,R)\) defina \(\theta_{i+1} = \theta_i\), isso é chamado de rejeição de \(\theta'\).
  5. Vá para a etapa (2) e continue até que amostras suficientes sejam retiradas.

Aqui não será demonstrado que as amostras assim extraídas podem ser consideradas amostras da distribuição a posteriori \(p(\cdot \, | \, e)\). No entanto, espera-se que fique claro a partir do algoritmo que ele gastará mais tempo onde \(p(e \, | \, \theta´)p(\theta´)\) for grande e menos tempo onde for pequeno.

Isso ocorre porque a taxa de aceitação será sempre 1 quando o candidato nos mover para uma área de maior probabilidade. Quando a amostra candidata nos move para uma área de menor probabilidade, a taxa de aceitação será menor que 1 e podemos ficar onde estamos.

Isso ocorre porque a taxa de aceitação será sempre 1 quando o candidato nos mover para uma área de maior probabilidade. Quando a amostra candidata nos move para uma área de menor probabilidade, a taxa de aceitação será menor que 1 e podemos ficar onde estamos.

Porque \(p(e \, | \, \theta´)p(\theta´)\) é grande exatamente nos mesmos lugares que a distribuição a posteriori \[ p(\theta´ \, | \, e) = \dfrac{p(e \, | \, \theta´)p(\theta´)}{\displaystyle \int p(e \, | \, \phi)p(\phi)\mbox{d}\phi}, \]

for grande, mais amostras serão retiradas das áreas do espaço de probabilidade onde a densidade a posteriori é grande.

Uma vez que temos amostras da distribuição a posteriori, podemos estimar a resposta para a maioria das questões que surgem em problemas estatísticos Bayesianos. Por exemplo, suponha que temos \(n\) amostras de \(\theta\) e queremos saber o valor esperado do primeiro parâmetro \(\theta_0\). Então podemos usar a estatística clássica para aproximá-lo: \[ \mbox{E}(\theta_0)\approx \frac{1}{n}\sum_{i=1}^n \theta_{0,i}\cdot \]

É irônico, mas muitas ou a maioria das técnicas Bayesianas no final dependem da estatística clássica para entender o comportamento produzido pelas cadeias MCMC.


3.4. Exemplo de MCMC no R


Vamos examinar um problema Bayesiano simples que não pode ser resolvido usando conjugados Bayesianos, mas pode ser rapidamente analisado com uma técnica MCMC. Lembre-se do exemplo na Seção 2.

O objetivo foi encontrar a média posterior e o desvio padrão (risco de parâmetro) da sinistralidade esperada. Assumimos que o risco do processo era distribuído log-normal com \(\sigma=0.25\), ignorando assim a oportunidade de aprender sobre o risco do nosso processo a partir dos dados.

Desta vez, vamos usar os mesmos números, mas modelar a incerteza sobre \(\sigma\). Assumiremos novamente que o risco de processo dos índices de sinistralidade reais é distribuído log-normal, com parâmetros \(\mu\) e \(\sigma\). A distribuição a priori de \(\mu\) será novamente normalmente distribuída com os mesmos parâmetros \(\mu_0\) e \(\sigma_0\). No entanto, agora nossa distribuição a priori de \(\sigma\) será gama inversa com parâmetros \(k_0\) e \(\theta_0\), \(\sigma\) não pode ser negativa, então a distribuição Gama inversa é comumente usada para modelar uma distribuição da variância a priori.

Podemos definir \(k_0\) e \(\theta_0\) usando o método dos momentos. 

prior.mean.mean <- 0.7
prior.mean.sd <- 0.1
prior.sd.mean <- 0.25
prior.sd.sd <- 0.1


Usando o fato de que para uma distribuição lognormal, o coeficiente de variação é \(\sqrt{e^{s^2}-1}\), 


sigma0 <- sqrt(log((prior.mean.sd / prior.mean.mean)^2 + 1))
mu0 <- log(prior.mean.mean) - sigma0^2 / 2


e para a distribuição Gama inversa \[ \mbox{E}(X)=\dfrac{\theta_0}{k-1} \] e o coeficiente de variação é \[ \dfrac{1}{\sqrt{k-1}}, \] 

k0 <- 2 + (prior.sd.mean / prior.sd.sd)^2
theta0 <- (k0 - 1) * prior.sd.mean


As observações registradas foram 

r <- c(0.958, 0.614, 0.977, 0.921, 0.756)
r.log <- log(r)
n <- length(r)


Para usar a função dinvgamma(), use o pacote de funções MCMCpack. O código para executar simulações MCMC está aqui: 

library(MCMCpack)
RunSim <- function(M, delta, mu, sigma) {
  output.df <- data.frame(mu=rep(NA, M), sigma=NA)
  set.seed(0)
  cur.prior.log <- (dnorm(mu, mu0, sigma0, log=TRUE)
                    + log(dinvgamma(sigma, k0, theta0)))
  cur.like.log <- sum(dnorm(r.log, mu, sigma, log=TRUE))
  for (i in seq_len(M)) {
    mu.cand <- mu + rnorm(1, sd=delta)
    sigma.cand <- max(1e-5, sigma + rnorm(1, sd=delta/2))
    cand.prior.log <- (dnorm(mu.cand, mu0, sigma0, log=TRUE)
                    + log(dinvgamma(sigma.cand, k0, theta0)))
    cand.like.log <- sum(dnorm(r.log, mu.cand, sigma.cand, log=TRUE))
    cand.ratio <- exp(cand.prior.log + cand.like.log
                    - cur.prior.log - cur.like.log)
  if (runif(1) < cand.ratio) {
    mu <- mu.cand
    sigma <- sigma.cand
    cur.prior.log <- cand.prior.log
    cur.like.log <- cand.like.log
  }
  output.df[i, ] <- c(mu, sigma) # escrevendo uma amostra para a saída
  }
  return(output.df)
  }
delta005.df <- RunSim(6000, .005, log(1.2), 0.5)
delta20.df <- RunSim(6000, .20, log(1.2), 0.5)
delta80.df <- RunSim(6000, .80, log(1.2), 0.5)

Os comentários no código explicam amplamente o que está acontecendo. As linhas aqui calculam a taxa de aceitação \[ R = \dfrac{p(e \, | \, \theta')p(\theta')}{p(e \, | \, \theta_i)p(\theta_i)} \]

e atualizam os valores para \(\mu\) e \(\sigma\) movendo-se apenas se a taxa de aceitação for alta o suficiente.

Observe que o tamanho do passo delta \(\delta\), o número de amostras para executar M e os valores iniciais são decisões puramente computacionais, eles não são determinados pela estrutura probabilística do problema, portanto, acima os tornamos entradas funcionais.

As últimas linhas executam três simulações, cada uma com 6.000 amostras. Os três usam as mesmas entradas escolhidas arbitrariamente; a diferença entre os três é o tamanho do passo \(\delta\). Examinando os três arquivos de saída, descobrimos que a taxa de aceitação quando \(\delta= 0.005\) é de cerca de 95%; quando \(\delta= 0.10\) é 26%, e quando \(\delta= 0.80\), é apenas 2.4%. Os gráficos a seguir mostram as primeiras amostras de cada simulação e ilustram o fenômeno. Estes são às vezes chamados de traceplots.

Eles podem ser obtidos usando 

library(coda)
traceplot(mcmc(delta20.df[1001:6000, ]))


por exemplo, ou 

lr <- function(m) exp(m[,1]+.5*m[,2]^2)
lrsd <- function(m) sqrt((exp(m[,2]^2)-1)*exp(2*m[,1]+m[,2]^2))
plot(100*lr(delta005.df[1001:6000,]), type="l", ylim=c(62,145), 
     ylab="Mean of Loss Ratio (in %)")
abline(h=100,col="grey",lty=2)

plot(100*lrsd(delta005.df[1001:6000,]),type="l",ylim=c(0,75),
     ylab="Std. Dev. of Loss Ratio (in %)")


O objetivo de nossa simulação é rastejar pelo espaço de probabilidade de modo que alcancemos todas as partes importantes do espaço e que o tempo que passamos lá seja proporcional a \(p(e \, | \, \theta)p(\theta)\). Quando \(\delta= 0.005\), estamos nos movendo muito lentamente pelo espaço e levamos muito tempo, ou seja, muitos ciclos de processador para viajar. Quando \(\delta= 0.8\), nossos saltos são muito grandes e continuam sendo rejeitados, então novamente não viajamos com eficiência. \(\delta= 0.2\) atinge uma taxa de aceitação de 26%, o que parece razoável. De fato, os resultados teóricos sugerem que as taxas de aceitação na faixa de 20%-50% são ideais. Diz-se que uma simulação MCMC que está amostrando eficientemente a partir da distribuição a posteriori tem boa mistura.

Outro detalhe a ser considerado antes de usar a saída é o valor inicial de \(mu\) e \(\sigma\). Idealmente, os valores iniciais teriam alta probabilidade a posteriori, mas normalmente essa informação não está disponível, pois aprender a probabilidade a posteriori é o ponto principal do cálculo do MCMC. Em vez disso, os valores iniciais geralmente têm probabilidade muito baixa, como no caso deste exemplo. Se incluíssemos as amostras iniciais, as probabilidades iniciais teriam muito peso. Podemos corrigir isso excluindo a primeira parte da cadeia. Esses resultados descartados geralmente são chamados de amostras de burn-in (amostras queimadas). Nos gráficos, podemos ver que a simulação com \(\delta= 0.2\), foi queimada após 1.000 amostras, portanto, usaremos apenas as amostras 1.001-6.000 em nossa análise.

Um detalhe final para se preocupar ao usar técnicas MCMC é a convergência. Como estamos usando uma técnica de Monte Carlo, quaisquer respostas fornecidas serão apenas assintoticamente corretas e o erro pode ser muito grande se não usarmos amostras suficientes. Como as amostras MCMC serão autocorrelacionadas, uma abordagem é usar uma série temporal simples para os resultados da simulação e estimar o erro usando resultados comprovados sobre séries temporais. Isso pode ser feito usando o pacote coda no R.

Abaixo estão os resultados da coda. Observe que a estimativa de erro de série temporal é quase três vezes maior do que a estimativa de erro ingênua, devido à autocorrelação. São possíveis proporções muito mais altas. É importante ter isso em mente. Não presuma apenas que seu erro é baixo porque você tem 100.000 amostras, porque elas podem ser equivalentes a 100 amostras verdadeiramente independentes! 

summary(mcmc(delta20.df[1001:6000, ]))
## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##         Mean      SD  Naive SE Time-series SE
## mu    -0.245 0.08387 0.0011861       0.003176
## sigma  0.222 0.05754 0.0008138       0.002783
## 
## 2. Quantiles for each variable:
## 
##          2.5%     25%     50%     75%    97.5%
## mu    -0.4160 -0.3013 -0.2452 -0.1848 -0.08963
## sigma  0.1354  0.1801  0.2116  0.2527  0.35980

Embora mixagem, burn-in sejam considerações importantes de convergência ao executar análises MCMC, elas ainda são apenas problemas práticos. Em teoria, todas essas simulações forneceram amostras corretas da distribuição a posteriori sem burn-in, mesmo se começarmos com uma baixa, ou seja, diferente de zero região de probabilidade e mesmo que nossa mixagem seja lenta, porém, na prática, um problema com qualquer um desses problemas pode tornar os resultados do MCMC inutilizáveis.

Além de usar as estatísticas resumidas produzidas pelo coda, podemos usar as amostras de \(\mu\) e \(\sigma\) diretamente. Por exemplo, se quiséssemos simular índices de perda futuros de uma forma que refletisse tanto o risco de processo quanto o risco de parâmetro, poderíamos fazer isso primeiro por amostragem de \(\mu\) e \(\sigma\) a partir de delta20.df e, em seguida, extraindo de uma distribuição lognormal com esses parâmetros.


4. Regressão Bayesiana


Para ilustrar a aplicação potencial da metodologia Bayesiana, vamos descrever a regressão bayesiana começando com o modelo linear e discutindo, brevemente, os modelos lineares generalizados.


4.1. Modelo Linear de uma perspectiva Bayesiana


No modelo linear gaussiano padrão, assumimos que \[ Y=\pmb{X}\beta+\epsilon, \] onde \(\epsilon\) é um ruído gaussiano i.i.d., centrado, com variância \(\sigma^2\) e \(\pmb{X} = (1,X_1,\cdots,X_k)\).

A verossimilhança deste modelo, dada uma amostra \(\{(Y_i,\pmb{X}_i\}\), para \(i = 1,\cdots,n\) é então \[ \mathcal{L}(\beta,\sigma^2 \, | \, Y,\pmb{X})=\left(\dfrac{1}{\sqrt{2\pi\sigma^2}}\right)^n\exp\left( -\dfrac{1}{2\sigma^2} (Y-\pmb{X}\beta)^\top (Y-\pmb{X}\beta)\right)\cdot \]

Conforme sabido, o estimador de máxima verossimilhança de \(\beta\) é \[ \widehat{\beta} = (\pmb{X}^\top\pmb{X})^{-1}\pmb{X}Y \]

e um estimador imparcial do parâmetro de variância \(\sigma^2\) é

\[ \widehat{\sigma}^2=\dfrac{(Y-\pmb{X}\beta)^\top (Y-\pmb{X}\beta)}{n-k-1}=\dfrac{s^2}{n-k-1}\cdot \]

Observe que a probabilidade pode ser escrita como \[ \mathcal{L}(\beta,\sigma^2 \, | \, Y,\pmb{X})=\left(\dfrac{1}{\sqrt{2\pi\sigma^2}}\right)^n\exp\left( -\dfrac{1}{2\sigma^2} (Y-\pmb{X}\beta)^\top (Y-\pmb{X}\beta)-\dfrac{1}{2\sigma^2}(\beta-\widehat{\beta})^\top (\pmb{X}^\top\pmb{X})^{-1}(\beta-\widehat{\beta}) \right), \] ou \[ \mathcal{L}(\beta,\sigma^2 \, | \, Y,\pmb{X})=\pi(\sigma^2)\times \pi(\beta \, | \, \sigma^2), \]

onde a primeira distribuição é uma distribuição Gama inversa e a segunda é uma normal multivariada com matriz de covariância proporcional a \(\sigma^2\).

Para facilidade computacional, a priori conjugada pode ser construída com essas famílias. Mais especificamente, suponha que as distribuições a priori sejam \[ \beta \, | \, \sigma^2 \sim N(\beta_0 , \sigma^2 M_0^{-1}), \]

para alguma matriz simétrica \((k + 1)\times (k + 1)\) definida positiva \(M_0\), enquanto \[ \sigma^2 \sim IG(a_1,b_0)\cdot \]

Nesse caso, a distribuição a posteriori condicional de \(\beta\) dados \(\sigma^2\) e \(Y\) é \[ \beta \, | \, \sigma^2, Y \sim N\Big( \big(M_0+(\pmb{X}^\top\pmb{X})\big)^{-1}\big(M_0\beta_0+(\pmb{X}^\top\pmb{X})\widehat{\beta} \big) \, , \, \sigma^2\big(M_0+(\pmb{X}^\top\pmb{X})\big)^{-1}\Big), \] enquanto a distribuição a posteriori para \(\sigma^2\) é \[ \sigma^2 \, | \, Y \sim IG\left( \dfrac{n}{2}a_0 \, , \, b_0+\dfrac{(\beta_0-\widehat{\beta})^\top\big(M_0+(\pmb{X}^\top\pmb{X})\big)^{-1}(\beta_0-\widehat{\beta})}{2}\right)\cdot \]

Com esta estrutura, as distribuições a priori e a posteriori estão na mesma família. Cálculos mostram que \[ \mbox{E}(\beta \, | \, Y,\pmb{X}) = \big(M_0+(\pmb{X}^\top\pmb{X})\big)^{-1}\big(M_0\beta_0+(\pmb{X}^\top\pmb{X})\widehat{\beta} \big) \] e \[ \mbox{E}(\sigma^2 \, | \, Y,\pmb{X})=\dfrac{2b_0+s^2+(\beta_0-\widehat{\beta})^\top\big(M_0^{-1}+(\pmb{X}^\top\pmb{X})^{-1}\big)^{-1}(\beta_0-\widehat{\beta})}{n+2(a_0-1)}\cdot \]

O interessante é que, como temos distribuições a posteriori explícitas, é possível derivar intervalos de credibilidade para parâmetros, bem como previsões, \(\widehat{Y}\) para um dado \(x\). Se usarmos o Teorema Bayesiano do Limite Centra e \[ \mbox{E}(\widehat{Y} \, | \, Y,\pmb{X}) = x\big(M_0+(\pmb{X}^\top\pmb{X})\big)^{-1}\big(M_0\beta_0+(\pmb{X}^\top\pmb{X})\widehat{\beta} \big) \] e \[ \mbox{Var}(\widehat{Y} \, | \, Y,\pmb{X}) = \sigma^2 \left(1+x\big(M_0+(\pmb{X}^\top\pmb{X})\big)^{-1}x^\top \right)\cdot \]

obtemos limites para um intervalo de credibilidade aproximado.

Para ilustrar esse modelo, considere o seguinte exemplo. 

# install.packages("CASdatasets", repos = "http://cas.uqam.ca/pub/", type="source")
library(CASdatasets)
data("Davis")
head(Davis)
##   sex weight height reportedWeight reportedHeight
## 1   M     77    182             77            180
## 2   F     58    161             51            159
## 3   F     53    161             54            158
## 4   M     68    177             70            175
## 5   F     59    157             59            155
## 6   M     76    170             76            165


A arquivo de dados contém 200 observações e 5 variáveis: sex um fator: M para o sexo masculino e F para o sexo feminino; weight peso, um numérico para o peso em Kg; height altura, um numérico para a altura em cm; reportedWeight peso reportado, um numérico para o peso em Kg e reportedHeight altura relatada, um numérico para a altura em cm. 

Y <- Davis$height
X1 <- Davis$sex
X2 <- Davis$weight
df <- data.frame(Y,X1,X2)
X <- cbind(rep(1,length(Y)),(X1=="M"),X2)
beta0 <- c(150,7,1/3)
beta <- as.numeric(lm(Y~X1+X2,data=df)$coefficients)
sigma2 <- summary(lm(Y~X1+X2,data=df))$sigma^2
M0 <- diag(c(100,100,1000^2))
Xi <- c(1,1,70)


então a previsão é 

(m <- Xi%*%solve(M0+t(X)%*%X) %*%(( t(X)%*%X)%*%beta +M0%*%beta0))
##          [,1]
## [1,] 176.7536


e um intervalo de credibilidade de 95% baseado na aproximação gaussiana é 

s2 <- sigma2*(1+ (t(Xi)%*% (solve(M0+t(X)%*%X)))%*%(Xi))
qnorm(c(.025,.975),mean=m,sd=sqrt(s2))
## [1] 166.6124 186.8947


Esses valores podem ser comparados com o limite inferior e superior obtidos usando a função predict(), usando uma abordagem estatística frequêntista. 

predict(lm(Y~X1+X2,data=df),newdata=data.frame(X1="M",X2=70),interval="prediction")
##       fit      lwr      upr
## 1 176.057 165.8188 186.2951

4.2. Extensão para Modelos Lineares Generalizados


Em Modelos Lineares Generalizados (McCullagh & Nelder, 1989), duas funções devem ser especificadas:

  1. A distribuição condicional de \(Y\) dado \(pmb{X}\), que pertence à família exponencial e seja \(f\) denota a densidade associada ou função de probabilidade, para variáveis discretas.

  2. Uma função de ligação \(g\) que relaciona a média condicional \(\mu(\pmb{x}) = \mbox{E}(Y \, | \, \pmb{X}=\pmb{x})\) e as covariáveis, como \(g\big(\mu(\pmb{x})\big) = \pmb{x}^\top\beta\).

Esta função de ligação deve ser bijetiva, por motivos de identificação e, então, \(\mu(\pmb{x})=g^{-1}(\pmb{x}^\top\beta)\). Dois modelos populares são o logit ou probit e os modelos de Poisson ou log-lineares. No modelo logit, \(Y\) é uma resposta binária, \(Y\, | \, \pmb{X}\sim B(\pi(\pmb{X)})\), onde a função de ligação é a transformação logit, \(g(x) = log\big(x/(1-x)\big)\). Uma alternativa é o modelo probit, onde a função de ligação é a função quantílica da distribuição \(N(0,1)\), \(g(x) = \Phi^{-1}(x)\). Para o modelo de Poisson, \(Y\) é uma resposta de contagem, com uma distribuição de Poisson \(Y \, | \, \pmb{X}\sim P(\lambda(\pmb{X}))\), onde a função de ligação é a transformação logarítmica, \(g(x) = \log(x)\).

O amostrador Probit Hastings-Metropolis, no contexto de um modelo Probit pode ser o seguinte: Para inicializar, comece com estimadores de máxima verossimilhança, \(\beta_{(0)} = \widehat{\beta}\), com matriz de função de variância assintótica \(\widehat{\Sigma}\). No estágio \(k\),

Uma ideia natural é usar a a priori, de modo que a distribuição a posteriori seja \[ \pi(\beta \, | \, Y) \propto \prod_{i=1}^n \Phi\big( \pmb{X}_i^\top \beta\big)^{Y_i}\Big(1-\Phi\big( \pmb{X}_i^\top \beta\big) \Big)^{1-Y_i}\cdot \]

Considere o conjunto de dados do advogado (Frees, 2009). Consideraremos aqui apenas cinquenta observações, sem valores ausentes: 

autobi = read.csv(file = "http://leg.ufpr.br/~lucambio/CM/AUTOBIsim.csv", header = TRUE)
head(autobi)
##   X ATTORNEY CLMSEX MARITAL CLMINSUR SEATBELT CLMAGE LOSS
## 1 1        1      2       2        2        1    1.8 6.38
## 2 2        2      2       1        2        1    5.6 0.95
## 3 3        1      1       2        2        1    1.4 1.88
## 4 4        2      1       1        2        1    5.4 0.11
## 5 5        1      2       2        2        1    0.9 6.95
## 6 6        2      2       2        2        1    2.3 0.15


index <- which((is.na(autobi$CLMAGE)==FALSE) & (is.na(autobi$CLMSEX)==FALSE))
set.seed(123)
index <- sample(index,size=50)


A variável de interesse é se o segurado tem advogado ou não: 

Y <- (autobi$ATTORNEY[index]-1)


Considere a seguinte regressão, onde as covariáveis são o custo do sinistro, o sexo do segurado e a idade do segurado. A densidade a posterior está aqui: 

X <- cbind(rep(1,nrow(autobi)), autobi$LOSS,autobi$CLMSEX==2,autobi$CLMAGE)
X <- X[index,]


A densidade a posteriori está aqui: 

posterior <- function(b) prod(dbinom(Y,size=1,prob=pnorm(X%*%b)))


O algoritmo Metropolis-Hasting será: 

model <- summary(glm(Y~0+X,family=binomial(link="probit")))
beta <- model$coeff[,1]
vbeta <- NULL
Sigma <- model$cov.unscaled
tau <- 1
nb.sim <- 10000
library(mnormt)
for(s in 1:nb.sim){
 sim.beta <- as.vector(rmnorm(n=1,mean=beta,varcov=tau^2*Sigma))
 prb <- min(1,posterior(sim.beta)/posterior(beta))
 beta <- cbind(beta,sim.beta)[,sample(1:2,size=1,prob=c(1-prb,prb))]
 vbeta=cbind(vbeta,beta)
 }


Para visualizar as sequências geradas, vamos retirar os primeiros 10% e mostrar a sequência de \(\beta_2\), parâmetro associado ao custo do sinistro: 

mcmc.beta <- vbeta[,seq(nb.sim/10,nb.sim)]
plot(mcmc.beta[2,],type="l",col="grey")
abline(h=model$coeff[2,1],lty=2)
grid()


Também podemos visualizar a distribuição do \(\beta_2\): 

hist(mcmc.beta[2,],proba=TRUE,col='grey',border='white')
vb <- seq(min(mcmc.beta[2,]),max(mcmc.beta[2,]),length=101)
lines(vb,dnorm(vb,model$coeff[2,1],model$coeff[2,2]),lwd=2)
lines(density(mcmc.beta[2,]),lty=2)

A linha plana é a distribuição gaussiana assintótica padrão, nos GLMs, a distribuição dos parâmetros é apenas assintótica e a linha pontilhada, o estimador baseado em kernel.


5. Exercícios