A Inferência Variacional é uma abordagem para se estimar distribuições a posteriori muito difíceis ou complexas.
Frequentemente, em Inferência Bayesiana estamos interessados em computar as medidas $p(z|x)$ para variáveis latentes $z_{1},\dots,z_{n}$ dado as observações $x_{1},\dots,x_{n}$. Quando essa distribuição é intratável podemos aproximar $p(z|x)$ por $q(z)$ onde a medida de proximidade é dada pela divergência de Kullback-Leibler.
Usualmente a Inferência Variacional é realizada através dos seguintes passos:
- Comece com um modelo para $p(z,x)$.
- Encontre a aproximação variacional $q(z|\nu)$ onde $\nu$ são os parâmetros variacionais.
- Encontre a Evidence Lower Bound (ELBO) dada por $L(\nu)=E_{Q}[\log(p(z,x))-\log(q(z|\nu))]$
- Minimize a quantia $L(\nu)$
- Volte para o passo 3 e faça isso até haver convergência.
Com base no exposto, vamos incialmente considerar um modelo de Regressão Probit construído da seguinte forma:
#Define a semente set.seed(8333) #Gera os dados x<-rnorm(1000) #Gera a latente z<- rnorm(1000,10+5*x,1) #Gera a variável dependente: y<-ifelse(z<0,0,1)O objetivo é estmar os parâmetros $a=10,b=5$ do modelo de Regressão Probit. Note que o Processo Gerador de Dados temos que a distribuição a posteriori é dada por: $p(a,b,z|y,x)= \prod_{i=1}^{n}y_{i}^{I(z_{i}>0)}(1-y_{i})^{1-I(z_{i}>0)}\times \frac{1}{\sqrt{2\pi}}\exp\left[-\frac{(z_{i}-(a+bx_{i}))^{2}}{2}\right]$ o qual também pode ser escrito como: $p(a,b,z|y,x)\propto \prod_{i=1}^{n}y_{i}^{I(z_{i}>0)}(1-y_{i})^{I(z_{i}\leq 0)}\times \exp\left[-\frac{(z_{i}-(a+bx_{i}))^{2}}{2}\right]$ para distribuições a priori uniformes para $a,b$. Dado o modelo precisamos definir a aproximação variacional. Nesse caso Blei et. al. (2017) sugere como regra ótima para a aproximação variacional fazer: $q(\theta_{j})\propto \exp\left[E_{\theta_{-j}}(\log(p(\theta|y,x)))\right]$ onde $\theta=(\theta_{1},\dots,\theta_{p})$ são os parâmetros e variáveis latentes do modelo. Nesse sentido, temos: $\log(p(a,b,z|y,x))=\sum_{i=1}^{n}I(z_{i}>0)\log(y_{i})+I(z_{i}\leq 0)\log(1-y_{i})-\sum_{i=1}^{n}\frac{(z_{i}-(a+bx_{i}))^{2}}{2}$ A distribuição varicional ótima é dada por: $q_{j}(z_{j})\propto \exp\left[E_{a,b,z_{-j}}(\log(p(a,b,z|y,x)))\right]$ Note que: $E_{a,b,z_{-j}}(\log(p(a,b,z|y,x)))=I(z_{j}>0)\log(y_{j})+I(z_{j}\leq 0)\log(1-y_{j})-\frac{E_{a,b}(z_{j}-(a+bx_{j}))^{2}}{2}$ Aplicando a exponencial temos que a distribuição ótima para $z_{j}$ é uma distribuição Normal Truncada, isso é $q_{j}(z_{j})\sim N^{+}(E(a)+E(b)x_{j},1)$ se $y_{j}=1$ e $q_{j}(z_{j})\sim N^{-}(E(a)+E(b)x_{j},1)$ se $y_{j}=0$. Similarmente, para $a$ temos $E_{b,z}[\log(p(a,b,z|y,x))]\propto E_{b,z}\left(-\frac{\sum_{i=1}^{n}(z_{i}-(a+bx_{i}))^{2}}{2}\right)$ , removendo os termos que não dependem de $a$ e complentando quadrados temos: $q_{a}(a)\sim N(\frac{\sum_{i=1}^{n}[E(z_{i})-E(b)x_{i}]}{n},\frac{1}{n})$. Finalmente para $b$ temos $E_{a,z}[\log(p(a,b,z|y,x))]\propto E_{a,z}\left(-\frac{\sum_{i=1}^{n}(z_{i}-(a+bx_{i}))^{2}}{2}\right)$, novamente removendo os termos que não dependem de $b$ e complentando quadrados temos: $q_{b}(b)\sim N(\frac{\sum_{i=1}^{n}[E(z_{i})-E(a)]x_{i}}{\sum_{i=1}^{n}x_{i}^{2}},\frac{1}{\sum_{i=1}^{n}x_{i}^{2}})$. Nesse caso, a otimização ocorre atribuindo os valores dos parâmetros a sua respectiva média da distribuição variacional. As equações para atualização dessas distribuições são dadas por:
#Equações de Update
update.newZj = function(newA,newB,j) {
mu = newA + newB*x[j]
if (y[j] == 1) {
return(mu + dnorm(-1*mu)/(1-pnorm(-1*mu)))
} else {
return(mu - dnorm(-1*mu)/(pnorm(-1*mu)))
}
}
update.newA = function(newZ,newB) {
return(sum(newZ-newB*x)/n)
}
update.newB = function(newZ,newA) {
return(sum(x*(newZ-newA))/sum(x^2))
}
Finalmente, podemos construir um sistema iterativo baseado nos valores otimizados:
#Total de observações
n<-length(z)
#Chute inicial
newA<-3
newB<-9
#Inicializa vetores
a.vec<-rep(NA,n)
b.vec<-rep(NA,n)
newZ<-rep(NA,n)
#Inicia o processo iterativo
for(it in 1:1000){
#Update Z
for (i in 1:n) {
newZ[i] = update.newZj(newA,newB,i)
}
#Correct Numerical
newZ[is.infinite(newZ)]<-100
#Update A
newA = update.newA(newZ,newB)
newB = update.newB(newZ,newA)
a.vec[it] = newA
b.vec[it] = newB
}
plot(a.vec,type="l")
summary(a.vec)
plot(b.vec,type="l")
summary(b.vec)
