<jupyter_start><jupyter_text>Bayesian Model SelectionAuthors: <NAME> and <NAME> - 12/10/2018 Last modified by <NAME> on 12/13/2018  Quick review of Bayesian model selectionThe problem of model selection consists in determining which model(s) best explain the available data $D$, given a set of candidate models $m_{1:M}$. Each model $m_{j}$ is parameterized by a set of parameters $\theta_{m_{j}} \in \Theta_{m_{j}}$, to be estimated based on data. In the Bayesian framework, model selection is perfomed by computing the posterior probability of each model $m_{j}$ using Bayes' theorem:$$P(m_{j} \vert D) = \frac{p(D \vert m_{j})P(m_{j})}{\sum_{j=1}^{M} P(D \vert m_{j})P(m_{j})}$$where $P(m_{j})$ is the prior assigned to model $m_{j}$ and $P(D \vert m_{j})$ is the model evidence, also called marginal likelihood.$$ p(D \vert m_{j}) = \int_{\Theta_{m_{j}}} p(D \vert m_{j}, \theta_{m_{j}}) p(\theta_{m_{j}} \vert m_{j}) d\theta_{m_{j}} $$where $p(\theta_{m_{j}} \vert m_{j})$ is the prior assigned to the parameter vector of model $m_{j}$.  Numerical exampleIn the following we present an example for which the posterior pdf of the parameters, evidences and model probabilities can be computed analytically. We drop the $m_{j}$ subscript when referring to model parameters for simplicity. Three models are considered (the domain $x$ is fixed and consists in 50 equally spaced points):\begin{align*}m_{linear}:& \quad y = \theta_{0} x + \epsilon \\m_{quadratic}:& \quad y = \theta_{0} x + \theta_{1} x^2 + \epsilon \\m_{cubic}:& \quad y = \theta_{0} x + \theta_{1} x^2+ \theta_{2} x^3 + \epsilon \\\end{align*}All three models can be written in a compact form as $y=X \theta + \epsilon$, where $X$ contains the necessary powers of $x$. For all three models, the prior is chosen to be Gaussian, $p(\theta) = N(\cdot, \theta_{prior}, \Sigma_{prior})  $, and so is the noise $\epsilon \sim N(\cdot; 0, \sigma_{n}^{2} I)$. Then the posterior of the parameters can be computed analytically as:\begin{align*}& p(\theta \vert D={x,y}) =  N(\cdot; \theta_{post}(D), \Sigma_{post}(D)) \\& \theta_{post}(D) = \left( \frac{1}{\sigma_{n}^{2}}X^{T}X + \Sigma_{prior}^{-1} \right)^{-1} \left(\frac{1}{\sigma_{n}^{2}}X^{T}y+\Sigma^{-1}\theta_{prior} \right) \\& \Sigma_{post}(D) = \left( \frac{1}{\sigma_{n}^{2}}X^{T}X + \Sigma_{prior}^{-1} \right)^{-1}\end{align*}Then the evidence of each model can be computed as $$ p(D) = \frac{p(D \vert \theta)p(\theta)}{p(\theta \vert D)} $$where $p(D \vert \theta) = N(\cdot; X\theta, \sigma_{n}^{2} I)$, $p(\theta) = N(\cdot, \theta_{prior}, \Sigma_{prior}) $ and $p(\theta \vert D) = N(\cdot, \theta_{post}(D), \Sigma_{post}(D))$. This formula can be computed at any point $\theta$. Generate data from the quadratic model<jupyter_code># import necessary packages
import numpy as np
import matplotlib.pyplot as plt
from UQpy.Inference import *
from UQpy.RunModel import RunModel # required to run the quadratic model
from sklearn.neighbors import KernelDensity # for the plots
from statsmodels.nonparametric.kde import KDEUnivariate
# Generate data from a quadratic function
import random
random.seed(8) #set seed for reproducibility
from scipy.stats import multivariate_normal
param_true = np.array([1.0, 2.0]).reshape(1, -1)
var_n = 1
error_covariance = var_n*np.eye(50)
print(param_true.shape)

z = RunModel(samples=param_true, model_script='pfn_models.py', model_object_name = 'model_quadratic', 
             var_names = ['theta_1', 'theta_2'])
data = z.qoi_list[0].reshape((-1,))+multivariate_normal.rvs(mean=None, cov=error_covariance, size=1)
print(data)<jupyter_output>(1, 2)
[ 1.55804518e+00  1.75925199e+00 -1.31493917e-01  9.26333803e-01
  2.42640783e+00  2.82549399e+00  2.20965222e+00  3.54850769e+00
  7.77924139e+00  8.70023087e+00  8.39858833e+00  1.20067656e+01
  1.51569944e+01  1.75151961e+01  1.92183358e+01  2.06834695e+01
  2.45151831e+01  2.92357643e+01  2.91794989e+01  3.39301161e+01
  3.69463783e+01  4.08398383e+01  4.83721570e+01  4.84018700e+01
  5.50772500e+01  5.58566672e+01  6.22728026e+01  6.83013460e+01
  6.99202590e+01  7.58944543e+01  8.10103745e+01  8.52979067e+01
  9.09147579e+01  9.73838282e+01  1.05668377e+02  1.08080734e+02
  1.15259336e+02  1.22235401e+02  1.28430556e+02  1.35515694e+02
  1.41475789e+02  1.47840578e+02  1.55672572e+02  1.64662612e+02
  1.71575766e+02  1.77402599e+02  1.85046286e+02  1.94755524e+02
  2.01294753e+02  2.10878518e+02]<jupyter_text>Define the models, compute the true values of the evidence.For all three models, a Gaussian prior is chosen for the parameters, with mean and covariance matrix of the appropriate dimensions. Each model is given prior probability $P(m_{j}) = 1/3$.<jupyter_code># Define the models
model_names = ['model_linear', 'model_quadratic', 'model_cubic']
model_n_params = [1, 2, 3]
model_prior_means = [[0], [0,0], [0,0,0]]
model_prior_stds = [[10], [1,1], [1,2,0.25]]
evidences = []
model_posterior_means = []
model_posterior_stds = []
for n, model in enumerate(model_names):
    # compute matrix X
    X = np.linspace(0, 10, 50).reshape((-1,1))
    if n == 1: # quadratic model
        X = np.concatenate([X, X**2], axis=1)
    if n == 2: # cubic model
        X = np.concatenate([X, X**2, X**3], axis=1)

    # compute posterior pdf
    m_prior = np.array(model_prior_means[n]).reshape((-1,1))
    S_prior = np.diag(np.array(model_prior_stds[n])**2)
    S_posterior = np.linalg.inv(1/var_n*np.matmul(X.T,X)+np.linalg.inv(S_prior))
    m_posterior = np.matmul(S_posterior, 
                            1/var_n*np.matmul(X.T, data.reshape((-1,1)))+np.matmul(np.linalg.inv(S_prior),m_prior))
    m_prior = m_prior.reshape((-1,))
    m_posterior = m_posterior.reshape((-1,))
    model_posterior_means.append(list(m_posterior))
    model_posterior_stds.append(list(np.sqrt(np.diag(S_posterior))))
    print('posterior mean and covariance for '+model)
    print(m_posterior, S_posterior)
    
    # compute evidence, evaluate the formula at the posterior mean
    like_theta = multivariate_normal.pdf(data, mean=np.matmul(X,m_posterior).reshape((-1,)), cov=error_covariance)
    prior_theta = multivariate_normal.pdf(m_posterior, mean=m_prior, cov=S_prior)
    posterior_theta = multivariate_normal.pdf(m_posterior, mean=m_posterior, cov=S_posterior)
    evidence = like_theta*prior_theta/posterior_theta
    evidences.append(evidence)
    print('evidence for '+model+'= {}\n'.format(evidence))
    
# compute the posterior probability of each model
tmp = [1/3*evidence for evidence in evidences]
model_posterior_probas = [p/sum(tmp) for p in tmp]

print('posterior probabilities of all three models')
print(model_posterior_probas)<jupyter_output>posterior mean and covariance for model_linear
[16.18513084] [[0.00059394]]
evidence for model_linear= 0.0

posterior mean and covariance for model_quadratic
[0.97602516 2.00753812] [[ 0.00941698 -0.00116538]
 [-0.00116538  0.00015392]]
evidence for model_quadratic= 4.771706227861077e-40

posterior mean and covariance for model_cubic
[ 0.76481102  2.07732049 -0.00518351] [[ 5.60383740e-02 -1.66457308e-02  1.15371992e-03]
 [-1.66457308e-02  5.29411877e-03 -3.83090725e-04]
 [ 1.15371992e-03 -3.83090725e-04  2.85512405e-05]]
evidence for model_cubic= 3.701094005456616e-41

posterior probabilities of all three models
[0.0, 0.928019706015724, 0.07198029398427594]<jupyter_text>Define the models for use in UQpy<jupyter_code># Define models
candidate_models = []
for n, model_name in enumerate(model_names):
    model = Model(n_params=model_n_params[n], model_type='python', 
                  model_script='pfn_models.py', model_object_name = model_name,
                  prior_name = ['normal']*model_n_params[n], 
                  prior_params = [[m, std] for (m,std) in zip(model_prior_means[n], model_prior_stds[n])],
                  error_covariance=error_covariance,
                  model_name=model_name)
    candidate_models.append(model)<jupyter_output><empty_output><jupyter_text>Run MCMC for all three models<jupyter_code># Linear model
from scipy.stats import norm

bayesMCMC = BayesParameterEstimation(data=data, model=candidate_models[0], sampling_method = 'MCMC', nsamples=1200,
                                     algorithm = 'MH', jump=10, nburn=100, pdf_proposal_type = 'Normal',
                                     pdf_proposal_scale = [0.1], 
                                     seed = model_prior_means[0]
                                     )
# plot prior, true posterior and estimated posterior
fig, ax = plt.subplots(figsize=(8,5))
domain_plot = np.linspace(0,20,200)
ax.plot(domain_plot, norm.pdf(domain_plot, loc=model_prior_means[0], scale=model_prior_stds[0]),
            label = 'prior', color='green', linestyle='--')
ax.plot(domain_plot, norm.pdf(domain_plot, loc=model_posterior_means[0], 
                                  scale=model_posterior_stds[0]),
            label = 'true posterior', color='red', linestyle='-')
ax.hist(bayesMCMC.samples[:,0], density=True, bins=100, label='estimated posterior MCMC')
plt.legend()
plt.title('MCMC for linear model')
plt.show(fig)
# Quadratic system
bayesMCMC = BayesParameterEstimation(data=data, model=candidate_models[1], sampling_method = 'MCMC', nsamples=3500,
                                     algorithm = 'MH', jump=10, nburn=100, pdf_proposal_type = 'Normal',
                                     pdf_proposal_scale = [0.1, 0.1], 
                                     seed = model_prior_means[1]
                                     )
# plot prior, true posterior and estimated posterior
fig, ax = plt.subplots(1,2,figsize=(16,5))
for n_p in range(2):
    domain_plot = np.linspace(-0.5,3,200)
    ax[n_p].plot(domain_plot, norm.pdf(domain_plot, loc=model_prior_means[1][n_p], scale=model_prior_stds[1][n_p]),
                label = 'prior', color='green', linestyle='--')
    ax[n_p].plot(domain_plot, norm.pdf(domain_plot, loc=model_posterior_means[1][n_p], 
                                      scale=model_posterior_stds[1][n_p]),
                label = 'true posterior', color='red', linestyle='-')
    ax[n_p].hist(bayesMCMC.samples[:,n_p], density=True, bins=30, label='estimated posterior MCMC')
    ax[n_p].legend()
    ax[n_p].set_title('MCMC for quadratic model')
plt.show(fig)
# Cubic system
domain_plot = [np.linspace(-0.5,3,200), np.linspace(0,3,200), np.linspace(-0.5,0.5,200)]

bayesMCMC = BayesParameterEstimation(data=data, model=candidate_models[2], sampling_method = 'MCMC', nsamples=12000,
                                     algorithm = 'MH', jump=30, nburn=500, pdf_proposal_type = 'Normal',
                                     pdf_proposal_scale = [0.15, 0.1, 0.05], 
                                     seed = model_prior_means[2]
                                     )
# plot prior, true posterior and estimated posterior
fig, ax = plt.subplots(1,3,figsize=(24,5))
for n_p in range(3):
    ax[n_p].plot(domain_plot[n_p], norm.pdf(domain_plot[n_p], loc=model_prior_means[2][n_p], 
                                            scale=model_prior_stds[2][n_p]),
                label = 'prior', color='green', linestyle='--')
    ax[n_p].plot(domain_plot[n_p], norm.pdf(domain_plot[n_p], loc=model_posterior_means[2][n_p], 
                                      scale=model_posterior_stds[2][n_p]),
                label = 'true posterior', color='red', linestyle='-')
    ax[n_p].hist(bayesMCMC.samples[:,n_p], density=True, bins=30, label='estimated posterior MCMC')
    ax[n_p].legend()
    ax[n_p].set_title('MCMC for cubic model')
plt.show(fig)<jupyter_output><empty_output><jupyter_text>Run Bayesian Model Selection for all three models<jupyter_code># Defines constants for the MCMC learning part, same as above
algos = ['MH']*3
proposal_types=['Normal']*3
scales = [[0.1],[0.1, 0.1],[0.15, 0.1, 0.05]]
nsamples = [1200,1500,12000]
nburn = [100,100,500]
jump = [10,10,30]
selection = BayesModelSelection(data=data, candidate_models=candidate_models, pdf_proposal_type=proposal_types, 
                                pdf_proposal_scale=scales, algorithm=algos, jump=jump, nsamples=nsamples, nburn=nburn,  
                                prior_probabilities=[1./3., 1./3., 1./3.],
                                seed=model_prior_means, verbose=True, sorted_outputs=True)
print('Sorted models:')
print(selection.model_names)
print('Evidence of sorted models:')
print(selection.evidences)
print('Posterior probabilities of sorted models:')
print(selection.probabilities)<jupyter_output>Sorted models:
['model_quadratic', 'model_cubic', 'model_linear']
Evidence of sorted models:
[3.947272595341599e-36, 2.154161353073316e-36, 0.0]
Posterior probabilities of sorted models:
[0.6469417892112159, 0.3530582107887841, 0.0]
{'fraction_non_alphanumeric': 0.142643822697265, 'fraction_numerical': 0.0915592580949387, 'mean_word_length': 3.758788332086761, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '127669', 'license': 'MIT', 'n_tokens_mistral': 5303, 'n_tokens_neox': 4514, 'n_words': 1001, 'path': 'example/Bayesian/Illustrations of Inference classes/Bayesian model selection.ipynb', 'repo_name': 'volpatto/UQpy', 'size': '12745'}