# Markov-Chain Monte Carlo (MCMC)

## Lotka-Volterra Model
    
$u:$ Prey

$v:$ Predator 

\begin{equation}
\begin{split}
\frac{du}{dt} =&  \alpha\cdot u-\beta \cdot v\cdot u \\
\frac{dv}{dt} =& -\gamma \cdot  v + \delta \cdot u\cdot v
\end{split}
\end{equation}


```python
import argparse
import os


import matplotlib
import matplotlib.pyplot as plt

from jax.experimental.ode import odeint
import jax.numpy as jnp
from jax.random import PRNGKey

import numpyro
import numpyro.distributions as dist
from numpyro.examples.datasets import LYNXHARE, load_dataset
from numpyro.infer import MCMC, NUTS, Predictive

```


```python

def dz_dt(z, t, theta):
    """
    Lotka–Volterra equations. Real positive parameters `alpha`, `beta`, `gamma`, `delta`
    describes the interaction of two species.
    """
    u = z[0]
    v = z[1]
    alpha, beta, gamma, delta = theta[..., 0], theta[..., 1], theta[..., 2], theta[..., 3]
    du_dt = (alpha - beta * v) * u
    dv_dt = (-gamma + delta * u) * v

    return jnp.stack([du_dt, dv_dt])


def model(N, y=None):
    """
    :param int N: number of measurement times
    :param numpy.ndarray y: measured populations with shape (N, 2)
    """
    # initial population
    z_init = numpyro.sample("z_init", dist.LogNormal(jnp.log(10), 1).expand([2]))
    # measurement times
    ts = jnp.arange(float(N))
    # parameters alpha, beta, gamma, delta of dz_dt
    theta = numpyro.sample(
        "theta",
        dist.TruncatedNormal(low=0., loc=jnp.array([1.0, 0.05, 1.0, 0.05]),
                             scale=jnp.array([0.5, 0.05, 0.5, 0.05])))
    # integrate dz/dt, the result will have shape N x 2
    z = odeint(dz_dt, z_init, ts, theta, rtol=1e-6, atol=1e-5, mxstep=1000)
    # measurement errors
    sigma = numpyro.sample("sigma", dist.LogNormal(-1, 1).expand([2]))
    # measured populations
    numpyro.sample("y", dist.LogNormal(jnp.log(z), sigma), obs=y)

```


```python
device = 'cpu'

num_samples = 1000
num_chains = 2
num_warmup = 200

numpyro.set_platform(device)
numpyro.set_host_device_count(num_chains)

_, fetch = load_dataset(LYNXHARE, shuffle=False)
year, data = fetch()  # data is in hare -> lynx order

# use dense_mass for better mixing rate
mcmc = MCMC(NUTS(model, dense_mass=True),
            num_warmup, num_samples, num_chains=num_chains, progress_bar= True)
mcmc.run(PRNGKey(1), N=data.shape[0], y=jnp.log(data))
mcmc.print_summary()

```

    
                    mean       std    median      5.0%     95.0%     n_eff     r_hat
      sigma[0]      0.68      0.24      0.74      0.41      0.92      1.01     11.38
      sigma[1]      0.31      0.09      0.34      0.20      0.40      1.01      9.82
      theta[0]      0.64      0.06      0.67      0.52      0.69      2.78      1.21
      theta[1]      3.61      3.41      3.67      0.18      7.02      1.00    233.99
      theta[2]      2.38      1.68      2.55      0.63      4.06      1.00     35.01
      theta[3]      1.71      1.50      1.76      0.19      3.21      1.00     96.89
     z_init[0]      2.41      1.77      2.01      0.66      4.38      1.01     11.20
     z_init[1]      2.29      1.23      1.89      1.07      3.69      1.01      9.09
    
    Number of divergences: 0



```python
# predict populations
ax = plt.figure(figsize=(15.5, 7.2))
y_pred = Predictive(model, mcmc.get_samples())(PRNGKey(2), data.shape[0])["y"]
pop_pred = jnp.exp(y_pred)
mu, pi = jnp.mean(pop_pred, 0), jnp.percentile(pop_pred, (10, 90), 0)
plt.plot(year, data[:, 0], "ko", mfc="none", ms=4, label="true hare", alpha=0.67)
plt.plot(year, data[:, 1], "bx", label="true lynx")
plt.plot(year, mu[:, 0], "k-.", label="pred hare", lw=1, alpha=0.67)
plt.plot(year, mu[:, 1], "b--", label="pred lynx")
plt.fill_between(year, pi[0, :, 0], pi[1, :, 0], color="k", alpha=0.2)
plt.fill_between(year, pi[0, :, 1], pi[1, :, 1], color="b", alpha=0.3)
plt.gca().set(ylim=(0, 160), xlabel="year", ylabel="population (in thousands)")
plt.title("Posterior predictive (80% CI) with predator-prey pattern.")
plt.legend()
plt.tight_layout()
plt.show()

```


```python
params = mcmc.get_samples()['theta']

fig, ax = plt.subplots(4,2, figsize=(20, 40))

ax[0,0].plot(params[:num_samples,0]), ax[0,1].hist(params[:num_samples,0])
ax[1,0].plot(params[:num_samples,1]), ax[1,1].hist(params[:num_samples,1])
ax[2,0].plot(params[:num_samples,2]), ax[2,1].hist(params[:num_samples,2])
ax[3,0].plot(params[:num_samples,3]), ax[3,1].hist(params[:num_samples,3])

for i in range(4):
    ax[i,0].set_title('Trace plot')
    ax[i,1].set_title('Histogram')
    ax[i,0].set_xlabel('Chain Iter')
    ax[i,1].set_ylabel('Value')
    ax[i,1].set_ylabel('Freq')

plt.show()
```


```python

```

dcd052401fda90661fcb572b30b4c75a7c95a66f
ipynb
Jupyter Notebook
MCMC/example_numpyro.ipynb
ChaosDonkey06/MODSIM_IBIO
e2ce70538aaba3b01f22cf037db21a80316a802f
Apache-2.0
MCMC/example_numpyro.ipynb
ChaosDonkey06/MODSIM_IBIO
e2ce70538aaba3b01f22cf037db21a80316a802f
Apache-2.0
MCMC/example_numpyro.ipynb
ChaosDonkey06/MODSIM_IBIO
e2ce70538aaba3b01f22cf037db21a80316a802f
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES