Introdução à Inferência Variacional.

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:

1. Comece com um modelo para $p(z,x)$.
2. Encontre a aproximação variacional $q(z|\nu)$
onde $\nu$ são os parâmetros variacionais.
3. Encontre a Evidence Lower Bound (ELBO) dada por $L(\nu)=E_{Q}[\log(p(z,x))-\log(q(z|\nu))]$
4. Minimize a quantia $L(\nu)$
5. 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)