diff --git a/model_selection/Model Selection.ipynb b/model_selection/Model Selection.ipynb
new file mode 100644
index 0000000..2747012
--- /dev/null
+++ b/model_selection/Model Selection.ipynb
@@ -0,0 +1,644 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Model Selection\n",
+ " This notebook shows how to implement model selection using ELFI for the models defined in the following paper: \"Semi-automatic selection of summary statistics for ABC model choice\" by Prangle et al, 2013. https://arxiv.org/abs/1302.5624"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "**NOTE:** It's required to have *ELFI* of version >= 0.7.1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import elfi\n",
+ "import scipy.stats as ss\n",
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "from elfi.methods.model_selection import compare_models"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Helper method for comparing different models\n",
+ "# Compare different models by using rejection sampling\n",
+ "def compare_rejection(distances, n_samples=100, threshold=1e5, batch_size=1):\n",
+ " samples = []\n",
+ " for dist in distances:\n",
+ " rej = elfi.Rejection(dist, batch_size=batch_size)\n",
+ " samples.append(rej.sample(n_samples, threshold=threshold))\n",
+ " \n",
+ " return compare_models(samples)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Case 1. 3 models:\n",
+ "### A1: Poisson($\\theta$), where $\\theta \\sim Exp(1)$ .\n",
+ "### A2: Geometric($\\theta$), where $\\theta \\sim Uniform(0, 1)$.\n",
+ "### A3: Binomial(10, $\\theta$), where $\\theta \\sim Beta(1, 9)$.\n",
+ "\n",
+ "As a summary statistics use $S_{10}(x) = \\{x^{(5)}, x^{(15)}, \\ldots, x^{(95)}\\}$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "elfi.new_model(name='Model 1')\n",
+ "\n",
+ "# Define priors for A1, A2, A3 respectively\n",
+ "a1_theta = elfi.Prior('expon', 1, name='a1_theta')\n",
+ "a2_theta = elfi.Prior('uniform', 0, 1, name='a2_theta')\n",
+ "a3_theta = elfi.Prior('beta', 1, 9, name='a3_theta')\n",
+ "\n",
+ "# Define models A1, A2, A3 respectively\n",
+ "def a1_model(theta, N=100, batch_size=1, random_state=None):\n",
+ " return ss.poisson.rvs(theta[:, None], size=(batch_size, N), random_state=random_state)\n",
+ "\n",
+ "def a2_model(theta, N=100, batch_size=1, random_state=None):\n",
+ " return ss.geom.rvs(theta[:, None], size=(batch_size, N), random_state=random_state)\n",
+ "\n",
+ "def a3_model(theta, N=100, batch_size=1, random_state=None):\n",
+ " return ss.binom.rvs(10, theta[:, None], size=(batch_size, N), random_state=random_state)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define summary statistics for all the models A1, A2, A3.\n",
+ "def a_summary(data):\n",
+ " percentiles = np.linspace(5, 95, num=10)\n",
+ " return np.percentile(data, percentiles, axis=1).T"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Let's generate observed data a_y0 from model A1. \n",
+ "# You can change it to any value you want and generate it from any model. \n",
+ "a_y0 = a1_model(np.array([1]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define simulators for each model.\n",
+ "a1_sim = elfi.Simulator(a1_model, a1_theta, observed=a_y0, name='a1_sim')\n",
+ "a2_sim = elfi.Simulator(a2_model, a2_theta, observed=a_y0, name='a2_sim')\n",
+ "a3_sim = elfi.Simulator(a3_model, a3_theta, observed=a_y0, name='a3_sim')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define summary statistics for each model.\n",
+ "a1_summary = elfi.Summary(a_summary, a1_sim, name='a1_summary')\n",
+ "a2_summary = elfi.Summary(a_summary, a2_sim, name='a2_summary')\n",
+ "a3_summary = elfi.Summary(a_summary, a3_sim, name='a3_summary')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define euclidean distances for each model.\n",
+ "a1_dist = elfi.Distance('euclidean', a1_summary, name='a1_dist')\n",
+ "a2_dist = elfi.Distance('euclidean', a2_summary, name='a2_dist')\n",
+ "a3_dist = elfi.Distance('euclidean', a3_summary, name='a3_dist')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# Draw graphs for models to make sure everything is correctly connected.\n",
+ "elfi.draw(a1_dist)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N_EXP = 30\n",
+ "# Perform experiment N_EXP times and plot the results\n",
+ "p1_probs = []\n",
+ "p2_probs = []\n",
+ "p3_probs = []\n",
+ "for n in range(N_EXP):\n",
+ " p1, p2, p3 = compare_rejection([a1_dist, a2_dist, a3_dist], batch_size=100)\n",
+ " p1_probs.append(p1)\n",
+ " p2_probs.append(p2)\n",
+ " p3_probs.append(p3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "{'boxes': [,\n",
+ " ,\n",
+ " ],\n",
+ " 'caps': [,\n",
+ " ,\n",
+ " ,\n",
+ " ,\n",
+ " ,\n",
+ " ],\n",
+ " 'fliers': [,\n",
+ " ,\n",
+ " ],\n",
+ " 'means': [],\n",
+ " 'medians': [,\n",
+ " ,\n",
+ " ],\n",
+ " 'whiskers': [,\n",
+ " ,\n",
+ " ,\n",
+ " ,\n",
+ " ,\n",
+ " ]}"
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADBtJREFUeJzt3X+o3fddx/HnaykRtU251xbZml8FIzPomHrpRJSNWaFFSYQJNkNYpZINDA6mfxSUJq0otMr+C7IghTnQbFaEqHH9Y3N/CG70VktnVrLG6ki6P5bt3llBbFZ8+0dvx9nl3pxzb86535x3ng8Ived7Ppz7br/NM4dP7vd7UlVIknp529ADSJKmz7hLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrotqG+8V133VUHDx4c6ttL0lx6/vnnv1lVd49bN1jcDx48yPLy8lDfXpLmUpKvTbLObRlJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0NdhHTzSLJ0CNsmZ97K2mcWz7uswplEiMsaTBuy0hSQ8Zdkhoy7pLUkHGXpIaMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8ZdkhqaKO5JHkhyMcmlJI9u8PzDSa4meWHt129Of1RJ0qRuG7cgyS7gNPCLwBXguSTnquor65Z+uqpOzGBGSdIWTfLO/T7gUlW9UlXXgLPA0dmOJUm6EZPE/R7g8sjjK2vH1vtAkheTPJNk31SmkyRty9htmQn9LfCXVfV6kg8DnwTev35RkuPAcYD9+/dP/OKLi4usrq5OadSdk2ToEbZkYWGBlZWVoceQNAWTxP1VYPSd+N61Y99VVd8aefhnwFMbvVBVnQHOACwtLdWkQ66urlI18XJt07z9YSRpc5NsyzwHHEpyb5LdwEPAudEFSd4+8vAI8NL0RpQkbdXYd+5V9UaSE8CzwC7g6aq6kOQJYLmqzgG/neQI8AawAjw8w5klSWNkqO2OpaWlWl5enmhtErdldoD/naWbX5Lnq2pp3DqvUJWkhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNTStW/7OVJ3cA6fuHHqM9urknqFHkDQlcxH3PP6a9zzZAUmoU0NPIWka3JaRpIaMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDxl2SGjLuktSQcZekhubiClV48+pJzdbCwsLQI0iakrmI+zzeeiDJXM4tqQe3ZSSpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ3Nx4zBJ/czjnV7n6WaAE71zT/JAkotJLiV59DrrPpCkkixNb0RJHVXVTH7N+rXnxdi4J9kFnAYeBA4Dx5Ic3mDdHcBHgS9Ne0hJ0tZM8s79PuBSVb1SVdeAs8DRDdb9AfAk8L9TnE+StA2TxP0e4PLI4ytrx74ryU8B+6rq76c4myRpm274p2WSvA34OPA7E6w9nmQ5yfLVq1dv9FtLkjYxSdxfBfaNPN67duwtdwA/DnwhyX8CPwOc2+gvVavqTFUtVdXS3Xffvf2pJUnXNUncnwMOJbk3yW7gIeDcW09W1X9V1V1VdbCqDgJfBI5U1fJMJpYkjTU27lX1BnACeBZ4CfhMVV1I8kSSI7MeUJK0dRNdxFRV54Hz6449tsna9934WDtnlhdSzOq15+3nbSXtvFv+ClVDKakj7y0jSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrIuEtSQ8Zdkhoy7pLUkHGXpIaMuyQ1ZNwlqSHjLkkN3fKfxCRpc4uLi6yurg49xpbN8uMzZ2FhYYGVlZWpvqZxl7Sp1dVVP4pyB8ziDyO3ZSSpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNeRFTJI2VSf3wKk7hx6jvTq5Z+qvadwlbSqPv+YVqjsgCXVquq/ptowkNWTcJakh4y5JDRl3SWrIuEtSQ8ZdkhqaKO5JHkhyMcmlJI9u8PxHknw5yQtJ/inJ4emPKkma1Ni4J9kFnAYeBA4DxzaI919U1U9U1buBp4CPT31SSdLEJnnnfh9wqapeqaprwFng6OiCqnpt5OEPAl71IEkDmuQK1XuAyyOPrwDvWb8oyW8BHwN2A++fynSSpG2Z2u0Hquo0cDrJB4HfBz60fk2S48BxgP3790/rW0uaoVl8eLO+18LCwtRfc5K4vwrsG3m8d+3YZs4Cf7rRE1V1BjgDsLS05NaNdJObx/vKJJnLuadtkj3354BDSe5Nsht4CDg3uiDJoZGHvwS8PL0RJUlbNfade1W9keQE8CywC3i6qi4keQJYrqpzwIkk9wPfAVbZYEtGkrRzJtpzr6rzwPl1xx4b+fqjU55LknQDvEJVkhoy7pLUkHGXpIaMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDU7txmLTT5vWGVt73RDvBuGtuzSqS3nhKHbgtI0kNGXdJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPeW0bSIGZ547dZvfY83XPIuEsaxDyFch65LSNJDRl3SWrIuEtSQ8ZdkhryL1Q1c4uLi6yurg49xpbM20f4LSwssLKyMvQYuokYd83c6uqqPxkxY/P2h5Fmz20ZSWrIuEtSQ8Zdkhoy7pLUkHGXpIaMuyQ1ZNwlqaGJ4p7kgSQXk1xK8ugGz38syVeSvJjkc0kOTH9USdKkxsY9yS7gNPAgcBg4luTwumX/CixV1buAZ4Cnpj2oJGlyk7xzvw+4VFWvVNU14CxwdHRBVf1jVf3P2sMvAnunO6YkaSsmuf3APcDlkcdXgPdcZ/0jwD9s9ESS48BxgP379084ouZdndwDp+4ceozW6uSeoUfQTWaq95ZJ8uvAEvDejZ6vqjPAGYClpSVvNnKLyOOveW+ZGUtCnRp6Ct1MJon7q8C+kcd71459jyT3A78HvLeqXp/OeJKk7Zhkz/054FCSe5PsBh4Czo0uSPKTwCeAI1X1jemPKUnairFxr6o3gBPAs8BLwGeq6kKSJ5IcWVv2x8DtwF8leSHJuU1eTpK0Aybac6+q88D5dcceG/n6/inPJUm6AV6hKkkNGXdJasiP2dOO8GPgZmthYWHoEXSTMe6auXn7GfckczeztJ7bMpLUkHGXpIaMuyQ1ZNwlqSHjLkkNGXdJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ0Zd0lqyLhLUkPGXZIaMu6S1JBxl6SGjLskNWTcJakh4y5JDRl3SWrotqEHkLYryVy+dlXN7LWltxh3zS0jKW3ObRlJasi4S1JDxl2SGjLuktSQcZekhoy7JDVk3CWpIeMuSQ1lqAtBklwFvjbIN98ZdwHfHHoIbYvnbr51P38HqurucYsGi3t3SZaramnoObR1nrv55vl7k9syktSQcZekhoz77JwZegBtm+duvnn+cM9dklrynbskNWTcpyDJrySpJO8cOfbZJN9O8ndDzqbx1p+/JO9O8s9JLiR5McmvDT2jNrfB+TuQ5F+SvLB2Dj8y9IxDcFtmCpJ8GngH8PmqOrl27BeAHwA+XFW/POR8ur715y/JjwJVVS8neQfwPPBjVfXtQQfVhjY4f7t5s22vJ7kd+DfgZ6vq64MOusN8536D1v7n+TngEeCht45X1eeA/x5qLk1mo/NXVV+tqpfXvv468A1g7EUj2nmbnL9rVfX62pLv4xbt3C35Lz1lR4HPVtVXgW8l+emhB9KWXPf8JbkP2A38+xDDaawNz1+SfUleBC4DT95q79rBuE/DMeDs2tdn1x5rfmx6/pK8HfgU8BtV9X8DzKbxNjx/VXW5qt4F/AjwoSQ/PNB8g3HP/QYkWQSuAFeBAnat/fNAVVWS9wG/6577zel65w+4A/gC8EdV9cxQM2pz437/jax7Gjh/q51H37nfmF8FPlVVB6rqYFXtA/4D+PmB59Jkrnf+/gb481stCHNm0/OX5PsBkizw5p78xQHnHMRtQw8w544BT6479tfAsSR/CLwTuD3JFeCRqnp2pwfUdW12/j4J7AV+KMnDa8cfrqoXdnA2jbfZ+fsE8J0kBQT4k6r68k4PNzS3ZSSpIbdlJKkh4y5JDRl3SWrIuEtSQ8Zdkhoy7pLUkHGXpIaMuyQ19P+onGdct82/iwAAAABJRU5ErkJggg==\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.boxplot((p1_probs, p2_probs, p3_probs), labels=['A1', 'A2', 'A3'], widths=0.8)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As we can see model $A_1$ is a bit more probable than model $A_3$, but much less improbable than $A_2$, when we're using Rejection Sampling."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Case 2. 2 models:\n",
+ "### B1: Laplace($\\theta$, $\\frac{1}{\\sqrt{2}}$), where $\\theta \\sim Normal(0, 2^2)$ .\n",
+ "### B2: Normal($\\theta$, $1$), where $\\theta \\sim Normal(0, 2^2)$.\n",
+ "\n",
+ "As a summary statistics use $S_{10}(x) = \\{x^{(5)}, x^{(15)}, \\ldots, x^{(95)}\\}$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "elfi.new_model(name='Model 2')\n",
+ "\n",
+ "# Define priors for B1 and B2 respectively\n",
+ "b1_theta = elfi.Prior('norm', 0, 2, name='b1_theta')\n",
+ "b2_theta = elfi.Prior('norm', 0, 2, name='b2_theta')\n",
+ "\n",
+ "# Define models A1, A2, A3 respectively\n",
+ "def b1_model(theta, N=100, batch_size=1, random_state=None):\n",
+ " return ss.laplace.rvs(theta[:, None], 1.0/np.sqrt(2), size=(batch_size, N), random_state=random_state)\n",
+ "\n",
+ "def b2_model(theta, N=100, batch_size=1, random_state=None):\n",
+ " return ss.norm.rvs(theta[:, None], 1, size=(batch_size, N), random_state=random_state)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define summary statistics for each model\n",
+ "def b_summary(data):\n",
+ " percentiles = np.linspace(5, 95, num=10)\n",
+ " return np.percentile(data, percentiles, axis=1).T"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Let's generate observed data y0 from model B1. \n",
+ "# You can change it to any value you want and generate it from any model. \n",
+ "b_y0 = b1_model(np.array([0]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define simulators for each model.\n",
+ "b1_sim = elfi.Simulator(b1_model, b1_theta, observed=b_y0, name='b1_sim')\n",
+ "b2_sim = elfi.Simulator(b2_model, b2_theta, observed=b_y0, name='b2_sim')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define summary statistics for each model.\n",
+ "b1_summary = elfi.Summary(b_summary, b1_sim, name='b1_summary')\n",
+ "b2_summary = elfi.Summary(b_summary, b2_sim, name='b2_summary')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Define euclidean distances for each model.\n",
+ "b1_dist = elfi.Distance('euclidean', b1_summary, name='b1_dist')\n",
+ "b2_dist = elfi.Distance('euclidean', b2_summary, name='b2_dist')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 19,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# Draw graphs for models to make sure everything is correctly connected.\n",
+ "elfi.draw(b1_dist)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Perform experiment N_EXP times and plot the results\n",
+ "p1_probs = []\n",
+ "p2_probs = []\n",
+ "for n in range(N_EXP):\n",
+ " p1, p2 = compare_rejection([b1_dist, b2_dist], batch_size=100)\n",
+ " p1_probs.append(p1)\n",
+ " p2_probs.append(p2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "{'boxes': [,\n",
+ " ],\n",
+ " 'caps': [,\n",
+ " ,\n",
+ " ,\n",
+ " ],\n",
+ " 'fliers': [,\n",
+ " ],\n",
+ " 'means': [],\n",
+ " 'medians': [,\n",
+ " ],\n",
+ " 'whiskers': [,\n",
+ " ,\n",
+ " ,\n",
+ " ]}"
+ ]
+ },
+ "execution_count": 20,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADdxJREFUeJzt3V+MXPV5h/Hni6mbP+DIYKvlj93lgiiNCEIwtZKqTSQqiCtV5gK1IalaVioyUmSRm6pymwuIuUrU5qq+sRAVpFIhshRpKVFcihQlqiDyGBlHxiU4VlOWRMJgE5rSxDF5e+FZZxjW3vHOeGfY3/ORVtpz5pzZd/HRs2fPMnNSVUiS2nDJpAeQJK0coy9JDTH6ktQQoy9JDTH6ktQQoy9JDTH6ktQQoy9JDTH6ktSQSyc9wKANGzbUzMzMpMeQpPeUAwcOvFZVG5fabuqiPzMzQ7fbnfQYkvSekuRHw2zn5R1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGTN2Ls1aTJJMeYap4P+bp4bH5Ti0dm0b/IpqGAynJVMyh6TItx4TH58rz8o4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDhop+kq1JXkxyNMnORR6fTXI8ycHexz0Dj69LMp/kH8c1uCTpwi35fvpJ1gC7gduAeWB/krmqemFg08erasc5nuZB4DsjTSpJGtkwZ/pbgKNVdayqTgGPAXcM+wWS3AL8FvBvyxtRkjQuw0T/GuDlvuX53rpBdyY5lGRvkk0ASS4B/gH465EnlSSNbFx/yH0CmKmqG4GngEd66z8PfLOq5s+3c5LtSbpJusePHx/TSJKkQcPcI/cVYFPf8rW9dWdV1et9iw8BX+l9/gngD5N8HrgMWJvkZ1W1c2D/PcAegE6n4w0zJekiGSb6+4Hrk1zHmdjfBXyuf4MkV1XVT3qL24AjAFX1533bzAKdweBLklbOktGvqtNJdgD7gDXAw1V1OMkuoFtVc8B9SbYBp4ETwOxFnFmStEypmq6rKZ1Op7rd7qTHWDWSMG3/xtICj8/xSXKgqjpLbecrciWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIcPcOes954orruDkyZOTHmNqJJn0CFNh/fr1nDhxYtJjSBO1KqN/8uRJb8ygd/GHn+TlHUlqitGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYMFf0kW5O8mORokp2LPD6b5HiSg72Pe3rrb0ryTJLDSQ4l+cy4vwFJ0vCWfD/9JGuA3cBtwDywP8lcVb0wsOnjVbVjYN1bwF9W1UtJrgYOJNlXVW+MY3hJ0oUZ5kx/C3C0qo5V1SngMeCOYZ68qn5QVS/1Pv8x8CqwcbnDSpJGM0z0rwFe7lue760bdGfvEs7eJJsGH0yyBVgL/HBZk0qSRjau2yU+AfxLVf0iyb3AI8CtCw8muQr4GnB3Vf1qcOck24HtAJs3bx7TSNJ08h7O7+RtLM9YqXs4DxP9V4D+M/dre+vOqqrX+xYfAr6ysJBkHfAk8MWqenaxL1BVe4A9AJ1Ox5vbalXzHs5azEr98Bvm8s5+4Pok1yVZC9wFzPVv0DuTX7ANONJbvxb4BvBoVe0dz8iSpOVa8ky/qk4n2QHsA9YAD1fV4SS7gG5VzQH3JdkGnAZOALO93f8M+CRwZZKFdbNVdXC834YkaRiZtl8zO51OdbvdkZ4jib8+612m5biYljk0XUY9LpIcqKrOUtv5ilxJaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGGH1JaojRl6SGjOt2iVOl7l8HD3xo0mNoytT96yY9gjRxqzL6+dKbvl+53iUJ9cCkp5Amy8s7ktQQoy9JDTH6ktQQoy9JDTH6ktQQoy9JDTH6ktQQoy9JDTH6ktQQoy9JDTH6ktQQoy9JDTH6ktQQoy9JDRkq+km2JnkxydEkOxd5fDbJ8SQHex/39D12d5KXeh93j3N4SdKFWfL99JOsAXYDtwHzwP4kc1X1wsCmj1fVjoF9rwDuBzpAAQd6+54cy/SSpAsyzJn+FuBoVR2rqlPAY8AdQz7/p4GnqupEL/RPAVuXN6okaVTDRP8a4OW+5fneukF3JjmUZG+STRe4ryRpBYzrD7lPADNVdSNnzuYfuZCdk2xP0k3SPX78+JhGkiQNGib6rwCb+pav7a07q6per6pf9BYfAm4Zdt/e/nuqqlNVnY0bNw47uyTpAg0T/f3A9UmuS7IWuAuY698gyVV9i9uAI73P9wG3J1mfZD1we2+dJGkClvy/d6rqdJIdnIn1GuDhqjqcZBfQrao54L4k24DTwAlgtrfviSQPcuYHB8CuqjpxEb4PSdIQUlWTnuEdOp1OdbvdkZ4jCdP2fWnypuW4mJY5NF1GPS6SHKiqzlLb+YpcSWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWrIkjdRkTRedf86eOBDkx5DU6buX7ciX8foSyssX3rTm6joXZJQD1z8r+PlHUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYYfUlqiNGXpIYMFf0kW5O8mORokp3n2e7OJJWk01v+jSSPJPl+kiNJ/nZcg0uSLtyS0U+yBtgN/DHwUeCzST66yHaXA18Avte3+k+B36yqjwG3APcmmRl9bEnScgxzpr8FOFpVx6rqFPAYcMci2z0IfBn4ed+6Aj6Y5FLg/cAp4M3RRpYkLdcw0b8GeLlveb637qwkNwObqurJgX33Av8L/AT4b+Dvq+rE8seVJI1i5NslJrkE+Cowu8jDW4C3gauB9cB3k/x7VR0beI7twHaAzZs3jzrSwnOO5Xm0eqxfv37SI0gTN0z0XwE29S1f21u34HLgBuDbvdD+NjCXZBvwOeBbVfVL4NUk/wF0gHdEv6r2AHsAOp3OyDcP9f6jv5bE/x6Szhrm8s5+4Pok1yVZC9wFzC08WFU/raoNVTVTVTPAs8C2qupy5pLOrQBJPgh8HPjPMX8PkqQhLRn9qjoN7AD2AUeAr1fV4SS7emfz57MbuCzJYc788Pinqjo06tCSpOXJtP3q3+l0qtvtTnqMVcPLO9PHfxMtZtTjIsmBquostZ2vyJWkhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWqI0Zekhhh9SWrIyLdLlHThvJ2nBq3U7TyNvrTCfC/9X/PeAivPyzuS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kNGSr6SbYmeTHJ0SQ7z7PdnUkqSadv3Y1JnklyOMn3k7xvHINLki7cku+nn2QNsBu4DZgH9ieZq6oXBra7HPgC8L2+dZcC/wz8RVU9n+RK4JdjnF+SdAGGOdPfAhytqmNVdQp4DLhjke0eBL4M/Lxv3e3Aoap6HqCqXq+qt0ecWZK0TMNE/xrg5b7l+d66s5LcDGyqqicH9v0wUEn2JXkuyd+MNK0kaSQj3y4xySXAV4HZczz/HwC/B7wFPJ3kQFU9PfAc24HtAJs3bx51JEnSOQxzpv8KsKlv+dreugWXAzcA307yX8DHgbneH3Pnge9U1WtV9RbwTeDmwS9QVXuqqlNVnY0bNy7vO5EkLWmY6O8Hrk9yXZK1wF3A3MKDVfXTqtpQVTNVNQM8C2yrqi6wD/hYkg/0/qj7KeCFd38JSdJKWDL6VXUa2MGZgB8Bvl5Vh5PsSrJtiX1PcubSz37gIPDcItf9JUkrJFU16RneodPpVLfbnfQYq0YSpu3fWFrg8Tk+vb+XdpbazlfkSlJDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDRr5Hrs4tyaRHAKZnDt83fXpMyzEB0zFLS8em0b+IWjqQ9N7isdkuL+9IUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1JNP2Io0kx4EfTXqOVWQD8Nqkh5DOweNzfH6nqjYutdHURV/jlaRbVZ1JzyEtxuNz5Xl5R5IaYvQlqSFGf/XbM+kBpPPw+FxhXtOXpIZ4pi9JDTH6q0iSt5McTPJ8kueS/H7fY99K8kaSf53kjGrTuY7NJDcleSbJ4SSHknxm0rOudl7eWUWS/KyqLut9/mng76rqU73lPwI+ANxbVX8ywTHVoHMdm0k+DFRVvZTkauAA8LtV9cYk513NPNNfvdYBJxcWqupp4H8mN4501tljs6p+UFUv9T7/MfAqsOQLjLR83i5xdXl/koPA+4CrgFsnPI+0YMljM8kWYC3wwxWerSlGf3X5v6q6CSDJJ4BHk9xQXsPT5J332ExyFfA14O6q+tUE51z1vLyzSlXVM5x5XxN/VdZUGTw2k6wDngS+WFXPTnK2Fhj9VSrJR4A1wOuTnkXq139sJlkLfAN4tKr2TnayNnh5Z3VZuG4KEM78qvw2QJLvAh8BLksyD/xVVe2b0Jxqz6LHZpLPAp8Erkwy23t8tqoOLvYkGp3/y6YkNcTLO5LUEKMvSQ0x+pLUEKMvSQ0x+pLUEKMvSQ0x+pLUEKMvSQ35f+D2s+bJwmudAAAAAElFTkSuQmCC\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "plt.boxplot((p1_probs, p2_probs), labels=['B1', 'B2'], widths=0.8)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As we can see model $B_1$ is a bit more probable than model $B_2$, when we're using Rejection Sampling."
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "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.5.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}