# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Section 5.1 Model Comparison Methods
# *How can we determine which model better fits our needs?*

# +
import os
import time

import arviz as az
import numpy as np
import pymc3 as pm
import scipy.stats as stats


np.random.seed(0)

# Path Constants
if os.path.split(os.getcwd())[-1] != "notebooks":
    os.chdir(os.path.join(".."))

from utils import metropolis_hastings

NETCDF_DIR = "inference_data"
# -

az.style.use('arviz-white')

# ## Learning Objectives
# * Understanding of how to interpret WAIC, PSIS-LOO numerical metrics
# * Understanding of how to interpret plot_compare

# ## Infinite parameters and Infinite Models
# As Bayesian modelers not only do we have to handle infinite model parameters, we also have to handle infinite model definitions.
#
# Take our example estimating water on a planet's surface from Section 1.3:
# $$
# \theta \sim \operatorname{Uniform}(0,1) \\
# p_{water} \sim \operatorname{Binom}(\theta)
# $$
#
# In this model we're evaulating not just one possible proportion of water on a planet, but an infinite number of proportions from 0 to 1.
#
# But this begs the question why 0 to 1, or why even this model? Here is another valid model:
#
# $$
# \theta \sim \operatorname{Beta}(1,5) \\
# p_{water} \sim \operatorname{Binom}(\theta)
# $$
#
# And another:
#
# $$ 
# \theta \sim \operatorname{Uniform}(0,1) \\
# \sigma \sim \operatorname{Uniform}(0,100) \\
# p_{water} = Norm(\theta, \sigma)
# $$
#
# ### How do we pick?
# The flexibility to design whatever model you like is simultaneously wonderful and challenging. While this philosophy is nice, for the practioner there still is the question of which model to choose. The reason this question is particularly challenging in Bayesian statistics because we don't get just one prediction, we get distribution of predictions. However, this does not mean that Bayesian analysis is not suspect to phenomena such as overfitting or excess complexity.
#
# Luckily there are tools that help. In particular we'll cover *Widely Applicable Information Criterion* and how it's used in cojunction with the `plot_compare` function.

# ### A history of Information Theory
# When you hear Information Criterion you usually hear the words entropy, divergence, and deviance.  Let's briefly cover how these ideas build upon each other by walking through information theory.
#
#
# \*Note:* This section heavily draws from Richard McElreath's examples. While we'll be covering the material at a high level, to get a full understanding we recommend his  [lectures](https://www.youtube.com/watch?v=gjrsYDJbRh0) and [books](https://www.amazon.com/Statistical-Rethinking-Bayesian-Examples-Chapman/dp/1482253445/). 
#
# To restate our prior question "If have  3, or 5, or 10 models, which most usefully explains the truth, given the truth is also uncertain?"
#
# Over the years a number of folks have contributed mathematical tools that help answer this question:
#
# 1. How do we quantify uncertainty? (Information Entropy)
# 2. How do we quantify the additional uncertainty added from an estimate when we know the truth? (KL Divergence)
# 3. How do we estimate the additional uncertainty added from an estimate when we don't know the truth? (Deviance)
# 4. How do we make sure we're not underestimating uncertainty with complex models? (Information Criteria)
#
# This section is quite heavy of math. However it is not critical to understand the fine details, but more so the path of thinking.
#
# #### *Information Entropy*: Quantifying Uncertainty
#
# Weather is uncertain, and dice rolls are uncertain. But weather in certain areas is delightfully consistent, like my hometown of Orange County, and "loaded dice" are designed to be consistent.
#
# Information Entropy is a summarization of our uncertainty through the following formula  
#
# $$ H_p(\mathbf{x}) = -\sum_{j=1}^{N} p(x_j)\log(p(x_{j}))$$
#
# For some intuition, this is the minimum number of yes/no questions you need to determine $\mathbf{x} = \{x_1, x_2, \ldots, x_N\}$. 
#
# **Example 1.** if it is sunny every day in Orange County, and I have a collection of 365 observations of whether it is sunny in Orange County, you already know $\mathbf{p} = \{1, 1, \ldots, 1\}$. To confirm with the formula that $H(\mathbf{p}) = 0$, notice that $p(\text{sunny}) = 1$, so $\log(p(\text{sunny})) = 0$, and we just have a sum of $0 \cdot 1$, which is 0.
#
# **Example 2.** If we have a fair coin, and $N$ flips, then intuitively, I will have to go through and ask whether each flip was a head or a tail to determine the exact results, so we expect the entropy to be $N$. To confirm, notice that $p(\text{heads}) = p(\text{tails}) = 0.5$. By convention, we use base 2 logarithms, so 
# $$
# \log(p(\text{heads})) = \log(p(\text{tails})) = \log_2(0.5) = -1
# $$
#
# Then the entropy is $H(\mathbf{x}) = -\sum_{j=1}^N -1 = N$.
#
#
# #### *Divergence*: Quantifying additional uncertainty added from estimations of the truth
# Now that we can quantify the uncertainty in the distribution of reality, we can compare the uncertainty in the distribution our model to see how far off we are. The formula is as follows, where $p$ is the truth and $q$ is the estimate
#
# $$KL_{p, q}(\mathbf{x}) = -\sum_{j=1}^{N} p(x_{j})(\log(q(x_{j})) - \log(p(x_{j}))) = H_p(\mathbf{x}) - \sum p(x_j)\log(q(x_j))$$
#
# This is the idea of **divergence**, a measurement of the "distance" between two distributions. One important note this distance is not symmetric. For some quick intuition It was way more surprising for me to go from always sunny Los Angeles and experience snow in Wisconsin, than it is for a native Wisconsite to see sun in Los Angeles. People in Los Angeles only see sun, people in Wisconsin see sun and snow, therefore a Los Angeles is much farther removed from the reality of a mixed sunny/snowy climate, than a Wisconsinite is to an always sunny climate.
#
# #### *Deviance*: Estimating additional uncertainty because we don't know the truth using deviance
# If we knew the truth then there wouldn't be anything to do. Since we don't know the truth this unfortunately means we can't use the formula above, but through some clever thinking it turns out you don't need to know the truth to compare two models, just how far off one model is relative from another. By using true data (The weather in the past), and a model estimates, its possible to calculate **deviance** which is an estimate of the divergence.
#
# Taken from <NAME>'s book [Bayesian Data Analysis in Python](https://www.packtpub.com/big-data-and-business-intelligence/bayesian-analysis-python-second-edition)
#
# $$D_{kl}(p \mid \mid q) - D_{kl}(q \mid \mid r) \approx \sum_{n=1}^{N} log(q_{i}) - \sum_{n=1}^{N} log(r_{i})$$
# Remember this works because we're not trying to see how far one model is from the truth, we're just trying to see how far two models are from the truth *relative to each other*
#
# #### *Information Criterion*: Making sure we're not overfitting to observed data using complex models
# Deviance uses data that has been observed to make an estimation. As it turns out it's impossible to use data you haven't seen because you haven't seen it. As it also turns out if you make models complex you can explain anything perfectly. *Information Criterion* uses both estimates of deviance and a penalizing for excess model complexity to "score" a model fit.
#

# ## Information criterion in detail
# Information Criterion calcuations are composed of two terms
#
# 1. Estimate of a the deviance
# 2. Number of parameters in the model
#
# Below is the summarized formula for Widely Applicable Information Criterion (Also Watanabe–Akaike information criterion) 
#
# $$WAIC = -2 lppd + p_{WAIC}$$ 
#
# *lppd* stands for log pointwise predictive density, and $p_{WAIC}$ is an effective number of the parameters. Between the two we get an estimate of model "accuracy" balanced by model complexity
#
# ### What about AIC, BIC, DIC?
# There are numerous variations of Information Criterion (IC) estimates. However due to assumptions in the calculations some of the IC's estimations are less general than others. As the name suggests Widely Applicable IC, is well, the most widely applicable. <NAME> goes through the full explanation in his book.
#
#

# ## Compare
# Compare Dataframe is a handy way of comparing a set of models using Information Criteria.

data1 = az.load_arviz_data("non_centered_eight")
data2 = az.load_arviz_data("centered_eight")
compare_dict = {"non centered": data1, "centered": data2}
az.compare(compare_dict)

az.compare(compare_dict, scale="log")

# Here we load two versions 
# In order
# 1. WAIC is the WAIC 
# 2. p_waic is the estimate of number of parameters
# 3. Weight is the Akaike weight which can be used for model average
# 4. se is the standard error of the WAIC estimate (Remember WAIC estimates are also distributions)
# 5. dse is the standard error of the top ranked WAIC and each model
# 6. Warning that WAIC may be failing
# 7. Scale is the "multiplier" on the in sample deviation portion of WAIC

# ## Plot Compare
# ArviZ provides plot compare as a way to summarize the results of WAIC, and in particular when comparing multiple models. 

model_compare = az.compare({'Centered 8 schools': az.load_arviz_data('centered_eight'),
                           'Non-centered 8 schools': az.load_arviz_data('non_centered_eight')})
az.plot_compare(model_compare, plot_ic_diff=True, plot_standard_error=True, insample_dev=True)

# * The open circles are the WAIC values of each model. This includes both the in sample deviance and the model complexity penalization.
# * The dark circles are in sample deviance. These measure how well the model fits the data.
# * The standard error of the differences ofplotted in by the gray triangles and gray lines. This is the standard error of the different 
# * The dark gray lines are the standard error of WAIC
#
# All except WAIC can be turned on or off using the boolean argument shown above.

# ## Enough with the math, what's the practical advice
# * Models with lower Information Criterion tend to be better
# * Lowest Information Criterion does not always mean the best model
#
# That being said the interpretation of Information Criterion is not always straightforward. Interpretation depends on your data, how your models choice and parametrization, priors etc. Remember it's just a tool to help you understand your models, and not a strict rule.

# ## Additional Methods Pareto Smoothed Importance Sample Leave one out

# Leave One Out validation is another method to estimate a model's ability to fit unseen data.
#
# Leave One Out validation is a method where the model is trained on all the data points but one, then the posterior is used to estimate the likelihood of the last point. This method is an effective way of testing a models ability to fit unseeen data, but unfortunately requires the Inference run to be repeated N times for N data points, this is usually takes too long.
#
# PSIS-Loo, developed by <NAME>, does some magic to pareto sample stuff and make this easier.



StarcoderdataJupyterScriptsDedupFiltered