Last updated: 2018-09-05
workflowr checks: (Click a bullet for more information) ✔ R Markdown file: up-to-date
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
✔ Environment: empty
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
✔ Seed:
set.seed(20180529)
The command set.seed(20180529) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
✔ Session information: recorded
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
✔ Repository version: ecdacb8
wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.Rhistory
Ignored: docs/.DS_Store
Ignored: docs/figure/
Unstaged changes:
Modified: analysis/_site.yml
Variational Bayes is a particular variational method which aims to find some approximate joint distribution \(q(\beta; \theta)\) over hidden variables \(\beta\) to approximate the true joint \(p(\beta)\), and define the ‘closeness’ as the KL divergence \(KL\left[q(\beta; \theta)|| p(\beta)\right]\). The mean field form of VB assumes that q factorises into single-variable factors, \(q(\beta) = \prod_{i} q_{i}(\beta_{i}|\theta_{i})\).
We wish to find a set of distributions \(\{q_{i}(\beta_{i};\theta_{i})\}\) to minimize the KL divergence: \[ KL\left[ q(\beta) \| p(\beta|D) \right] = \int q(\beta) \log \frac{q(\beta)}{p(\beta|D)} d\beta \]
\[ q(\beta) = \prod_{i} q_{i}(\beta_{i}|\theta_{i}) \] subject to the normalization constraints: \[ \forall i \int q(\beta_{i}) d\beta_{i} = 1 \]
Note: q approximates the joint, but the individual \(q_{i}(\beta_{i})\) are poor approximations to the true marginals \(p_{i}(\beta_{i})\).
\[ KL\left[q(\beta)\|p(\beta|D)\right] = \int q(\beta)\log \frac{q(\beta)}{p(\beta|D)} d\beta = -\int q(\beta)\log \frac{p(\beta|D)}{q(\beta)} d\beta \]
\[ p(\beta|D) = \frac{p(\beta, D)}{p(D)} \quad \rightarrow \log p(\beta|D) = \log p(\beta, D) - \log p(D) \]
\[ KL\left[q(\beta)\|p(\beta|D)\right] = -\int q(\beta)\log \frac{p(\beta,D)}{q(\beta)} d\beta + \log p(D) \] We can ignore the \(\log p(D)\) term, since it does not involve Q.
Define \[ L\left[q(\beta) \right] = \int q(\beta) \log \frac{p(\beta,D)}{q(\beta)} d\beta \]
So \[ KL\left[q(\beta)\|p(\beta|D)\right] = -L + \log p(D) \] To minimize the KL divergence, we must maximize L.
\(L(q(\beta))\) are lower bounds on the model log-likelihood, \(p(D) = p(D|M)\). \[ \log p(D) = L\left[q(\beta)\right] + KL\left[q(\beta)\|p(\beta|D)\right] \] When the KL-divergence is zero (a perfect fit), L is equal to the model log-likelihood. The KL-divergence is always positive.
\[ \begin{align*} L\left[q(\beta)\right] &= \int q(\beta) \log \frac{p(\beta,D)}{q(\beta)} d\beta \\ &= \int q(\beta) \log p(\beta,D) d\beta - \int q(\beta) \log q(\beta) d\beta \\ &= \left<E(\beta, D)\right>_{q(\beta)} + H\left[q(\beta)\right] \end{align*} \] where we define energy as \(E = \log p\) and entropy as \(H\left[q(\beta)\right] = - \int q(\beta) \log q(\beta) d\beta\).
By the mean-field assumption: \[ L\left[q(\beta)\right] = \int \left(\prod_{i}q_{i}(\beta_{i})\right) E(\beta,D) d\beta - \int \left(\prod_{i}q_{i}(\beta_{i})\right) \sum_{i} \log q_{i}(\beta_{i}) d\beta \]
Consider the partitions \(\beta = \{\beta_{i}, \bar{\beta}_{i}\}\), \(\bar{\beta}_{i} = \beta\setminus\beta_{i}\)
The entropy term: \[ \begin{align*} \int \left(\prod_{i}q_{i}(\beta_{i})\right) \sum_{i} \log q_{i}(\beta_{i}) d\beta &= \sum_{i} \int \left(\prod_{i}q_{i}(\beta_{i})\right) \log q_{i}(\beta_{i}) d\beta \\ &= \sum_{i} \int q_{i}(\beta_{i}) q(\bar{\beta}_{i}) \log q_{i}(\beta_{i}) d\beta_{i} d\bar{\beta}_{i} \\ &= \sum_{i} \int q_{i}(\beta_{i}) \log q_{i}(\beta_{i}) d\beta_{i} \end{align*} \]
The energy term: \[ \begin{align*} \int \left(\prod_{i}q_{i}(\beta_{i})\right) E(\beta,D) d\beta &= \int q_{i}(\beta_{i}) \int q(\bar{\beta}_{i}) E(\beta,D) d \bar{\beta}_{i} d\beta_{i} \\ &= \int q_{i}(\beta_{i}) \left<E(\beta,D)\right>_{q(\bar{\beta}_{i})} d\beta_{i} \\ &= \int q_{i}(\beta_{i}) \log \exp \left<E(\beta,D)\right>_{q(\bar{\beta}_{i})} d\beta_{i} \\ &= \int q_{i}(\beta_{i})\log q_{i}^{\star}(\beta_{i}) d\beta_{i} + \log Z \\ \end{align*} \] where we define \(q_{i}^{\star} = \frac{1}{Z} \exp\left<E(\beta,D)\right>_{q(\bar{\beta}_{i})}\), and Z normalized \(q_{i}^{\star}(\beta_{i})\).
So \[ \begin{align*} L\left[q(\beta)\right] &= \int q_{i}(\beta_{i})\log q_{i}^{\star}(\beta_{i}) d\beta_{i} - \sum_{i} \int q_{i}(\beta_{i}) \log q_{i}(\beta_{i}) d\beta_{i} + \log Z \\ &= \{\int q_{i}(\beta_{i})\log q_{i}^{\star}(\beta_{i}) d\beta_{i} - \int q_{i}(\beta_{i}) \log q_{i}(\beta_{i}) d\beta_{i}\} + H\left[q(\bar{\beta}_{i}) \right] + \log Z \\ &= \int q_{i}(\beta_{i}) \log \frac{q_{i}^{\star}(\beta_{i})}{q_{i}(\beta_{i})} d\beta_{i} + H\left[q(\bar{\beta}_{i}) \right] + \log Z \\ &= -KL\left[q_{i}(\beta_{i}) \| q_{i}^{\star}(\beta_{i}) \right] + H\left[q(\bar{\beta}_{i}) \right] + \log Z \end{align*} \] We convert the problem from minimizing KL divergence between large joint distributions (hard) to minimizing KL divergence between individual 1D distributions (easier).
We wish to maximise L with respect to each \(q_{i}\), subject to the constraint that all \(q_{i}\) are normalizes to unity. L will be maximized when the KL divergence is zero \[ q_{i}(\beta_{i}) = q_{i}^{\star}(\beta_{i}) \] The optimal \(q_{i}(\beta_{i})\) is \[ q_{i}(\beta_{i}) = \frac{1}{Z}\exp\left<E(\beta_{i}, \bar{\beta}_{i},D)\right>_{q(\bar{\beta}_{i})} = \frac{1}{Z}\exp \int \log p(\beta_{i}, \bar{\beta}_{i},D) q(\bar{\beta}_{i}) d \bar{\beta}_{i} \]
\[ \log q_{i}(\beta_{i}) \leftarrow \int \log p(\beta_{i}, \bar{\beta}_{i},D) q(\bar{\beta}_{i}) d \bar{\beta}_{i} \]
Suppose we have M observations $D = {D_{i}}_{i=1}^{M} $, \[ D_{i}|\mu, \gamma \sim N(\mu, \gamma^{-1}) \quad \mu \sim N(m,\beta^{-1}) \quad \gamma \sim Gamma(a,b) \] The task is to infer the posterior distribution.
Assuming the conjugate exponential priors, we have \[ q(x) = q(\mu)q(\gamma) = N(\mu; m, \beta^{-1}) \Gamma(\gamma; a,b) \]
\[ \begin{align*} \log q(\mu) &\leftarrow \int q(\gamma) \log p(\mu, \gamma,D) d \gamma \\ &= \int \Gamma(\gamma; a,b)\log N(\mu; m, \beta^{-1}) \Gamma(\gamma; a,b)\prod_{i}N(D_{i}|\mu,\gamma^{-1}) d \gamma \\ &= \int \Gamma(\gamma; a,b)\{\log N(\mu; m, \beta^{-1}) + \log \Gamma(\gamma; a,b) + \sum_{i} \log N(D_{i}|\mu,\gamma^{-1})\} d \gamma \\ &= \int \Gamma(\gamma; a,b)\log N(\mu; m, \beta^{-1})d \gamma + \int \Gamma(\gamma; a,b)\log \Gamma(\gamma; a,b) d \gamma + \sum_{i} \int \Gamma(\gamma; a,b) \log N(D_{i}|\mu,\gamma^{-1}) d \gamma \\ &= \log N(\mu; m, \beta^{-1}) + \sum_{i} \int \Gamma(\gamma; a,b) \log N(D_{i}|\mu,\gamma^{-1}) d \gamma + \log \frac{1}{Z}\\ \end{align*} \]
\[ \begin{align*} \int \Gamma(\gamma; a,b) \log N(D_{i}|\mu,\gamma^{-1}) d \gamma &= \int \Gamma(\gamma; a,b) \left[c + (D_{i}-\mu)\gamma(D_{i}-\mu) \right] d \gamma \\ &\propto (D_{i}-\mu)\left[\int \gamma \Gamma(\gamma; a,b) d\gamma\right] (D_{i}-\mu) \\ &= (D_{i}-\mu)\left<\gamma\right>(D_{i}-\mu) \\ \end{align*} \]
Note: \(\left<\gamma\right> = \left[\int \gamma \Gamma(\gamma; a,b) d\gamma\right] = \mathbb{E}(\Gamma(\gamma; a,b)) = ab^{-1}\)
\[ \begin{align*} \log q(\mu) &= \log N(\mu; m, \beta^{-1}) + \sum_{i} (D_{i}-\mu)\left<\gamma\right>(D_{i}-\mu) + \log \frac{1}{Z}\\ q(\mu) &= \frac{1}{Z}N(\mu; m, \beta^{-1})\exp{\sum_{i} (D_{i}-\mu)\left<\gamma\right>(D_{i}-\mu)} \\ &= \frac{1}{Z}N(\mu; m, \beta^{-1})\prod_{i}N(\mu | D_{i}, \left<\gamma\right>^{-1}) \end{align*} \]
Therefore, using conjugate prior knowledge \[ q(\mu) = N(\mu| m', \beta'^{-1}) \quad \beta' = \beta + M\left<\gamma\right> \quad m' = \beta'^{-1}(\beta m + \left<\gamma\right>\sum_{i=1}^{M} D_{i}) \]
\[ \mathbf{y} = X_{n\times p}\mathbf{\beta}_{p\times 1} + \mathbf{\epsilon} \quad \mathbf{\epsilon} \sim N_{n}(0, \sigma^2 I) \]
\[ \beta_{j} \sim \pi_{j}(\beta_{j}) \] Only one of the \(\beta\) is non zero. The task is to find the posterior distribution of \(\beta\).
\[ q(\beta) = \prod_{j} q(\beta_{j}) \]
\[ \begin{align*} \log q(\beta_{j}) &\leftarrow \int q(\bar{\beta}_{j}) \log p(\mathbf{y}, \beta) d \bar{\beta}_{j} = \mathbb{E}_{q_{j'}, j'\neq j}\left[\log p(\mathbf{y}, \beta)\right]\\ &= \mathbb{E}_{q_{j'}, j'\neq j}\left[\log \left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^{n}\exp\left(-\frac{\|\mathbf{y} - X\beta\|^2}{2\sigma^2}\right)\prod_{j=1}^{p}\pi_{j}(\beta_{j}) \right]\\ &= \mathbb{E}_{q_{j'}, j'\neq j}\left[-\frac{\|\mathbf{y} - X\beta\|^2}{2\sigma^2} + \sum_{j=1}^{p}\log \pi_{j}(\beta_{j}) \right]\\ &= -\frac{\|\mathbf{y} - X_{-j}\bar{\beta}_{-j} - \mathbf{x}_{j}\beta_{j}\|^2}{2\sigma^2} + \log \pi_{j}(\beta_{j}) + C \\ &= -\frac{\|\mathbf{y}_{resid, -j} - \mathbf{x}_{j}\beta_{j}\|^2}{2\sigma^2} + \log \pi_{j}(\beta_{j}) + C \\ q(\beta_{j}) &\propto \exp \left(-\frac{\|\mathbf{y}_{resid, -j} - \mathbf{x}_{j}\beta_{j}\|^2}{2\sigma^2}\right)\pi_{j}(\beta_{j}) \end{align*} \]
\(\bar{\beta}_{-j}\) is the vector of estimated posterior mean other than j.
\[ \beta_{j} \sim \begin{cases} N(0, \sigma^2 \sigma_{\beta}^2) & \text{w.p } \ p \\ 0 & \text{w.p } \ 1-p \end{cases} \] \[ \pi_{j}(\beta_{j}) = p\frac{1}{\sqrt{2\pi\sigma^2\sigma_{\beta}^2}}\exp \left(-\frac{\beta_{j}^2}{2\sigma^2 \sigma_{\beta}^2}\right) + (1-p) \delta_{0} \] Therefore, \[ \beta_{j} \sim_{q_{j}} \begin{cases} N(\tilde{\mu}_{j}, \tilde{\sigma}_{j}^{2}) & \text{w.p } \alpha_{j} \\ 0 & \text{w.p } 1-\alpha_{j} \end{cases} \]
\[ \tilde{\sigma}_{j}^{2} = \frac{\sigma^2\sigma_{\beta}^2}{\mathbf{x}_{j}^{T}\mathbf{x}_{j}\sigma_{\beta}^2 + 1} = \left( \frac{\mathbf{x}_{j}^{T}\mathbf{x}_{j}}{\sigma^2} + \frac{1}{\sigma^2\sigma_{\beta}^2} \right)^{-1} \]
\[ \tilde{\mu}_{j} = \frac{\mathbf{x}_{j}^{T}\mathbf{y}_{res,-j}}{\mathbf{x}_{j}^{T}\mathbf{x}_{j} + 1/\sigma_{\beta}^2} = \frac{\tilde{\sigma}_{j}^{2}}{\sigma^2} \mathbf{x}_{j}^{T}\mathbf{y}_{res,-j} \]
\[ \frac{\alpha_{j}}{1-\alpha_{j}} = \frac{\pi}{1-\pi} \frac{\tilde{\sigma}_{j}}{\sigma_{\beta}\sigma} \exp(\frac{\tilde{\mu}_{j}^{2}}{2\tilde{\sigma}_{j}^{2}}) \quad \quad \alpha_{j} \propto \pi N(\frac{\mathbf{x}_{j}^{T}\mathbf{y}_{res,-j}}{\mathbf{x}_{j}^{T}\mathbf{x}_{j}}; 0, \frac{\sigma^2}{\mathbf{x}_{j}^{T}\mathbf{x}_{j}} + \sigma^2\sigma_{\beta}^2 ) \]
sessionInfo()R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] workflowr_1.1.1 Rcpp_0.12.18 digest_0.6.15
[4] rprojroot_1.3-2 R.methodsS3_1.7.1 backports_1.1.2
[7] git2r_0.23.0 magrittr_1.5 evaluate_0.11
[10] stringi_1.2.4 whisker_0.3-2 R.oo_1.22.0
[13] R.utils_2.6.0 rmarkdown_1.10 tools_3.5.1
[16] stringr_1.3.1 yaml_2.2.0 compiler_3.5.1
[19] htmltools_0.3.6 knitr_1.20 This reproducible R Markdown analysis was created with workflowr 1.1.1