```python
%pylab inline
%config InlineBackend.figure_format = 'retina'
from ipywidgets import interact
import scipy.stats as stats
import pandas as pd
# import pymc3 as pm
import arviz as az
import seaborn as sns

import tensorflow as tf
import tensorflow_probability as tfp
tfd = tfp.distributions
tfb = tfp.bijectors

import warnings # the warning spam is pointless and annoying
warnings.simplefilter(action="ignore", category=FutureWarning)
```

    Populating the interactive namespace from numpy and matplotlib


# Example
Suppose we are modeling animal movement, specifically polar bears. We are planning to utilize data that comes from attaching a tracking collar to our polar bear. But first we are exploring possible models. 

Let's start by modeling the movement of a bear as a sequence of random steps.

## Suppose for model 1 (m1) we pick normal random steps $x,y \sim \text{Normal}(0, 1)$


```python
Nsteps = 10000

## normally distributed steps
X_m1 = normal(0, 1, [Nsteps, 2])
```


```python
## random path of a hypothetical polar bear
Z_m1 = X_m1.cumsum(axis=0) ## this is the path that results from adding up steps
plot(Z_m1[:, 0], Z_m1[:, 1])
```


```python
H, xb, yb = histogram2d(X_m1[:, 0], X_m1[:, 1])
pcolor(xb, yb, H);
# sns.kdeplot(X[:, 0], X[:, 1])
```


```python
sns.distplot(X_m1[:, 1])
```

## Suppose for model 2 (m2) we pick a random step length $r\sim\text{Normal}(0, 1)$ and a random direction $\theta \sim \text{Uniform}(0, 2\pi)$


```python
## Random step length and a random direction
R_m2 = normal(0, 1, Nsteps) ## normal step length
theta_m2 = 2*pi*rand(Nsteps) ## uniform direction
```


```python
X_m2 = stack([R_m2*cos(theta_m2), R_m2*sin(theta_m2)], 1)
Z_m2 = X_m2.cumsum(axis=0) ## this is the path that results from adding up steps
plot(Z_m2[:, 0], Z_m2[:, 1]);
```


```python
H, xb, yb = histogram2d(X_m2[:, 0], X_m2[:, 1])
pcolor(xb, yb, H);
# sns.kdeplot(X2[:, 0], X2[:, 1])
```


```python
sns.distplot(R_m2);
```

## Are the two models equivalent?


```python
sns.distplot(X_m1[:, 0], label='model 1')
sns.distplot(X_m2[:, 0], label='model 2')# a slice of our 2D density does not look Normal
legend()
```

## what happens if we look at step lengths from the first model where $x,y\sim\text{Normal}(0, 1)$?


```python
R_m1 = sqrt(X_m1[:, 0]**2 + X_m1[:, 1]**2) ## length of the step
sns.distplot(R_m1, label='model 1')
sns.distplot(absolute(R_m2), label='model 2')
xlabel(r'$r$', fontsize=24)
legend();
```

## Conclusion: these are different models and not equivalent
The fundamental reason for the confusion in our intuition is that we failed to account for how the *area* (volume in general) of an infinitesimal region changes when we change variables.

Consider the region $S_{dr}$ to be the small circle of radius $dr$ around $x,y=0,0$. The area of $S_{dr}$ is $A=\pi dr^2$. In the first model, the center of $S_{dr}$ is the most likely spot, yet it has very small area, so samples will be very unlikely to be found in $S_{dr}$. This is kind of like the bulls eye of a dart board. The probability of hitting the bulls eye is $P(x \in S_{dr}) \approx \text{Area}(S_{dr}) p(0)$. **We have the competition of two factors: the small area $\text{Area}(S_{dr})$ and the large density $p(0)$.** The combination yields a vanishing density in the variable $r$ as $r\to 0^+$.

In the second model we pick normal random step lengths (in model 2) which will result in the density for $r$ being Gaussian and therefore maximal at $r=0$.

### We need a conversion factor that tells us how $\text{Area}(x, y; dx, dy)$ is related to $\text{Area}(r, \theta; dr, d\theta)$

We have 
$$ x = r\cos(\theta), \quad y = r\sin(\theta) $$
The Jacobian 
$$ J = \begin{bmatrix}
\cos(\theta) & -r\sin(\theta) \\
\sin(\theta) & r\cos(\theta)
\end{bmatrix}$$
tells us
$$
\begin{bmatrix}
dx \\ dy
\end{bmatrix}
= 
\begin{bmatrix}
\cos(\theta) & -r\sin(\theta) \\
\sin(\theta) & r\cos(\theta)
\end{bmatrix}
\begin{bmatrix}
dr \\ d\theta
\end{bmatrix}
$$
The area of an infinitesimal region satisfies
$$ dxdy = \det(J)dr d\theta.$$
The determinant of the Jacobian is 
$$\det(J) = r.$$

--------------------------------------------- 

# Change of variables for a random variable
Suppose that we have a RV $X\in \mathbb{R}$ with density function $p(x)$. Suppose that $f(\cdot)$ is a bijection so that 
$$ y = f(x) \Leftarrow\Rightarrow f^{-1}(y) = x.$$

Consider the new RV $$Y = f(X).$$ What is the density function for $Y$? Is it $$ g(y) = p(f^{-1}(y))?$$

Consider the integral
$$\int_{-\infty}^{\infty}p(x)dx = 1$$
We must also have
$$\int_{-\infty}^{\infty}g(y)dy = 1$$
So then
$$\int_{-\infty}^{\infty}p(f^{-1}(y))dy = ?$$

## Recall Calculus
How do we do a change of variables in an integral?
$$ y = f(x) \Rightarrow dy = f'(x)dx \Rightarrow dx = \frac{dy}{f'(x)}$$

$$ 1 = \int_{-\infty}^{\infty}p(x)dx 
= \int_{-\infty}^{\infty}p(f^{-1}(y)) dx 
= \int_{-\infty}^{\infty}p(f^{-1}(y)) \frac{dy}{f'(x)}
= \int_{-\infty}^{\infty}\frac{p(f^{-1}(y))}{f'(f^{-1}(y))} dy$$

Alternatively, 
$$ x = f^{-1}(y) \Rightarrow dx = (f^{-1})'(y)dy,$$
$$ 1 = \int_{-\infty}^{\infty} (f^{-1})'(y) p(f^{-1}(y)) dy$$

## Conclusion
The probability density for $Y = f(X)$ is
$$g(y) = (f^{-1})'(y) p(f^{-1}(y)) $$

## In higher dimensions
If $X\in \mathbb{R}^n$, $Y \in\mathbb{R}^n$, and $Y = f(X)$ then the density function for $Y$ is
$$ g(y) = \det(D_y [f^{-1}](y)) p(f^{-1}(y)), $$
where $D_y$ denotes the operator that yields the Jacobian matrix of partial derivatives.

## The volume conversion factor we were looking for is given by
$$\det(D_y [f^{-1}](y))$$

## Measure theory generalization
The volume conversion factor is called the [Radon-Nikodym derivative](https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem).

----------------------

# Bijectors
A [bijector](https://en.wikipedia.org/wiki/Bijection) is a function that can be inverted. That is, if $y = f(x)$ then there exists an inverse function such that $x = f^{-1}(y)$.

## Experiment
The following is an experiment with the kernel density estimator (KDE). The standard KDE is meant to approximate the density function for an unbounded RV $x\in\mathbb{R}$ (or higher dimensions). The idea in this experiment is to impose constraints on the RV so that, for example, the density is restricted to $x\geq 0$ or $x \in [a, b]$. This is a reasonable idea, but **fair warning, I do not know how well behaved the resulting approximation will be.**


```python
Nsamples = 5000
u = rand(Nsamples)
lam = 5.
tau = -log(u)/lam
f = tfb.Log()
y = f.forward(tau)
py = stats.kde.gaussian_kde(y)
# plot(py(linspace(-5, 3)))
def p_tau(tau):
    y = array(f.forward(tau))
    logJy = array(f.forward_log_det_jacobian(tau, 0))
    return exp(logJy + py.logpdf(y))
tau_plot = linspace(1e-2, 2, 100)
density_plot = p_tau(tau_plot)
plot(tau_plot, density_plot, 'r')
plot(tau_plot, lam*exp(-lam*tau_plot), 'k')
trapz(density_plot, tau_plot)
sns.distplot(tau)
```


```python
Nsamples = 5000
u = -1 + 2*rand(Nsamples)
# f = tfb.Invert(tfb.Tanh())
f = tfb.Invert(tfb.SoftClip(low=-1-1e-7, high=1+1e-7, hinge_softness=4))
y = f.forward(u)
py = stats.kde.gaussian_kde(y)
# plot(py(linspace(-5, 3)))
def p_u(u):
    y = array(f.forward(u))
    logJy = array(f.forward_log_det_jacobian(u, 0))
    return exp(logJy + py.logpdf(y))
u_plot = linspace(-1, 1, 100)
density_plot = p_u(u_plot)
plot(u_plot, density_plot, 'r')
ylim(0, density_plot.max()*1.2)
trapz(density_plot, u_plot)
sns.distplot(u)
```

# Bijectors with MCMC
We will recreate our globe tossing model. Recall the probability parameter in a Binomial is $0\leq p \leq 1$. How do we sample the posterior $\pi(p)$ when Hamiltonian Monte Carlo is designed to sample distributions in $\mathbb{R}^n$?


```python
w = 6 ## from Chapter 3 and Week 3 lecture
n = 9

# data = 1*(rand(100) < 0.7)
# w = data.sum()
# n = data.size
```


```python
Nsteps = 2000
N_burn_in = Nsteps

Nchains = 4

def log_target_density(p):
    return tfd.Binomial(n, probs=p).unnormalized_log_prob(w)

@tf.function()
def sampler():
    def trace_everything(states, previous_kernel_results):
        return previous_kernel_results
    Xinit = rand(Nchains)
    kernel = tfp.mcmc.NoUTurnSampler(
        target_log_prob_fn=log_target_density,
        max_tree_depth=8,
        max_energy_diff=1000,
        step_size=tf.cast(0.1, tf.float64))

    # kernel = tfp.mcmc.HamiltonianMonteCarlo(
    #     target_log_prob_fn=log_target_density,
    #     num_leapfrog_steps=10,
    #     step_size=tf.cast(0.1, tf.float64))
    
    bijectors = [tfb.NormalCDF()]
    
    transformed_kernel = tfp.mcmc.TransformedTransitionKernel(
        inner_kernel=kernel,
        bijector=bijectors)

    adaptive_kernel = tfp.mcmc.DualAveragingStepSizeAdaptation(
        inner_kernel=transformed_kernel,
        num_adaptation_steps=int(0.8 * N_burn_in),
        target_accept_prob=tf.cast(0.7, tf.float64))
    
    Xchain, trace = tfp.mcmc.sample_chain(
        num_results=Nsteps,
        num_burnin_steps=N_burn_in,
        current_state=Xinit,
        trace_fn=trace_everything,
        kernel=adaptive_kernel,
        parallel_iterations=4)
    return Xchain, trace
Xchain, trace = sampler()
Xchain = array(Xchain)
```

    WARNING:tensorflow:@custom_gradient grad_fn has 'variables' in signature, but no ResourceVariables were used on the forward pass.



```python
tfp.mcmc.effective_sample_size(Xchain)
```




    <tf.Tensor: shape=(4,), dtype=float64, numpy=array([1339.05453023, 1264.62717084, 1414.23880154, 1536.30757449])>




```python
sns.distplot(Xchain.flatten())
x_to_plot = linspace(0, 1, 100)
plot(x_to_plot, stats.beta.pdf(x_to_plot, w + 1, n - w + 1), '--k', label="True posterior")
legend()
xlabel(r'$p$', fontsize=24)
ylabel('probability density', fontsize=20);
```

# A full working Tensorflow + Bijector version of the random time HMC from Week 11


```python
w = 6 ## from Chapter 3 and Week 3 lecture
n = 9
```


```python
Nsteps = 2000
N_burn_in = Nsteps
Nchains = 8

lam = 0.7
dt = 0.1

bijector = tfb.NormalCDF() ## our bijector

x0 = normal(0, 1, [Nchains, 1])

@tf.function()
def sample(x0):
    def L(x):
        log_likelihood = tfd.Binomial(n, probs=bijector.forward(x)).unnormalized_log_prob(w)
        log_prior = 0. ## we could change this to an informative prior
        log_J = bijector.forward_log_det_jacobian(x, 0)
        return log_likelihood + log_prior + log_J
    def gradL(x): 
        with tf.GradientTape() as tape:
            tape.watch(x)
            L_at_x = L(x)
        return tape.gradient(L_at_x, x)
    def H(x, p):
        return -L(x) + 0.5*p**2

    x = tf.identity(x0) ## make a copy of the values in the Variable x0
    
    ## Burn in steps
    for j in tf.range(N_burn_in): ## notice that this is not `arange` or `range`
        p0 = tf.random.normal(mean=0, stddev=1, shape=[Nchains, 1], dtype=tf.float64)
        p = tf.identity(p0) ## make a copy of p0
        tau = -tf.math.log(tf.random.uniform(shape=[1]))/lam ## generate an exponential RV for time
        K = tf.cast(tau/dt, tf.int64) ## cast (convert) a float to an integer
        for k in tf.range(K[0]):
            _p = p + 0.5*dt*gradL(x)
            x = x + dt*_p
            p = _p + 0.5*dt*gradL(x)

    Xchain = tf.TensorArray( ## We have to use this to store values in a tensorflow loop
        tf.float64,
        size=Nsteps,
        element_shape=[Nchains, 1])
    Xchain = Xchain.write(0, x) ## The write method saves a time step (the first in this case)

    for j in tf.range(Nsteps-1):
        p0 = tf.random.normal(mean=0, stddev=1, shape=[Nchains, 1], dtype=tf.float64)
        x0 = tf.identity(x) ## make a copy of the values in the Variable x
        p = tf.identity(p0)
        tau = -tf.math.log(tf.random.uniform(shape=[1]))/lam ## generate an exponential RV for time
        K = tf.cast(tau/dt, tf.int64)
        for k in tf.range(K[0]):
            _p = p + 0.5*dt*gradL(x)
            x = x + dt*_p
            p = _p + 0.5*dt*gradL(x)
        log_p_MH = H(x0, p0) - H(x, p) ## if our numerical integrator is perfect, this should always be one
        u = tf.random.uniform(shape=[Nchains, 1], dtype=tf.float64)
        chi = tf.cast(tf.math.log(u) < log_p_MH, tf.float64) ## convert a boolean (True/False) to a float
        x = x0*(1 - chi) + x*chi
        Xchain = Xchain.write(j+1, x)
    return bijector.forward(Xchain.stack()) ## note the final application of the bijector
Xchain = sample(x0)
Xchain = array(Xchain)
```


```python
tfp.mcmc.effective_sample_size(Xchain)
```




    <tf.Tensor: shape=(8, 1), dtype=float64, numpy=
    array([[1597.48126007],
           [1221.27567889],
           [1412.30453425],
           [1440.66529007],
           [1630.61013096],
           [1382.13776459],
           [1635.34584095],
           [1465.03662131]])>




```python
sns.distplot(Xchain.flatten())
x_to_plot = linspace(0, 1, 100)
plot(x_to_plot, stats.beta.pdf(x_to_plot, w + 1, n - w + 1), '--k', label="True posterior")
legend()
xlabel(r'$p$', fontsize=24)
ylabel('probability density', fontsize=20);
```


```python

```

# State space models
Let us continue with our example about tracking polar bears. Suppose we now have data from a tracking collar of one of our bears. The GPS measurements that come from the collar are samples at equal time intervals of 1 hour. Our collar measurements will be given by
$$ z_{k} \sim \text{Normal}(x_{k}, \sigma_{\rm obs}),$$
where $x_{k}$ is the actual position of the bear at times $t_k = k$. We will model the actual position of the bear as a simple random walk
$$ x_k \sim \text{Normal}(x_{k-1}, \sigma)$$

The full generative model is then
\begin{align}
  z_{k} &\sim \text{Normal}(x_{k}, \sigma_{\rm obs}) \\
   x_k &\sim \text{Normal}(x_{k-1}, \sigma) \\
   x_0 & \sim  \text{Normal}(0, 25)\\
   \sigma_{\rm obs} &\sim ??? \\
   \sigma &\sim ???
\end{align}

## The variables $x_i$ are unobserved
Things we might be interested in knowing:
  1. The full joint posterior over the unobserved positions
    $$\pi(x_0, \ldots, x_T | z_1, \ldots, z_T; \theta) = \frac{\prod_{k=1}^{T}g(z_k| x_k)f(x_k | x_{k-1})p_0(x_0)}{\int_{R}\cdots \int_{R}\prod_{k=1}^{T}g(z_k| x_k)f(x_k | x_{k-1})p_0(x_0) \,dx_0\cdots dx_T}$$
  2. The marginals
    $$\pi(x_k | z_1, \ldots, z_T)$$


# Expectation Maximization Algorithm
  1. E-step: given fixed parameters $\theta_j$ compute 
  $$ X_{1:T} | z_{1:T}, \theta_j$$
  \begin{align}
  Q(\theta | \theta_j) &= E[\log(\prod_{k=1}^{T}g(z_k| X_k, \theta)f(X_k | X_{k-1}, \theta)p_0(X_0|\theta)) | \theta_j]\\
  &= E[\sum_{k=1}^{T}\log(g(z_k| X_k, \theta)) + \sum_{k=1}^{T}\log(f(X_k | X_{k-1}, \theta)) + \log(p_0(X_0|\theta)) | \theta_j]
  \end{align}
  2. M-step: compute $$\theta_{j+1} = \text{arg max}_{\theta} Q(\theta | \theta_j)  $$
  
Repeat these steps until convergence is reached.

## Simulate collar tracking GPS data


```python
Nt = 10
sigma_obs = 0.2
sigma = 2.

X0true = -25. + 50.*rand(1, 2)
Xtrue = X0true + normal(0, 2., [Nt, 2]).cumsum(axis=0)
Z = Xtrue + normal(0, sigma_obs, [Nt, 2])
plot(Z[:, 0], Z[:, 1]);
```

## State space model in Tensorflow
First, we will sample the posterior
$$\pi(x_0, x_1, \ldots, x_T | z_1, \ldots, z_T; \theta) \propto \prod_{k=1}^{T}g(z_k| x_k)f(x_k | x_{k-1})p_0(x_0)$$
The generative model is
$$p(x_0, x_1, \ldots, x_T, z_1, \ldots, z_T | \theta) = \prod_{k=1}^{T}g(z_k| x_k)f(x_k | x_{k-1})p_0(x_0)$$



```python
Nsteps = 50
N_burn_in = Nsteps
Nchains = 2


Sigma_obs = sigma_obs*tf.ones([1, 1, 2], dtype=tf.float64)
Sigma = sigma*tf.ones([1, 1, 2], dtype=tf.float64)
Sigma_x0 = 100*tf.ones([1, 2], dtype=tf.float64)
def log_target_density(Xfull):
    ## X : [parallel chains, time, space]
    X0 = Xfull[:, 0, :]
    X = Xfull[:, 1:, :] ## same shape as `Z`
    ## observation likelihood: density for Z given X
    G = tfd.MultivariateNormalDiag(
        loc=X, 
        scale_diag=Sigma_obs)
    logG = tf.reduce_sum(
        G.unnormalized_log_prob(Z),
        axis=1)
    ## motion model: density for X
    F = tfd.MultivariateNormalDiag(
        loc=Xfull[:, :-1, :], 
        scale_diag=Sigma)
    logF = tf.reduce_sum(
        F.unnormalized_log_prob(X),
        axis=1)
    ## initial position prior: density for X0
    F0 = tfd.MultivariateNormalDiag(
        loc=tf.zeros([1, 2], dtype=tf.float64), 
        scale_diag=Sigma_x0)
    logF0 = F0.unnormalized_log_prob(X0)
    return logG + logF + logF0

@tf.function()
def track_sampler(Xinit):
    kernel = tfp.mcmc.NoUTurnSampler(
        target_log_prob_fn=log_target_density,
        max_tree_depth=8,
        max_energy_diff=1000,
        step_size=tf.cast(0.1, tf.float64))

    # kernel = tfp.mcmc.HamiltonianMonteCarlo(
    #     target_log_prob_fn=log_target_density,
    #     num_leapfrog_steps=10,
    #     step_size=tf.cast(0.1, tf.float64))
    

    adaptive_kernel = tfp.mcmc.DualAveragingStepSizeAdaptation(
        inner_kernel=kernel,
        num_adaptation_steps=int(0.8 * N_burn_in),
        target_accept_prob=tf.cast(0.7, tf.float64))
    
    Xchain = tfp.mcmc.sample_chain(
        num_results=Nsteps,
        num_burnin_steps=N_burn_in,
        current_state=Xinit,
        #trace_fn=trace_everything,
        trace_fn=None,
        kernel=adaptive_kernel,
        parallel_iterations=4)
    return Xchain

Xinit = zeros([Nchains, Nt+1, 2])
Xinit[:, 0, :] = normal(0, 25., [Nchains, 2])
Xinit[:, 1:, :] = Z + normal(0, 5., [Nchains, Nt, 2])
Xchain = track_sampler(Xinit)
Xchain = array(Xchain)
```


```python
mean(tfp.mcmc.effective_sample_size(Xchain))
```




    44.1706914348532




```python
## plot results
```


```python
figure(1, [7, 7])
## Xchain: [mcmc steps, parallel chains, time, space]
plot(Xchain[..., 0].reshape(-1, Nt+1), Xchain[..., 1].reshape(-1, Nt+1), '.')
plot(Z[:, 0], Z[:, 1], 'r')
xlabel('x', fontsize=20)
ylabel('y', fontsize=20);
```

# Suppose we want to sample $X$ at times that are in between the observation times


```python
Nt = 100
sigma_subsample = sigma/sqrt(10)

Zpadded = zeros([Nt, 2])
Zpadded[::10] = Z
### where Z is zero, we will modify our model so that these values of Z have no effect

## this is an indicator variable for when we have measured positions
CHI_observed = zeros([Nt]) 
CHI_observed[::10] = 1.
```


```python
Nsteps = 50
N_burn_in = Nsteps
Nchains = 5


Sigma_obs = sigma_obs*tf.ones([1, 1, 2], dtype=tf.float64)
Sigma = sigma_subsample*tf.ones([1, 1, 2], dtype=tf.float64)
Sigma_x0 = 100*tf.ones([1, 2], dtype=tf.float64)
def log_target_density(Xfull):
    ## X : [parallel chains, time, space]
    X0 = Xfull[:, 0, :]
    X = Xfull[:, 1:, :] ## same shape as `Z`
    ## observation likelihood: density for Z given X
    G = tfd.MultivariateNormalDiag(
        loc=X, 
        scale_diag=Sigma_obs)
    logG = tf.reduce_sum(
        CHI_observed[None, :]*G.unnormalized_log_prob(Zpadded),
        axis=1)
    ## motion model: density for X
    F = tfd.MultivariateNormalDiag(
        loc=Xfull[:, :-1, :], 
        scale_diag=Sigma)
    logF = tf.reduce_sum(
        F.unnormalized_log_prob(X),
        axis=1)
    ## initial position prior: density for X0
    F0 = tfd.MultivariateNormalDiag(
        loc=tf.zeros([1, 2], dtype=tf.float64), 
        scale_diag=Sigma_x0)
    logF0 = F0.unnormalized_log_prob(X0)
    return logG + logF + logF0

@tf.function()
def track_sampler(Xinit):
    ## We are sampling a distribution over 200 dimensions
    
    kernel = tfp.mcmc.NoUTurnSampler(
        target_log_prob_fn=log_target_density,
        max_tree_depth=8,
        max_energy_diff=1000,
        step_size=tf.cast(0.1, tf.float64))

    # kernel = tfp.mcmc.HamiltonianMonteCarlo(
    #     target_log_prob_fn=log_target_density,
    #     num_leapfrog_steps=10,
    #     step_size=tf.cast(0.1, tf.float64))
    

    adaptive_kernel = tfp.mcmc.DualAveragingStepSizeAdaptation(
        inner_kernel=kernel,
        num_adaptation_steps=int(0.8 * N_burn_in),
        target_accept_prob=tf.cast(0.7, tf.float64))
    
    Xchain = tfp.mcmc.sample_chain(
        num_results=Nsteps,
        num_burnin_steps=N_burn_in,
        current_state=Xinit,
        #trace_fn=trace_everything,
        trace_fn=None,
        kernel=adaptive_kernel,
        parallel_iterations=4)
    return Xchain

Xinit = zeros([Nchains, Nt+1, 2])
Xinit[:, 0, :] = normal(0, 25., [Nchains, 2])
Xinit[:, 1:, :] = Zpadded + normal(0, 5., [Nchains, Nt, 2])
Xchain = track_sampler(Xinit)
Xchain = array(Xchain)
```


```python
mean(tfp.mcmc.effective_sample_size(Xchain))
```




    46.864604441106906



## This plot shows the MCMC result as a time series of marginal distributions for the position at time $t$


```python
xmin = Xchain.min()
xmax = Xchain.max()
@interact(t=(0, Nt, 1))
def pfun(t=0):
    fig = figure(1, [7, 7])
    plot(Xchain[:, :, t, 0], Xchain[:, :, t, 1], 'r.', ms=2)
    if t > 0 and CHI_observed[t-1] == 1:
        plot(Zpadded[t-1, 0], Zpadded[t-1, 1], 'ko', mfc='none', ms=25)
    xlim(xmin, xmax)
    ylim(xmin, xmax)
    xlabel(r'$x$', fontsize=24)
    ylabel(r'$y$', fontsize=24)
```


    interactive(children=(IntSlider(value=0, description='t'), Output()), _dom_classes=('widget-interact',))


## This shows the MCMC result as a collection of sample trajectories
Each parallel chain (at a given step) represents an independent path $X_{0:T}$. Here each is shown as a different color.


```python
xmin = Xchain.min()
xmax = Xchain.max()
@interact(t=(0, Nt, 1))
def pfun(t=0):
    fig = figure(1, [7, 7])
    for i in arange(5):
        plot(Xchain[-1, i, t, 0], Xchain[-1, i, t, 1], '.', ms=10)
    if t > 0 and CHI_observed[t-1] == 1:
        plot(Zpadded[t-1, 0], Zpadded[t-1, 1], 'ko', mfc='none', ms=25)
    xlim(xmin, xmax)
    ylim(xmin, xmax)
    xlabel(r'$x$', fontsize=24)
    ylabel(r'$y$', fontsize=24)
```


    interactive(children=(IntSlider(value=0, description='t'), Output()), _dom_classes=('widget-interact',))



```python

```

5faeec2aa977c650949cdbf31b2c025f3e5682c2
ipynb
Jupyter Notebook
Week 12.ipynb
newby-jay/MATH509-Fall2021-JupyterNotebooks
29d43b0b5f46b55149efb9836f2190fc957c2e3f
Apache-2.0
Week 12.ipynb
newby-jay/MATH509-Fall2021-JupyterNotebooks
29d43b0b5f46b55149efb9836f2190fc957c2e3f
Apache-2.0
Week 12.ipynb
newby-jay/MATH509-Fall2021-JupyterNotebooks
29d43b0b5f46b55149efb9836f2190fc957c2e3f
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES