# Lecture 24 - Sequential Monte Carlo in `PyMC3`


```python
import numpy as np
import pymc3 as pm
import theano as T
from theano import shared, function, tensor as tt
from sample_smc import sample_smc 

try:
    import sympy
except:
    _=!pip install sympy
    import sympy

import matplotlib.pyplot as plt 
import seaborn as sns
sns.set()

import logging

%matplotlib inline 
```

## Objectives

+ Compute the model evidence using `PyMC3`.
+ Do model selection with `PyMC3`.

## Sanity check - Does the calculation of the evidence with PySMC work?

Let
$$
p(\theta) = \mathcal{N}(\theta|0, 1),
$$
and
$$
p(y|\theta) = \mathcal{N}(y|\theta,0).
$$
The posterior of $\theta$ given $y$ is:
$$
p(\theta|y) = \frac{p(y|\theta)p(\theta)}{Z},
$$
where
$$
Z = \int_{-\infty}^{\infty} p(y|\theta)p(\theta)d\theta.
$$
Let's first calculate $Z$ analytically.


```python
import sympy.stats
sympy.init_printing()
y, t = sympy.symbols('y \\theta')
q = 1. / sympy.sqrt(2. * sympy.pi) * sympy.exp(-0.5 * (y - t) ** 2) * \
    1. / sympy.sqrt(2. * sympy.pi) * sympy.exp(-0.5 * t ** 2)
sympy.simplify(sympy.integrate(q, (t, -sympy.oo, sympy.oo)))
```

So, if the observed $y$ was zero, then the Z should be:
$$
Z = \frac{1}{2\sqrt{\pi}}.
$$


```python
Z = 1 /  2. / np.sqrt(np.pi)
print('log Z = {0:.3f}'.format(np.log(Z)))
```

    log Z = -1.266


All, right. Now, let's program this thing in pysmc and compare the results.

We start with the model:


```python
model = pm.Model()
yobs = 0.
with model:
    # prior over theta 
    theta = pm.Normal('theta', mu=0., sigma=1.,testval=0.)
    
    # log likelihood 
    llk = pm.Potential('llk', pm.Normal.dist(theta, 1.).logp(yobs))
    
    trace, smcres = sample_smc(1000)
```

    Sample initial stage: ...
    Stage:   0 Beta: 1.000 Steps:  25 Acce: 1.000



```python
# get the model evidence 
log_evidence_smc = np.log(smcres.model.marginal_likelihood)
print('True log evidence: %.4f \nSMC log evidence: %.4f'%(np.log(Z), log_evidence_smc))
```

    True log evidence: -1.2655 
    SMC log evidence: -1.2880


Which is close to the truth.

### Questions

+ Repeat the calculations above for a varying number of SMC particles. Start from 10 and go up to 10,000.

## Polynomial Regression


```python
def compute_design_matrix(X, phi):
    """
    Arguments:
    
    X   -  The observed inputs (1D array)
    phi -  The basis functions.
    """
    num_observations = X.shape[0]
    num_basis = phi.num_basis
    Phi = np.ndarray((num_observations, num_basis))
    for i in range(num_observations):
        Phi[i, :] = phi(X[i, :])
    return Phi

class PolynomialBasis(object):
    """
    A set of linear basis functions.
    
    Arguments:
    degree  -  The degree of the polynomial.
    """
    def __init__(self, degree):
        self.degree = degree
        self.num_basis = degree + 1
    def __call__(self, x):
        return np.array([x[0] ** i for i in range(self.degree + 1)])
    

class FourierBasis(object):
    """
    A set of linear basis functions.
    
    Arguments:
    num_terms  -  The number of Fourier terms.
    L          -  The period of the function.
    """
    def __init__(self, num_terms, L):
        self.num_terms = num_terms
        self.L = L
        self.num_basis = 2 * num_terms
    def __call__(self, x):
        res = np.ndarray((self.num_basis,))
        for i in range(num_terms):
            res[2 * i] = np.cos(2 * i * np.pi / self.L * x[0])
            res[2 * i + 1] = np.sin(2 * (i+1) * np.pi / self.L * x[0])
        return res
    

class RadialBasisFunctions(object):
    """
    A set of linear basis functions.
    
    Arguments:
    X   -  The centers of the radial basis functions.
    ell -  The assumed lengthscale.
    """
    def __init__(self, X, ell):
        self.X = X
        self.ell = ell
        self.num_basis = X.shape[0]
    def __call__(self, x):
        return np.exp(-.5 * (x - self.X) ** 2 / self.ell ** 2).flatten()
```

Let's generate some fake data.


```python
np.random.seed(12345)
def getdata(N, sigma2):
    X = 2 * np.random.rand(N) - 1.
    y = 0.5 * X  ** 3 - 0.3 * X ** 2 + np.sqrt(sigma2) * np.random.rand(N)
    return X, y

num_samples = 50
sigma2 = 1e-3
X, y = getdata(num_samples, sigma2)
plt.figure(figsize=(8, 6))
plt.plot(X, y, 'o', markeredgewidth=2, label='Data')
plt.xlabel('$x$', fontsize=20)
plt.ylabel('$y$', fontsize=20)
plt.legend(loc='best', fontsize=20)
plt.tight_layout()
```

We are going to implement a standard Bayesian linear regression and train it with `PyMC3`.
We will compute the evidence in order to select the best class of basis functions.
The model is as follows:

The output $y$ conditioned on the input $x$, the weights of the basis functions $w$ and
the noise variance $\sigma^2$ has likelihood:
$$
p(y|x,w,\sigma, \mathcal{M}) = \mathcal{N}(y|w^T\phi_{\mathcal{M}}(x), \sigma^2),
$$
where $\phi_{\mathcal{M},1}(\cdot), \dots, \phi_{\mathcal{M},m_{\mathcal{M}}}(\cdot)$ are the
$m_{\mathcal{M}}$ basis functions of the model $\mathcal{M}$.
We put a normal prior on the weights:
$$
p(w|\alpha) = \mathcal{N}(w|0, \alpha I_{m_{\mathcal{M}}}),
$$
and an inverse Gamma prior for $\sigma$ and $\alpha$:
$$
p(\sigma^2) = \mathrm{IG}(\sigma^2|1, 1),
$$
and
$$
p(\alpha) = \mathrm{IG}(\alpha|1,1).
$$

Assume that the data we have observed are:
$$
x_{1:n} = \{x_1,\dots,x_n\},\;\mathrm{and}\;y_{1:n} = \{y_1,\dots,y_n\}.
$$
Consider the design matrix $\Phi_{\mathcal{M}}\in\mathbb{R}^{n\times m}$:
$$
\Phi_{\mathcal{M},ij} = \phi_{\mathcal{M},j}(x_i).
$$
The likelihood of the data is:
$$
p(y_{1:n} | x_{1:n}, w, \sigma, \mathcal{M}) = \mathcal{N}(y_{1:n}|\Phi_{\mathcal{M}}w, \sigma^2I_n).
$$
Let's turn this into `PyMC3` code.


```python
def make_model(Phi, y):
    """
    INPUTS:
        Phi -> Design matrix.
        y   -> Target vector. 
        
    RETURNS:
        model -> `pymc3.model` context.
    """
    num_data, num_features = Phi.shape
    
    # define the model 
    with pm.Model() as model:
        # prior on the weights 
        alpha = pm.InverseGamma('alpha', alpha=1., beta=1.)
        w = pm.Normal('w', mu=0., tau=alpha, shape=num_features)
        
        # prior on the likelihood noise variance 
        sigma2 = pm.InverseGamma('sigma2', alpha=5., beta=0.1)
        
        # the data likelihood mean 
        ymean = pm.Deterministic('ymean', tt.dot(Phi, w))
        
        # likelihood 
        y = pm.Normal('y', ymean, sigma2, shape=num_data, observed=y)
        #llk = pm.Potential('llk', pm.Normal.dist(ymean, tt.sqrt(sigma2)).logp_sum(y))
    return model
```

Now, let's create a function that trains the model using pysmc for a polynomial basis with a given order.


```python
def fit_poly(phi, X, y, num_particles=100):
    """
    
    RETURNS:
        1. An instance of pymc3.Model for the SMC model.
        2. The SMC trace.
        3. An instance of pymc3.smc.SMC containing sampling information.
    """
    Phi = compute_design_matrix(X[:, None], phi)
    smcmodel = make_model(Phi, y)
    trace, res = sample_smc(draws=num_particles, 
                         model=smcmodel, 
                         progressbar=True, 
                         threshold=0.8)
    return smcmodel, trace, res

phi = PolynomialBasis(3)
model, trace, res = fit_poly(phi, X, y)
```

    Sample initial stage: ...
    Stage:   0 Beta: 0.000 Steps:  25 Acce: 1.000
    
      0%|          | 0/100 [00:00<?, ?it/s][A
    100%|██████████| 100/100 [00:00<00:00, 933.04it/s][A
    Stage:   1 Beta: 0.000 Steps:  25 Acce: 0.155
    
    100%|██████████| 100/100 [00:00<00:00, 1061.16it/s]
    Stage:   2 Beta: 0.000 Steps:  25 Acce: 0.159
    
    100%|██████████| 100/100 [00:00<00:00, 1048.84it/s]
    Stage:   3 Beta: 0.000 Steps:  25 Acce: 0.197
    
    100%|██████████| 100/100 [00:00<00:00, 1242.81it/s]
    Stage:   4 Beta: 0.001 Steps:  21 Acce: 0.257
    
    100%|██████████| 100/100 [00:00<00:00, 1725.92it/s]
    Stage:   5 Beta: 0.001 Steps:  15 Acce: 0.289
    
    100%|██████████| 100/100 [00:00<00:00, 1986.61it/s]
    Stage:   6 Beta: 0.002 Steps:  13 Acce: 0.276
    
    100%|██████████| 100/100 [00:00<00:00, 1770.82it/s]
    Stage:   7 Beta: 0.004 Steps:  14 Acce: 0.211
    
    100%|██████████| 100/100 [00:00<00:00, 1358.40it/s]
    Stage:   8 Beta: 0.006 Steps:  19 Acce: 0.228
    
    100%|██████████| 100/100 [00:00<00:00, 1510.53it/s]
    Stage:   9 Beta: 0.011 Steps:  17 Acce: 0.246
    
    100%|██████████| 100/100 [00:00<00:00, 1636.34it/s]
    Stage:  10 Beta: 0.017 Steps:  16 Acce: 0.196
    
    100%|██████████| 100/100 [00:00<00:00, 1238.94it/s]
    Stage:  11 Beta: 0.025 Steps:  21 Acce: 0.250
    
    100%|██████████| 100/100 [00:00<00:00, 1612.00it/s]
    Stage:  12 Beta: 0.035 Steps:  16 Acce: 0.244
    
    100%|██████████| 100/100 [00:00<00:00, 1618.49it/s]
    Stage:  13 Beta: 0.050 Steps:  16 Acce: 0.196
    
    100%|██████████| 100/100 [00:00<00:00, 1229.92it/s]
    Stage:  14 Beta: 0.068 Steps:  21 Acce: 0.285
    
    100%|██████████| 100/100 [00:00<00:00, 1966.60it/s]
    Stage:  15 Beta: 0.090 Steps:  13 Acce: 0.313
    
    100%|██████████| 100/100 [00:00<00:00, 2164.50it/s]
    Stage:  16 Beta: 0.116 Steps:  12 Acce: 0.338
    
    100%|██████████| 100/100 [00:00<00:00, 2285.18it/s]
    Stage:  17 Beta: 0.151 Steps:  11 Acce: 0.289
    
    100%|██████████| 100/100 [00:00<00:00, 1935.20it/s]
    Stage:  18 Beta: 0.192 Steps:  13 Acce: 0.278
    
    100%|██████████| 100/100 [00:00<00:00, 1840.10it/s]
    Stage:  19 Beta: 0.240 Steps:  14 Acce: 0.272
    
    100%|██████████| 100/100 [00:00<00:00, 1826.72it/s]
    Stage:  20 Beta: 0.298 Steps:  14 Acce: 0.264
    
    100%|██████████| 100/100 [00:00<00:00, 1681.72it/s]
    Stage:  21 Beta: 0.378 Steps:  15 Acce: 0.276
    
    100%|██████████| 100/100 [00:00<00:00, 1877.26it/s]
    Stage:  22 Beta: 0.479 Steps:  14 Acce: 0.289
    
    100%|██████████| 100/100 [00:00<00:00, 1988.47it/s]
    Stage:  23 Beta: 0.610 Steps:  13 Acce: 0.272
    
    100%|██████████| 100/100 [00:00<00:00, 1843.80it/s]
    Stage:  24 Beta: 0.816 Steps:  14 Acce: 0.252
    
    100%|██████████| 100/100 [00:00<00:00, 1749.14it/s]
    Stage:  25 Beta: 1.000 Steps:  15 Acce: 0.261
    
    100%|██████████| 100/100 [00:00<00:00, 1668.52it/s]


### Postprocessing 

Once you have the `trace` object for the SMC simulation you can apply all the standard postprocessing tools from `PyMC3` as usual. 

Here's the posterior distribution over the weights precision and the the likelihood noise, $\alpha$ and $\sigma^2$ respectively:


```python
_=pm.plot_posterior(trace, var_names=['alpha', 'sigma2'])
```

Here's the posterior predictive mean of the output $y$, i.e., $\mathbb{E}[y|x, w, \sigma]$:


```python
ppsamples = pm.sample_posterior_predictive(model=model, 
                               trace=trace, var_names=['ymean'])['ymean']
```

    
    100%|██████████| 100/100 [00:00<00:00, 3874.47it/s]



```python
idx = np.argsort(X)
plt.figure(figsize=(10, 8))
plt.plot(X[idx], ppsamples.mean(0)[idx], linewidth=2.5, label='Posterior Predictive Mean' )
plt.plot(X, y, 'x', markeredgewidth=2.5, markersize=10, label='Observed data')

plt.legend(loc='best', fontsize=20)
plt.tight_layout()
```

SMC does a particle approximation of the posterior distribution. The particles themselves can be obtained from the `trace` object and the particle weights can be obtained from the  `res` object. 

Recall that the approximate posterior distribution is of the form $p(\theta|\mathcal{D}) = \sum_{j=1}^{N} w_j \delta(\theta - \theta_j)$.


```python
particles_w = trace.w
particles_alpha = trace.alpha
particle_weights = res.weights   # <- these are the ws from the above equation
```

## Model comparison

Since SMC can approximate the model evidence it provides a principled way of comparing models. Let's compare 5 different polynomial regression models where we change the degree of the polynomial from 1 to 5. 


```python
# Evaluate the evidence for the various degrees
log_Zs = []
D = [1, 2, 3, 4, 5]
for d in D:
    phi = PolynomialBasis(d)
    _, _, res = fit_poly(phi, X, y, num_particles=500)
    log_Z = np.log(res.model.marginal_likelihood)
    log_Zs.append(log_Z)
```

    Sample initial stage: ...
    Stage:   0 Beta: 0.000 Steps:  25 Acce: 1.000
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     21%|██▏       | 107/500 [00:00<00:00, 1069.53it/s][A
     43%|████▎     | 214/500 [00:00<00:00, 1067.47it/s][A
     64%|██████▍   | 319/500 [00:00<00:00, 1061.39it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1060.87it/s][A
    Stage:   1 Beta: 0.000 Steps:  25 Acce: 0.288
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     41%|████      | 203/500 [00:00<00:00, 2027.15it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2003.37it/s][A
    Stage:   2 Beta: 0.000 Steps:  13 Acce: 0.259
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 180/500 [00:00<00:00, 1799.16it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1763.08it/s][A
    Stage:   3 Beta: 0.000 Steps:  15 Acce: 0.271
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     37%|███▋      | 187/500 [00:00<00:00, 1869.01it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1860.91it/s][A
    Stage:   4 Beta: 0.001 Steps:  14 Acce: 0.256
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 179/500 [00:00<00:00, 1787.86it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1784.59it/s][A
    Stage:   5 Beta: 0.001 Steps:  15 Acce: 0.259
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 177/500 [00:00<00:00, 1761.36it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1753.73it/s][A
    Stage:   6 Beta: 0.002 Steps:  15 Acce: 0.223
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     30%|███       | 151/500 [00:00<00:00, 1508.31it/s][A
     60%|█████▉    | 298/500 [00:00<00:00, 1493.93it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1473.24it/s][A
    Stage:   7 Beta: 0.003 Steps:  18 Acce: 0.242
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▎      | 168/500 [00:00<00:00, 1676.20it/s][A
     67%|██████▋   | 334/500 [00:00<00:00, 1670.55it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1644.80it/s][A
    Stage:   8 Beta: 0.006 Steps:  16 Acce: 0.228
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     32%|███▏      | 159/500 [00:00<00:00, 1584.49it/s][A
     63%|██████▎   | 315/500 [00:00<00:00, 1575.19it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1561.86it/s][A
    Stage:   9 Beta: 0.009 Steps:  17 Acce: 0.235
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     32%|███▏      | 158/500 [00:00<00:00, 1576.07it/s][A
     63%|██████▎   | 314/500 [00:00<00:00, 1570.26it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1550.86it/s][A
    Stage:  10 Beta: 0.015 Steps:  17 Acce: 0.234
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███▏      | 157/500 [00:00<00:00, 1567.88it/s][A
     62%|██████▏   | 311/500 [00:00<00:00, 1556.50it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1535.76it/s][A
    Stage:  11 Beta: 0.024 Steps:  17 Acce: 0.246
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 166/500 [00:00<00:00, 1654.25it/s][A
     65%|██████▌   | 327/500 [00:00<00:00, 1638.17it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1618.74it/s][A
    Stage:  12 Beta: 0.038 Steps:  16 Acce: 0.235
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 154/500 [00:00<00:00, 1537.05it/s][A
     62%|██████▏   | 309/500 [00:00<00:00, 1539.67it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1536.17it/s][A
    Stage:  13 Beta: 0.061 Steps:  17 Acce: 0.260
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 177/500 [00:00<00:00, 1762.24it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1745.45it/s][A
    Stage:  14 Beta: 0.095 Steps:  15 Acce: 0.237
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 154/500 [00:00<00:00, 1538.01it/s][A
     60%|██████    | 302/500 [00:00<00:00, 1518.98it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1519.69it/s][A
    Stage:  15 Beta: 0.149 Steps:  17 Acce: 0.251
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 176/500 [00:00<00:00, 1751.92it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1734.64it/s][A
    Stage:  16 Beta: 0.239 Steps:  15 Acce: 0.243
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 167/500 [00:00<00:00, 1661.72it/s][A
     66%|██████▌   | 330/500 [00:00<00:00, 1650.53it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1639.98it/s][A
    Stage:  17 Beta: 0.381 Steps:  16 Acce: 0.258
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 176/500 [00:00<00:00, 1756.82it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1755.52it/s][A
    Stage:  18 Beta: 0.600 Steps:  15 Acce: 0.246
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 167/500 [00:00<00:00, 1661.42it/s][A
     66%|██████▌   | 328/500 [00:00<00:00, 1643.73it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1624.14it/s][A
    Stage:  19 Beta: 0.947 Steps:  16 Acce: 0.291
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     40%|████      | 201/500 [00:00<00:00, 2009.38it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1991.53it/s][A
    Stage:  20 Beta: 1.000 Steps:  13 Acce: 0.230
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 156/500 [00:00<00:00, 1552.94it/s][A
     61%|██████    | 306/500 [00:00<00:00, 1535.65it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1505.63it/s][A
    Sample initial stage: ...
    Stage:   0 Beta: 0.000 Steps:  25 Acce: 1.000
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     22%|██▏       | 111/500 [00:00<00:00, 1102.46it/s][A
     44%|████▍     | 219/500 [00:00<00:00, 1094.13it/s][A
     65%|██████▌   | 327/500 [00:00<00:00, 1087.08it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1077.47it/s][A
    Stage:   1 Beta: 0.000 Steps:  25 Acce: 0.222
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     30%|███       | 151/500 [00:00<00:00, 1501.80it/s][A
     59%|█████▉    | 296/500 [00:00<00:00, 1482.92it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1464.64it/s][A
    Stage:   2 Beta: 0.000 Steps:  18 Acce: 0.250
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▍      | 170/500 [00:00<00:00, 1691.09it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1677.64it/s][A
    Stage:   3 Beta: 0.000 Steps:  16 Acce: 0.237
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     32%|███▏      | 159/500 [00:00<00:00, 1586.73it/s][A
     63%|██████▎   | 314/500 [00:00<00:00, 1573.64it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1559.97it/s][A
    Stage:   4 Beta: 0.000 Steps:  17 Acce: 0.224
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     30%|███       | 150/500 [00:00<00:00, 1491.62it/s][A
     59%|█████▉    | 297/500 [00:00<00:00, 1478.53it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1462.24it/s][A
    Stage:   5 Beta: 0.001 Steps:  18 Acce: 0.225
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     30%|██▉       | 149/500 [00:00<00:00, 1482.93it/s][A
     60%|█████▉    | 299/500 [00:00<00:00, 1485.79it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1482.70it/s][A
    Stage:   6 Beta: 0.002 Steps:  18 Acce: 0.208
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     28%|██▊       | 142/500 [00:00<00:00, 1413.95it/s][A
     56%|█████▋    | 282/500 [00:00<00:00, 1407.08it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1399.51it/s][A
    Stage:   7 Beta: 0.003 Steps:  19 Acce: 0.201
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     26%|██▋       | 132/500 [00:00<00:00, 1315.07it/s][A
     54%|█████▎    | 268/500 [00:00<00:00, 1325.73it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1335.19it/s][A
    Stage:   8 Beta: 0.005 Steps:  20 Acce: 0.222
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     30%|██▉       | 149/500 [00:00<00:00, 1482.87it/s][A
     60%|█████▉    | 298/500 [00:00<00:00, 1484.76it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1484.78it/s][A
    Stage:   9 Beta: 0.008 Steps:  18 Acce: 0.233
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     32%|███▏      | 159/500 [00:00<00:00, 1588.49it/s][A
     63%|██████▎   | 315/500 [00:00<00:00, 1576.88it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1554.79it/s][A
    Stage:  10 Beta: 0.013 Steps:  17 Acce: 0.241
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▎      | 168/500 [00:00<00:00, 1674.23it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1667.77it/s][A
    Stage:  11 Beta: 0.021 Steps:  16 Acce: 0.243
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▍      | 170/500 [00:00<00:00, 1694.05it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1670.03it/s][A
    Stage:  12 Beta: 0.033 Steps:  16 Acce: 0.262
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 179/500 [00:00<00:00, 1785.59it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1778.25it/s][A
    Stage:  13 Beta: 0.052 Steps:  15 Acce: 0.257
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 181/500 [00:00<00:00, 1802.70it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1776.71it/s][A
    Stage:  14 Beta: 0.079 Steps:  15 Acce: 0.243
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▍      | 169/500 [00:00<00:00, 1682.21it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1667.50it/s][A
    Stage:  15 Beta: 0.120 Steps:  16 Acce: 0.275
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 190/500 [00:00<00:00, 1893.38it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1881.68it/s][A
    Stage:  16 Beta: 0.181 Steps:  14 Acce: 0.264
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 178/500 [00:00<00:00, 1777.73it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1775.35it/s][A
    Stage:  17 Beta: 0.264 Steps:  15 Acce: 0.276
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     39%|███▊      | 193/500 [00:00<00:00, 1927.86it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1895.63it/s][A
    Stage:  18 Beta: 0.391 Steps:  14 Acce: 0.266
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 190/500 [00:00<00:00, 1894.60it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1893.48it/s][A
    Stage:  19 Beta: 0.597 Steps:  14 Acce: 0.259
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▋      | 182/500 [00:00<00:00, 1814.00it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1784.45it/s][A
    Stage:  20 Beta: 0.861 Steps:  15 Acce: 0.260
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 178/500 [00:00<00:00, 1772.80it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1772.45it/s][A
    Stage:  21 Beta: 1.000 Steps:  15 Acce: 0.269
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     37%|███▋      | 185/500 [00:00<00:00, 1848.93it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1843.60it/s][A
    Sample initial stage: ...
    Stage:   0 Beta: 0.000 Steps:  25 Acce: 1.000
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     21%|██        | 106/500 [00:00<00:00, 1054.14it/s][A
     42%|████▏     | 212/500 [00:00<00:00, 1054.90it/s][A
     64%|██████▍   | 320/500 [00:00<00:00, 1059.56it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1063.03it/s][A
    Stage:   1 Beta: 0.000 Steps:  25 Acce: 0.179
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     23%|██▎       | 116/500 [00:00<00:00, 1156.89it/s][A
     46%|████▋     | 232/500 [00:00<00:00, 1157.22it/s][A
     70%|██████▉   | 348/500 [00:00<00:00, 1155.89it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1150.37it/s][A
    Stage:   2 Beta: 0.000 Steps:  23 Acce: 0.183
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     24%|██▍       | 121/500 [00:00<00:00, 1203.27it/s][A
     48%|████▊     | 239/500 [00:00<00:00, 1193.19it/s][A
     72%|███████▏  | 359/500 [00:00<00:00, 1194.53it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1187.77it/s][A
    Stage:   3 Beta: 0.000 Steps:  22 Acce: 0.229
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███▏      | 157/500 [00:00<00:00, 1569.56it/s][A
     62%|██████▏   | 310/500 [00:00<00:00, 1554.74it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1535.60it/s][A
    Stage:   4 Beta: 0.000 Steps:  17 Acce: 0.205
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     27%|██▋       | 134/500 [00:00<00:00, 1331.19it/s][A
     53%|█████▎    | 263/500 [00:00<00:00, 1317.43it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1313.46it/s][A
    Stage:   5 Beta: 0.001 Steps:  20 Acce: 0.220
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     29%|██▉       | 145/500 [00:00<00:00, 1448.67it/s][A
     57%|█████▋    | 287/500 [00:00<00:00, 1439.91it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1442.99it/s][A
    Stage:   6 Beta: 0.002 Steps:  18 Acce: 0.209
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     28%|██▊       | 139/500 [00:00<00:00, 1384.92it/s][A
     55%|█████▌    | 275/500 [00:00<00:00, 1375.67it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1364.39it/s][A
    Stage:   7 Beta: 0.003 Steps:  19 Acce: 0.213
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     27%|██▋       | 136/500 [00:00<00:00, 1355.78it/s][A
     55%|█████▍    | 274/500 [00:00<00:00, 1361.41it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1357.11it/s][A
    Stage:   8 Beta: 0.004 Steps:  19 Acce: 0.244
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     32%|███▏      | 159/500 [00:00<00:00, 1587.49it/s][A
     64%|██████▍   | 321/500 [00:00<00:00, 1594.92it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1603.08it/s][A
    Stage:   9 Beta: 0.008 Steps:  16 Acce: 0.255
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▍      | 174/500 [00:00<00:00, 1733.17it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1730.94it/s][A
    Stage:  10 Beta: 0.012 Steps:  15 Acce: 0.250
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 176/500 [00:00<00:00, 1754.75it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1729.70it/s][A
    Stage:  11 Beta: 0.019 Steps:  15 Acce: 0.248
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 164/500 [00:00<00:00, 1634.20it/s][A
     66%|██████▌   | 329/500 [00:00<00:00, 1638.88it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1621.96it/s][A
    Stage:  12 Beta: 0.027 Steps:  16 Acce: 0.228
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 154/500 [00:00<00:00, 1532.29it/s][A
     61%|██████    | 306/500 [00:00<00:00, 1525.70it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1520.49it/s][A
    Stage:  13 Beta: 0.037 Steps:  17 Acce: 0.260
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 176/500 [00:00<00:00, 1753.18it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1744.77it/s][A
    Stage:  14 Beta: 0.049 Steps:  15 Acce: 0.269
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 188/500 [00:00<00:00, 1873.02it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1856.05it/s][A
    Stage:  15 Beta: 0.064 Steps:  14 Acce: 0.259
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 176/500 [00:00<00:00, 1750.41it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1744.21it/s][A
    Stage:  16 Beta: 0.082 Steps:  15 Acce: 0.301
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▍     | 222/500 [00:00<00:00, 2216.88it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2160.81it/s][A
    Stage:  17 Beta: 0.104 Steps:  12 Acce: 0.325
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     48%|████▊     | 238/500 [00:00<00:00, 2379.86it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2340.54it/s][A
    Stage:  18 Beta: 0.130 Steps:  11 Acce: 0.306
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▎     | 218/500 [00:00<00:00, 2169.89it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2133.08it/s][A
    Stage:  19 Beta: 0.166 Steps:  12 Acce: 0.339
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     47%|████▋     | 237/500 [00:00<00:00, 2368.87it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2349.21it/s][A
    Stage:  20 Beta: 0.209 Steps:  11 Acce: 0.346
    
      0%|          | 0/500 [00:00<?, ?it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2540.49it/s][A
    Stage:  21 Beta: 0.264 Steps:  10 Acce: 0.299
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▎     | 218/500 [00:00<00:00, 2172.21it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2135.84it/s][A
    Stage:  22 Beta: 0.337 Steps:  12 Acce: 0.277
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     37%|███▋      | 186/500 [00:00<00:00, 1851.76it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1859.78it/s][A
    Stage:  23 Beta: 0.433 Steps:  14 Acce: 0.283
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     40%|███▉      | 199/500 [00:00<00:00, 1983.40it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1983.67it/s][A
    Stage:  24 Beta: 0.556 Steps:  13 Acce: 0.271
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     37%|███▋      | 186/500 [00:00<00:00, 1857.38it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1838.67it/s][A
    Stage:  25 Beta: 0.716 Steps:  14 Acce: 0.246
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 165/500 [00:00<00:00, 1646.21it/s][A
     65%|██████▌   | 327/500 [00:00<00:00, 1638.23it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1622.28it/s][A
    Stage:  26 Beta: 0.938 Steps:  16 Acce: 0.251
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▍      | 171/500 [00:00<00:00, 1700.48it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1721.75it/s][A
    Stage:  27 Beta: 1.000 Steps:  15 Acce: 0.245
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     32%|███▏      | 161/500 [00:00<00:00, 1609.55it/s][A
     64%|██████▍   | 319/500 [00:00<00:00, 1598.80it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1577.64it/s][A
    Sample initial stage: ...
    Stage:   0 Beta: 0.000 Steps:  25 Acce: 1.000
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     22%|██▏       | 108/500 [00:00<00:00, 1072.85it/s][A
     42%|████▏     | 212/500 [00:00<00:00, 1060.37it/s][A
     63%|██████▎   | 317/500 [00:00<00:00, 1056.47it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1054.29it/s][A
    Stage:   1 Beta: 0.000 Steps:  25 Acce: 0.172
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     22%|██▏       | 112/500 [00:00<00:00, 1118.55it/s][A
     44%|████▍     | 222/500 [00:00<00:00, 1111.77it/s][A
     67%|██████▋   | 333/500 [00:00<00:00, 1111.16it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1102.35it/s][A
    Stage:   2 Beta: 0.000 Steps:  24 Acce: 0.200
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     26%|██▋       | 132/500 [00:00<00:00, 1318.88it/s][A
     52%|█████▏    | 262/500 [00:00<00:00, 1310.00it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1310.18it/s][A
    Stage:   3 Beta: 0.000 Steps:  20 Acce: 0.225
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     29%|██▉       | 145/500 [00:00<00:00, 1449.41it/s][A
     58%|█████▊    | 288/500 [00:00<00:00, 1440.46it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1432.79it/s][A
    Stage:   4 Beta: 0.000 Steps:  18 Acce: 0.237
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     32%|███▏      | 162/500 [00:00<00:00, 1619.89it/s][A
     64%|██████▍   | 322/500 [00:00<00:00, 1611.07it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1602.81it/s][A
    Stage:   5 Beta: 0.001 Steps:  16 Acce: 0.239
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 163/500 [00:00<00:00, 1629.53it/s][A
     65%|██████▍   | 323/500 [00:00<00:00, 1619.02it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1607.04it/s][A
    Stage:   6 Beta: 0.001 Steps:  16 Acce: 0.261
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▍      | 172/500 [00:00<00:00, 1715.83it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1730.58it/s][A
    Stage:   7 Beta: 0.002 Steps:  15 Acce: 0.272
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     37%|███▋      | 185/500 [00:00<00:00, 1847.38it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1831.29it/s][A
    Stage:   8 Beta: 0.004 Steps:  14 Acce: 0.235
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 155/500 [00:00<00:00, 1547.42it/s][A
     62%|██████▏   | 308/500 [00:00<00:00, 1542.07it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1522.64it/s][A
    Stage:   9 Beta: 0.007 Steps:  17 Acce: 0.234
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 155/500 [00:00<00:00, 1548.53it/s][A
     62%|██████▏   | 310/500 [00:00<00:00, 1546.09it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1531.79it/s][A
    Stage:  10 Beta: 0.011 Steps:  17 Acce: 0.248
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 163/500 [00:00<00:00, 1624.51it/s][A
     65%|██████▌   | 325/500 [00:00<00:00, 1622.72it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1618.81it/s][A
    Stage:  11 Beta: 0.017 Steps:  16 Acce: 0.226
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 154/500 [00:00<00:00, 1536.28it/s][A
     61%|██████    | 303/500 [00:00<00:00, 1521.00it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1511.12it/s][A
    Stage:  12 Beta: 0.026 Steps:  17 Acce: 0.247
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 164/500 [00:00<00:00, 1635.40it/s][A
     65%|██████▌   | 325/500 [00:00<00:00, 1626.33it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1612.26it/s][A
    Stage:  13 Beta: 0.037 Steps:  16 Acce: 0.277
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 189/500 [00:00<00:00, 1885.08it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1849.35it/s][A
    Stage:  14 Beta: 0.051 Steps:  14 Acce: 0.272
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     37%|███▋      | 185/500 [00:00<00:00, 1848.78it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1839.42it/s][A
    Stage:  15 Beta: 0.066 Steps:  14 Acce: 0.262
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 156/500 [00:00<00:00, 1552.66it/s][A
     65%|██████▌   | 326/500 [00:00<00:00, 1593.78it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1647.52it/s][A
    Stage:  16 Beta: 0.084 Steps:  15 Acce: 0.266
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 191/500 [00:00<00:00, 1908.80it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1870.24it/s][A
    Stage:  17 Beta: 0.105 Steps:  14 Acce: 0.279
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     37%|███▋      | 184/500 [00:00<00:00, 1839.96it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1839.83it/s][A
    Stage:  18 Beta: 0.130 Steps:  14 Acce: 0.292
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     40%|████      | 200/500 [00:00<00:00, 1994.91it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1979.62it/s][A
    Stage:  19 Beta: 0.159 Steps:  13 Acce: 0.289
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     41%|████      | 203/500 [00:00<00:00, 2024.84it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1999.22it/s][A
    Stage:  20 Beta: 0.198 Steps:  13 Acce: 0.307
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▍     | 219/500 [00:00<00:00, 2189.91it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2167.67it/s][A
    Stage:  21 Beta: 0.244 Steps:  12 Acce: 0.317
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     43%|████▎     | 216/500 [00:00<00:00, 2152.30it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2141.15it/s][A
    Stage:  22 Beta: 0.301 Steps:  12 Acce: 0.306
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▎     | 218/500 [00:00<00:00, 2172.72it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2155.81it/s][A
    Stage:  23 Beta: 0.379 Steps:  12 Acce: 0.285
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     39%|███▉      | 197/500 [00:00<00:00, 1962.11it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1962.52it/s][A
    Stage:  24 Beta: 0.473 Steps:  13 Acce: 0.299
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▍     | 221/500 [00:00<00:00, 2204.88it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2136.37it/s][A
    Stage:  25 Beta: 0.609 Steps:  12 Acce: 0.262
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 175/500 [00:00<00:00, 1742.79it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1725.43it/s][A
    Stage:  26 Beta: 0.774 Steps:  15 Acce: 0.260
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▍      | 174/500 [00:00<00:00, 1737.90it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1744.96it/s][A
    Stage:  27 Beta: 1.000 Steps:  15 Acce: 0.255
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▍      | 170/500 [00:00<00:00, 1696.58it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1673.65it/s][A
    Sample initial stage: ...
    Stage:   0 Beta: 0.000 Steps:  25 Acce: 1.000
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     22%|██▏       | 109/500 [00:00<00:00, 1086.97it/s][A
     43%|████▎     | 217/500 [00:00<00:00, 1083.80it/s][A
     65%|██████▌   | 325/500 [00:00<00:00, 1081.61it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1071.84it/s][A
    Stage:   1 Beta: 0.000 Steps:  25 Acce: 0.185
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     25%|██▌       | 125/500 [00:00<00:00, 1242.61it/s][A
     49%|████▉     | 246/500 [00:00<00:00, 1232.08it/s][A
     74%|███████▎  | 368/500 [00:00<00:00, 1226.11it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1214.16it/s][A
    Stage:   2 Beta: 0.000 Steps:  22 Acce: 0.199
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     27%|██▋       | 133/500 [00:00<00:00, 1327.45it/s][A
     53%|█████▎    | 266/500 [00:00<00:00, 1325.57it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1321.80it/s][A
    Stage:   3 Beta: 0.000 Steps:  20 Acce: 0.202
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     27%|██▋       | 134/500 [00:00<00:00, 1338.78it/s][A
     53%|█████▎    | 267/500 [00:00<00:00, 1333.01it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1330.21it/s][A
    Stage:   4 Beta: 0.000 Steps:  20 Acce: 0.227
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███▏      | 157/500 [00:00<00:00, 1569.09it/s][A
     63%|██████▎   | 316/500 [00:00<00:00, 1572.47it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1561.78it/s][A
    Stage:   5 Beta: 0.001 Steps:  17 Acce: 0.248
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▍      | 169/500 [00:00<00:00, 1681.45it/s][A
     67%|██████▋   | 333/500 [00:00<00:00, 1666.36it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1641.69it/s][A
    Stage:   6 Beta: 0.001 Steps:  16 Acce: 0.252
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 178/500 [00:00<00:00, 1775.74it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1777.11it/s][A
    Stage:   7 Beta: 0.002 Steps:  15 Acce: 0.272
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 190/500 [00:00<00:00, 1898.56it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1889.00it/s][A
    Stage:   8 Beta: 0.004 Steps:  14 Acce: 0.256
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 181/500 [00:00<00:00, 1806.21it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1776.93it/s][A
    Stage:   9 Beta: 0.006 Steps:  15 Acce: 0.254
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 179/500 [00:00<00:00, 1789.59it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1762.80it/s][A
    Stage:  10 Beta: 0.010 Steps:  15 Acce: 0.245
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 166/500 [00:00<00:00, 1656.20it/s][A
     67%|██████▋   | 333/500 [00:00<00:00, 1657.09it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1653.33it/s][A
    Stage:  11 Beta: 0.016 Steps:  16 Acce: 0.232
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███       | 155/500 [00:00<00:00, 1547.77it/s][A
     62%|██████▏   | 312/500 [00:00<00:00, 1552.81it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1543.57it/s][A
    Stage:  12 Beta: 0.025 Steps:  17 Acce: 0.226
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     31%|███▏      | 157/500 [00:00<00:00, 1561.83it/s][A
     63%|██████▎   | 313/500 [00:00<00:00, 1558.66it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1547.65it/s][A
    Stage:  13 Beta: 0.037 Steps:  17 Acce: 0.257
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▌      | 177/500 [00:00<00:00, 1765.20it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1759.77it/s][A
    Stage:  14 Beta: 0.050 Steps:  15 Acce: 0.247
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     33%|███▎      | 166/500 [00:00<00:00, 1656.54it/s][A
     66%|██████▌   | 330/500 [00:00<00:00, 1650.65it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1647.94it/s][A
    Stage:  15 Beta: 0.065 Steps:  16 Acce: 0.268
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 188/500 [00:00<00:00, 1876.47it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1866.43it/s][A
    Stage:  16 Beta: 0.081 Steps:  14 Acce: 0.243
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     34%|███▎      | 168/500 [00:00<00:00, 1671.31it/s][A
     66%|██████▌   | 331/500 [00:00<00:00, 1658.09it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1636.48it/s][A
    Stage:  17 Beta: 0.100 Steps:  16 Acce: 0.276
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 192/500 [00:00<00:00, 1917.25it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1884.04it/s][A
    Stage:  18 Beta: 0.124 Steps:  14 Acce: 0.298
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     41%|████      | 203/500 [00:00<00:00, 2026.37it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2000.13it/s][A
    Stage:  19 Beta: 0.153 Steps:  13 Acce: 0.300
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▍     | 220/500 [00:00<00:00, 2191.91it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2165.52it/s][A
    Stage:  20 Beta: 0.185 Steps:  12 Acce: 0.314
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▍     | 221/500 [00:00<00:00, 2201.72it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2173.88it/s][A
    Stage:  21 Beta: 0.221 Steps:  12 Acce: 0.298
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     44%|████▍     | 219/500 [00:00<00:00, 2184.89it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2168.14it/s][A
    Stage:  22 Beta: 0.264 Steps:  12 Acce: 0.288
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     41%|████      | 203/500 [00:00<00:00, 2027.77it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2002.81it/s][A
    Stage:  23 Beta: 0.317 Steps:  13 Acce: 0.297
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     40%|████      | 202/500 [00:00<00:00, 2013.44it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1997.53it/s][A
    Stage:  24 Beta: 0.387 Steps:  13 Acce: 0.292
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     41%|████      | 205/500 [00:00<00:00, 2044.67it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2013.66it/s][A
    Stage:  25 Beta: 0.472 Steps:  13 Acce: 0.268
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 190/500 [00:00<00:00, 1895.87it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1874.79it/s][A
    Stage:  26 Beta: 0.576 Steps:  14 Acce: 0.267
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 189/500 [00:00<00:00, 1880.88it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1860.04it/s][A
    Stage:  27 Beta: 0.723 Steps:  14 Acce: 0.274
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 192/500 [00:00<00:00, 1911.78it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1880.73it/s][A
    Stage:  28 Beta: 0.910 Steps:  14 Acce: 0.262
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     36%|███▌      | 179/500 [00:00<00:00, 1783.10it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1760.96it/s][A
    Stage:  29 Beta: 1.000 Steps:  15 Acce: 0.251
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     35%|███▍      | 173/500 [00:00<00:00, 1727.79it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1705.16it/s][A



```python
for d, log_Z in zip(D, log_Zs):
    print('degree %d gives %.4f'%(d, log_Z))
```

    degree 1 gives 25.0212
    degree 2 gives 41.0753
    degree 3 gives 142.4574
    degree 4 gives 139.3597
    degree 5 gives 137.5729



```python
plt.figure(figsize=(10, 5))
_=plt.bar(D, log_Zs, width=0.3)
_=plt.xticks(D)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel('Polynomial degree', fontsize=20)
plt.ylabel('Model Evidence', fontsize=20)
plt.tight_layout()
```

## Questions
+ The model with degree 3 polynomials has the gradest evidence. However, the degree 4 and 5 seem also very plausible. Is this a problem for the theory of Bayesian model selection? What complicates things here, is that model 3 is included in model 4 which is included in model 5. This requires us to design special priors for the models being right. They have to be consistent in some sense. For example, if model 3 is right then model 4 must be right, etc.

+ Revisit the motorcycle dataset problem. Evaluate the model evidence for a 1) Polynomial basis; 2) a Fourier basis; and 3) a Radial basis function basis.

## Revisiting Challenger Disaster Problem (Model Selection)


```python
challenger_data = np.genfromtxt("challenger_data.csv", skip_header=1,
                                usecols=[1, 2], missing_values="NA",
                                delimiter=",")
# drop the NA values
challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])]

# plot it, as a function of temperature (the first column)
print("Temp (F), O-Ring failure?")
print(challenger_data)
```

    Temp (F), O-Ring failure?
    [[66.  0.]
     [70.  1.]
     [69.  0.]
     [68.  0.]
     [67.  0.]
     [72.  0.]
     [73.  0.]
     [70.  0.]
     [57.  1.]
     [63.  1.]
     [70.  1.]
     [78.  0.]
     [67.  0.]
     [53.  1.]
     [67.  0.]
     [75.  0.]
     [70.  0.]
     [81.  0.]
     [76.  0.]
     [79.  0.]
     [75.  1.]
     [76.  0.]
     [58.  1.]]



```python
# plot it, as a function of temperature (the first column)
plt.figure(figsize=(12, 5))
plt.plot(challenger_data[:, 0], challenger_data[:, 1], 'ro', 
         markersize=15)
plt.ylabel("Damage Incident?",fontsize=20)
plt.xlabel("Outside temperature (Fahrenheit)",fontsize=20)
plt.title("Defects of the Space Shuttle O-Rings vs temperature",
          fontsize=20)
plt.yticks([0, 1], fontsize=15)
plt.xticks(fontsize=15)
plt.tight_layout()
```


```python
# gather the data and apply preprocessing if any 
temp = challenger_data[:, 0]
temp_scaled = (temp - np.mean(temp))/np.std(temp)
data = challenger_data[:, 1]

# instantiate the pymc3 model
challenger_model = pm.Model()

# define the graph 
with challenger_model:
    # define the prior
    alpha = pm.Normal('alpha', mu=0., sigma=10.)
    beta = pm.Normal('beta', mu=0., sigma=10.)
    
    # get the probabilities of failure at each observed temp 
    p = pm.Deterministic('p', 1./(1. + tt.exp(alpha + beta*temp_scaled)))
    
    # define the likelihood 
    x = pm.Bernoulli('x', p=p, observed=data)
print("Challenger space shuttle disaster model:")
challenger_model
```

    Challenger space shuttle disaster model:





$$
            \begin{array}{rcl}
            \text{alpha} &\sim & \text{Normal}(\mathit{mu}=0.0,~\mathit{sigma}=10.0)\\\text{beta} &\sim & \text{Normal}(\mathit{mu}=0.0,~\mathit{sigma}=10.0)\\\text{p} &\sim & \text{Deterministic}(\text{Constant},~\text{Constant},~\text{alpha},~\text{beta},~\text{Constant})\\\text{x} &\sim & \text{Bernoulli}(\mathit{p}=\text{p})
            \end{array}
            $$




```python
num_particles = 500
trace, smc = sample_smc(model=challenger_model, 
                           draws=num_particles,
                          threshold=0.8, 
                          progressbar=True)
```

    Sample initial stage: ...
    Stage:   0 Beta: 0.000 Steps:  25 Acce: 1.000
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     24%|██▍       | 121/500 [00:00<00:00, 1204.65it/s][A
     48%|████▊     | 239/500 [00:00<00:00, 1195.81it/s][A
     72%|███████▏  | 360/500 [00:00<00:00, 1199.98it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1192.56it/s][A
    Stage:   1 Beta: 0.011 Steps:  25 Acce: 0.549
    
    100%|██████████| 500/500 [00:00<00:00, 5654.91it/s]
    Stage:   2 Beta: 0.028 Steps:   5 Acce: 0.406
    
      0%|          | 0/500 [00:00<?, ?it/s][A
    100%|██████████| 500/500 [00:00<00:00, 3532.90it/s][A
    Stage:   3 Beta: 0.056 Steps:   8 Acce: 0.346
    
      0%|          | 0/500 [00:00<?, ?it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2808.64it/s][A
    Stage:   4 Beta: 0.106 Steps:  10 Acce: 0.301
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     48%|████▊     | 240/500 [00:00<00:00, 2394.73it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2361.99it/s][A
    Stage:   5 Beta: 0.182 Steps:  12 Acce: 0.275
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     41%|████      | 205/500 [00:00<00:00, 2048.77it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2041.58it/s][A
    Stage:   6 Beta: 0.313 Steps:  14 Acce: 0.266
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     42%|████▏     | 212/500 [00:00<00:00, 2117.50it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2076.29it/s][A
    Stage:   7 Beta: 0.529 Steps:  14 Acce: 0.250
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     39%|███▉      | 196/500 [00:00<00:00, 1959.82it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1943.86it/s][A
    Stage:   8 Beta: 0.897 Steps:  15 Acce: 0.270
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     42%|████▏     | 209/500 [00:00<00:00, 2085.22it/s][A
    100%|██████████| 500/500 [00:00<00:00, 2060.88it/s][A
    Stage:   9 Beta: 1.000 Steps:  14 Acce: 0.261
    
      0%|          | 0/500 [00:00<?, ?it/s][A
     38%|███▊      | 192/500 [00:00<00:00, 1912.39it/s][A
    100%|██████████| 500/500 [00:00<00:00, 1860.62it/s][A



```python
ppsamples = pm.sample_posterior_predictive(model=challenger_model, 
                                           trace=trace, 
                                           var_names=['p'])['p']
```

    
      0%|          | 0/500 [00:00<?, ?it/s][A
    100%|██████████| 500/500 [00:00<00:00, 3612.98it/s][A



```python
ppmean = ppsamples.mean(0)
pp_lower, pp_upper = np.percentile(ppsamples, axis=0, q=[2.5, 97.5])

plt.figure(figsize=(15, 8))
plt.plot(temp, data, 'ro', markersize=12, label='Observed data')
idx=np.argsort(temp)
plt.plot(temp[idx], ppmean[idx], linestyle='--', linewidth=2.5, 
         label='Post. pred. mean prob.')
plt.fill_between(temp[idx], pp_lower[idx], pp_upper[idx], 
                color='purple', alpha=0.25, label='95% Confidence')
plt.ylabel("Probability estimate",fontsize=20)
plt.xlabel("Outside temperature (Fahrenheit)",fontsize=20)
plt.title("Defects of the Space Shuttle O-Rings vs temperature",
          fontsize=20)
plt.yticks(np.arange(0., 1.01, 0.2), fontsize=20)
plt.xticks(fontsize=20)
plt.legend(loc='best', fontsize=20)
plt.tight_layout()
```


```python
logZ_temp = np.log(smc.model.marginal_likelihood)
```


```python

```

a6d773aefd28776fa3a851cb7f12a7b25ace083f
ipynb
Jupyter Notebook
lectures/lecture_24.ipynb
PredictiveScienceLab/uq-course
ddbe0865c9f91c4bd9b12e9b85d4293168306438
MIT
2016-01-04T15:31:44.000Z
2022-03-23T20:09:27.000Z
lectures/lecture_24.ipynb
ragusa/uq-course
ddbe0865c9f91c4bd9b12e9b85d4293168306438
MIT
2019-02-22T08:13:54.000Z
2020-02-08T19:25:16.000Z
lectures/lecture_24.ipynb
ragusa/uq-course
ddbe0865c9f91c4bd9b12e9b85d4293168306438
MIT
2016-01-05T18:50:34.000Z
2022-03-15T04:33:28.000Z
Qwen/Qwen-72B
1. YES
2. YES