Algorítmo de Gibbs: Modelo Poisson com prioris (conjugadas) em dois estágios

Modelo e expressões

Seja um modelo com estrutura hierárquica no qual assume-se a distribuição \(Poisson(\lambda)\) para as observações, priori \(Gamma(a, \beta)\) com \(a\) considerado aqui constante e hiperpriori \(GammaInversa(c, d)\) para \(\beta\). Neste modelo não é possível amostrar diretamente da posteriori e obter as marginais em forma fechada. Implementar o algoritmo de Gibbs.

\[\begin{align*} Y|\lambda &\sim {\rm Po}(\lambda) \\ \lambda|\beta &\sim {\rm Ga}(a, \beta) \\ \beta &\sim {\rm IGa}(c,d) \end{align*}\]

Neste modelo os parâmetros são \(\theta = (\lambda, \beta)^{\prime}\) e \((a, c, d)\) são constantes conhecidas/fixadas.

A função de probabilidade para uma observação \(y_i\) é \[f(y; \lambda) = \frac{\exp\{-\lambda\} \lambda^{y_i}}{y_i!}\] e a verossimilhança para as observações, consideradas condicionalmente independentes dado \(\lambda\): \[L(\lambda; y) = \frac{\exp\{-n\lambda\} \lambda^{\sum{y_i}}}{\prod y!}\] As prioris são: \[\begin{align*} f(\lambda|\beta) &= \frac{1}{\beta^a \Gamma(a)} \lambda^{a-1} \exp\{-\lambda/\beta\}\\ f(\beta) &= \frac{1}{d^c \Gamma(c)} \beta^{-c-1} \exp\{-1/d\beta\} \end{align*}\]

A posteriori é obtida por \[f(\lambda,\beta|y) = \frac{f(\lambda,\beta,y)}{f(y)} \propto f(\lambda,\beta,y) = f(y|\lambda,\beta)f(\lambda|\beta)f(\beta) = f(y|\lambda)f(\lambda|\beta)f(\beta) \] Na expressão anterior, a estrutura de independência condicional do modelo implica que \(f(y|\lambda,\beta) = f(y|\lambda)\).

Tem-se então que \[\begin{align*} f(\lambda,\beta|y) &= \frac{\exp\{-n\lambda\} \lambda^{\sum{y_i}}}{\prod y!} \frac{1}{\beta^a \Gamma(a)} \lambda^{a-1} \exp\{-\lambda/\beta\} \frac{1}{d^c \Gamma(c)} \beta^{-c-1} \exp\{-1/d\beta\} \\ &\propto \lambda^{a-1} \lambda^{\sum y_i} \beta^{-a}\beta^{-c-1} \exp\{-n\lambda\}\exp\{-\lambda/\beta\}\exp\{-1/d\beta\} \\ & = \lambda^{a+\sum y_i} \beta^{-(a+c)-1} \exp\{-n\lambda -\frac{\lambda}{\beta} -\frac{1}{d\beta}\} \end{align*}\] Neste exemplo é possível obter a posteriori \(f(\lambda,\beta|y)\) analiticamente que pode ser escrita como \[f(\lambda,\beta|y) = f(\lambda|\beta,y) f(\beta|y). \] As marginais são obtidas analiticamente integrando-se a conjunta \[\begin{align*} f(\lambda|y) &= \int f(\lambda,\beta|y) {\rm d} \beta \\ \mbox{e} & \\ f(\beta|y) &= \int f(\lambda,\beta|y) {\rm d} \lambda . \end{align*}\] Portanto a inferência baseada na posteriori é disponível e, pode-se amostrar diretamente desta distribuição com os passos: Entretanto, o objetivo aqui é ilustrar a implementação do algoritmo de Gibbs que requer apenas as condicionais completas, que são facilmente obtidas conforme a seguir. \[\begin{align*} f(\lambda|y, \beta) & \propto \lambda^{a+(\sum{y_i}) -1} \exp\{-n\lambda -\frac{\lambda}{\beta}\} = \lambda^{a+(\sum{y_i}) -1} \exp\{-\left(\frac{\beta}{n\beta+1}\right)^{-1} \lambda\} \\ [\lambda|y, \beta] &\sim {\rm G}(a+\sum y_i, \frac{\beta}{n\beta + 1}) \\ f(\beta|y, \lambda) &= \beta^{-(a+c)-1} \exp\{-\frac{\lambda}{\beta} -\frac{1}{d\beta}\} = \beta^{-(a+c)-1} \exp\{-\left(\frac{d}{d \lambda + 1}\right)^{-1}\frac{1}{\beta}\}\\ [\beta|y, \lambda] &\sim {\rm IG}(a+c, \frac{d}{d\lambda + 1}) \end{align*}\]

Códigos em R para o Exemplo

Inicialmente vamos definir duas funções auxiliares, uma para densidade da distribuição \(Gama-Inversa(a,b)\) e outra que retorna a esperança e variância desta distribuição.

dinvgamma <- function(x, shape, scale) dgamma(1/x, shape=shape, scale=scale)/x^2
##
EVIG <- function(a,b){
    E <- ifelse(a > 1, 1/(b * (a-1)), "não pode ser calculada para este valor de a")
    V  <- ifelse(a > 2, 1/((b^2) * ((a-1)^2) * (a-2)), "não pode ser calculada para este valor de a")
    return(c(E,V))
}

Consideramos neste exemplo \(a=4,5\), \(c=2,2\) e \(d=0,6\). Na próxima figura pode-se visualizar a priori para \(\beta\). Note que para \(\lambda\) o modelos define apenas a condicional \(\lambda|\beta\) mas pode-se obter a marginal. \[\begin{align*} f(\lambda) &= \int f(\lambda, \beta) {\rm d}\beta = \int f(\lambda|\beta) f(\beta) {\rm d}\beta. \end{align*}\]

par(mar=c(3,3,0,0), mgp=c(2,1,0))
curve(dinvgamma(x, shape=2.2, scale=0.6),
      from=0, to=5, xlab=expression(beta), ylab=expression(f(beta)))

A seguir será simulado um conjunto de dados com 50 observações. Para isto são inicialmente fixadas as constantes no modelo. Na sequencia simulam-se valores de \(\beta\) e para cada um deles simula-se um valor de \(\lambda\). Finalmente simula-se, respectivamente, uma observação da variável resposta \(Y\).

set.seed(2018)
ctes <- list(a=4.5, c=2.2, d=0.6, n=50)
with(ctes, EVIG(c, d))

## [1] 1.388889 9.645062

betas <- with(ctes, 1/rgamma(n, shape=c, scale=d))
c(mean(betas),var(betas))

## [1] 1.385987 1.247596

lambdas <- with(ctes, rgamma(n, shape=a, scale=betas))
(ctes$y <- rpois(ctes$n, lambda=lambdas))

##  [1]  5 12  4  1  4  3  0 11  1  5  4  3 14  3  4  7 10  5  3  5  9  2 13
## [24] 18  5  2 10  2  4 11  2  4  3  2 13  6 35  0  3  1 13  2  8  4  1  0
## [47]  0  7  4  4

A seguir temos estatísticas dscritivas do dados. Notamos que a média é bem inferior à variância, caracterizando uma sobredispersão. No gráfico dos dados (linhas escura) vê-se claramente qua a distribuição Poisson (linhas vermelhas) com parâmetro constante não se ajusta bem aos dados.

with(ctes, c(media=mean(y), var=var(y)))

##    media      var 
##  5.84000 35.97388

with(ctes, plot(prop.table(table(y)), type="h", ylab="P[Y=y]"))
with(ctes,lines((0:max(y))+0.1, dpois(0:max(y), lambda=mean(y)), type="h", col=2))

Vamos agora implementar um amostrador de Gibbs para obter amostras da posteriori. Serão tomadas 11.000 simulações.

ctes$sumY <- sum(ctes$y)
N <- 11000
beta.sam <- lambda.sam <- numeric(N) 
beta.sam[1] <- lambda.sam[1] <- 10 ## valores iniciais
{
for(i in 2:N){
    beta.sam[i] <- with(ctes, 1/rgamma(1, shape=a+c, scale=d/(d*lambda.sam[i-1]+1)))
    lambda.sam[i] <- with(ctes, rgamma(1, shape=a+sumY, scale=beta.sam[i]/(n*beta.sam[i]+1)))
}}

A seguir temos a visualização das cadeias para uma inspeção visual inicial.

par(mfrow=c(2,1), mar=c(2.7,2.7,1,1), mgp=c(1.8,0.8,1))
plot(beta.sam, type="l")
plot(lambda.sam, type="l")

Entretanto é preciso atentar para a necessidade de aquecimento (//burn-in//) da cadeia. Ou seja, amostras iniciais podem não ser represetnativas da distribuição estacionária. Serão agora descartadas 1.000 como aquecimento da cadeia (//burn-in//).

Além disto iremos também proceder um raleamento (//thining//) armazenando apenas uma a cada dez amostras uma vez que as cedeias geradas tendem a ter muita autocorrelação entre as amostras o que implica em redundância de informação e uso desnecessário de memória.

B <- 1000
beta.sam <- beta.sam[-(1:B)]
lambda.sam <- lambda.sam[-(1:B)]
beta.sam <- beta.sam[seq(1,N-B, by=10)]
lambda.sam <- lambda.sam[seq(1,N-B, by=10)]

OBS: os códigos acima são ineficientes pois não seria necessário armazenar as amostras de aquecimento nem as que seriam descartadas no raleamento. Deixamos desta forma simplesmente para ilustrar os procedimentos e refazemos o código sem usar tal armazenamento no final destas notas.

A seguir temos a visualização das cadeias para uma nova inspeção visual.

par(mfrow=c(2,1), mar=c(2.7,2.7,0,0), mgp=c(1.8,0.8,0))
plot(beta.sam, type="l"); lines(lowess(beta.sam), col=2)
plot(lambda.sam, type="l"); lines(lowess(lambda.sam), col=2)

Não se percebe nenhum problema obvio com as cadeias mas nota-se a posteriori de \(\beta\) parece ter um comportamento assimétrico. Vamos então examinar a cadeia para escala logarítmica.

par(mfrow=c(1,1), mar=c(2.7,2.7,1,1), mgp=c(1.8,0.8,0))
plot(log(beta.sam), type="l"); lines(lowess(log(beta.sam)), col=2)

E vamos também verificar a autocorrelação após o raleamento.

par(mfrow=c(1,2), mar=c(2.7,2.7,0,0), mgp=c(1.8,0.8,0))
acf(beta.sam, xlab=expression(beta))
acf(lambda.sam, xlab=expression(lambda))

Estas análiss iniciais não apontam problemas óbvios e parece que temos uma boa amostra. Entretanto uma análise mais cuidadosa deveria incluir testes de diagnóstico de convergência, possivelmente com a obtenção de não apenas uma cadeia. Não vamos abordar estes testes aqui e seguiremos com a análise a partir da amostra obtida.

A seguir visualiza-se as amostras das distribuições marginais. Os valores indicados nos eixos são os gerados na simulação dos dados.

par(mfrow=c(1,2), mar=c(2.7,2.7,0,0), mgp=c(1.8,0.8,0))
plot(density(beta.sam), main=""); rug(betas)
plot(density(lambda.sam), main=""); rug(lambdas)

## Warning in rug(lambdas): some values will be clipped

A seguir visualiza-se as amostras da conjunta com \(\beta\) na escala original e logarítmica.

par(mfrow=c(1,2), mar=c(2.7,2.7,0,0), mgp=c(1.8,0.8,0))
plot(beta.sam, lambda.sam, pch=19, cex=0.4); rug(betas, side=1); rug(lambdas, side=2)

## Warning in rug(lambdas, side = 2): some values will be clipped

plot(log(beta.sam), lambda.sam, pch=19, cex=0.4); rug(log(betas), side=1); rug(lambdas, side=2)

## Warning in rug(log(betas), side = 1): some values will be clipped

## Warning in rug(log(betas), side = 1): some values will be clipped

Resumos dos dados e das posterioris dos parâmetros.

summary(ctes$y)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    2.00    4.00    5.84    7.75   35.00

summary(betas)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2341  0.7331  1.1373  1.3860  1.5907  5.6039

summary(beta.sam)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4102  0.9008  1.1640  1.2989  1.5469  6.2081

summary(lambdas)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1716  2.7200  4.0948  5.8029  6.4230 40.6322

summary(lambda.sam)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.634   5.569   5.809   5.814   6.053   6.835

dbetapost <- function(x, a, c, d, soma, n, log=FALSE){
res <- lgamma(a+soma+1) + (a+soma+1)*(log(x)-log(n*x+1))-(a+c+1)*log(x) - 1/(d*x)
if(log) return(res) else return(exp(res))
}
dlambdapost <- function(x, a, c, d, soma, n, log=FALSE){
res <- lgamma(a+c) + (a+c)*(log(d)-log(d*x-1))+(s+soma)*log(x) - n*x
if(log) return(res) else return(exp(res))
}