text:
```python
%matplotlib inline

import matplotlib.pyplot as plt

import math
import numpy as np

from scipy.integrate import solve_ivp
```

- https://docs.scipy.org/doc/scipy-1.3.0/reference/integrate.html#solving-initial-value-problems-for-ode-systems
- https://docs.scipy.org/doc/scipy-1.3.0/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp (the new API)

- https://docs.scipy.org/doc/scipy-1.3.0/reference/tutorial/integrate.html#ordinary-differential-equations-odeint (the old API)

## Conversion d'une ODE d'ordre $n$ en un système d'ODE d'ordre 1

C.f. "Toutes les mathématiques et les bases de l'informatique" de H.Stocker (Dunod) p.745

### Exemple 1

$$y'' + 5y'^2 + 7y^5 = 6x^3$$

On pose $y_1 = y'$ et donc $y_1' = y''$.

On obtient alors le système couplé:

$$
\left\{
    \begin{array}{ll}
        y_1' + 5 y_1^2 + 7 y^5 & = 6 x^3 \\
        y' - y_1 & = 0
    \end{array}
\right.
$$

### Exemple 2

$$\color{red}{y^{(6)}} + 8 x^2 \color{green}{y^{(4)}} + \ln y = 0$$

On pose :
- $y_1 = y'$
- $y_2 = y_1' = y''$
- $y_3 = y_2' = y_1'' = y'''$
- $\color{green}{y_4} = y_3' = \dots = \color{green}{y^{(4)}}$
- $y_5 = y_4' = \dots = y^{(5)}$
- $\color{red}{y_6} = y_5' = \dots = \color{red}{y^{(6)}}$

On obtient alors le système couplé:

$$
\left\{
    \begin{array}{ll}
        \color{red}{y_6} + 8x^2 \color{green}{y_4} + \ln y & = 0 \\
        y_1 - y' & = 0 \\
        y_2 - y_1' & = 0 \\
        y_3 - y_2' & = 0 \\
        y_4 - y_3' & = 0 \\
        y_5 - y_4' & = 0 \\
        y_6 - y_5' & = 0 \\
    \end{array}
\right.
$$

## Solving initial value problems (IVP) for ODE systems

- https://pundit.pratt.duke.edu/wiki/Python:Ordinary_Differential_Equations
- https://pundit.pratt.duke.edu/wiki/Python:Ordinary_Differential_Equations/Examples

### Example 1

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations&printable=yes#Simple_Example

Solve

$$y' + y = t$$

i.e.

$$\frac{\text{d}}{\text{d}t} y(t) = t - y(t)$$

for some times between 0 and 15, assuming the initial value for $y(0)$ is 2.

```python
integration_range = [0, 15]
initial_values = [2]

# Define the dependent function f
def f(t, y):
    dydt = t-y
    return dydt

sol = solve_ivp(f, integration_range, initial_values)
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

```python
integration_range = [0, 15]
initial_values = [2]

# Define the dependent function f
def f(t, y):
    dydt = t-y
    return dydt

sol = solve_ivp(f, integration_range, initial_values, dense_output=True)
print(sol)

t = np.linspace(0, 15, 50)
y = sol.sol(t)               # <- with dense_output=True, solve_ivp() returned a callable (function) solution instead of a sample (array)

plt.plot(t, y[0], 'o-');
```

```python
integration_range = [0, 15]
initial_values = [2]

# Define the dependent function f
def f(t, y):
    dydt = t-y
    #print(y.shape)
    return dydt

sol = solve_ivp(f, integration_range, initial_values,
                t_eval = np.linspace(0, 15, 30))         # <--
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 2

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations&printable=yes#Another_Example

Solve

$$\frac{\text{d}}{\text{d}t} y(t) + \frac{1}{4} y(t) - x(t) = 0$$

i.e.

$$\frac{\text{d}}{\text{d}t} y(t) = - \frac{1}{4} y(t) + x(t)$$

assuming

$$
\begin{align}
x(t) & = \cos(3t) \\
y(0) & = -1
\end{align}
$$

for some times between 0 and 15.

```python
integration_range = [0, 15]
initial_values = [-1]

def f(t, y):
    dydt = -y + np.cos(3. * t)
    return dydt

sol = solve_ivp(f, integration_range, initial_values,
                t_eval = np.linspace(0, 15, 100))
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 3 : constant rate of change

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Constant_Rate_of_Change

Solve

$$\frac{\text{d}}{\text{d}t} y(t) = c$$

assuming

$$
y(0) = 6
$$

for some times between 0 and 10.

```python
integration_range = [0, 10]
initial_values = [6]
c = 1.2

# Define the dependent function f
def f(t, y):
    dydt = [c]
    return dydt

sol = solve_ivp(f, integration_range, initial_values)
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 4 : time-dependent rate of change

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Time-dependent_Rate_of_Change

Solve

$$\frac{\text{d}}{\text{d}t} y(t) = 2 t^2 - 6 t + 3$$

assuming

$$
y(0) = 6
$$

for some times between 0 and 4.

```python
integration_range = [0, 4]
initial_values = [6]

# Define the dependent function f
def f(t, y):
    dydt = 2. * t**2 - 6. * t + 3.
    return dydt

sol = solve_ivp(f, integration_range, initial_values,
                t_eval=np.linspace(0, 4, 30))
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 5 : population growth

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Population_Growth

For population growth, the rate of change of population is dependent upon the number of people as well as some constant of proportionality:

$$\frac{\text{d}}{\text{d}t} y(t) = c y(t)$$

assuming

$$
\begin{align}
y(0) & = 10   \\
c    & = 1.02
\end{align}
$$

for some times between 0 and 3.

```python
integration_range = [0, 3]
initial_values = [10]
c = 1.02

# Define the dependent function f
def f(t, y):
    dydt = c * y
    return dydt

sol = solve_ivp(f, integration_range, initial_values,
                t_eval=np.linspace(0, 3, 30))
print(sol)

plt.plot(sol.t, sol.y[0], 'o-');
```

### Example 6 : multiple dependent variable models

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Multiple_Variable_Models

Solve:

$$
\left\{
\begin{align}
\frac{\text{d}}{\text{d}t} y_0(t) & = 4 \cos(3t) \\
\frac{\text{d}}{\text{d}t} y_1(t) & = -2 y_0(t) + 0.5 t
\end{align}
\right.
$$

assuming

$$
\begin{align}
y_1(0) & = 0   \\
y_2(0) & = -3
\end{align}
$$

for some times between 0 and 5.

```python
integration_range = [0, 5]
initial_values = [0, -3]

# Define the dependent function f
def f(t, y):
    dy1dt = 4. * math.cos(3. * t)
    dy2dt = -2. * y[0] + 0.5 * t
    return [dy1dt, dy2dt]

sol = solve_ivp(f, integration_range, initial_values,
                t_eval=np.linspace(0, 5, 50))
print(sol)

fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(16, 4))
ax1.plot(sol.t, sol.y[0], 'o-', label="y1")
ax1.legend()
ax2.plot(sol.t, sol.y[1], 'o-', label="y2")
ax2.legend();
```

### Example 7 : higher Order Differential Equations

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations/Examples&printable=yes#Higher_Order_Differential_Equations

Solve:

$$
\frac{\text{d}^3}{\text{d}t^3} y(t) = c
$$

This third-order equation can be rewritten as a system of three first-order differential equations:

$$
\left\{
\begin{align}
y_0 = y \\
\frac{\text{d}}{\text{d}t}y_0(t) & = \frac{\text{d}}{\text{d}t}y(t) = y_1 \\
\frac{\text{d}}{\text{d}t}y_1(t) & = \frac{\text{d}^2}{\text{d}t^2} y_0(t) = \frac{\text{d}^2}{\text{d}t^2} y(t) = y_2 \\
\frac{\text{d}}{\text{d}t}y_2(t) & = \frac{\text{d}^2}{\text{d}t^2} y_1(t) = \frac{\text{d}^3}{\text{d}t^3} y_0(t) = \frac{\text{d}^3}{\text{d}t^3} y(t) = c \\
\end{align}
\right.
$$

$$
\Leftrightarrow
\left\{
\begin{align}
y_0 = y \\
\frac{\text{d}}{\text{d}t} y_0(t) & = y_1 \\
\frac{\text{d}}{\text{d}t} y_1(t) & = y_2 \\
\frac{\text{d}}{\text{d}t} y_2(t) & = c \\
\end{align}
\right.
$$

assuming

$$
\begin{align}
y_1(0) & = 6   \\
y_2(0) & = 2   \\
y_3(0) & = -4
\end{align}
$$

i.e.

- an initial position of 6
- an initial velocity of 2
- an initial acceleration of -4

and a constant acceleration of 1.3
over a period of 8 seconds.

```python
integration_range = [0, 8]
initial_values = [6, 2, -4]
c = 1.3

# Define the dependent function f
def f(t, y):
    dy1dt = y[1]
    dy2dt = y[2]
    dy3dt = c
    return [dy1dt, dy2dt, dy3dt]

sol = solve_ivp(f, integration_range, initial_values,
                t_eval=np.linspace(0, 8, 50))
print(sol)

fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(14, 4))
ax1.plot(sol.t, sol.y[0], 'o-', label="y1")
ax1.legend()
ax2.plot(sol.t, sol.y[1], 'o-', label="y2")
ax2.legend();
ax3.plot(sol.t, sol.y[2], 'o-', label="y3")
ax3.legend();
```

### Remark about the initial value

Taken from https://pundit.pratt.duke.edu/piki/index.php?title=Python:Ordinary_Differential_Equations&printable=yes#Notes_.2F_Troubleshooting

Lets imagine we want to solve

$$\frac{\text{d}}{\text{d}t} y(t) = c$$

assuming

$$
y(0) = 1
$$

for some times between **5** and 10.

```python
integration_range = [5, 10]
initial_values = [1]
c = 1.2

# Define the dependent function f
def f(t, y):
    dydt = [c]
    return dydt

res = solve_ivp(f, integration_range, initial_values, t_eval=np.linspace(5, 10, 10))
print(res)

plt.plot(res.t, res.y[0], 'o-');
```

Here, we can see that scipy considers that the actual initial value $y(0) = 1$ is not for $t=0$ but for the start of the time span !

ext: ipynb

max_stars_repo_name: jdhp-docs/python-notebooks

max_issues_repo_name: jdhp-docs/python-notebooks

max_forks_repo_name: jdhp-docs/python-notebooks


text:
## Portfolio Optimization Example

In portfolio optimization [1], [2], we are interested in finding an optimal portfolio $w\in \mathbb{R}^n$ of $n$ assets. The return $r$ has expected value $\alpha$ and covariance $\Sigma \in \mathbb{S}_{++}^n$. The risk aversion factor is denoted by $\gamma \geq 0$, and transaction (short-selling) cost is denoted by $\kappa_\mathrm{tc(sh)} \geq 0$. We solve the optimization

\begin{equation}
\begin{array}{ll}
\text{maximize} \quad &\alpha^T w - \gamma w^T \Sigma w - \kappa_\mathrm{tc}^T |w-w^\mathrm{prev}| + \kappa_\mathrm{sh}^T (w)_- \\
\text{subject to} \quad &\mathbb{1}^T w = 1 \\
&\Vert w \Vert_1 \leq L,
\end{array}
\end{equation}

where $w \in \mathbb{R}^n$ is the variable. The previous portfolio is $w^\mathrm{prev}$, $(\cdot)_-$ represents the argument's negative part, and $L \geq 1$ is the maximum value of total long positions plus the total magnitude of short positions.

With the factor model $\Sigma = F \Sigma^f F^T + D$ [3], we can rewrite an equivalent and [DPP-compliant](https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming) problem, i.e.,

\begin{equation}
\begin{array}{ll}
\text{maximize} \quad &a^T w - \Vert \left(\Sigma^f\right)^{1/2} f \Vert_2^2 - \Vert D^{1/2} w\Vert_2^2 - k_\mathrm{tc}^T |\Delta w| + k_\mathrm{sh}^T (w)_- \\
\text{subject to} &f = F^T w \\
&\mathbb{1}^T w = 1 \\
&\Vert w \Vert_1 \leq L \\
&\Delta w = w-w^\mathrm{prev}, \\
\end{array}
\end{equation}

where $w \in \mathbb{R}^n$ and $f \in \mathbb{R}^m$ are variables. To summarize, the parameters are:

\begin{equation}
\begin{array}{ll}
a = \frac{\alpha}{\gamma} \\
F \\
\left(\Sigma^f\right)^{1/2} \\
D^{1/2} \\
k_\mathrm{tc} = \frac{\kappa_\mathrm{tc}}{\gamma} \\
k_\mathrm{sh} = \frac{\kappa_\mathrm{sh}}{\gamma} \\
w^\mathrm{prev} \\
L.
\end{array}
\end{equation}

Note that we divided the objective function by the risk aversion factor $\gamma$. This way, updating the value of $\gamma$ only affects the linear part of the objective function, avoiding to compute a matrix factorization when solving the problem repeatedly.

Let's define the corresponding CVXPY problem. Note that we define the diagonal of $D^{1/2}$ as a vector.

```python
import cvxpy as cp
import numpy as np

# define dimensions
n, m = 100, 10

# define variables
w = cp.Variable(n, name='w')
delta_w = cp.Variable(n, name='delta_w')
f = cp.Variable(m, name='f')

# define parameters
a = cp.Parameter(n, name='a')
F = cp.Parameter((n, m), name='F')
Sig_f_sqrt = cp.Parameter((m, m), name='Sig_f_sqrt')
d_sqrt = cp.Parameter(n, name='d_sqrt')
k_tc = cp.Parameter(n, nonneg=True, name='k_tc')
k_sh = cp.Parameter(n, nonneg=True, name='k_sh')
w_prev = cp.Parameter(n, name='w_prev')
L = cp.Parameter(nonneg=True, name='L')

# define objective
objective = cp.Maximize(a@w
                        -cp.sum_squares(Sig_f_sqrt@f)
                        -cp.sum_squares(cp.multiply(d_sqrt, w))
                        -k_tc@cp.abs(delta_w)
                        +k_sh@cp.minimum(0, w))

# define constraints
constraints = [f == F.T@w,
               np.ones(n)@w == 1,
               cp.norm(w, 1) <= L,
               delta_w == w-w_prev]

# define problem
problem = cp.Problem(objective, constraints)
```

Assign parameter values and solve the problem.

```python
np.random.seed(0)
gamma = 1
alpha = np.random.randn(n)
kappa_tc = 0.01*np.ones(n)
kappa_sh = 0.05*np.ones(n)

a.value = alpha/gamma
F.value = np.random.randn(n, m)
Sig_f_sqrt.value = np.random.rand(m, m)
d_sqrt.value = np.random.rand(n)
k_tc.value = kappa_tc/gamma
k_sh.value = kappa_sh/gamma
w_prev_unnormalized = np.random.rand(n)
w_prev.value = w_prev_unnormalized/np.linalg.norm(w_prev_unnormalized)
L.value = 1.6

val = problem.solve()
```

Generating C source for the problem is as easy as:

```python
from cvxpygen import cpg

cpg.generate_code(problem, code_dir='portfolio_code')
```

Now, you can use a python wrapper around the generated code as a custom CVXPY solve method.

```python
from portfolio_code.cpg_solver import cpg_solve
import numpy as np
import pickle
import time

# load the serialized problem formulation
with open('portfolio_code/problem.pickle', 'rb') as f:
    prob = pickle.load(f)

# assign parameter values
n, m = 100, 10
np.random.seed(1)
gamma = 1
alpha = np.random.randn(n)
kappa_tc = 0.01*np.ones(n)
kappa_sh = 0.05*np.ones(n)
prob.param_dict['a'].value = alpha/gamma
prob.param_dict['F'].value = np.round(np.random.randn(n, m))
prob.param_dict['Sig_f_sqrt'].value = np.diag(np.random.rand(m))
prob.param_dict['d_sqrt'].value = np.random.rand(n)
prob.param_dict['k_tc'].value = kappa_tc/gamma
prob.param_dict['k_sh'].value = kappa_sh/gamma
prob.param_dict['w_prev'].value = np.zeros(n)
prob.param_dict['L'].value = 1.6

# solve problem conventionally
t0 = time.time()
# CVXPY chooses eps_abs=eps_rel=1e-5, max_iter=10000, polish=True by default,
# however, we choose the OSQP default values here, as they are used for code generation as well
val = prob.solve(eps_abs=1e-3, eps_rel=1e-3, max_iter=4000, polish=False)
t1 = time.time()
print('\nCVXPY\nSolve time: %.3f ms' % (1000 * (t1 - t0)))
print('Objective function value: %.6f\n' % val)

# solve problem with C code via python wrapper
prob.register_solve('CPG', cpg_solve)
t0 = time.time()
val = prob.solve(method='CPG')
t1 = time.time()
print('\nCVXPYgen\nSolve time: %.3f ms' % (1000 * (t1 - t0)))
print('Objective function value: %.6f\n' % val)
```

### References

[1] Lobo, Miguel Sousa, Maryam Fazel, and Stephen Boyd. "Portfolio optimization with linear and fixed transaction costs." Annals of Operations Research 152.1 (2007): 341-365.

[2] Moehle, N., Kochenderfer, M. J., Boyd, S., and Ang, A. "Tax-Aware Portfolio Construction via Convex Optimization." Journal of Optimization Theory and Applications 189.2 (2021): 364-383.

[3] Ng, Victor, Robert F. Engle, and Michael Rothschild. "A multi-dynamic-factor model for stock returns." Journal of Econometrics 52.1-2 (1992): 245-266.

max_stars_repo_name: cvxgrp/cvxpygen

max_issues_repo_name: cvxgrp/cvxpygen

max_forks_repo_name: cvxgrp/cvxpygen


text:
# Courte introduction à notebook

Les documents **ipython notebook** sont constitués de cellules (**cell**) contenant soit du texte, soit du code exécutable.

```
print "Hello world!"
```

    Hello world!

La sortie des cellules de code est imprimé juste en dessous, et permet d'imprimer autant du texte que des images. (Il n'est pas important de comprendre le code pour l'instant.

```
%pylab inline
from pylab import *

x = linspace(0, 5, 10)
y = x ** 2

figure()
plot(x, y, 'r')
xlabel('x')
ylabel('y')
title('title')
show()
```

Lorsqu'une cellule de code est exécuté, les variables contenu dans celles-ci sont accèssibles dans les cellules subséquentes

```
print x
```

    [ 0.          0.55555556  1.11111111  1.66666667  2.22222222  2.77777778
      3.33333333  3.88888889  4.44444444  5.        ]

Remarquez que la variable `x` a été défini à la cellule prédécente **In [4]**.

Les cellules contenant du texte permettent une mise en page plus flexible que de simples commentaires en ligne. Il est donc possible de rédiger un rapport (propre) contenant le code source. Cliquez sur cette cellules pour voir comment modifier le texte.

###On peut écrire des titres

* des
* listes

1. à points
2. ou numéroté.

**Éditer du texte en gras** ou *en italique*

> **Citer du texte**
> Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum

ou intégrer du code (non exécutable) pour le commenter

    def function_simple(x):
        return x+1

Et aussi (surtout?), on peut intégrer des fonctions mathématiques suivant la même syntaxe qu'en LaTeX

\begin{equation}f(x) = \sum_{i=0}^{m} w_ix_i + b\end{equation}

Libre à vous d'expérimenter avec l'interface notebook et de découvrir ses nombreux avantages.

## Exercices python

Maintenant, créez une nouvelle cellule pour le code et commencez votre lecture du tutoriel python! [version française](http://www.afpy.org/doc/python/2.7/tutorial/) [version anglaise](http://docs.python.org/tutorial/)

Lisez et reproduisez les exemples du tutoriel des parties suivantes: 1,3,4,5,6,7,8,10; en sautant les sous-parties: 3.1.3, 4.7.6, 5.8, 8.4, 8.5, 8.6, 10.5, 10.7, 10.10, 10.11, 10.12.

```

```


text:
```python
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.style.use('classic')
%matplotlib inline
```

# Class 10: Stochastic Time Series Processes

Most models of the business cycle are *stochastic* time series models. That is, the models incorporate randomness as a determinant of equilibrium. Randomness is necessary for two reasons. First, from a practical point of view, any model will be incomplete and will not fit the data perfectly and so randomness in a macroeconomic model is analagous to an error term in an OLS model. Second, from a philosophical perspective, randomness in a macroeconomic model reflects economists' consensus that fundamentally unpredictable forces cause the business cycle.

## Simulating Normal Random Variables with Numpy

Recall that the `numpy.random` module has functions for generating (pseudo) random variables. Learn more about the module by reading the documentation: https://docs.scipy.org/doc/numpy/reference/routines.random.html

We use the `numpy.random.normal()` function to crate arrays of random draws from the normal distribution. The function takes three arguments:
* `loc`: the mean of the distribution (default=0)
* `scale`: the standard deviation of the distribution (default=1)
* `size`: how many to numbers to draw (default = `None`)

The default is to draw numbers from the *standard normal* distribution.

### Example

Draw 500 values each from the $\mathcal{N}(0,1)$ and $\mathcal{N}(0,2^2)$ distributions. Plot.

```python
# Set the seed for the random number generator to 126
np.random.seed(126)

# Create two arrays:
#     x: 500 draws from the normal(0,1) distribution
#     y: 500 draws from the normal(0,4) distribution
x = np.random.normal(loc=0,scale=1,size=500)
y = np.random.normal(loc=0,scale=2,size=500)

# Plot x and y
plt.plot(x,lw=2,alpha = 0.6,label='$\sigma=1$')
plt.plot(y,lw=2,alpha = 0.6,label='$\sigma=2$')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid()
```

## The White Noise Process

In the previous example, we created two variables that stored draws from normal distrbutions with means of zero but with different standard deviations. Both of the variables were simulations of *white noise processes*. A white noise process is a random variable $\epsilon_t$ with constant mean and constant variance. We are concerned only with the zero-mean white noise process and we'll denote that a variable is a zero-mean white noise process with the following shorthand notation:

\begin{align}
\epsilon_t & \sim \text{WN}(0,\sigma^2),
\end{align}

where $\sigma^2$ is the variance of the processes. Strictly speaking, a white noise process can follow any distribution as long as the mean and variance are constant, but we'll concentrate exclusively white noise process drawn from the normal distribution.

## The AR(1) Process

A random variable $y_t$ is an *autoregressive process of order 1*  or AR(1) process if it can be written in the following form:

\begin{align}
y_t & = \rho y_{t-1} + \epsilon_t,
\end{align}

where $\rho$ is a constant and $\epsilon \sim \text{WN}(0,\sigma^2)$. The AR(1) process is the stochastic analog of the first-order difference equation where the random variable $\epsilon_t$ replaces the exogenous variable $w_t$.

### Example

Simulate an AR(1) process for 101 periods ($t = 0,\ldots, 100$) using the following parameter values:

\begin{align}
\rho & = 0.5\\
\sigma &  = 1\\
y_0 & = 0
\end{align}

Plot the simulated values for $y$.

```python
# Set the seed for the random number generator to 126
np.random.seed(126)

# Initialize values for T, y0, rho, and sigma
T = 101
y0 = 0
rho = 0.5
sigma = 1

# Initialize an array of zeros for y
y = np.zeros(T)

# Set the first value of y equal to y0
y[0] = y0

# Create a variable called 'epsilon' equal to an array containing T draws from the normal(0,sigma^2) process
eps= np.random.normal(loc=0,scale=sigma,size=T)

# Iterate over t in range(T-1) to compute y
for t in range(T-1):
    y[t+1] = rho*y[t] + eps[t+1]

# Plot y
plt.plot(y,lw=2)
plt.grid(linestyle=':')
```

The AR(1) process wtih $\rho = 0.5$ seems to fluctuate around the value $y=0$ and so the process appears to be stable

### Example

Simulate an AR(1) process for 101 periods ($t = 0,\ldots, 100$) using the following parameter values:

\begin{align}
\rho & = 1.5\\
\sigma &  = 1\\
y_0 & = 0
\end{align}

Plot the simulated values for $y$.

```python
# Set the seed for the random number generator to 126
np.random.seed(126)

# Initialize values for T, y0, rho, and sigma
T = 101
y0 = 0
rho = 1.5
sigma = 1

# Initialize an array of zeros for y
y = np.zeros(T)

# Set the first value of y equal to y0
y[0] = y0

# Create a variable called 'epsilon' equal to an array containing T draws from the normal(0,sigma^2) process
eps= np.random.normal(loc=0,scale=sigma,size=T)

# Iterate over t in range(T-1) to compute y
for t in range(T-1):
    y[t+1] = rho*y[t] + eps[t+1]

# Plot y
plt.plot(y,lw=2)
plt.grid(linestyle=':')
```

The AR(1) process wtih $\rho = 1.5$ seems to be approaching infinity in magnitude and so the process appears to be explosive.

In general, like the first-order difference equation, if $|\rho| < 1$, the the AR(1) process is stable. If $|\rho|>1$ then the process is explosive. A special case is the *random walk process* which is an AR(1) process with $\rho = 1$. The random walk process has many important applications including asset pricing theory.

The function in the next cell computes a simulation of an AR(1) process.

```python
# Define a function for simulating an AR(1) process. CELL PROVIDED
def ar1_sim(rho=0,sigma=1,y0=0,T=25):
    '''Funciton for simulating an AR(1) process for T periods

            Args:
            rho (float):   autoregressive parameter
            sigma (float): standard deviation of the white noise process
            y0 (float):    initial value of the process
            T (int):       number of periods to simulate

        Returns:
            numpy array
    '''

    # initialize y array
    y = np.zeros(T)
    y[0] = y0

    # draw random numbers for white noise process
    eps= np.random.normal(loc=0,scale=sigma,size=T-1)
    for t in range(T-1):
        y[t+1] = rho*y[t] + eps[t]

    return y
```

### Example:

Use the `ar1_sim()` function to simulate an AR(1) process for 201 periods ($t = 0,\ldots, 200$) using the following parameter values:

\begin{align}
\rho & = -0.99\\
\sigma &  = 0.5\\
y_0 & = 0
\end{align}

Set the seed for the NumPy random number generator to 126. Plot the simulated values for $y$.

```python
# Set the seed for the random number generator to 126
np.random.seed(126)

# Simulate y and plot
plt.plot(ar1_sim(rho=-0.99,sigma=0.5,y0=0,T=201))
plt.grid()
```

## Application: Estimate TFP as an AR(1) process

It is routine to model quarterly fulctuations in TFP as an AR(1) process and we will encounter this in our business cycle models. Here we will fit the following AR(1) model to US data:

\begin{align}
\log\left(A_t/A^{trend}_t\right) & = \rho \log\left(A_{t-1}/A^{trend}_{t-1}\right) + \epsilon_t
\end{align}

where $\log\left(A_t/A^{trend}_t\right) = \log A_t - \log A^{trend}$ is the log-deviation of TFP from its trend and $\epsilon_t$ is a white noize process with mean 0 and variance $\sigma^2$.

### The Data
The file `rbc_data_actual_trend.csv`, available at https://github.com/letsgoexploring/econ126/raw/master/Data/Csv/rbc_data_actual_trend.csv, contains actual and trend data for real GDP per capita, real consumption per capita, real investment per capita, real physical capital per capita, TFP, and hours per capita at quarterly frequency. The GDP, consumption, investment, and capital data are in terms of 2012 dollars. Hours is measured as an index with the value in October 2012 set to 100. All of the data are *real* quantities. That is, there are no *nominal* quantities like money or inflation or a nominal interest rate. The reason is that the first theory that we will encounter is called *real business cycle* or RBC theory and, in that theory, there is no place for nominal quantities. RBC theory seeks to explain fluctuations in real quantities as being primarily due to TFP shocks; i.e., shocks to the production function.

### Objectives

1. Import TFP data (trend and actual) and compute cyclical component as log-deviation from trend.
2. Construct a scatter plot of log-deviation of TFP from trend against *lagged* log-deviation of TFP from trend.
3. Estimate the AR(1) model of log-deviation of TFP from trend to obtain estimates of $\rho$ and $\sigma$.

```python
# Read business_cycle_data_actual_trend.csv into a Pandas DataFrame called 'df' with the first column set as the index and parse_dates=True
df = pd.read_csv('https://github.com/letsgoexploring/econ126/raw/master/Data/Csv/rbc_data_actual_trend.csv',index_col=0,parse_dates=True)
```

```python
# Construct a plot of TFP with its trend with:
#    1. Actual line: blue with lw=1, alpha=0.7, label = 'actual'
#    2. Trend line:  red with lw=3, alpha=0.7, label = 'trend'
plt.plot(df['tfp'],lw=2,alpha=0.6,label='Actual')
plt.plot(df['tfp_trend'],'r',lw=2,alpha=0.6,label='Trend')
plt.title('TFP for US from '+df.index[0].strftime('%B %Y')+' to '+df.index[-1].strftime('%B %Y'))
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid()
```

Now we need to compute the log-deviation of TFP from its trend. Note that since $\log(a/b) = \log(a) - \log(b)$, we can write:

\begin{align}
\log A_t - \log A^{trend} & = \log\left(A_t/A^{trend}_t\right)
\end{align}

The term on the right will make the AR(1) model easier to read.

\begin{align}
\log\left(A_t/A^{trend}_t\right) & = \rho \log\left(A_{t-1}/A^{trend}_{t-1}\right) + \epsilon_t
\end{align}

```python
# Create a new column to df called tfp_cycle equal to the log difference between actual TFP and it's trend:
df['tfp_cycle'] = np.log(df['tfp']/df['tfp_trend'])
```

```python
# Plot the log deviation of TFP from its trend (times 100)
plt.plot(df['tfp_cycle'],lw=2,alpha=0.6)
plt.title('TFP for US from '+df.index[0].strftime('%B %Y')+' to '+df.index[-1].strftime('%B %Y'))
plt.grid()
```

```python
# Create a new column to df called tfp_cycle_lag that, for each date, contains values
# in 'tfp_cycle' at the previous date
df['tfp_cycle_lag'] = df['tfp_cycle'].shift()

# Print the first five rows of only the 'tfp_cycle' and 'tfp_cycle_lag' columns. PROVIDED
print(df[['tfp_cycle','tfp_cycle_lag']].head())
```

                tfp_cycle  tfp_cycle_lag
    1948-01-01   0.001157            NaN
    1948-04-01   0.007843       0.001157
    1948-07-01  -0.001230       0.007843
    1948-10-01  -0.004375      -0.001230
    1949-01-01  -0.017405      -0.004375

Notice that there is a missing value inthe `tfp_cycle_lag` column for the first date in the index because there is no prior observation in the `tfp_cycle` column. To proceed with the estimation, we need to get rid of the row with missing values.

```python
# Use dropna() method on df to remove the row iwth the missing value
df= df.dropna()

# Print the first five rows of only the 'tfp_cycle' and 'tfp_cycle_lag' columns. PROVIDED
print(df[['tfp_cycle','tfp_cycle_lag']].head())
```

                tfp_cycle  tfp_cycle_lag
    1948-04-01   0.007843       0.001157
    1948-07-01  -0.001230       0.007843
    1948-10-01  -0.004375      -0.001230
    1949-01-01  -0.017405      -0.004375
    1949-04-01  -0.019563      -0.017405

Next, we should make a scatter plot of log-deviation of TFP from trend against the one-period lag of log-deviation of TFP from trend to see if out AR(1) model is a good idea.

```python
# Construct a scatter plot of log-deviation of TFP from trend against the one-period lag of log-deviation
# of TFP from trend with:
#    1. scatter points sized at least 50, opacity (alpha) no greater than 0.25
#    2. x- and y-axis limits: [-0.04,0.04]
plt.scatter(df['tfp_cycle_lag'],df['tfp_cycle'],s=50,alpha=0.25)
plt.xlabel('Lag log-deviation from trend')
plt.ylabel('Current log-deviation from trend')
plt.title('TFP for US from '+df.index[0].strftime('%B %Y')+' to '+df.index[-1].strftime('%B %Y'))
plt.xlim([-0.04,0.04])
plt.ylim([-0.04,0.04])
plt.grid()
```

Now we use StatsModels to fit the model. First we estimate the AR(1) model to obtain an estimate of $\rho$. Then we find the standard deviation of the residuals of the regression to estimate $\sigma$.

```python
# Create a variable 'X' to be the independent variable of the OLS model. Do not add a constant to X
X = df['tfp_cycle_lag']

# Create a variable 'Y' to be the dependent variable of the OLS model.
Y = df['tfp_cycle']

# Create a variable called 'model' that initializes the OLS model
model = sm.OLS(Y,X)

# Create a variable called 'results' that stores the results of the estimated OLS model
results = model.fit()

# Print output of summary2() method of results
print(results.summary2())
```

                            Results: Ordinary least squares
    ================================================================================
    Model:                  OLS              Adj. R-squared (uncentered): 0.540
    Dependent Variable:     tfp_cycle        AIC:                         -2124.8056
    Date:                   2020-02-12 12:52 BIC:                         -2121.1497
    No. Observations:       286              Log-Likelihood:              1063.4
    Df Model:               1                F-statistic:                 336.5
    Df Residuals:           285              Prob (F-statistic):          3.57e-50
    R-squared (uncentered): 0.541            Scale:                       3.4631e-05
    ------------------------------------------------------------------------------------
                        Coef.      Std.Err.        t        P>|t|      [0.025     0.975]
    ------------------------------------------------------------------------------------
    tfp_cycle_lag       0.7359       0.0401     18.3450     0.0000     0.6569     0.8148
    --------------------------------------------------------------------------------
    Omnibus:                    7.555             Durbin-Watson:               1.880
    Prob(Omnibus):              0.023             Jarque-Bera (JB):            9.064
    Skew:                       -0.243            Prob(JB):                    0.011
    Kurtosis:                   3.724             Condition No.:               1
    ================================================================================

Estimated coefficients are stored as a Pandas `Series` in the `params` attribute of `results`. Index values correspond to names of columns in the dependent variable `X`

```python
# Print the contents of the params method of results
print(results.params)
```

    tfp_cycle_lag    0.735887
    dtype: float64

```python
# Create a variable called 'rho' that equals the value of the estimated coefficient on df['tfp_cycle_lag']
rho = results.params['tfp_cycle_lag']
```

The residuals from the regression are stored in an attribute of `results` called `resid`

```python
# Create a variable called 'sigma' that equals the standard deviation of the residuals of the regression
sigma = results.resid.std()

# Print the value of sigma
print(sigma)
```

    0.00588482113457377

And that's how you estimate the parameters of an AR(1) model of the cyclical component of TFP for the US.

## Application: A Stochastic Solow Growth Model (Optional)

Consider the Solow growth model with written in "per worker" terms:

\begin{align}
y_t & = A_tk_t^{\alpha}\\
k_{t+1} & = i_t + (1-\delta)k_{t}\\
y_t & = c_t + i_t\\
c_t & = (1-s)y_t,
\end{align}

where $y_t$ is output per worker, $k_t$ is capital per worker, $c_t$ is consumption per worker, $i_t$ is investment per worker, and $A_t$ is time-varying TFP. Assume:

\begin{align}
\log A_{t+1} & = \rho \log A_{t} + \epsilon_{t+1}
\end{align}

where $\epsilon_t$ is a white noise process with mean 0 and variance $\sigma^2$. Since capital and TFP are the only two variables, that depend on past value of themselves, they are the two state variables of the model. Therefore, we can simulate simulate the model in two steps. First, simulate $k_t$ and $A_t$ using the following two laws of motion:

\begin{align}
k_{t+1} & = sAk_t^{\alpha} + (1-\delta)k_{t} \label{eqn:capital_solution}\\
\log A_{t+1} & = \rho \log A_{t} + \epsilon_{t+1}. \label{eqn:tfp_solution}
\end{align}

Second, compute simulated values for $y_t$, $c_t$, and $i_t$ using the following static relationships:

\begin{align}
y_t & = A_tk_t^{\alpha}\\
i_t & = sA_tk_t^{\alpha}\\
c_t & = (1-s)A_tk_t^{\alpha},
\end{align}

In the previous example, we estimated $\rho$ and $\sigma$ for the US so let's use those values for this simulation. For the other parameter, use the following values for the simulation:

| $A_0$ | $k_0$ | $s$ | $\alpha$ | $\delta $ | $T$  |
|-------|-------|-----|----------|-----------|------|
| 1     | 8.43  | 0.1 | 0.35     |  0.025    | 201  |

Where $T$ is the total number of simulation periods (i.e., $t$ will range from $0$ to $200$).

### Function

The function in the next cell simulates a stochastic Solow growth model. It returns a `DataFrame` with columns equal to the log-deviations of the simulated variables relative to the trend (i.e., nonstochastic steady state) implied by the model.

```python
# Define a function that returns a DataFrame of simulated value from the Solow model with exogenous labor and TFP growth. CELL PROVIDED
def solow_stochastic(s,alpha,delta,k0,A0,rho,sigma,T):

    '''Function for computing a simulation of the Solow growth model with exogenous labor and TFP growth.

            y[t]   = A[t]*K[t]^alpha
            k[t+1] = i[t] + (1-delta)*k[t]
            y[t]   = c[t] + i[t]
            c[t]   = (1-s)*y[t]
            A[t+1] = exp(rho*logA[t] + epsilon[t+1])

    Args:
        s (float):      Saving rate
        alpha (float):  Capital share in Cobb-Douglas production function
        delta (float):  Capital depreciation rate
        T (int):        Number of periods to simulate
        k0 (float):     Initial value of capital per worker
        A0 (float):     Initial TFP
        rho (float):    AR coeficient for log A[t]
        sigma (float):  Standard deviation of shock to log A[t]

     Returns:
         Pandas DataFrame
     '''
    # Create epsilon values
    epsilon = np.random.normal(scale=sigma,size=T)

    # Initialize capital values
    capital = np.zeros(T)

    # Set first value of capital equal to k0
    capital[0] = k0

    # Initialize TFP values
    tfp = np.zeros(T)

    # Set first value of TFP equal to A0
    tfp[0] = A0

    # Iterate over t in range(T-1) to update subsequent values in the capital and tfp arrays
    for t in range(T-1):
        capital[t+1] = s*tfp[t]*capital[t]**alpha + (1-delta)*capital[t]
        tfp[t+1] = np.exp(rho*np.log(tfp[t]) + epsilon[t])

    # Compute the values of the other aggregate variables
    output = tfp*capital**alpha
    consumption = (1-s)*output
    investment = s*output

    # Compute steady state (or trend) of endogenous vars
    capital_ss = (s/delta)**(1/(1-alpha))
    tfp_ss = 1
    output_ss = tfp_ss*capital_ss**alpha
    consumption_ss = (1-s)*output_ss
    investment_ss = s*output_ss

    # Put simulated data into a DataFrame
    df = pd.DataFrame({'output_log_dev':np.log(output/output_ss),
                       'consumption_log_dev':np.log(consumption/consumption_ss),
                       'investment_log_dev':np.log(investment/investment_ss),
                       'capital_log_dev':np.log(capital/capital_ss),
                       'tfp_log_dev':np.log(tfp/tfp_ss)})

    # Return the simulated data
    return df
```

### Simulate the Stochastic Solow Model

```python
# CELL PROVIDED
# Set parameters for simulation
alpha=0.35
s = 0.1
delta = 0.025
k0 = 8.43
A0=1
T = 201

np.random.seed(126)
# Simulate the model and store output in a variable called 'solow_df'
solow_df = solow_stochastic(s,alpha,delta,k0,A0,rho,sigma,T)

# Plot
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
solow_df.plot(ax=ax,grid=True)
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
```

### Simulation Results

Compute the standard deviation of the simulated data and the correlation coefficients of the simulated data.

```python
# Standard deviations. CELL PROVIDED
solow_df.std()*100
```

    output_log_dev         0.791095
    consumption_log_dev    0.791095
    investment_log_dev     0.791095
    capital_log_dev        0.173041
    tfp_log_dev            0.783475
    dtype: float64

```python
# Correlation coefficients. CELL PROVIDED
solow_df.corr()
```

<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>output_log_dev</th>
      <th>consumption_log_dev</th>
      <th>investment_log_dev</th>
      <th>capital_log_dev</th>
      <th>tfp_log_dev</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>output_log_dev</th>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>0.163488</td>
      <td>0.997088</td>
    </tr>
    <tr>
      <th>consumption_log_dev</th>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>0.163488</td>
      <td>0.997088</td>
    </tr>
    <tr>
      <th>investment_log_dev</th>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>1.000000</td>
      <td>0.163488</td>
      <td>0.997088</td>
    </tr>
    <tr>
      <th>capital_log_dev</th>
      <td>0.163488</td>
      <td>0.163488</td>
      <td>0.163488</td>
      <td>1.000000</td>
      <td>0.087776</td>
    </tr>
    <tr>
      <th>tfp_log_dev</th>
      <td>0.997088</td>
      <td>0.997088</td>
      <td>0.997088</td>
      <td>0.087776</td>
      <td>1.000000</td>
    </tr>
  </tbody>
</table>
</div>

It looks like the stochastic Solow model does a good job replicating the volatility (standard deviation) of ouput and consumption relative to the data. However the simulated investment volatility is too small by about a factor of 6. The stochastic Solow model does capture the correlation between output, consumption, and investment, but it implies perfect correlation which is too much.

## The Random Walk Process (Optional)

The *random walk process* is an AR(1) process with $\rho=1$:

\begin{align}
y_t = y_{t-1} + \epsilon_t
\end{align}

The random walk process has an important place in finance since the evidence suggests that stock prices follow a random walk process.

### Example

Simulate 7 random walk processes for 501 periods. Set $\sigma = 1$. Plot all 7 simulated processes on the same axes.

```python
# CELL PROVIDED
np.random.seed(126)
for i in range(7):
    plt.plot(ar1_sim(rho=1,T=501))

plt.title('Seven random walk processes')
plt.grid()
```


text:
# Table of Contents

1  TF intro

2  MLP on MNIST

3  VAEs!

# TF intro

This guide assumes you can read through basic Python code or use your google skills to catch up on that as needed. We begin by understanding how tensorflow works. The key point to remember is that all the tensorflow computation happens in a graph, and all that you get to do in python is to manipulate and run that graph. This creates a programming paradigm that looks a lot like python but is actually quite different. We begin with a simple logistic regression example.

```python
import tensorflow as tf
import numpy as np

n_obs = 1000
n_features = 5

x_ph = tf.placeholder(shape=(n_obs, n_features), name="x_ph", dtype=tf.float32)
beta_init = np.random.normal(size=(n_features, 1))
beta_hat = tf.Variable(beta_init, dtype=tf.float32, name="beta_hat")
y_hat = tf.nn.sigmoid(tf.matmul(x_ph, beta_hat), name="yhat")
```

We can visualize thie graph, and then explain what the different parts are:

This code is a modification from the DeepDream notebook. There is more visualization / exploration that can be done using tensorboard, which is the tool this uses.

```python
# TensorFlow Graph visualizer code
# https://stackoverflow.com/questions/41388673/visualizing-a-tensorflow-graph-in-jupyter-doesnt-work
import numpy as np
from IPython.display import clear_output, Image, display, HTML

def strip_consts(graph_def, max_const_size=32):
    """Strip large constant values from graph_def."""
    strip_def = tf.GraphDef()
    for n0 in graph_def.node:
        n = strip_def.node.add()
        n.MergeFrom(n0)
        if n.op == 'Const':
            tensor = n.attr['value'].tensor
            size = len(tensor.tensor_content)
            if size > max_const_size:
                tensor.tensor_content = bytes("<stripped %d bytes>"%size, 'utf-8')
    return strip_def

def show_graph(graph_def, max_const_size=32):
    """Visualize TensorFlow graph."""
    if hasattr(graph_def, 'as_graph_def'):
        graph_def = graph_def.as_graph_def()
    strip_def = strip_consts(graph_def, max_const_size=max_const_size)
    code = """

        <link rel="import" href="https://tensorboard.appspot.com/tf-graph-basic.build.html" onload=load()>
        <div style="height:800px">
          <tf-graph-basic id="{id}"></tf-graph-basic>
        </div>
    """.format(data=repr(str(strip_def)), id='graph'+str(np.random.rand()))

    iframe = """

    """.format(code.replace('"', '"'))
    display(HTML(iframe))
```

```python
show_graph(tf.get_default_graph().as_graph_def())
```

What we do in tensorflow is construct graphs like this, and then evaluate nodes. Each graph node is associated with some code. When we evaluate a node like `y_hat`, tensorflow figures out what nodes it depends on, evaluates all of those nodes, and then evaluates `y_hat`. In this graph, there are three types of nodes (`Tensor`s), a `Variable`, a `placeholder`, and a vanilla `Tensor` that is none of the above. Tensors have no state: they are computable from the rest of the graph. Variables have state (they're the only thing we can optimize / save / load). Placeholders are not computable and don't have state: we must feed values into them. We can also feed values into other tensors, but TF will explicitly complain if we fail to feed a value into a placeholder.

Our python objects like `beta_hat` are references to the TF graph nodes, not the nodes themselves (i.e. copying the python object does not dupliate the graph node).

To evaluate a graph we need to associate it with a session, and then either a tensor's `eval` method or the session's `run` method. The difference is that we can `run` multiple tensors together, which might be useful if they share dependencies.

Now we define the loss, generate synth data, and optimize:

```python
from scipy.special import expit as logistic

true_beta = np.random.normal(size=(n_features, 1))
x = np.random.normal(size=(n_obs, n_features))
y = np.random.binomial(n=1, p=logistic(x @ true_beta))

y_ph = tf.placeholder(shape=(n_obs, 1), name="y_ph", dtype=tf.float32)

logistic_loss = -tf.reduce_sum(y_ph * tf.log(1e-10 + y_hat)+ (1-y_ph) * tf.log(1e-10 + 1 - y_hat))

# if we needed the gradients for some reason, e.g. to pass to an external optimizier or to plot
# grads = tf.gradients(logistic_loss, beta_hat)

# create the optimizer
optimizer = tf.train.GradientDescentOptimizer(learning_rate=0.0005)

# this is the long way of doing it
# grads_and_vars = optimizer.compute_gradients(logistic_loss, beta_hat)
# optionally, modify gradients here (e.g. threshold)
# minimize_op = optimizer.apply_gradients(grads_and_vars)

# the short way
minimize_op = optimizer.minimize(logistic_loss) # by default, minimize w.r.t. all variables

sess = tf.Session()
sess.run(tf.global_variables_initializer())
for i in range(2500):
    if i % 100 == 0:
        l = logistic_loss.eval(session=sess, feed_dict={x_ph:x, y_ph:y})
        print("Iter %i, loss %f" % (i, l))
    sess.run(minimize_op, feed_dict={x_ph:x, y_ph:y})

accuracy = np.average(tf.equal(tf.round(y_hat), y_ph).eval(session=sess, feed_dict={x_ph:x, y_ph:y}))
print(accuracy)
```

    Iter 0, loss 1257.209106
    Iter 100, loss 425.954163
    Iter 200, loss 425.796875
    Iter 300, loss 425.795349
    Iter 400, loss 425.795319
    Iter 500, loss 425.795227
    Iter 600, loss 425.795227
    Iter 700, loss 425.795227
    Iter 800, loss 425.795227
    Iter 900, loss 425.795227
    Iter 1000, loss 425.795227
    Iter 1100, loss 425.795227
    Iter 1200, loss 425.795227
    Iter 1300, loss 425.795227
    Iter 1400, loss 425.795227
    Iter 1500, loss 425.795227
    Iter 1600, loss 425.795227
    Iter 1700, loss 425.795227
    Iter 1800, loss 425.795227
    Iter 1900, loss 425.795227
    Iter 2000, loss 425.795227
    Iter 2100, loss 425.795227
    Iter 2200, loss 425.795227
    Iter 2300, loss 425.795227
    Iter 2400, loss 425.795227
    0.814

```python
# another optimizer
optimizer = tf.train.AdamOptimizer(learning_rate=0.01)
minimize_op = optimizer.minimize(logistic_loss) # by default, minimize w.r.t. all variables

sess.run(tf.global_variables_initializer())
for i in range(1000):
    if i % 100 == 0:
        l = logistic_loss.eval(session=sess, feed_dict={x_ph:x, y_ph:y})
        print("Iter %i, loss %f" % (i, l))
    sess.run(minimize_op, feed_dict={x_ph:x, y_ph:y})

accuracy = np.average(tf.equal(tf.round(y_hat), y_ph).eval(session=sess, feed_dict={x_ph:x, y_ph:y}))
print(accuracy)
```

    Iter 0, loss 1257.209106
    Iter 100, loss 621.341553
    Iter 200, loss 453.863708
    Iter 300, loss 430.841217
    Iter 400, loss 426.745361
    Iter 500, loss 425.957336
    Iter 600, loss 425.818573
    Iter 700, loss 425.798065
    Iter 800, loss 425.795532
    Iter 900, loss 425.795319
    0.814

One more tutorial point is on how to print things. Since a tensor only has a value when the graph is executed, inspecting things is trickier than usual. The `Print` op returns the same op as its input, but prints as a side-effect. This means we need to inject the op into the graph. Unfortunately, the print happens on the C++ end, so we will see the logging messages in the jupyter server log (or in a shell).

```python
logistic_loss_with_print = tf.Print(input_=logistic_loss, data=[x, logistic_loss])

_ = logistic_loss_with_print.eval(session=sess, feed_dict={x_ph:x, y_ph:y})
```

# MLP on MNIST

Now we build some neural network building blocks we will reuse for VAEs.

```python
from tensorflow.examples.tutorials.mnist import input_data

global_dtype = tf.float32
mnist = input_data.read_data_sets('MNIST_data', one_hot=True)
input_size = mnist.train.images.shape[1]
```

    Successfully downloaded train-images-idx3-ubyte.gz 9912422 bytes.
    Extracting MNIST_data/train-images-idx3-ubyte.gz
    Successfully downloaded train-labels-idx1-ubyte.gz 28881 bytes.
    Extracting MNIST_data/train-labels-idx1-ubyte.gz
    Successfully downloaded t10k-images-idx3-ubyte.gz 1648877 bytes.
    Extracting MNIST_data/t10k-images-idx3-ubyte.gz
    Successfully downloaded t10k-labels-idx1-ubyte.gz 4542 bytes.
    Extracting MNIST_data/t10k-labels-idx1-ubyte.gz

Define our neural network building blocks.

```python
def _dense_mlp_layer(x, input_size, out_size, nonlinearity=tf.nn.softmax, name_prefix=""):
    w_init = tf.truncated_normal(shape=[input_size, out_size], stddev=0.001)
    b_init = tf.ones(shape=[out_size]) * 0.1
    W = tf.Variable(w_init, name="%s_W" % name_prefix)
    b = tf.Variable(b_init, name="%s_b" % name_prefix)
    out = nonlinearity(tf.matmul(x, W) + b)
    return out, [W, b]

def _mlp(x, n_layers, units_per_layer, input_size, out_size, nonlinearity=tf.tanh):
    train_vars = []

    x, v = _dense_mlp_layer(x, input_size, units_per_layer, nonlinearity, name_prefix="into_hidden")
    train_vars.extend(v)
    # exploit the fact that repeatedly calling the same TF function creates multiple ops.
    # no need to hang onto the intermediate layer handles (though we can get them back if we need them)
    for l in range(n_layers-1):
        x, v = _dense_mlp_layer(x, units_per_layer, units_per_layer, nonlinearity, name_prefix="hidden")
        train_vars.extend(v)

    x, v = _dense_mlp_layer(x, units_per_layer, out_size, nonlinearity, name_prefix="readout")
    train_vars.extend(v)
    return x, train_vars
```

Now we construct the graph. The graph and scope boilerplate makes our life easier as far as visualization and debugging is concerned. We can visualize/run only this graph and not the graph for logistic regression (above).

```python
mlp_graph = tf.Graph()

with mlp_graph.as_default():
    with tf.name_scope("Feedforward_Net"):
        x = tf.placeholder(shape=[None, input_size], dtype=global_dtype, name='x')
        y = tf.placeholder(shape=[None, 10], dtype=global_dtype, name='y')
        y_hat, mlp_test_vars = _mlp(x, n_layers=2, units_per_layer=30, input_size=784, out_size=10)

    with tf.name_scope("Opt_and_loss"):
        cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(labels=y, logits=y_hat))
        # learning rates are much smaller for optimizers like Adam and RMSProp
        train_step_mlptest = tf.train.AdamOptimizer(0.01).minimize(cross_entropy)

    with tf.name_scope("Support_stuff"):
        init = tf.global_variables_initializer()
        correct_prediction = tf.equal(tf.argmax(y,1), tf.argmax(y_hat,1))
        mlp_acc = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
```

We visualize our graph:

```python
show_graph(mlp_graph.as_graph_def())
```

Next we create a session, initialize our variables, and train the network:

```python
sess = tf.Session(graph=mlp_graph)
sess.run(init)
train_steps = 2500

acc = np.zeros(train_steps)

# create this op outside of the loop so we don't create it 5000 times
for i in range(train_steps):
    batch_xs, batch_ys = mnist.train.next_batch(100)
    acc[i] = mlp_acc.eval(session = sess, feed_dict = {x: batch_xs, y: batch_ys})
    sess.run(train_step_mlptest, feed_dict={x: batch_xs, y: batch_ys})

test_acc = mlp_acc.eval(session=sess, feed_dict={x: mnist.test.images, y: mnist.test.labels})
print("Test accuracy: %f" % test_acc)

import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(acc)
```

# VAEs!
$$
\DeclareMathOperator{\Tr}{Tr}
\newcommand{\trp}{{^\top}} % transpose
\newcommand{\trace}{\text{Trace}} % trace
\newcommand{\inv}{^{-1}}
\newcommand{\mb}{\mathbf{b}}
\newcommand{\M}{\mathbf{M}}
\newcommand{\G}{\mathbf{G}}
\newcommand{\A}{\mathbf{A}}
\newcommand{\R}{\mathbf{R}}
\renewcommand{\S}{\mathbf{S}}
\newcommand{\B}{\mathbf{B}}
\newcommand{\Q}{\mathbf{Q}}
\newcommand{\mH}{\mathbf{H}}
\newcommand{\U}{\mathbf{U}}
\newcommand{\mL}{\mathbf{L}}
\newcommand{\diag}{\mathrm{diag}}
\newcommand{\etr}{\mathrm{etr}}
\renewcommand{\H}{\mathbf{H}}
\newcommand{\vecop}{\mathrm{vec}}
\newcommand{\I}{\mathbf{I}}
\newcommand{\X}{\mathbf{X}_{ij}}
\newcommand{\Y}{\mathbf{Y}_{jk}}
\newcommand{\Z}{\mathbf{Z}_{ik}}
$$

In our generative model, we would like to estimate the density of some complicated probability density $\log P_{\theta}(X)$. This is slightly odd notation but seems standard in these papers: it says that $X$ is a random variable but $\theta$ are parameters. To do this, we will write it as follows:

$$
\log P_{\theta}(X) =  \log P_{\theta}(X,Z) - \log P_{\theta}(Z\mid X)
$$

This is just the definition of marginal probability. Next, we add/subtract $\log Q_{\phi}(Z\mid X)$ which sums to 0:

$$
\log P_{\theta}(X) =  \log P_{\theta}(X,Z) - \log P_{\theta}(Z\mid X) - \log Q_{\phi}(Z\mid X) + \log Q_{\phi}(Z\mid X)
$$

We take expectation of both sides. Note that this expectation has to be w.r.t. the conditional distribution $Q(Z\mid X)$ for the rest of this to work properly:

$$
\mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X) =  \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X,Z) - \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(Z\mid X) - \mathbb{E}_{Q(Z\mid X)}\log Q_{\phi}(Z\mid X) + \mathbb{E}_{Q(Z\mid X)}\log Q_{\phi}(Z\mid X)
$$

Since the LHS is independent of Z, that expectation just goes away. We rearrange the terms and realize we ended up with the evidence lower bound (ELBO) and a KL divergence:

$$
\begin{align}
\log P_{\theta}(X) =&  \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X,Z) - \mathbb{E}_{Q(Z\mid X)}\log Q_{\phi}(Z\mid X) + \mathbb{E}_{Q(Z\mid X)}\log Q_{\phi}(Z\mid X)- \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(Z\mid X) \\
  =&  \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X,Z) - \mathbb{E}_{Q(Z\mid X)}\log Q_{\phi}(Z\mid X) + \int_Z Q_{\phi}(Z\mid X)\log Q_{\phi}(Z\mid X)- \int_Z Q_{\phi}(Z\mid X)\log P_{\theta}(Z\mid X) \\
=&  \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X,Z) +\mathbb{E}_{Q(Z\mid X)}\log Q_{\phi}(Z\mid X) - \int_Z Q_{\phi}(Z\mid X)\log \frac{Q_{\phi}(Z\mid X)}{P_{\theta}(Z\mid X)} \\
=&  \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X,Z) + \mathbb{E}_{Q(Z\mid X)}\log Q_{\phi}(Z\mid X) -\mathcal{D}_{KL}(Q_{\theta}(Z\mid X)||  P_{\phi}(Z\mid X))
\end{align}
$$

The ELBO is a lower bound because KL divergence is greater than equal to 0. So if we want to maximize the LHS, we can choose to either minimize the KL divergence, or maximize the ELBO.  We would rather do the former, because the KL divergence contains the posterior $P(Z\mid X)$ and if we knew how to do that, we wouldn't be going through this hassle. The nice thing is, since this holds for any $Q$ and any $Z$ we can define both distributions to be as nice as we like. So we're going to say that the likelihood $P_{\theta}(X,Z)$, prior $P_{\theta}(Z)$ and approximate posterior $Q_{\theta}(X,Z)$ are all gaussian. We pick the easiest possible marginal distribution over $Q(Z)$, an identity-covariance gaussian:

$$
\begin{align}
P_{\theta}(X\mid Z) :=& \mathcal{N}(a(Z,\theta), b(Z,\theta)b(Z,\theta)\trp) \\
Q_{\phi}(Z\mid X):=&\mathcal{N}(f(X,\phi),g(X,\phi)g(X,\phi)^{\intercal})\\
P_{\theta}(Z) :=& \mathcal{N}(0,\mathbf{I})
\end{align}
$$

The distributions $P_{\theta}(X\mid Z)$ and $Q_{\phi}(Z\mid X)$ are parameterized by mean and covariance-square-root functions $a,b,f,g$ which we leave unspecified for now. In practice for VAEs people use the equivalent of a mean-field assumption, which means that the covariance functions will just return SDs/variances, but I'd like to write the reparameterization trick in general form. We can additionally rewrite the expression to get a second KL divergence:

$$
\begin{align}
\log P_{\theta}(X) =& \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X\mid Z) + \mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(Z) + \mathbb{E}_{Q(Z\mid X)}\log Q_{\phi}(Z\mid X)-\mathcal{D}_{KL}(Q_{\theta}(Z\mid X)||  P_{\phi}(Z\mid X)) \\
=&\mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X\mid Z) + \mathcal{D}_{KL}( Q_{\phi}(Z\mid X)||P_{\theta}(Z)) - \mathcal{D}_{KL}(Q_{\theta}(Z\mid X)||  P_{\phi}(Z\mid X))
\end{align}
$$

Conveniently, the KL divergence between two gaussians (the prior and approximate posterior) is analytic. What remains is the expectation $\mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X\mid Z)$, which we can compute from its empirical, sample-based mean:

$$
\mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X\mid Z) = \int Q_{\phi}(Z\mid X)\log P_{\theta}(X|Z)d Q \approx \frac{1}{N} \sum \log P_{\theta}(X|Z) Z_i, \\\
Z_i\sim\mathcal{N}(f(X,\phi),g(X,\phi)g(X,\phi)\trp)
$$

The naive gradient estimator here has very high variance according to the VAE paper, though it is used in Blei, Jordan and Paisley 2012 (ICML):

$$
\nabla_{\phi}\mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X\mid Z) \approx \frac{1}{N} \sum \nabla_{\phi}\log P_{\theta}(X|Z) Z_i
$$

What the VAE paper does instead is apply the reparameterization trick:

\begin{align}
\epsilon&\sim\mathcal{N}(0, \I)\\
Z &= f(X,\phi) + g(X,\phi)\epsilon\\
\mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X\mid Z)  &\approx \frac{1}{N} \sum \log P_{\theta}(X|Z)(f(X,\phi) + g(X,\phi)\epsilon)\\
\nabla_{\phi}\mathbb{E}_{Q(Z\mid X)}\log P_{\theta}(X\mid Z)  &\approx \frac{1}{N} \sum \nabla_{\phi}\log P_{\theta}(X|Z)(f(X,\phi) + g(X,\phi)\epsilon)
\end{align}

Now let's see if we can implement it using tensorflow. We begin by sanity-checking a basic MLP and then go to VAEs. First, import some things we'll need and download MNIST:

```python
from tensorflow.examples.tutorials.mnist import input_data

encoder_depth = 2
decoder_depth = 2
encoder_units = 500
decoder_units = 500
latent_size = 10
global_dtype = tf.float32
minibatch_size = 100
input_size = mnist.train.images.shape[1]
train_steps = mnist.train.num_examples // minibatch_size
encoder_nonlinearity = tf.nn.sigmoid
decoder_nonlinearity = tf.nn.sigmoid
n_epochs = 10
```

Construct the graph:

```python
vae_graph = tf.Graph()

with vae_graph.as_default():

    with tf.name_scope("Encoder_Q"):

        x = tf.placeholder(shape=[None, input_size], dtype=global_dtype, name='x')
        q_network, q_mu_vars = _mlp(x, n_layers=encoder_depth, units_per_layer=encoder_units, input_size=input_size, out_size=encoder_units, nonlinearity=encoder_nonlinearity)
        w_mu = tf.Variable(tf.truncated_normal(shape=[encoder_units, latent_size], stddev=0.1), name="w_mu")
        w_logsig = tf.Variable(tf.truncated_normal(shape=[encoder_units, latent_size], stddev=0.1), name="w_logsig")
        b_mu = tf.Variable(tf.truncated_normal(shape=[latent_size], stddev=0.1), name="b_mu")
        b_logsig = tf.Variable(tf.truncated_normal(shape=[latent_size], stddev=0.1), name="b_logsig")
        q_mu = tf.matmul(q_network, w_mu) + b_mu
        q_logsigma = tf.matmul(q_network, w_logsig) + b_logsig
        epsilon = tf.random_normal([minibatch_size, latent_size])
        z = q_mu + tf.sqrt(tf.exp(q_logsigma)) * epsilon

    with tf.name_scope("Decoder_P"):

        p_mu, p_mu_vars = _mlp(z, n_layers=decoder_depth, units_per_layer=decoder_units, input_size=latent_size,  out_size=input_size, nonlinearity=decoder_nonlinearity)

    with tf.name_scope("Opt_and_loss"):
        kld = 0.5 * tf.reduce_sum(1 + q_logsigma - tf.square(q_mu) - tf.exp(q_logsigma), 1)
        ll = tf.reduce_sum(x * tf.log(1e-10 + p_mu)+ (1-x) * tf.log(1e-10 + 1 - p_mu), 1)
        elbo = tf.reduce_mean(ll + kld)
        minimize_op = tf.train.AdamOptimizer(0.001).minimize(-elbo)

    with tf.name_scope("Support_stuff"):
        init = tf.global_variables_initializer()

```

```python
show_graph(vae_graph.as_graph_def())
```

Now we run and visualize:

```python
sess = tf.Session(graph=vae_graph)
sess.run(init)

elbo_log = np.zeros(n_epochs * train_steps)

for i in range(n_epochs):
    for j in range(train_steps):
        batch_xs, batch_ys = mnist.train.next_batch(minibatch_size)
        sess.run(minimize_op, feed_dict={x: batch_xs})
        elbo_log[i*train_steps + j] = elbo.eval(session=sess, feed_dict={x: batch_xs})
        if j % 10 == 0:
            print("Epoch %i, step %i: average elbo=%f" % (i, j, elbo_log[i*train_steps + j]))
```

    Epoch 0, step 0: average elbo=-498.922028
    Epoch 0, step 10: average elbo=-214.074295
    Epoch 0, step 20: average elbo=-214.997269
    Epoch 0, step 30: average elbo=-212.165344
    Epoch 0, step 40: average elbo=-210.181015
    Epoch 0, step 50: average elbo=-213.656189
    Epoch 0, step 60: average elbo=-207.370346
    Epoch 0, step 70: average elbo=-203.687546
    Epoch 0, step 80: average elbo=-210.040588
    Epoch 0, step 90: average elbo=-207.408844
    Epoch 0, step 100: average elbo=-200.401810
    Epoch 0, step 110: average elbo=-211.652115
    Epoch 0, step 120: average elbo=-203.729889
    Epoch 0, step 130: average elbo=-198.393951
    Epoch 0, step 140: average elbo=-206.218460
    Epoch 0, step 150: average elbo=-207.235901
    Epoch 0, step 160: average elbo=-209.720154
    Epoch 0, step 170: average elbo=-217.805222
    Epoch 0, step 180: average elbo=-207.267563
    Epoch 0, step 190: average elbo=-205.444199
    Epoch 0, step 200: average elbo=-203.058441
    Epoch 0, step 210: average elbo=-211.681366
    Epoch 0, step 220: average elbo=-213.850220
    Epoch 0, step 230: average elbo=-210.270935
    Epoch 0, step 240: average elbo=-204.285355
    Epoch 0, step 250: average elbo=-202.927597
    Epoch 0, step 260: average elbo=-205.129791
    Epoch 0, step 270: average elbo=-197.652893
    Epoch 0, step 280: average elbo=-205.535339
    Epoch 0, step 290: average elbo=-192.104752
    Epoch 0, step 300: average elbo=-202.904068
    Epoch 0, step 310: average elbo=-202.045944
    Epoch 0, step 320: average elbo=-206.704285
    Epoch 0, step 330: average elbo=-219.215546
    Epoch 0, step 340: average elbo=-215.538666
    Epoch 0, step 350: average elbo=-205.286438
    Epoch 0, step 360: average elbo=-201.292419
    Epoch 0, step 370: average elbo=-206.491440
    Epoch 0, step 380: average elbo=-207.749573
    Epoch 0, step 390: average elbo=-206.005707
    Epoch 0, step 400: average elbo=-203.475937
    Epoch 0, step 410: average elbo=-195.294220
    Epoch 0, step 420: average elbo=-200.652267
    Epoch 0, step 430: average elbo=-208.812149
    Epoch 0, step 440: average elbo=-200.045319
    Epoch 0, step 450: average elbo=-199.878784
    Epoch 0, step 460: average elbo=-206.345490
    Epoch 0, step 470: average elbo=-207.140121
    Epoch 0, step 480: average elbo=-204.130676
    Epoch 0, step 490: average elbo=-203.469315
    Epoch 0, step 500: average elbo=-205.958557
    Epoch 0, step 510: average elbo=-203.972565
    Epoch 0, step 520: average elbo=-201.650574
    Epoch 0, step 530: average elbo=-204.986053
    Epoch 0, step 540: average elbo=-198.621933
    Epoch 1, step 0: average elbo=-202.771149
    Epoch 1, step 10: average elbo=-201.271362
    Epoch 1, step 20: average elbo=-195.510757
    Epoch 1, step 30: average elbo=-209.701096
    Epoch 1, step 40: average elbo=-198.979568
    Epoch 1, step 50: average elbo=-198.095337
    Epoch 1, step 60: average elbo=-210.112228
    Epoch 1, step 70: average elbo=-201.188522
    Epoch 1, step 80: average elbo=-205.822937
    Epoch 1, step 90: average elbo=-198.407654
    Epoch 1, step 100: average elbo=-193.962448
    Epoch 1, step 110: average elbo=-194.385345
    Epoch 1, step 120: average elbo=-204.923477
    Epoch 1, step 130: average elbo=-208.576111
    Epoch 1, step 140: average elbo=-204.448181
    Epoch 1, step 150: average elbo=-204.872742
    Epoch 1, step 160: average elbo=-200.421188
    Epoch 1, step 170: average elbo=-198.152557
    Epoch 1, step 180: average elbo=-200.047562
    Epoch 1, step 190: average elbo=-206.334961
    Epoch 1, step 200: average elbo=-202.665192
    Epoch 1, step 210: average elbo=-204.729080
    Epoch 1, step 220: average elbo=-197.487152
    Epoch 1, step 230: average elbo=-196.357513
    Epoch 1, step 240: average elbo=-207.793549
    Epoch 1, step 250: average elbo=-209.511246
    Epoch 1, step 260: average elbo=-202.061935
    Epoch 1, step 270: average elbo=-197.077194
    Epoch 1, step 280: average elbo=-195.298264
    Epoch 1, step 290: average elbo=-194.038712
    Epoch 1, step 300: average elbo=-199.171677
    Epoch 1, step 310: average elbo=-188.342758
    Epoch 1, step 320: average elbo=-203.318146
    Epoch 1, step 330: average elbo=-196.126175
    Epoch 1, step 340: average elbo=-191.844864
    Epoch 1, step 350: average elbo=-193.601959
    Epoch 1, step 360: average elbo=-201.827423
    Epoch 1, step 370: average elbo=-188.135818
    Epoch 1, step 380: average elbo=-197.268051
    Epoch 1, step 390: average elbo=-192.352890
    Epoch 1, step 400: average elbo=-202.524338
    Epoch 1, step 410: average elbo=-195.991409
    Epoch 1, step 420: average elbo=-198.012756
    Epoch 1, step 430: average elbo=-197.807190
    Epoch 1, step 440: average elbo=-203.761795
    Epoch 1, step 450: average elbo=-195.336060
    Epoch 1, step 460: average elbo=-195.639923
    Epoch 1, step 470: average elbo=-197.546234
    Epoch 1, step 480: average elbo=-194.885941
    Epoch 1, step 490: average elbo=-186.946213
    Epoch 1, step 500: average elbo=-200.764420
    Epoch 1, step 510: average elbo=-201.621597
    Epoch 1, step 520: average elbo=-197.619278
    Epoch 1, step 530: average elbo=-190.513000
    Epoch 1, step 540: average elbo=-202.596405
    Epoch 2, step 0: average elbo=-196.629211
    Epoch 2, step 10: average elbo=-201.397049
    Epoch 2, step 20: average elbo=-196.755905
    Epoch 2, step 30: average elbo=-190.376953
    Epoch 2, step 40: average elbo=-191.373810
    Epoch 2, step 50: average elbo=-197.041595
    Epoch 2, step 60: average elbo=-196.976074
    Epoch 2, step 70: average elbo=-196.638824
    Epoch 2, step 80: average elbo=-193.106171
    Epoch 2, step 90: average elbo=-194.680786
    Epoch 2, step 100: average elbo=-188.465073
    Epoch 2, step 110: average elbo=-194.780670
    Epoch 2, step 120: average elbo=-191.621445
    Epoch 2, step 130: average elbo=-188.459686
    Epoch 2, step 140: average elbo=-205.755707
    Epoch 2, step 150: average elbo=-189.846405
    Epoch 2, step 160: average elbo=-200.592209
    Epoch 2, step 170: average elbo=-195.979645
    Epoch 2, step 180: average elbo=-187.060059
    Epoch 2, step 190: average elbo=-190.503784
    Epoch 2, step 200: average elbo=-188.892029
    Epoch 2, step 210: average elbo=-194.118698
    Epoch 2, step 220: average elbo=-193.117111
    Epoch 2, step 230: average elbo=-188.780701
    Epoch 2, step 240: average elbo=-195.566757
    Epoch 2, step 250: average elbo=-188.729858
    Epoch 2, step 260: average elbo=-197.187729
    Epoch 2, step 270: average elbo=-188.502106
    Epoch 2, step 280: average elbo=-194.004059
    Epoch 2, step 290: average elbo=-191.728455
    Epoch 2, step 300: average elbo=-195.860443
    Epoch 2, step 310: average elbo=-197.077072
    Epoch 2, step 320: average elbo=-188.211090
    Epoch 2, step 330: average elbo=-184.821442
    Epoch 2, step 340: average elbo=-184.030563
    Epoch 2, step 350: average elbo=-194.977005
    Epoch 2, step 360: average elbo=-202.008591
    Epoch 2, step 370: average elbo=-191.488541
    Epoch 2, step 380: average elbo=-188.002853
    Epoch 2, step 390: average elbo=-193.981567
    Epoch 2, step 400: average elbo=-184.338333
    Epoch 2, step 410: average elbo=-184.117584
    Epoch 2, step 420: average elbo=-190.043304
    Epoch 2, step 430: average elbo=-192.296600
    Epoch 2, step 440: average elbo=-189.440826
    Epoch 2, step 450: average elbo=-194.788101
    Epoch 2, step 460: average elbo=-193.212387
    Epoch 2, step 470: average elbo=-189.541321
    Epoch 2, step 480: average elbo=-189.744705
    Epoch 2, step 490: average elbo=-181.076370
    Epoch 2, step 500: average elbo=-192.771469
    Epoch 2, step 510: average elbo=-191.602188
    Epoch 2, step 520: average elbo=-189.750305
    Epoch 2, step 530: average elbo=-195.219528
    Epoch 2, step 540: average elbo=-186.222946
    Epoch 3, step 0: average elbo=-199.638092
    Epoch 3, step 10: average elbo=-188.312546
    Epoch 3, step 20: average elbo=-192.338470
    Epoch 3, step 30: average elbo=-185.960022
    Epoch 3, step 40: average elbo=-192.097534
    Epoch 3, step 50: average elbo=-192.343826
    Epoch 3, step 60: average elbo=-183.228088
    Epoch 3, step 70: average elbo=-200.525116
    Epoch 3, step 80: average elbo=-183.065903
    Epoch 3, step 90: average elbo=-192.388733
    Epoch 3, step 100: average elbo=-197.787506
    Epoch 3, step 110: average elbo=-183.872620
    Epoch 3, step 120: average elbo=-183.707382
    Epoch 3, step 130: average elbo=-186.933151
    Epoch 3, step 140: average elbo=-184.939850
    Epoch 3, step 150: average elbo=-190.825287
    Epoch 3, step 160: average elbo=-184.756760
    Epoch 3, step 170: average elbo=-188.678238
    Epoch 3, step 180: average elbo=-191.503235
    Epoch 3, step 190: average elbo=-178.095306
    Epoch 3, step 200: average elbo=-183.801208
    Epoch 3, step 210: average elbo=-179.513519
    Epoch 3, step 220: average elbo=-177.637024
    Epoch 3, step 230: average elbo=-186.102478
    Epoch 3, step 240: average elbo=-187.571564
    Epoch 3, step 250: average elbo=-183.121918
    Epoch 3, step 260: average elbo=-185.857697
    Epoch 3, step 270: average elbo=-172.899399
    Epoch 3, step 280: average elbo=-172.112228
    Epoch 3, step 290: average elbo=-179.014175
    Epoch 3, step 300: average elbo=-180.813889
    Epoch 3, step 310: average elbo=-173.349045
    Epoch 3, step 320: average elbo=-182.294907
    Epoch 3, step 330: average elbo=-180.487030
    Epoch 3, step 340: average elbo=-176.551819
    Epoch 3, step 350: average elbo=-177.762909
    Epoch 3, step 360: average elbo=-174.468552
    Epoch 3, step 370: average elbo=-177.856094
    Epoch 3, step 380: average elbo=-175.555847
    Epoch 3, step 390: average elbo=-182.225510
    Epoch 3, step 400: average elbo=-173.955200
    Epoch 3, step 410: average elbo=-178.312073
    Epoch 3, step 420: average elbo=-169.398071
    Epoch 3, step 430: average elbo=-173.134094
    Epoch 3, step 440: average elbo=-181.403122
    Epoch 3, step 450: average elbo=-184.885941
    Epoch 3, step 460: average elbo=-173.984009
    Epoch 3, step 470: average elbo=-175.292831
    Epoch 3, step 480: average elbo=-171.830017
    Epoch 3, step 490: average elbo=-168.667892
    Epoch 3, step 500: average elbo=-179.198853
    Epoch 3, step 510: average elbo=-175.162155
    Epoch 3, step 520: average elbo=-183.335052
    Epoch 3, step 530: average elbo=-174.732819
    Epoch 3, step 540: average elbo=-172.761948
    Epoch 4, step 0: average elbo=-169.987076
    Epoch 4, step 10: average elbo=-179.686371
    Epoch 4, step 20: average elbo=-181.981567
    Epoch 4, step 30: average elbo=-170.259186
    Epoch 4, step 40: average elbo=-175.775681
    Epoch 4, step 50: average elbo=-167.100113
    Epoch 4, step 60: average elbo=-173.629333
    Epoch 4, step 70: average elbo=-172.521881
    Epoch 4, step 80: average elbo=-165.267502
    Epoch 4, step 90: average elbo=-171.267242
    Epoch 4, step 100: average elbo=-174.288712
    Epoch 4, step 110: average elbo=-166.835739
    Epoch 4, step 120: average elbo=-166.294037
    Epoch 4, step 130: average elbo=-169.909088
    Epoch 4, step 140: average elbo=-177.945663
    Epoch 4, step 150: average elbo=-176.419724
    Epoch 4, step 160: average elbo=-167.970230
    Epoch 4, step 170: average elbo=-175.322250
    Epoch 4, step 180: average elbo=-174.118790
    Epoch 4, step 190: average elbo=-182.828949
    Epoch 4, step 200: average elbo=-165.007339
    Epoch 4, step 210: average elbo=-178.020233
    Epoch 4, step 220: average elbo=-170.333313
    Epoch 4, step 230: average elbo=-170.001465
    Epoch 4, step 240: average elbo=-171.660828
    Epoch 4, step 250: average elbo=-165.533676
    Epoch 4, step 260: average elbo=-175.753433
    Epoch 4, step 270: average elbo=-159.924728
    Epoch 4, step 280: average elbo=-167.253693
    Epoch 4, step 290: average elbo=-168.621872
    Epoch 4, step 300: average elbo=-166.031387
    Epoch 4, step 310: average elbo=-166.500565
    Epoch 4, step 320: average elbo=-167.965042
    Epoch 4, step 330: average elbo=-168.473206
    Epoch 4, step 340: average elbo=-169.340530
    Epoch 4, step 350: average elbo=-169.905426
    Epoch 4, step 360: average elbo=-169.395905
    Epoch 4, step 370: average elbo=-172.168900
    Epoch 4, step 380: average elbo=-164.849762
    Epoch 4, step 390: average elbo=-166.733353
    Epoch 4, step 400: average elbo=-166.215134
    Epoch 4, step 410: average elbo=-162.192719
    Epoch 4, step 420: average elbo=-162.633209
    Epoch 4, step 430: average elbo=-164.450668
    Epoch 4, step 440: average elbo=-171.088303
    Epoch 4, step 450: average elbo=-153.585678
    Epoch 4, step 460: average elbo=-174.374100
    Epoch 4, step 470: average elbo=-167.306763
    Epoch 4, step 480: average elbo=-164.465836
    Epoch 4, step 490: average elbo=-164.340775
    Epoch 4, step 500: average elbo=-166.516174
    Epoch 4, step 510: average elbo=-165.074493
    Epoch 4, step 520: average elbo=-161.835327
    Epoch 4, step 530: average elbo=-159.028549
    Epoch 4, step 540: average elbo=-158.972260
    Epoch 5, step 0: average elbo=-161.663269
    Epoch 5, step 10: average elbo=-155.588303
    Epoch 5, step 20: average elbo=-158.082321
    Epoch 5, step 30: average elbo=-155.623718
    Epoch 5, step 40: average elbo=-156.327423
    Epoch 5, step 50: average elbo=-171.485825
    Epoch 5, step 60: average elbo=-161.134216
    Epoch 5, step 70: average elbo=-161.478607
    Epoch 5, step 80: average elbo=-154.923416
    Epoch 5, step 90: average elbo=-166.341919
    Epoch 5, step 100: average elbo=-161.477417
    Epoch 5, step 110: average elbo=-157.946991
    Epoch 5, step 120: average elbo=-164.460999
    Epoch 5, step 130: average elbo=-156.922211
    Epoch 5, step 140: average elbo=-155.834213
    Epoch 5, step 150: average elbo=-159.178268
    Epoch 5, step 160: average elbo=-163.428238
    Epoch 5, step 170: average elbo=-157.463638
    Epoch 5, step 180: average elbo=-154.159225
    Epoch 5, step 190: average elbo=-153.533615
    Epoch 5, step 200: average elbo=-154.449234
    Epoch 5, step 210: average elbo=-154.205353
    Epoch 5, step 220: average elbo=-157.540344
    Epoch 5, step 230: average elbo=-151.772644
    Epoch 5, step 240: average elbo=-155.326782
    Epoch 5, step 250: average elbo=-168.326431
    Epoch 5, step 260: average elbo=-158.993423
    Epoch 5, step 270: average elbo=-155.198044
    Epoch 5, step 280: average elbo=-160.985458
    Epoch 5, step 290: average elbo=-153.521072
    Epoch 5, step 300: average elbo=-149.906235
    Epoch 5, step 310: average elbo=-154.694489
    Epoch 5, step 320: average elbo=-147.911850
    Epoch 5, step 330: average elbo=-157.380417
    Epoch 5, step 340: average elbo=-154.107620
    Epoch 5, step 350: average elbo=-151.641586
    Epoch 5, step 360: average elbo=-152.034393
    Epoch 5, step 370: average elbo=-162.538422
    Epoch 5, step 380: average elbo=-150.565140
    Epoch 5, step 390: average elbo=-155.165543
    Epoch 5, step 400: average elbo=-156.571915
    Epoch 5, step 410: average elbo=-161.094864
    Epoch 5, step 420: average elbo=-170.060822
    Epoch 5, step 430: average elbo=-151.670639
    Epoch 5, step 440: average elbo=-155.040726
    Epoch 5, step 450: average elbo=-157.455261
    Epoch 5, step 460: average elbo=-154.977341
    Epoch 5, step 470: average elbo=-160.902618
    Epoch 5, step 480: average elbo=-153.935730
    Epoch 5, step 490: average elbo=-153.944031
    Epoch 5, step 500: average elbo=-149.662338
    Epoch 5, step 510: average elbo=-157.111954
    Epoch 5, step 520: average elbo=-162.256409
    Epoch 5, step 530: average elbo=-150.760132
    Epoch 5, step 540: average elbo=-152.029526
    Epoch 6, step 0: average elbo=-155.732666
    Epoch 6, step 10: average elbo=-156.560745
    Epoch 6, step 20: average elbo=-152.059937
    Epoch 6, step 30: average elbo=-152.967682
    Epoch 6, step 40: average elbo=-158.282028
    Epoch 6, step 50: average elbo=-145.727631
    Epoch 6, step 60: average elbo=-152.622086
    Epoch 6, step 70: average elbo=-149.752701
    Epoch 6, step 80: average elbo=-152.749710
    Epoch 6, step 90: average elbo=-153.806732
    Epoch 6, step 100: average elbo=-147.049759
    Epoch 6, step 110: average elbo=-150.624130
    Epoch 6, step 120: average elbo=-152.360565
    Epoch 6, step 130: average elbo=-149.232956
    Epoch 6, step 140: average elbo=-162.980728
    Epoch 6, step 150: average elbo=-154.330383
    Epoch 6, step 160: average elbo=-155.444077
    Epoch 6, step 170: average elbo=-150.516769
    Epoch 6, step 180: average elbo=-151.843826
    Epoch 6, step 190: average elbo=-148.542816
    Epoch 6, step 200: average elbo=-151.447815
    Epoch 6, step 210: average elbo=-155.992691
    Epoch 6, step 220: average elbo=-151.121979
    Epoch 6, step 230: average elbo=-149.608627
    Epoch 6, step 240: average elbo=-152.166275
    Epoch 6, step 250: average elbo=-153.272659
    Epoch 6, step 260: average elbo=-148.642929
    Epoch 6, step 270: average elbo=-147.143890
    Epoch 6, step 280: average elbo=-155.657715
    Epoch 6, step 290: average elbo=-139.734222
    Epoch 6, step 300: average elbo=-148.644791
    Epoch 6, step 310: average elbo=-142.534683
    Epoch 6, step 320: average elbo=-156.466125
    Epoch 6, step 330: average elbo=-150.992905
    Epoch 6, step 340: average elbo=-157.106628
    Epoch 6, step 350: average elbo=-153.045273
    Epoch 6, step 360: average elbo=-150.769318
    Epoch 6, step 370: average elbo=-147.912750
    Epoch 6, step 380: average elbo=-151.458710
    Epoch 6, step 390: average elbo=-152.496918
    Epoch 6, step 400: average elbo=-152.984879
    Epoch 6, step 410: average elbo=-151.655121
    Epoch 6, step 420: average elbo=-151.046295
    Epoch 6, step 430: average elbo=-144.511719
    Epoch 6, step 440: average elbo=-152.132370
    Epoch 6, step 450: average elbo=-149.324203
    Epoch 6, step 460: average elbo=-145.197098
    Epoch 6, step 470: average elbo=-151.296692
    Epoch 6, step 480: average elbo=-158.355743
    Epoch 6, step 490: average elbo=-149.179382
    Epoch 6, step 500: average elbo=-158.373520
    Epoch 6, step 510: average elbo=-153.022156
    Epoch 6, step 520: average elbo=-152.593201
    Epoch 6, step 530: average elbo=-153.738419
    Epoch 6, step 540: average elbo=-146.962357
    Epoch 7, step 0: average elbo=-151.509415
    Epoch 7, step 10: average elbo=-149.951126
    Epoch 7, step 20: average elbo=-160.073105
    Epoch 7, step 30: average elbo=-152.959381
    Epoch 7, step 40: average elbo=-149.960999
    Epoch 7, step 50: average elbo=-152.841812
    Epoch 7, step 60: average elbo=-143.995834
    Epoch 7, step 70: average elbo=-148.298065
    Epoch 7, step 80: average elbo=-161.453522
    Epoch 7, step 90: average elbo=-144.098984
    Epoch 7, step 100: average elbo=-148.892197
    Epoch 7, step 110: average elbo=-152.635895
    Epoch 7, step 120: average elbo=-150.777634
    Epoch 7, step 130: average elbo=-146.769760
    Epoch 7, step 140: average elbo=-150.729996
    Epoch 7, step 150: average elbo=-154.418198
    Epoch 7, step 160: average elbo=-140.354202
    Epoch 7, step 170: average elbo=-145.693726
    Epoch 7, step 180: average elbo=-149.687698
    Epoch 7, step 190: average elbo=-143.482346
    Epoch 7, step 200: average elbo=-149.290619
    Epoch 7, step 210: average elbo=-150.292618
    Epoch 7, step 220: average elbo=-145.138580
    Epoch 7, step 230: average elbo=-137.252075
    Epoch 7, step 240: average elbo=-148.015137
    Epoch 7, step 250: average elbo=-147.578522
    Epoch 7, step 260: average elbo=-142.746353
    Epoch 7, step 270: average elbo=-146.899185
    Epoch 7, step 280: average elbo=-149.724579
    Epoch 7, step 290: average elbo=-147.446030
    Epoch 7, step 300: average elbo=-148.224670
    Epoch 7, step 310: average elbo=-146.019089
    Epoch 7, step 320: average elbo=-151.039597
    Epoch 7, step 330: average elbo=-143.525162
    Epoch 7, step 340: average elbo=-149.361511
    Epoch 7, step 350: average elbo=-149.946106
    Epoch 7, step 360: average elbo=-145.240402
    Epoch 7, step 370: average elbo=-138.103790
    Epoch 7, step 380: average elbo=-141.549301
    Epoch 7, step 390: average elbo=-142.104797
    Epoch 7, step 400: average elbo=-149.644714
    Epoch 7, step 410: average elbo=-142.766296
    Epoch 7, step 420: average elbo=-151.002380
    Epoch 7, step 430: average elbo=-145.616943
    Epoch 7, step 440: average elbo=-135.332092
    Epoch 7, step 450: average elbo=-144.044830
    Epoch 7, step 460: average elbo=-141.058121
    Epoch 7, step 470: average elbo=-148.847992
    Epoch 7, step 480: average elbo=-137.432938
    Epoch 7, step 490: average elbo=-144.435593
    Epoch 7, step 500: average elbo=-149.157562
    Epoch 7, step 510: average elbo=-143.406174
    Epoch 7, step 520: average elbo=-153.542984
    Epoch 7, step 530: average elbo=-149.728409
    Epoch 7, step 540: average elbo=-148.322662
    Epoch 8, step 0: average elbo=-152.375198
    Epoch 8, step 10: average elbo=-140.552338
    Epoch 8, step 20: average elbo=-144.333649
    Epoch 8, step 30: average elbo=-143.684143
    Epoch 8, step 40: average elbo=-143.386276
    Epoch 8, step 50: average elbo=-140.615875
    Epoch 8, step 60: average elbo=-147.880035
    Epoch 8, step 70: average elbo=-151.309311
    Epoch 8, step 80: average elbo=-147.549103
    Epoch 8, step 90: average elbo=-147.350769
    Epoch 8, step 100: average elbo=-146.331955
    Epoch 8, step 110: average elbo=-151.270279
    Epoch 8, step 120: average elbo=-146.124786
    Epoch 8, step 130: average elbo=-144.839905
    Epoch 8, step 140: average elbo=-149.929428
    Epoch 8, step 150: average elbo=-149.361938
    Epoch 8, step 160: average elbo=-142.470749
    Epoch 8, step 170: average elbo=-142.187073
    Epoch 8, step 180: average elbo=-140.513519
    Epoch 8, step 190: average elbo=-144.875626
    Epoch 8, step 200: average elbo=-149.019836
    Epoch 8, step 210: average elbo=-140.572418
    Epoch 8, step 220: average elbo=-143.931229
    Epoch 8, step 230: average elbo=-150.187775
    Epoch 8, step 240: average elbo=-146.801804
    Epoch 8, step 250: average elbo=-142.372375
    Epoch 8, step 260: average elbo=-152.931595
    Epoch 8, step 270: average elbo=-135.314545
    Epoch 8, step 280: average elbo=-148.277603
    Epoch 8, step 290: average elbo=-142.926910
    Epoch 8, step 300: average elbo=-144.686462
    Epoch 8, step 310: average elbo=-149.680130
    Epoch 8, step 320: average elbo=-146.330063
    Epoch 8, step 330: average elbo=-151.952103
    Epoch 8, step 340: average elbo=-151.465546
    Epoch 8, step 350: average elbo=-149.537109
    Epoch 8, step 360: average elbo=-140.898880
    Epoch 8, step 370: average elbo=-147.172318
    Epoch 8, step 380: average elbo=-138.434204
    Epoch 8, step 390: average elbo=-146.433090
    Epoch 8, step 400: average elbo=-150.392761
    Epoch 8, step 410: average elbo=-140.727234
    Epoch 8, step 420: average elbo=-141.578613
    Epoch 8, step 430: average elbo=-148.095291
    Epoch 8, step 440: average elbo=-145.144806
    Epoch 8, step 450: average elbo=-143.880737
    Epoch 8, step 460: average elbo=-137.393417
    Epoch 8, step 470: average elbo=-140.097153
    Epoch 8, step 480: average elbo=-146.022690
    Epoch 8, step 490: average elbo=-144.845795
    Epoch 8, step 500: average elbo=-143.871811
    Epoch 8, step 510: average elbo=-149.038940
    Epoch 8, step 520: average elbo=-139.470383
    Epoch 8, step 530: average elbo=-148.597366
    Epoch 8, step 540: average elbo=-144.198151
    Epoch 9, step 0: average elbo=-140.357361
    Epoch 9, step 10: average elbo=-141.146988
    Epoch 9, step 20: average elbo=-141.967148
    Epoch 9, step 30: average elbo=-146.998611
    Epoch 9, step 40: average elbo=-144.164032
    Epoch 9, step 50: average elbo=-148.779556
    Epoch 9, step 60: average elbo=-151.063293
    Epoch 9, step 70: average elbo=-137.717789
    Epoch 9, step 80: average elbo=-152.716034
    Epoch 9, step 90: average elbo=-144.909103
    Epoch 9, step 100: average elbo=-143.754822
    Epoch 9, step 110: average elbo=-140.322479
    Epoch 9, step 120: average elbo=-134.213715
    Epoch 9, step 130: average elbo=-147.544861
    Epoch 9, step 140: average elbo=-143.714844
    Epoch 9, step 150: average elbo=-138.135941
    Epoch 9, step 160: average elbo=-153.722473
    Epoch 9, step 170: average elbo=-150.919556
    Epoch 9, step 180: average elbo=-140.218307
    Epoch 9, step 190: average elbo=-142.421478
    Epoch 9, step 200: average elbo=-134.036087
    Epoch 9, step 210: average elbo=-149.913147
    Epoch 9, step 220: average elbo=-145.744049
    Epoch 9, step 230: average elbo=-136.944687
    Epoch 9, step 240: average elbo=-137.228867
    Epoch 9, step 250: average elbo=-135.546371
    Epoch 9, step 260: average elbo=-141.177231
    Epoch 9, step 270: average elbo=-142.965271
    Epoch 9, step 280: average elbo=-138.939972
    Epoch 9, step 290: average elbo=-136.309570
    Epoch 9, step 300: average elbo=-138.670303
    Epoch 9, step 310: average elbo=-148.282288
    Epoch 9, step 320: average elbo=-131.374817
    Epoch 9, step 330: average elbo=-144.719391
    Epoch 9, step 340: average elbo=-141.071304
    Epoch 9, step 350: average elbo=-144.354645
    Epoch 9, step 360: average elbo=-144.292847
    Epoch 9, step 370: average elbo=-139.482513
    Epoch 9, step 380: average elbo=-143.472778
    Epoch 9, step 390: average elbo=-146.002213
    Epoch 9, step 400: average elbo=-144.582214
    Epoch 9, step 410: average elbo=-148.937836
    Epoch 9, step 420: average elbo=-136.593323
    Epoch 9, step 430: average elbo=-146.098541
    Epoch 9, step 440: average elbo=-137.195114
    Epoch 9, step 450: average elbo=-150.506149
    Epoch 9, step 460: average elbo=-136.688766
    Epoch 9, step 470: average elbo=-140.248764
    Epoch 9, step 480: average elbo=-130.937668
    Epoch 9, step 490: average elbo=-141.994659
    Epoch 9, step 500: average elbo=-145.443604
    Epoch 9, step 510: average elbo=-140.903854
    Epoch 9, step 520: average elbo=-138.184250
    Epoch 9, step 530: average elbo=-141.018204
    Epoch 9, step 540: average elbo=-134.956665

```python
plt.plot(elbo_log)
```

```python
x_sample = mnist.test.next_batch(minibatch_size)[0]
x_reconstruct = p_mu.eval(session=sess, feed_dict={x:x_sample})[0]
plt.subplot(1, 2, 1)
plt.imshow(x_sample[0].reshape(28, 28), vmin=0, vmax=1, cmap="gray")
plt.title("Test input")
plt.colorbar()
plt.subplot(1, 2, 2)
plt.imshow(x_reconstruct.reshape(28, 28), vmin=0, vmax=1, cmap="gray")
plt.title("Reconstruction")
plt.colorbar()
plt.tight_layout()
```

max_stars_repo_name: LiwenxuanNJU/TVpgGLM

max_issues_repo_name: LiwenxuanNJU/TVpgGLM

max_forks_repo_name: LiwenxuanNJU/TVpgGLM


text:
```python
import numpy as np
import matplotlib.pyplot as plt
from sympy import Symbol, integrate
%matplotlib inline
```

### Smooth local paths
We will use cubic spirals to generate smooth local paths. Without loss of generality, as $\theta$ smoothly changes from 0 to 1, we impose a condition on the curvature as follows

$\kappa = f'(x) = K(x(1-x))^n $

This ensures curvature vanishes at the beginning and end of the path. Integrating, the yaw changes as
$\theta = \int_0^x f'(x')dx'$

With $n = 1$ we get a cubic spiral, $n=2$ we get a quintic spiral and so on. Let us use the sympy package to find the family of spirals

1. Declare $x$ a Symbol

2. You want to find Integral of $f'(x)$

3. You can choose $K$ so that all coefficients are integers

Verify if $\theta(0) = 0$ and $\theta(1) = 1$

```python
K =  6 #choose for cubic/quintic
n =  1 #choose for cubic/ quintic
x =  Symbol('x') #declare as Symbol
print(integrate(K*(x*(1-x)), x)) # complete the expression
```

    -2*x**3 + 3*x**2

```python
K =  30 #choose for cubic/quintic
n =  2 #choose for cubic/ quintic
x =  Symbol('x') #declare as Symbol
print(integrate(K*(x*(1-x))**n, x)) # complete the expression
```

    6*x**5 - 15*x**4 + 10*x**3

```python
x = np.linspace(0, 1, num=100)
thetas = -2*x**3 + 3*x**2
plt.figure()
plt.plot(x, thetas,'.')
thetas = 6*x**5 - 15*x**4 + 10*x**3
plt.plot(x, thetas,'.')
```

```python
#write function to compute a cubic spiral
#input can be any theta_i and theta_f (not just 0 and 1)
def cubic_spiral(theta_i, theta_f, n=10):
    x = np.linspace(0, 1, num=n)
    thetas = -2*x**3 + 3*x**2
    return (theta_f - theta_i)*thetas + theta_i

def quintic_spiral(theta_i, theta_f, n=10):
    x = np.linspace(0, 1, num=n)
    thetas = 6*x**5 - 15*x**4 + 10*x**3
    return (theta_f - theta_i)*thetas + theta_i
```

### Plotting
Plot cubic, quintic spirals along with how $\theta$ will change when moving in a circular arc. Remember circular arc is when  $\omega $ is constant

```python
num_pts=100
plt.figure()
plt.plot(np.pi/2*(1-np.linspace(0,1,num_pts)), label='Circular')
plt.plot(cubic_spiral(np.pi/2,0,num_pts), label='Cubic')
plt.plot(quintic_spiral(np.pi/2,0,num_pts),label='Quintic')
plt.grid()
plt.legend()
```

## Trajectory

Using the spirals, convert them to trajectories $\{(x_i,y_i,\theta_i)\}$. Remember the unicycle model

$dx = v\cos \theta dt$

$dy = v\sin \theta dt$

$\theta$ is given by the spiral functions you just wrote. Use cumsum() in numpy to calculate {}

What happens when you change $v$?

```python
num_pts = 100
v = 1
dt = 0.1

dx = v*np.cos(cubic_spiral(np.pi/2,0,num_pts))*dt
dy = v*np.sin(cubic_spiral(np.pi/2,0,num_pts))*dt
xc = np.cumsum(dx)
yc = np.cumsum(dy)

dx = v*np.cos(quintic_spiral(np.pi/2,0,num_pts))*dt
dy = v*np.sin(quintic_spiral(np.pi/2,0,num_pts))*dt
xq = np.cumsum(dx)
yq = np.cumsum(dy)

dx = v*np.cos(np.pi/2*(1-np.linspace(0,1,num_pts)))*dt
dy = v*np.sin(np.pi/2*(1-np.linspace(0,1,num_pts)))*dt
x = np.cumsum(dx)
y = np.cumsum(dy)

# plot trajectories for circular/ cubic/ quintic
plt.figure()
plt.plot(xc, yc, label='Cubic')
plt.plot(xq, yq, label='Quintic')
plt.plot(x, y, label='Circular')
plt.legend()
plt.grid()
```

```python
num_pts = 100
v = 1
dt = 0.1

dx = v*np.cos(cubic_spiral(np.pi/2,np.pi,num_pts))*dt
dy = v*np.sin(cubic_spiral(np.pi/2,np.pi,num_pts))*dt
xc = np.cumsum(dx)
yc = np.cumsum(dy)

dx = v*np.cos(quintic_spiral(np.pi/2,np.pi,num_pts))*dt
dy = v*np.sin(quintic_spiral(np.pi/2,np.pi,num_pts))*dt
xq = np.cumsum(dx)
yq = np.cumsum(dy)

dx = v*np.cos(np.pi/2*(1+np.linspace(0,1,num_pts)))*dt
dy = v*np.sin(np.pi/2*(1+np.linspace(0,1,num_pts)))*dt
x = np.cumsum(dx)
y = np.cumsum(dy)

# plot trajectories for circular/ cubic/ quintic
plt.figure()
plt.plot(xc, yc, label='Cubic')
plt.plot(xq, yq, label='Quintic')
plt.plot(x, y, label='Circular')
plt.legend()
plt.grid()
```

## Symmetric poses

We have been doing only examples with $|\theta_i - \theta_f| = \pi/2$.

What about other orientation changes? Given below is an array of terminal angles (they are in degrees!). Start from 0 deg and plot the family of trajectories

```python
dt = 0.1
thetas = [15, 30, 45, 60, 90, 120, 150, 180] #convert to radians
plt.figure()
for tf in thetas:
    t = cubic_spiral(0, np.deg2rad(tf),50)
    x = np.cumsum(np.cos(t)*dt)
    y = np.cumsum(np.sin(t)*dt)
    plt.plot(x, y)

# On the same plot, move from 180 to 180 - theta
thetas = np.deg2rad([15, 30, 45, 60, 90, 120, 150, 180])
#plt.figure()
for tf in thetas:
    t = cubic_spiral(np.pi, np.pi-tf,50)
    x = np.cumsum(np.cos(t)*dt)
    y = np.cumsum(np.sin(t)*dt)
    plt.plot(x, y)

plt.grid()
```

Modify your code to print the following for the positive terminal angles $\{\theta_f\}$
1. Final x, y position in corresponding trajectory: $x_f, y_f$
2. $\frac{y_f}{x_f}$ and $\tan \frac{\theta_f}{2}$

What do you notice?
What happens when $v$ is doubled?

```python
dt = 0.1
thetas = [15, 30, 45, 60, 90, 120, 150, 180] #convert to radians
plt.figure()
for tf in thetas:
    t = cubic_spiral(0, np.deg2rad(tf),50)
    x = np.cumsum(np.cos(t)*dt)
    y = np.cumsum(np.sin(t)*dt)
    print(f"tf:{tf:0.1f}--- xf:{x[-1]:0.3f}, yf:{y[-1]:0.3f}, yf/xf:{y[-1]/x[-1]:0.3f}, tan(theta/2):{np.tan(np.deg2rad(tf)/2):0.3f}")
```

    tf:15.0--- xf:4.936, yf:0.650, yf/xf:0.132, tan(theta/2):0.132
    tf:30.0--- xf:4.748, yf:1.272, yf/xf:0.268, tan(theta/2):0.268
    tf:45.0--- xf:4.444, yf:1.841, yf/xf:0.414, tan(theta/2):0.414
    tf:60.0--- xf:4.041, yf:2.333, yf/xf:0.577, tan(theta/2):0.577
    tf:90.0--- xf:3.015, yf:3.015, yf/xf:1.000, tan(theta/2):1.000
    tf:120.0--- xf:1.865, yf:3.231, yf/xf:1.732, tan(theta/2):1.732
    tf:150.0--- xf:0.800, yf:2.987, yf/xf:3.732, tan(theta/2):3.732
    tf:180.0--- xf:0.000, yf:2.382, yf/xf:42906192676138176.000, tan(theta/2):16331239353195370.000

    <Figure size 432x288 with 0 Axes>

These are called *symmetric poses*. With this spiral-fitting approach, only symmetric poses can be reached.

In order to move between any 2 arbitrary poses, you will have to find an intermediate pose that is pair-wise symmetric to the start and the end pose.

What should be the intermediate pose? There are infinite possibilities. We would have to formulate it as an optimization problem. As they say, that has to be left for another time!

```python

```