notebooks/6. Model Checking.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Checking\n",
    "\n",
    "After running an MCMC simulation, `sample` returns a `MutliTrace` object containing the samples for all the stochastic and deterministic random variables. The final step in Bayesian computation is model checking, in order to ensure that inferences derived from your sample are valid. There are two components to model checking:\n",
    "\n",
    "1. Convergence diagnostics\n",
    "2. Goodness of fit\n",
    "\n",
    "Convergence diagnostics are intended to detect lack of convergence in the Markov chain Monte Carlo sample; it is used to ensure that you have not halted your sampling too early. However, a converged model is not guaranteed to be a good model. The second component of model checking, goodness of fit, is used to check the internal validity of the model, by comparing predictions from the model to the data used to fit the model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convergence Diagnostics\n",
    "\n",
    "Valid inferences from sequences of MCMC samples are based on the\n",
    "assumption that the samples are derived from the true posterior\n",
    "distribution of interest. Theory guarantees this condition as the number\n",
    "of iterations approaches infinity. It is important, therefore, to\n",
    "determine the **minimum number of samples** required to ensure a reasonable\n",
    "approximation to the target posterior density. Unfortunately, no\n",
    "universal threshold exists across all problems, so convergence must be\n",
    "assessed independently each time MCMC estimation is performed. The\n",
    "procedures for verifying convergence are collectively known as\n",
    "*convergence diagnostics*.\n",
    "\n",
    "One approach to analyzing convergence is **analytical**, whereby the\n",
    "variance of the sample at different sections of the chain are compared\n",
    "to that of the limiting distribution. These methods use distance metrics\n",
    "to analyze convergence, or place theoretical bounds on the sample\n",
    "variance, and though they are promising, they are generally difficult to\n",
    "use and are not prominent in the MCMC literature. More common is a\n",
    "**statistical** approach to assessing convergence. With this approach,\n",
    "rather than considering the properties of the theoretical target\n",
    "distribution, only the statistical properties of the observed chain are\n",
    "analyzed. Reliance on the sample alone restricts such convergence\n",
    "criteria to **heuristics**. As a result, convergence cannot be guaranteed.\n",
    "Although evidence for lack of convergence using statistical convergence\n",
    "diagnostics will correctly imply lack of convergence in the chain, the\n",
    "absence of such evidence will not *guarantee* convergence in the chain.\n",
    "Nevertheless, negative results for one or more criteria may provide some\n",
    "measure of assurance to users that their sample will provide valid\n",
    "inferences.\n",
    "\n",
    "For most simple models, convergence will occur quickly, sometimes within\n",
    "a the first several hundred iterations, after which all remaining\n",
    "samples of the chain may be used to calculate posterior quantities. For\n",
    "more complex models, convergence requires a significantly longer burn-in\n",
    "period; sometimes orders of magnitude more samples are needed.\n",
    "Frequently, lack of convergence will be caused by **poor mixing**. \n",
    "Recall that *mixing* refers to the degree to which the Markov\n",
    "chain explores the support of the posterior distribution. Poor mixing\n",
    "may stem from inappropriate proposals (if one is using the\n",
    "Metropolis-Hastings sampler) or from attempting to estimate models with\n",
    "highly correlated variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import seaborn as sns; sns.set_context('notebook')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import exp, Normal, Binomial, sample, Model\n",
    "\n",
    "# Samples for each dose level\n",
    "n = 5 * np.ones(4, dtype=int)\n",
    "# Log-dose\n",
    "dose = np.array([-.86, -.3, -.05, .73])\n",
    "deaths = np.array([0, 1, 3, 5])\n",
    "\n",
    "def invlogit(x):\n",
    "    return exp(x) / (1 + exp(x))\n",
    "\n",
    "with Model() as bioassay_model:\n",
    "\n",
    "    # Logit-linear model parameters\n",
    "    alpha = Normal('alpha', 0, 0.01)\n",
    "    beta = Normal('beta', 0, 0.01)\n",
    "\n",
    "    # Calculate probabilities of death\n",
    "    theta = invlogit(alpha + beta * dose)\n",
    "\n",
    "    # Data likelihood\n",
    "    deaths = Binomial('deaths', n=n, p=theta, observed=deaths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import Metropolis\n",
    "\n",
    "with bioassay_model:\n",
    "    step = Metropolis(scaling=0.0001)\n",
    "    bioassay_trace = sample(1000, step=step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import traceplot\n",
    "\n",
    "traceplot(bioassay_trace[500:], varnames=['alpha'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Informal Methods\n",
    "\n",
    "The most straightforward approach for assessing convergence is based on\n",
    "simply **plotting and inspecting traces and histograms** of the observed\n",
    "MCMC sample. If the trace of values for each of the stochastics exhibits\n",
    "asymptotic behavior over the last $m$ iterations, this may be\n",
    "satisfactory evidence for convergence. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with bioassay_model:\n",
    "    bioassay_trace = sample(10000)\n",
    "    \n",
    "traceplot(bioassay_trace[9000:], varnames=['beta'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A similar approach involves\n",
    "plotting a histogram for every set of $k$ iterations (perhaps 50-100)\n",
    "beyond some burn in threshold $n$; if the histograms are not visibly\n",
    "different among the sample intervals, this may be considered some evidence for\n",
    "convergence. Note that such diagnostics should be carried out for each\n",
    "stochastic estimated by the MCMC algorithm, because convergent behavior\n",
    "by one variable does not imply evidence for convergence for other\n",
    "variables in the analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "beta_trace = bioassay_trace['beta']\n",
    "\n",
    "fig, axes = plt.subplots(2, 5, figsize=(14,6))\n",
    "axes = axes.ravel()\n",
    "for i in range(10):\n",
    "    axes[i].hist(beta_trace[500*i:500*(i+1)])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An extension of this approach can be taken\n",
    "when multiple parallel chains are run, rather than just a single, long\n",
    "chain. In this case, the final values of $c$ chains run for $n$\n",
    "iterations are plotted in a histogram; just as above, this is repeated\n",
    "every $k$ iterations thereafter, and the histograms of the endpoints are\n",
    "plotted again and compared to the previous histogram. This is repeated\n",
    "until consecutive histograms are indistinguishable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another *ad hoc* method for detecting lack of convergence is to examine\n",
    "the traces of several MCMC chains initialized with different starting\n",
    "values. Overlaying these traces on the same set of axes should (if\n",
    "convergence has occurred) show each chain tending toward the same\n",
    "equilibrium value, with approximately the same variance. Recall that the\n",
    "tendency for some Markov chains to converge to the true (unknown) value\n",
    "from diverse initial values is called *ergodicity*. This property is\n",
    "guaranteed by the reversible chains constructed using MCMC, and should\n",
    "be observable using this technique. Again, however, this approach is\n",
    "only a heuristic method, and cannot always detect lack of convergence,\n",
    "even though chains may appear ergodic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with bioassay_model:\n",
    "    \n",
    "    bioassay_trace = sample(1000, njobs=2, start=[{'alpha':0.5}, {'alpha':5}])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bioassay_trace.get_values('alpha', chains=0)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(bioassay_trace.get_values('alpha', chains=0)[:200], 'r--')\n",
    "plt.plot(bioassay_trace.get_values('alpha', chains=1)[:200], 'k--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A principal reason that evidence from informal techniques cannot\n",
    "guarantee convergence is a phenomenon called ***metastability***. Chains may\n",
    "appear to have converged to the true equilibrium value, displaying\n",
    "excellent qualities by any of the methods described above. However,\n",
    "after some period of stability around this value, the chain may suddenly\n",
    "move to another region of the parameter space. This period\n",
    "of metastability can sometimes be very long, and therefore escape\n",
    "detection by these convergence diagnostics. Unfortunately, there is no\n",
    "statistical technique available for detecting metastability.\n",
    "\n",
    "### Formal Methods\n",
    "\n",
    "Along with the *ad hoc* techniques described above, a number of more\n",
    "formal methods exist which are prevalent in the literature. These are\n",
    "considered more formal because they are based on existing statistical\n",
    "methods, such as time series analysis.\n",
    "\n",
    "PyMC currently includes three formal convergence diagnostic methods. The\n",
    "first, proposed by [Geweke (1992)](http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.ss/1177011446), is a time-series approach that\n",
    "compares the mean and variance of segments from the beginning and end of\n",
    "a single chain.\n",
    "\n",
    "$$z = \\frac{\\bar{\\theta}_a - \\bar{\\theta}_b}{\\sqrt{S_a(0) + S_b(0)}}$$\n",
    "\n",
    "where $a$ is the early interval and $b$ the late interval, and $S_i(0)$ is the spectral density estimate at zero frequency for chain segment $i$. If the\n",
    "z-scores (theoretically distributed as standard normal variates) of\n",
    "these two segments are similar, it can provide evidence for convergence.\n",
    "PyMC calculates z-scores of the difference between various initial\n",
    "segments along the chain, and the last 50% of the remaining chain. If\n",
    "the chain has converged, the majority of points should fall within 2\n",
    "standard deviations of zero.\n",
    "\n",
    "In PyMC, diagnostic z-scores can be obtained by calling the `geweke` function. It\n",
    "accepts either (1) a single trace, (2) a Node or Stochastic object, or\n",
    "(4) an entire Model object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import geweke\n",
    "\n",
    "with bioassay_model:\n",
    "    tr = sample(2000)\n",
    "    \n",
    "z = geweke(tr, intervals=15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(*z['alpha'].T)\n",
    "plt.hlines([-1,1], 0, 1000, linestyles='dotted')\n",
    "plt.xlim(0, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The arguments expected are the following:\n",
    "\n",
    "- `x` : The trace of a variable.\n",
    "- `first` : The fraction of series at the beginning of the trace.\n",
    "- `last` : The fraction of series at the end to be compared with the section at the beginning.\n",
    "- `intervals` : The number of segments.\n",
    "\n",
    "Plotting the output displays the scores in series, making it is easy to\n",
    "see departures from the standard normal assumption."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A second convergence diagnostic provided by PyMC is the Gelman-Rubin\n",
    "statistic [Gelman and Rubin (1992)](http://projecteuclid.org/DPubS?service=UI&version=1.0&verb=Display&handle=euclid.ss/1177011136). This diagnostic uses multiple chains to\n",
    "check for lack of convergence, and is based on the notion that if\n",
    "multiple chains have converged, by definition they should appear very\n",
    "similar to one another; if not, one or more of the chains has failed to\n",
    "converge.\n",
    "\n",
    "The Gelman-Rubin diagnostic uses an analysis of variance approach to\n",
    "assessing convergence. That is, it calculates both the between-chain\n",
    "varaince (B) and within-chain varaince (W), and assesses whether they\n",
    "are different enough to worry about convergence. Assuming $m$ chains,\n",
    "each of length $n$, quantities are calculated by:\n",
    "\n",
    "$$\\begin{align}B &= \\frac{n}{m-1} \\sum_{j=1}^m (\\bar{\\theta}_{.j} - \\bar{\\theta}_{..})^2 \\\\\n",
    "W &= \\frac{1}{m} \\sum_{j=1}^m \\left[ \\frac{1}{n-1} \\sum_{i=1}^n (\\theta_{ij} - \\bar{\\theta}_{.j})^2 \\right]\n",
    "\\end{align}$$\n",
    "\n",
    "for each scalar estimand $\\theta$. Using these values, an estimate of\n",
    "the marginal posterior variance of $\\theta$ can be calculated:\n",
    "\n",
    "$$\\hat{\\text{Var}}(\\theta | y) = \\frac{n-1}{n} W + \\frac{1}{n} B$$\n",
    "\n",
    "Assuming $\\theta$ was initialized to arbitrary starting points in each\n",
    "chain, this quantity will overestimate the true marginal posterior\n",
    "variance. At the same time, $W$ will tend to underestimate the\n",
    "within-chain variance early in the sampling run. However, in the limit\n",
    "as $n \\rightarrow \n",
    "\\infty$, both quantities will converge to the true variance of $\\theta$.\n",
    "In light of this, the Gelman-Rubin statistic monitors convergence using\n",
    "the ratio:\n",
    "\n",
    "$$\\hat{R} = \\sqrt{\\frac{\\hat{\\text{Var}}(\\theta | y)}{W}}$$\n",
    "\n",
    "This is called the potential scale reduction, since it is an estimate of\n",
    "the potential reduction in the scale of $\\theta$ as the number of\n",
    "simulations tends to infinity. In practice, we look for values of\n",
    "$\\hat{R}$ close to one (say, less than 1.1) to be confident that a\n",
    "particular estimand has converged. In PyMC, the function\n",
    "`gelman_rubin` will calculate $\\hat{R}$ for each stochastic node in\n",
    "the passed model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import gelman_rubin\n",
    "\n",
    "gelman_rubin(bioassay_trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the best results, each chain should be initialized to highly\n",
    "dispersed starting values for each stochastic node.\n",
    "\n",
    "By default, when calling the `forestplot` function using nodes with\n",
    "multiple chains, the $\\hat{R}$ values will be plotted alongside the\n",
    "posterior intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import forestplot\n",
    "\n",
    "forestplot(bioassay_trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Autocorrelation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import autocorrplot\n",
    "\n",
    "autocorrplot(tr);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bioassay_trace['alpha'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import effective_n\n",
    "\n",
    "effective_n(bioassay_trace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goodness of Fit\n",
    "\n",
    "Checking for model convergence is only the first step in the evaluation\n",
    "of MCMC model outputs. It is possible for an entirely unsuitable model\n",
    "to converge, so additional steps are needed to ensure that the estimated\n",
    "model adequately fits the data. One intuitive way of evaluating model\n",
    "fit is to compare model predictions with the observations used to fit\n",
    "the model. In other words, the fitted model can be used to simulate\n",
    "data, and the distribution of the simulated data should resemble the\n",
    "distribution of the actual data.\n",
    "\n",
    "Fortunately, simulating data from the model is a natural component of\n",
    "the Bayesian modelling framework. Recall, from the discussion on\n",
    "imputation of missing data, the posterior predictive distribution:\n",
    "\n",
    "$$p(\\tilde{y}|y) = \\int p(\\tilde{y}|\\theta) f(\\theta|y) d\\theta$$\n",
    "\n",
    "Here, $\\tilde{y}$ represents some hypothetical new data that would be\n",
    "expected, taking into account the posterior uncertainty in the model\n",
    "parameters. Sampling from the posterior predictive distribution is easy\n",
    "in PyMC. The code looks identical to the corresponding data stochastic,\n",
    "with two modifications: (1) the node should be specified as\n",
    "deterministic and (2) the statistical likelihoods should be replaced by\n",
    "random number generators. Consider the `gelman_bioassay` example, \n",
    "where deaths are modeled as a binomial random variable for which\n",
    "the probability of death is a logit-linear function of the dose of a\n",
    "particular drug."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import Normal, Binomial, Deterministic, invlogit\n",
    "\n",
    "# Samples for each dose level\n",
    "n = 5 * np.ones(4, dtype=int)\n",
    "# Log-dose\n",
    "dose = np.array([-.86, -.3, -.05, .73])\n",
    "\n",
    "with Model() as model:\n",
    "\n",
    "    # Logit-linear model parameters\n",
    "    alpha = Normal('alpha', 0, 0.01)\n",
    "    beta = Normal('beta', 0, 0.01)\n",
    "\n",
    "    # Calculate probabilities of death\n",
    "    theta = Deterministic('theta', invlogit(alpha + beta * dose))\n",
    "\n",
    "    # Data likelihood\n",
    "    deaths = Binomial('deaths', n=n, p=theta, observed=[0, 1, 3, 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior predictive distribution of deaths uses the same functional\n",
    "form as the data likelihood, in this case a binomial stochastic. Here is\n",
    "the corresponding sample from the posterior predictive distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    \n",
    "    deaths_sim = Binomial('deaths_sim', n=n, p=theta, shape=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the observed stochastic `Binomial` has been replaced with a stochastic node that is identical in every respect to `deaths`, except that its values are not fixed to be the observed data -- they are left to vary according to the values of the fitted parameters.\n",
    "\n",
    "The degree to which simulated data correspond to observations can be evaluated in at least two ways. First, these quantities can simply be compared visually. This allows for a qualitative comparison of model-based replicates and observations. If there is poor fit, the true value of the data may appear in the tails of the histogram of replicated data, while a good fit will tend to show the true data in high-probability regions of the posterior predictive distribution. The Matplot package in PyMC provides an easy way of producing such plots, via the `gof_plot` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    \n",
    "    gof_trace = sample(2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import forestplot\n",
    "\n",
    "forestplot(gof_trace, varnames=['deaths_sim'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Meta-analysis of beta blocker effectiveness\n",
    "\n",
    "Carlin (1992) considers a Bayesian approach to meta-analysis, and includes the following examples of 22 trials of beta-blockers to prevent mortality after myocardial infarction.\n",
    "\n",
    "In a random effects meta-analysis we assume the true effect (on a log-odds scale) $d_i$ in a trial $i$\n",
    "is drawn from some population distribution. Let $r^C_i$ denote number of events in the control group in trial $i$,\n",
    "and $r^T_i$ denote events under active treatment in trial $i$. Our model is:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "r^C_i &\\sim \\text{Binomial}\\left(p^C_i, n^C_i\\right) \\\\\n",
    "r^T_i &\\sim \\text{Binomial}\\left(p^T_i, n^T_i\\right) \\\\\n",
    "\\text{logit}\\left(p^C_i\\right) &= \\mu_i \\\\\n",
    "\\text{logit}\\left(p^T_i\\right) &= \\mu_i + \\delta_i \\\\\n",
    "\\delta_i &\\sim \\text{Normal}(d, t) \\\\\n",
    "\\mu_i &\\sim \\text{Normal}(m, s)\n",
    "\\end{aligned}$$\n",
    "\n",
    "We want to make inferences about the population effect $d$, and the predictive distribution for the effect $\\delta_{\\text{new}}$ in a new trial. Build a model to estimate these quantities in PyMC, and (1) use convergence diagnostics to check for convergence and (2) use posterior predictive checks to assess goodness-of-fit.\n",
    "\n",
    "Here are the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_t_obs = [3, 7, 5, 102, 28, 4, 98, 60, 25, 138, 64, 45, 9, 57, 25, 33, 28, 8, 6, 32, 27, 22]\n",
    "n_t_obs = [38, 114, 69, 1533, 355, 59, 945, 632, 278,1916, 873, 263, 291, 858, 154, 207, 251, 151, 174, 209, 391, 680]\n",
    "r_c_obs = [3, 14, 11, 127, 27, 6, 152, 48, 37, 188, 52, 47, 16, 45, 31, 38, 12, 6, 3, 40, 43, 39]\n",
    "n_c_obs = [39, 116, 93, 1520, 365, 52, 939, 471, 282, 1921, 583, 266, 293, 883, 147, 213, 122, 154, 134, 218, 364, 674]\n",
    "N = len(n_c_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your answer here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science. A Review Journal of the Institute of Mathematical Statistics, 457–472.\n",
    "\n",
    "Geweke, J., Berger, J. O., & Dawid, A. P. (1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In Bayesian Statistics 4.\n",
    "\n",
    "Brooks, S. P., Catchpole, E. A., & Morgan, B. J. T. (2000). Bayesian Animal Survival Estimation. Statistical Science. A Review Journal of the Institute of Mathematical Statistics, 15(4), 357–376. doi:10.1214/ss/1177010123\n",
    "\n",
    "Gelman, A., Meng, X., & Stern, H. (1996). Posterior predicitive assessment of model fitness via realized discrepencies with discussion. Statistica Sinica, 6, 733–807.\n",
    "\n",
    "Raftery, A., & Lewis, S. (1992). One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo. Statistical Science. A Review Journal of the Institute of Mathematical Statistics, 7, 493–497.\n",
    "\n",
    "[CrossValidated: How to use scikit-learn's cross validation functions on multi-label classifiers](http://stats.stackexchange.com/questions/65828/how-to-use-scikit-learns-cross-validation-functions-on-multi-label-classifiers)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}