# Bayesian Logistic Regression

Contact: mihaelacr@google.com

Aim of this colab:
Logistic regression with a twist! We will use variational inference to learn a distribution over the learned parameters.

Througout this work we will be focusing on a binary classification task, on 
a small dataset (the UCD breast cancer dataset). 

Check out the slides for this colab [here](https://github.com/eeml2019/PracticalSessions/tree/master/unsupervised).

## Logistic regression

Learn the weights which best classify the $(x, y)$ pairs in the dataset. 
The weights are learned via the following maximization problem:

\begin{equation}
\mathbb{E}_{p^*(x, y)} \log p(y|x, w) = \\
 \mathbb{E}_{p^*(x, y)} \left[y \log \sigma(w x + b) + (1 - y) \log(1- \sigma(w x + b)\right] 
\end{equation}

Here $\sigma$ denotes the sigmoid function.

## Bayesian logistic regression

Learn a distribution over the weights which best classify the $(x, y)$ pairs in the dataset.

We now want to learn the distributional parameters of $q(w)$ in order to maximize:

\begin{equation}
\mathbb{E}_{p^*(x, y)} \log p_w(y| x) = \\
\mathbb{E}_{p^*(x, y)} \log \int p(y|x, w) p(w) \delta w = \\
\mathbb{E}_{p^*(x, y)} \log \int p(y|x, w) p(w) \frac{q(w)}{q(w)} \delta w \ge \\
\mathbb{E}_{p^*(x, y)} \mathbb{E}_{q(w)} \log \left[p(y| x, w) \frac{p(w)}{q(w)} \right]= \\
\mathbb{E}_{q(w)} \mathbb{E}_{p^*(x, y)} \log p(y|x, w)  -  KL(q(w)||p(w)) = \\
 \mathbb{E}_{q(w)}\mathbb{E}_{p^*(x, y)} \left[y \log \sigma(w x + b) + (1 - y) \log(1- \sigma(w x + b)\right] -  KL(q(w)||p(w))
\end{equation}

In Bayesian Logistic Regression, we thus learn a distribution over parameters
$w$ which can explain the data, while staying close to a chosen prior $p(w)$.

Throughout this lab, we will use Gaussian Distributions for $p(w)$ and $q(w)$.

For more details, see 
[this paper](https://pdfs.semanticscholar.org/e407/ea7fda6d152d2186f4b5e27aa04ec2d32dcd.pdf) or 
[this tutorial](https://www.ece.rice.edu/~vc3/elec633/logistic.pdf).

### Gradient estimation

In Bayesian Logistic Regression, we aim to learn the parameters of the distribution $q(w)$. We will call these parameters $\theta$. In the case of a 
Gaussian distribution, these will be the mean and covariance matrix of a
multivariate distribution. Since we will use stochastic gradient descent to learn $\theta$, we need to be able to compute the gradients with respect to 
$\theta$ of our objective. Since our objective contains an expectation with respect to $q(w)$, this can be challenging. To illustrate this, from now on we will denote $q(w)$ as $q_{\theta}(w)$.

Specifically, we are interested in:
\begin{equation}
  \nabla_{\theta} \left[\mathbb{E}_{q_{\theta}(w)} \mathbb{E}_{p^*(x, y)} \log p(y|x, w)  -  KL(q_{\theta}(w)||p(w))\right]
\end{equation}

Since in we will be using Gaussian distributions for $q_{\theta}(w)$ and $p(w)$ 
$KL(q_{\theta}(w)||p(w))$ can be computed in closed form and we can rely on 
TensorFlow's automatic differentiation to compute 
$\nabla_{\theta} KL(q_{\theta}(w)||p(w))$.

Let's now turn attention to computing $\nabla_{\theta} \mathbb{E}_{q_{\theta}(w)} \mathbb{E}_{p^*(x, y)} \log p(y| x, w)$. Unlike the KL, the integral 
$\mathbb{E}_{q_{\theta}(w)} \mathbb{E}_{p^*(x, y)} \log p(y|x, w)$
is not tractable, and thus cannot be computed in closed form. Thus, we are interested in other approaches to compute the gradient.

### Reinforce gradient estimation

A first approach to compute an estimate of the gradient is to use the
REINFORCE gradient estimator:

\begin{equation}
  \nabla_{\theta} \mathbb{E}_{q_{\theta}(w)} \mathbb{E}_{p^*(x, y)} \log p(y| x, w)  = \\
  \nabla_{\theta} \int q_{\theta}(w) \mathbb{E}_{p^*(x, y)} \log p(y| x, w) \delta w = \\
  \int \nabla_{\theta} \left[q_{\theta}(w) \mathbb{E}_{p^*(x, y)} \log p(y| x, w) \right] \delta w = \\
 \int \nabla_{\theta} \left[q_{\theta}(w)\right] \mathbb{E}_{p^*(x, y)} \log p(y| x, w) \delta w = \\
 \int q_{\theta}(w) \nabla_{\theta} \left[\log q_{\theta}(w)\right] \mathbb{E}_{p^*(x, y)} \log p(y| x, w)  \delta w = \\
 \mathbb{E}_{q_{\theta}(w)} \left[ \nabla_{\theta} \log q_{\theta}(w) \mathbb{E}_{p^*(x, y)} \log p(y| x, w) \right]
\end{equation}

We can use samples from $q_{\theta}(w)$ to compute the Monte Carlo estimate of
the last integral, to get an unbiased estimator of the true gradient. The more 
samples we use to estimate the integral, the more accurate the gradient estimate
will be.


### Pathwise gradient estimation (reparametrization)

In the REINFORCE estimator, we did not use any knowledge of the variational
distribution $q_{\theta}(w)$.
For a Gaussian distribution, we can use the reparametrization trick to obtain
an unbiased gradient estimator of  $\nabla_{\theta} \mathbb{E}_{q_{\theta}(w)} \mathbb{E}_{p^*(x, y)} \log p(y| x, w)$.

To do so, we will use the fact that 
\begin{equation}
    z \sim N(\mu, \sigma), z = \mu + \epsilon \sigma, {\text {with }} \epsilon \sim N(0, 1)
\end{equation}

Namely:
\begin{equation}
  \nabla_{\theta} \mathbb{E}_{q_{\theta}(w)} \mathbb{E}_{p^*(x, y)} \log p(y| x, w)  = \\
  \nabla_{\theta} \mathbb{E}_{p(\epsilon)} \mathbb{E}_{p^*(x, y)} \log p(y| x, \mu + \epsilon \sigma)  = \\
  \mathbb{E}_{p(\epsilon)} \nabla_{\theta} \mathbb{E}_{p^*(x, y)} \log p(y| x, \mu + \epsilon \sigma)
\end{equation}

Note that we were able to move the gradient inside the integral since 
$p(\epsilon)$ does not depend on $\theta$.
To estimate the last integral, we can use a Monte Carlo estimate using samples
of $p(\epsilon)$.



```
import enum
import numpy as np

import math
from sklearn import datasets
from sklearn import linear_model
import collections

from matplotlib import pyplot as plt
import tensorflow as tf
import tensorflow_probability as tfp
import seaborn as sns

tfd = tfp.distributions

```


```
sns.set(rc={"lines.linewidth": 2.8}, font_scale=2)
sns.set_style("whitegrid")
```

## Get the dataset and visualize it


```
def get_data(normalize=True):
    data = datasets.load_breast_cancer()

    print(data.target_names)
    print(data.feature_names)

    features = np.array(data.data, dtype=np.float32)
    targets = np.array(data.target, dtype=np.float32)

    if normalize:
        # Note: the data dimensions have very different scales.
        # We normalize by the mean and scale of the entire dataset.
        features = features - features.mean(axis=0)
        features = features / features.std(axis=0)

    # Add a dimension of 1 (bias).
    features = np.concatenate((features, np.ones((features.shape[0], 1))), axis=1)
    features = np.array(features, dtype=np.float32)
  
    return features, targets, data.feature_names
```

### Dataset

This datasets has meaningful features: we are trying to classify based on these features whether someone has breast cancer or not. We can look at the
features we are using to classify the data below.



```
features, targets, feature_names = get_data()
```

    ['malignant' 'benign']
    ['mean radius' 'mean texture' 'mean perimeter' 'mean area'
     'mean smoothness' 'mean compactness' 'mean concavity'
     'mean concave points' 'mean symmetry' 'mean fractal dimension'
     'radius error' 'texture error' 'perimeter error' 'area error'
     'smoothness error' 'compactness error' 'concavity error'
     'concave points error' 'symmetry error' 'fractal dimension error'
     'worst radius' 'worst texture' 'worst perimeter' 'worst area'
     'worst smoothness' 'worst compactness' 'worst concavity'
     'worst concave points' 'worst symmetry' 'worst fractal dimension']


### Trivial baseline

How well can we do with a constant classifier?

If we look at the average label, we can see how many 1s are in the data, so we know that our classifier has to get more than that accuracy


```
np.mean(targets)
```




    0.6274165




```
features.shape
```




    (569, 31)




```
dataset_size = features.shape[0]
data_dims = features.shape[-1]
```

## Logistic regression baseline - how well can we do on this dataset?

Let's look at a simple, non Bayesian classification approach to get an idea of how well we can do on the dataset. 


```
lr = linear_model.LogisticRegression()
```


```
lr.fit(features, targets)
lr.score(features, targets)
```




    0.9876977152899824



## Utilities


```
def get_shape_list(tensor):
  return tensor.shape.as_list()
```


```
def get_sklearn_data_as_tensors(
      features, targets, batch_size, dataset_name='breast_cancer'):
  """Read sklearn datasets as tf.Tensors.

  Args:
    batch_size: Integer or None. If None, the entire dataset is used.
    dataset_name: A string, the name of the dataset.

  Returns:
    A tuple of size two containing two tensors of rank 2 `[B, F]`, the data
    features and targets.
  """
  dataset = tf.data.Dataset.from_tensor_slices((features, targets))
  if batch_size:
    # Shuffle, repeat, and batch the examples.
    batched_dataset = dataset.shuffle(1000).repeat().batch(batch_size)
  else:
    batch_size = features.shape[0]
    batched_dataset = dataset.repeat().batch(batch_size)

  iterator = batched_dataset.make_one_shot_iterator()
  batch_features, batch_targets = iterator.get_next()

  data_dim = features.shape[1]
  batch_features.set_shape([batch_size, data_dim])
  batch_targets.set_shape([batch_size])

  return batch_features, batch_targets
```

## Bayesian Logistic Regression - define the prior and the posterior

We will learn a distribution over parameters, so we will first define prior distribution, and then a posterior distribution over parameters which we 
will learn.

Both of these distibutions will be Gaussian.


```
def multi_normal(loc, log_scale):
  scale = tf.exp(log_scale)
  return tfd.MultivariateNormalDiag(loc=loc, scale_diag=scale)

def diagonal_gaussian_posterior(data_dims):
  mean = tf.Variable(
      tf.zeros(shape=(data_dims), dtype=tf.float32), name='mean')
  log_scale = tf.Variable(
      tf.zeros(shape=(data_dims), dtype=tf.float32),
      name='log_scale')
  learned_vars = [mean, log_scale]
  return multi_normal(loc=mean, log_scale=log_scale), learned_vars
```

## Bayesian Logistic Regression - the model

### Hypers


```
BATCH_SIZE = 500
NUM_POSTERIOR_SAMPLES = 100
```


```
prior = multi_normal(loc=tf.zeros(data_dims), log_scale=tf.zeros(data_dims))
posterior, learned_vars  = diagonal_gaussian_posterior(data_dims)
```


```
features, targets = get_sklearn_data_as_tensors(features, targets, BATCH_SIZE)
```


```
def _predictions(logits):
  return tf.to_float(logits >= 0)

def _accuracy(targets, predictions):
  # `predictions` have rank 2: batch_size, num_posterior samples.
  # We expand dims the targets to compute the accuracy per sample.
  targets = tf.expand_dims(targets, axis=1)
  return tf.reduce_mean(tf.to_float(tf.equal(targets, predictions)))

def linear_model(data, targets, posterior_samples):
  num_posterior_samples = tf.shape(posterior_samples)[0]
  logits = tf.matmul(data, posterior_samples, transpose_b=True)
  # Make targets [B, 1] to use broadcasting.
  targets = tf.expand_dims(targets, axis=1)
  targets = targets * tf.ones([1, num_posterior_samples])

  log_probs = - tf.nn.sigmoid_cross_entropy_with_logits(
      labels=targets, logits=logits)
  return log_probs, logits
```


```
def run_model(features, targets, prior, posterior, num_samples):
  posterior_samples = posterior.sample(num_samples)

  # Compute the log probs of the data under all the models we have sampled.
  log_probs, logits = linear_model(features, targets, posterior_samples)
  # Compute the KL.
  kl = posterior.kl_divergence(prior)

  # Sum over data.
  param_log_probs = tf.reduce_mean(log_probs, axis=0) * dataset_size
  param_log_probs.shape.assert_is_compatible_with([num_samples])

  return param_log_probs, kl, posterior_samples, logits 
```


```
param_log_probs, kl, posterior_samples, logits = run_model(
      features, targets, prior, posterior, NUM_POSTERIOR_SAMPLES)
predictions = _predictions(logits)
accuracy = _accuracy(targets=targets, predictions=predictions)

accuracy = tf.reduce_mean(accuracy)
per_sample_elbo = param_log_probs - kl
elbo = tf.reduce_mean(per_sample_elbo)
```

## Training loop

## Training the model via reparametrization


```
# Stochastic loss depends on the output of a sampling operation  (posterior samples)
# but tensorflow implements reparametrization by default.
per_sample_reparametrization_loss = - per_sample_elbo
reparametrization_loss = tf.reduce_mean(per_sample_reparametrization_loss)
```


```
optimizer = tf.train.GradientDescentOptimizer(0.0001)
reparam_min_op = optimizer.minimize(reparametrization_loss)
```


```
NUM_ITERATIONS = 1000
```


```
sess = tf.Session()

# Initialize all variables
sess.run(tf.initialize_all_variables())

reparam_accuracies = []
reparam_kls = []
reparam_elbos = []

for i in xrange(NUM_ITERATIONS):
  sess.run(reparam_min_op)

  if i % 10 == 0:
    reparam_acc, reparam_kl, reparam_elbo = sess.run([accuracy, kl, elbo])
    reparam_accuracies += [reparam_acc]
    reparam_kls += [reparam_kl]
    reparam_elbos += [reparam_elbo]

    print('Iteration {}. Elbo {}. KL {}'.format(
        i, reparam_elbo, reparam_kl))

reparam_learned_mean, reparam_learned_log_scale = sess.run(learned_vars)
```

    Iteration 0. Elbo -1206.31604004. KL 0.00301633775234
    Iteration 10. Elbo -642.610412598. KL 0.293668031693
    Iteration 20. Elbo -577.772949219. KL 0.783629059792
    Iteration 30. Elbo -392.721923828. KL 1.2917201519
    Iteration 40. Elbo -351.349578857. KL 1.781883955
    Iteration 50. Elbo -306.663360596. KL 2.20695495605
    Iteration 60. Elbo -228.974685669. KL 2.57792901993
    Iteration 70. Elbo -234.82182312. KL 2.95425152779
    Iteration 80. Elbo -228.548690796. KL 3.29743432999
    Iteration 90. Elbo -196.769805908. KL 3.61361837387
    Iteration 100. Elbo -190.871459961. KL 3.9075319767
    Iteration 110. Elbo -220.950912476. KL 4.18836021423
    Iteration 120. Elbo -181.642272949. KL 4.45673751831
    Iteration 130. Elbo -163.075805664. KL 4.70559310913
    Iteration 140. Elbo -166.205795288. KL 4.96141529083
    Iteration 150. Elbo -161.878341675. KL 5.22007894516
    Iteration 160. Elbo -165.926605225. KL 5.44898462296
    Iteration 170. Elbo -158.645767212. KL 5.66910552979
    Iteration 180. Elbo -164.979400635. KL 5.88050985336
    Iteration 190. Elbo -144.492584229. KL 6.12092161179
    Iteration 200. Elbo -151.752105713. KL 6.33854484558
    Iteration 210. Elbo -165.249649048. KL 6.54486370087
    Iteration 220. Elbo -127.59765625. KL 6.73881149292
    Iteration 230. Elbo -119.381202698. KL 6.92593240738
    Iteration 240. Elbo -111.081077576. KL 7.11052560806
    Iteration 250. Elbo -115.483856201. KL 7.28692436218
    Iteration 260. Elbo -135.568496704. KL 7.45468378067
    Iteration 270. Elbo -126.830368042. KL 7.62879562378
    Iteration 280. Elbo -121.736724854. KL 7.79468870163
    Iteration 290. Elbo -120.828201294. KL 7.97476053238
    Iteration 300. Elbo -123.977989197. KL 8.13904094696
    Iteration 310. Elbo -127.896560669. KL 8.29942989349
    Iteration 320. Elbo -130.98664856. KL 8.4574213028
    Iteration 330. Elbo -107.724006653. KL 8.59646511078
    Iteration 340. Elbo -96.2202224731. KL 8.7516078949
    Iteration 350. Elbo -99.8892974854. KL 8.88774490356
    Iteration 360. Elbo -110.625091553. KL 9.02900981903
    Iteration 370. Elbo -110.705986023. KL 9.16001415253
    Iteration 380. Elbo -118.763916016. KL 9.29595279694
    Iteration 390. Elbo -101.959465027. KL 9.42527866364
    Iteration 400. Elbo -108.359840393. KL 9.55098152161
    Iteration 410. Elbo -94.4869918823. KL 9.68084907532
    Iteration 420. Elbo -88.2326431274. KL 9.8120470047
    Iteration 430. Elbo -104.121643066. KL 9.92966651917
    Iteration 440. Elbo -112.153770447. KL 10.0524768829
    Iteration 450. Elbo -102.999305725. KL 10.1651391983
    Iteration 460. Elbo -110.671699524. KL 10.2773513794
    Iteration 470. Elbo -99.833732605. KL 10.3971281052
    Iteration 480. Elbo -101.300125122. KL 10.5014190674
    Iteration 490. Elbo -106.583534241. KL 10.6122703552
    Iteration 500. Elbo -100.561035156. KL 10.7189235687
    Iteration 510. Elbo -99.5090408325. KL 10.8257017136
    Iteration 520. Elbo -94.6386489868. KL 10.9309864044
    Iteration 530. Elbo -92.851524353. KL 11.0329875946
    Iteration 540. Elbo -83.7906646729. KL 11.1409406662
    Iteration 550. Elbo -100.541503906. KL 11.242099762
    Iteration 560. Elbo -84.2009353638. KL 11.3463287354
    Iteration 570. Elbo -99.4818344116. KL 11.4444103241
    Iteration 580. Elbo -87.810798645. KL 11.5469789505
    Iteration 590. Elbo -90.5667190552. KL 11.6429700851
    Iteration 600. Elbo -89.6919326782. KL 11.7347326279
    Iteration 610. Elbo -89.1579208374. KL 11.8316688538
    Iteration 620. Elbo -95.3815689087. KL 11.924328804
    Iteration 630. Elbo -97.5506820679. KL 12.0179824829
    Iteration 640. Elbo -89.58984375. KL 12.1033706665
    Iteration 650. Elbo -87.0892181396. KL 12.1903495789
    Iteration 660. Elbo -83.8239746094. KL 12.2757358551
    Iteration 670. Elbo -98.5659866333. KL 12.3601484299
    Iteration 680. Elbo -85.6225814819. KL 12.454410553
    Iteration 690. Elbo -85.4044265747. KL 12.5378456116
    Iteration 700. Elbo -81.7902679443. KL 12.6191778183
    Iteration 710. Elbo -90.9003829956. KL 12.701210022
    Iteration 720. Elbo -86.0762786865. KL 12.7857275009
    Iteration 730. Elbo -81.9525756836. KL 12.8606281281
    Iteration 740. Elbo -82.8061523438. KL 12.9407939911
    Iteration 750. Elbo -75.7044754028. KL 13.014670372
    Iteration 760. Elbo -85.7485656738. KL 13.0864629745
    Iteration 770. Elbo -84.7195129395. KL 13.1637172699
    Iteration 780. Elbo -81.5038986206. KL 13.2388887405
    Iteration 790. Elbo -91.4418334961. KL 13.305138588
    Iteration 800. Elbo -91.3380737305. KL 13.378786087
    Iteration 810. Elbo -79.8998031616. KL 13.4527482986
    Iteration 820. Elbo -82.5501861572. KL 13.5232639313
    Iteration 830. Elbo -73.904045105. KL 13.5923223495
    Iteration 840. Elbo -94.8883132935. KL 13.6631288528
    Iteration 850. Elbo -82.671081543. KL 13.7367839813
    Iteration 860. Elbo -89.1314620972. KL 13.8053359985
    Iteration 870. Elbo -89.9466018677. KL 13.8741493225
    Iteration 880. Elbo -72.7582550049. KL 13.9450321198
    Iteration 890. Elbo -86.2754516602. KL 14.0125989914
    Iteration 900. Elbo -75.2065505981. KL 14.0771455765
    Iteration 910. Elbo -76.9281082153. KL 14.1364145279
    Iteration 920. Elbo -76.7940979004. KL 14.2024030685
    Iteration 930. Elbo -71.4757919312. KL 14.2663583755
    Iteration 940. Elbo -78.8379974365. KL 14.3307580948
    Iteration 950. Elbo -69.5632553101. KL 14.3894805908
    Iteration 960. Elbo -84.2538375854. KL 14.4518518448
    Iteration 970. Elbo -93.6064071655. KL 14.5128087997
    Iteration 980. Elbo -83.2115631104. KL 14.5737915039
    Iteration 990. Elbo -82.0588302612. KL 14.636089325



```
fig, axes = plt.subplots(1, 3, figsize=(3*8,5))

axes[0].plot(reparam_elbos, label='ELBO')
axes[0].set_title('Time', fontsize=15)
axes[0].set_ylim((-500, -50))
axes[0].legend()

axes[1].plot(reparam_kls, label='KL')
axes[1].set_title('Time', fontsize=15)
axes[1].legend()

axes[2].plot(reparam_accuracies, label='Accuracy')
axes[2].set_title('Time', fontsize=15)
axes[2].legend()
```

## Training the model via reinforce


```
# Important: the multiplication has to be done elementwise. 
# We have to train the model using a surrogate loss.
# Since TensorFlow does automatic differentation, we have to use 
# stop_gradient on the loss, as well as the samples!
per_sample_reinforce_loss =  - posterior.log_prob(tf.stop_gradient(posterior_samples)) * tf.stop_gradient(param_log_probs)
reinforce_loss = tf.reduce_mean(per_sample_reinforce_loss) + kl
```


```
optimizer = tf.train.GradientDescentOptimizer(0.0001)
reinforce_min_op = optimizer.minimize(reinforce_loss)
```


```
NUM_ITERATIONS = 1000
```


```
reinforce_accuracies = []
reinforce_kls = []
reinforce_elbos = []

sess = tf.Session()

# Initialize all variables
sess.run(tf.initialize_all_variables())

for i in xrange(NUM_ITERATIONS):
  sess.run(reinforce_min_op)

  if i % 10 == 0:
    reinforce_acc, reinforce_kl, reinforce_elbo = sess.run([accuracy, kl, elbo])
    reinforce_accuracies += [reinforce_acc]
    reinforce_kls += [reinforce_kl]
    reinforce_elbos += [reinforce_elbo]

    print('Iteration {}. Elbo {}. KL {}'.format(
        i, reinforce_elbo, reinforce_kl))

reinforce_learned_mean, reinforce_learned_log_scale = sess.run(learned_vars)
```

    Iteration 0. Elbo -1231.48754883. KL 0.0220593810081
    Iteration 10. Elbo -808.592407227. KL 0.338543713093
    Iteration 20. Elbo -524.581542969. KL 0.909575939178
    Iteration 30. Elbo -379.652740479. KL 1.41321325302
    Iteration 40. Elbo -301.671539307. KL 1.8704212904
    Iteration 50. Elbo -279.019744873. KL 2.25260615349
    Iteration 60. Elbo -287.911437988. KL 2.66793394089
    Iteration 70. Elbo -248.627319336. KL 3.03274250031
    Iteration 80. Elbo -274.209838867. KL 3.39950537682
    Iteration 90. Elbo -238.783615112. KL 3.7056684494
    Iteration 100. Elbo -180.144882202. KL 4.10226154327
    Iteration 110. Elbo -216.213302612. KL 4.40311670303
    Iteration 120. Elbo -183.531951904. KL 4.69419574738
    Iteration 130. Elbo -186.911499023. KL 4.96079969406
    Iteration 140. Elbo -173.515960693. KL 5.19095230103
    Iteration 150. Elbo -157.491165161. KL 5.43379211426
    Iteration 160. Elbo -158.108306885. KL 5.6262254715
    Iteration 170. Elbo -163.221206665. KL 5.8356795311
    Iteration 180. Elbo -146.565139771. KL 6.06200742722
    Iteration 190. Elbo -190.070983887. KL 6.17145061493
    Iteration 200. Elbo -153.055160522. KL 6.33472442627
    Iteration 210. Elbo -161.657592773. KL 6.56980705261
    Iteration 220. Elbo -116.549758911. KL 6.77828407288
    Iteration 230. Elbo -125.225143433. KL 7.03187656403
    Iteration 240. Elbo -147.91607666. KL 7.22355747223
    Iteration 250. Elbo -133.264282227. KL 7.41294574738
    Iteration 260. Elbo -126.578613281. KL 7.57341766357
    Iteration 270. Elbo -119.080780029. KL 7.75994729996
    Iteration 280. Elbo -109.886116028. KL 7.95444965363
    Iteration 290. Elbo -121.945098877. KL 8.08462715149
    Iteration 300. Elbo -124.018943787. KL 8.18129253387
    Iteration 310. Elbo -106.144706726. KL 8.37997722626
    Iteration 320. Elbo -121.718406677. KL 8.57028007507
    Iteration 330. Elbo -125.367019653. KL 8.72685909271
    Iteration 340. Elbo -118.693809509. KL 8.87324142456
    Iteration 350. Elbo -101.319114685. KL 9.01533222198
    Iteration 360. Elbo -114.103240967. KL 9.16397762299
    Iteration 370. Elbo -91.875869751. KL 9.27885913849
    Iteration 380. Elbo -114.098571777. KL 9.42788505554
    Iteration 390. Elbo -114.865745544. KL 9.54022693634
    Iteration 400. Elbo -94.1702957153. KL 9.65817832947
    Iteration 410. Elbo -109.974090576. KL 9.79271125793
    Iteration 420. Elbo -109.752044678. KL 9.85466575623
    Iteration 430. Elbo -110.549911499. KL 10.0302639008
    Iteration 440. Elbo -108.423713684. KL 10.163107872
    Iteration 450. Elbo -80.6597824097. KL 10.3107624054
    Iteration 460. Elbo -103.389984131. KL 10.4190883636
    Iteration 470. Elbo -101.75. KL 10.5393753052
    Iteration 480. Elbo -112.833908081. KL 10.6717405319
    Iteration 490. Elbo -94.2823638916. KL 10.7516059875
    Iteration 500. Elbo -102.636833191. KL 10.8633699417
    Iteration 510. Elbo -90.9130401611. KL 10.9743413925
    Iteration 520. Elbo -76.0304870605. KL 11.0320749283
    Iteration 530. Elbo -102.61164856. KL 11.130033493
    Iteration 540. Elbo -97.4822235107. KL 11.2614850998
    Iteration 550. Elbo -80.1745529175. KL 11.3975610733
    Iteration 560. Elbo -92.6554031372. KL 11.5362567902
    Iteration 570. Elbo -87.7685241699. KL 11.6449260712
    Iteration 580. Elbo -113.672859192. KL 11.7816791534
    Iteration 590. Elbo -89.4435577393. KL 11.870429039
    Iteration 600. Elbo -98.1239089966. KL 11.971988678
    Iteration 610. Elbo -86.9800796509. KL 12.0661430359
    Iteration 620. Elbo -105.8777771. KL 12.1732368469
    Iteration 630. Elbo -90.7897491455. KL 12.2917041779
    Iteration 640. Elbo -94.3220977783. KL 12.3896751404
    Iteration 650. Elbo -88.0863800049. KL 12.4850358963
    Iteration 660. Elbo -84.8598632812. KL 12.6095962524
    Iteration 670. Elbo -90.5580291748. KL 12.6718521118
    Iteration 680. Elbo -79.3591690063. KL 12.7729558945
    Iteration 690. Elbo -90.6020584106. KL 12.8668231964
    Iteration 700. Elbo -80.9606170654. KL 12.9421491623
    Iteration 710. Elbo -98.1525802612. KL 13.0128526688
    Iteration 720. Elbo -89.2859802246. KL 13.1471347809
    Iteration 730. Elbo -84.0707321167. KL 13.1970977783
    Iteration 740. Elbo -89.5597915649. KL 13.2572488785
    Iteration 750. Elbo -99.1416320801. KL 13.3330011368
    Iteration 760. Elbo -92.910118103. KL 13.4171209335
    Iteration 770. Elbo -76.833732605. KL 13.5153656006
    Iteration 780. Elbo -77.1132965088. KL 13.5576887131
    Iteration 790. Elbo -93.1914291382. KL 13.6402273178
    Iteration 800. Elbo -88.9732894897. KL 13.7351942062
    Iteration 810. Elbo -101.056121826. KL 13.8010444641
    Iteration 820. Elbo -99.5751037598. KL 13.8656616211
    Iteration 830. Elbo -91.9339675903. KL 13.9571743011
    Iteration 840. Elbo -70.7454147339. KL 14.0105953217
    Iteration 850. Elbo -91.7919006348. KL 14.0381126404
    Iteration 860. Elbo -82.0033187866. KL 14.1070785522
    Iteration 870. Elbo -81.9294433594. KL 14.1275062561
    Iteration 880. Elbo -103.577262878. KL 14.1831588745
    Iteration 890. Elbo -88.1100387573. KL 14.2206459045
    Iteration 900. Elbo -95.5690231323. KL 14.2572431564
    Iteration 910. Elbo -68.3534698486. KL 14.3253746033
    Iteration 920. Elbo -86.1409835815. KL 14.3638477325
    Iteration 930. Elbo -76.2172088623. KL 14.3877239227
    Iteration 940. Elbo -67.9566345215. KL 14.4778957367
    Iteration 950. Elbo -78.4613800049. KL 14.5488395691
    Iteration 960. Elbo -76.8335723877. KL 14.587480545
    Iteration 970. Elbo -75.7317123413. KL 14.6351766586
    Iteration 980. Elbo -85.0177459717. KL 14.7029514313
    Iteration 990. Elbo -96.061920166. KL 14.7556276321



```
fig, axes = plt.subplots(1, 3, figsize=(3*8,5))

axes[0].plot(reinforce_elbos, label='ELBO')
axes[0].set_title('Time', fontsize=15)
axes[0].set_ylim((-500, -50))
axes[0].legend()

axes[1].plot(reinforce_kls, label='KL')
axes[1].set_title('Time', fontsize=15)
axes[1].legend()

axes[2].plot(reinforce_accuracies, label='Accuracy')
axes[2].set_title('Time', fontsize=15)
axes[2].legend()
```

## Analysis

* How sensitive are the different estimators to the number of posterior samples used? Try 1, 10, 100.
* How sensitive are the different estimators to the learning rate used? Try 1e-4, 1e-3 and 1e-2 for both estimators.
* How sensitive are the different estimators to batch size? Try batch size 10 for both estimators.

# Answer
* Reparametriztion less sensitive to changes in posterior samples, learning rates, batch size.

# Visualize learned parameters

* Let's visualize the learned parameters, as well as check the different uncertanties around different parameters
* Can we now see which parameters are most important in the classification problem?


```
learned_mean = reparam_learned_mean
learned_scale = np.exp(reparam_learned_log_scale)
```


```
PARAM_INDEX = 19
```


```
def gaussian(x, mu, sig):
    return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
```


```
mean = learned_mean[PARAM_INDEX]
scale = learned_scale[PARAM_INDEX]
x_values = np.linspace(-3, 3, num=1000)
y_values = [gaussian(x, mean, scale) for x in x_values]

plt.figure(figsize=(7,5))
plt.plot(x_values, y_values)
plt.vlines(
      learned_mean[PARAM_INDEX], 
      ymin=0, 
      ymax=gaussian(mean, mean, scale),
      colors='g',
      alpha=0.8)

x_std_values = np.linspace(mean - scale, mean + scale, num=1000)
y_std_values = [gaussian(x, mean, scale) for x in x_std_values]
plt.fill_between(x_std_values, y_std_values, 0, alpha=0.3, color='b')

plt.title('Parameter name= {}'.format(feature_names[PARAM_INDEX]))
```

## Visualize gradient variance - which estimator has lower variance?

* Let's assess the gradient variance per of each estimator, by looking at how much gradients vary for each estimators. 


```
def per_batch_grads(ys, xs):
  """Computes the gradients of ys[i] wrt xs."""
  grads = []

  for y in tf.unstack(ys, axis=0):
    grads.append(tf.squeeze(tf.gradients(y, xs)[0]))
  return tf.stack(grads, axis=0)
```


```
def compute_grad_var(grads, param_index):
  return np.mean(grads[:, param_index]), np.std(grads[:, param_index])**2
```


```
per_batch_reinforce_grads = per_batch_grads(per_sample_reinforce_loss, learned_vars[0])
per_batch_reinforce_grads.shape
```




    TensorShape([Dimension(100), Dimension(31)])




```
per_batch_reparam_grads = per_batch_grads(per_sample_reparametrization_loss, learned_vars[0])
per_batch_reparam_grads.shape
```




    TensorShape([Dimension(100), Dimension(31)])




```
reinforce_grads = []

sess = tf.Session()

# Initialize all variables
sess.run(tf.initialize_all_variables())

for i in xrange(NUM_ITERATIONS):
  sess.run(reinforce_min_op)

  if i % 10 == 0:
    reinforce_grads.append(sess.run(per_batch_reinforce_grads))
```


```
reparam_grads = []

sess = tf.Session()

# Initialize all variables
sess.run(tf.initialize_all_variables())

for i in xrange(NUM_ITERATIONS):
  sess.run(reparam_min_op)

  if i % 10 == 0:
    reparam_grads.append(sess.run(per_batch_reparam_grads))
```


```
PARAM_INDEX = 19
```


```
reinforce_grad_vars = []
for timestep_grad in reinforce_grads:
  _, reinforce_grad_var = compute_grad_var(timestep_grad, PARAM_INDEX)
  reinforce_grad_vars.append(reinforce_grad_var)
```


```
reparam_grad_vars = []
for timestep_grad in reparam_grads:
  _, reparam_grad_var = compute_grad_var(timestep_grad, PARAM_INDEX)
  reparam_grad_vars.append(reparam_grad_var)
```


```
plt.figure(figsize=(7,5))
plt.plot(reinforce_grad_vars, label='reinforce')
plt.plot(reparam_grad_vars, label='reparametrization')
plt.yscale('log')
plt.xlabel('It')
plt.ylabel('Gradient variance')
plt.legend()
```

### Analysis

* What do you notice about gradient variance as training progresses?
* What do you notice about the variance of the different estimators?
* How do batch size and number of posterior samples affect the variance of the different estimators?
* How does gradient variance affect convergence?
* Which gradient estimator would you use in practice for this problem?


# Answer 
 * Gradient variance goes down as training progresses regardless of estimator.
 * Reinforce = bad. 
 * Higher batch size and number of samples = lower variance.
 * Learning rate: lower learning rate is a way to counter balance high gradient variance.
 * Unless you reduce learning rates in training, you might oscilate too much, which slows down learning.
 * Reparametrization. 

47bce430055e84368f4560ed8f63a5824344f5f4
ipynb
Jupyter Notebook
unsupervised/1_bayesian_logistic_regression.ipynb
hazratali/PracticalSessions
95c10a71267f6fb1f61a6bd6f8ea33041365ceb5
MIT
2020-02-11T23:34:05.000Z
2021-08-20T02:58:08.000Z
unsupervised/1_bayesian_logistic_regression.ipynb
hazratali/PracticalSessions
95c10a71267f6fb1f61a6bd6f8ea33041365ceb5
MIT
unsupervised/1_bayesian_logistic_regression.ipynb
hazratali/PracticalSessions
95c10a71267f6fb1f61a6bd6f8ea33041365ceb5
MIT
2019-06-27T12:15:38.000Z
2019-10-21T12:29:54.000Z
Qwen/Qwen-72B
1. YES
2. YES