{ "cells": [ { "cell_type": "markdown", "id": "c53eb8f5-80a7-4330-9f13-63220177fcc0", "metadata": { "tags": [] }, "source": [ "# Section 5.3 — Bayesian difference between means\n", "\n", "This notebook contains the code examples from [Section 5.3 Bayesian difference between means]() from the **No Bullshit Guide to Statistics**.\n", "\n", "See also [t-test.ipynb](./explorations/bambi/t-test.ipynb)\n" ] }, { "cell_type": "markdown", "id": "a2d8dda2-58a9-424e-9fb3-32ad6e8777d8", "metadata": { "tags": [] }, "source": [ "#### Notebook setup" ] }, { "cell_type": "code", "execution_count": 1, "id": "aea0223e-aee9-4875-a714-897b6646baaa", "metadata": {}, "outputs": [], "source": [ "# load Python modules\n", "import os\n", "import numpy as np\n", "import pandas as pd\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "efd86c5a-c9d2-4eab-b67d-a65e39b23ef2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Figures setup\n", "plt.clf() # needed otherwise `sns.set_theme` doesn\"t work\n", "from plot_helpers import RCPARAMS\n", "RCPARAMS.update({\"figure.figsize\": (5, 3)}) # good for screen\n", "# RCPARAMS.update({\"figure.figsize\": (5, 1.6)}) # good for print\n", "sns.set_theme(\n", " context=\"paper\",\n", " style=\"whitegrid\",\n", " palette=\"colorblind\",\n", " rc=RCPARAMS,\n", ")\n", "\n", "# High-resolution please\n", "%config InlineBackend.figure_format = \"retina\"\n", "\n", "# Where to store figures\n", "DESTDIR = \"figures/bayesian/dmeans\"" ] }, { "cell_type": "code", "execution_count": 3, "id": "df811a10-417d-4389-8bff-30e59b5f6aef", "metadata": {}, "outputs": [], "source": [ "# set random seed for repeatability\n", "np.random.seed(42)\n", "#######################################################" ] }, { "cell_type": "code", "execution_count": null, "id": "b473f755-7ed8-414e-819b-609cd9220b7d", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "7bdad75b-efdc-4a90-836a-62ce38d88dd2", "metadata": {}, "source": [ "## Example movies genres\n", "\n", "see https://www.andrewheiss.com/blog/2019/01/29/diff-means-half-dozen-ways/\n", "\n", "see also https://bookdown.org/content/3686/metric-predicted-variable-on-one-or-two-groups.html#two-groups\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "be3f6790-973b-4358-b21b-37538043d2f0", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
titleyearratinggenregenre_numeric
0Blowing Wild19535.6Action1
1No Way Back19955.2Action1
2New Jack City19916.1Action1
3Noigwon19834.2Action1
4Tarzan and the Jungle Boy19685.2Action1
\n", "
" ], "text/plain": [ " title year rating genre genre_numeric\n", "0 Blowing Wild 1953 5.6 Action 1\n", "1 No Way Back 1995 5.2 Action 1\n", "2 New Jack City 1991 6.1 Action 1\n", "3 Noigwon 1983 4.2 Action 1\n", "4 Tarzan and the Jungle Boy 1968 5.2 Action 1" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "movies = pd.read_csv(\"../datasets/movies.csv\")\n", "movies.head()" ] }, { "cell_type": "code", "execution_count": 5, "id": "a2d29469-678f-49a3-8347-43b06a7d968c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "genre\n", "Action 5.2845\n", "Comedy 5.9670\n", "Name: rating, dtype: float64" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "movies.groupby(\"genre\")[\"rating\"].mean()" ] }, { "cell_type": "code", "execution_count": 6, "id": "ccfc85e6-235a-41bb-b53e-1961225d2832", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "TtestResult(statistic=-4.47525173500199, pvalue=9.976981171112132e-06, df=398.0)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.stats import ttest_ind\n", "\n", "actions = movies[movies[\"genre\"]==\"Action\"][\"rating\"]\n", "comedies = movies[movies[\"genre\"]==\"Comedy\"][\"rating\"]\n", "\n", "ttest_ind(actions, comedies, equal_var=True)" ] }, { "cell_type": "code", "execution_count": 7, "id": "dc9e120d-2f46-42ff-bc1e-8cef602d29ea", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "TtestResult(statistic=-4.47525173500199, pvalue=9.978285839671782e-06, df=397.7995256063933)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ttest_ind(actions, comedies, equal_var=False)" ] }, { "cell_type": "code", "execution_count": 8, "id": "03d416c6-0fff-451d-a6ec-4313447db21d", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] } ], "source": [ "import bambi as bmb\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d7246a95-061c-4d12-a3d0-d8e5c0d9cef6", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 44, "id": "97ed3700-c5c1-4d28-b324-59fdee0b406e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Formula: rating ~ 1 + C(genre, levels=levels)\n", " Family: gaussian\n", " Link: mu = identity\n", " Observations: 400\n", " Priors: \n", " target = mu\n", " Common-level effects\n", " Intercept ~ Normal(mu: 0.0, sigma: 5.0)\n", " C(genre, levels=levels) ~ Normal(mu: 0.0, sigma: 1.0)\n", " \n", " Auxiliary parameters\n", " sigma ~ HalfStudentT(nu: 4.0, sigma: 1.559)\n" ] }, { "ename": "AttributeError", "evalue": "'NoneType' object has no attribute 'model'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[44], line 16\u001b[0m\n\u001b[1;32m 14\u001b[0m \u001b[38;5;66;03m# Get model description\u001b[39;00m\n\u001b[1;32m 15\u001b[0m \u001b[38;5;28mprint\u001b[39m(model_eq)\n\u001b[0;32m---> 16\u001b[0m \u001b[43mmodel_eq\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbackend\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mmodel\u001b[49m\n", "\u001b[0;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'model'" ] } ], "source": [ "# Model formula\n", "levels = [\"Comedy\", \"Action\"]\n", "formula = bmb.Formula(\"rating ~ 1 + C(genre, levels=levels)\")\n", "\n", "# Choose custom priors \n", "priors = {\n", " \"Intercept\": bmb.Prior(\"Normal\", mu=0, sigma=5),\n", " \"C(genre, levels=levels)\": bmb.Prior(\"Normal\", mu=0, sigma=1)\n", "}\n", "\n", "# Build model\n", "model_eq = bmb.Model(formula, priors=priors, data=movies)\n", "\n", "# Get model description\n", "print(model_eq)" ] }, { "cell_type": "code", "execution_count": 50, "id": "cc64f418-cebe-45fc-8a68-c9743a14b8cd", "metadata": {}, "outputs": [], "source": [ "# model_eq.build()\n", "# model_eq.graph()" ] }, { "cell_type": "code", "execution_count": 11, "id": "6e1a3f52-7ff6-4604-ae8c-3e7632067b05", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [rating_sigma, Intercept, C(genre, levels=levels)]\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [8000/8000 00:05<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 17 seconds.\n" ] } ], "source": [ "# Fit the model using 1000 on each chain\n", "idata_eq = model_eq.fit(draws=1000)" ] }, { "cell_type": "code", "execution_count": 12, "id": "ccbdcf94-ce9b-4672-8509-61dd04490f88", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
medianmadeti_3%eti_97%mcse_medianess_medianess_tailr_hat
Intercept5.9560.0725.7596.1610.0025919.2852712.01.0
C(genre, levels=levels)[Action]-0.6630.105-0.948-0.3770.0026063.0533142.01.0
rating_sigma1.5250.0371.4311.6330.0015276.6563283.01.0
\n", "
" ], "text/plain": [ " median mad eti_3% eti_97% mcse_median \\\n", "Intercept 5.956 0.072 5.759 6.161 0.002 \n", "C(genre, levels=levels)[Action] -0.663 0.105 -0.948 -0.377 0.002 \n", "rating_sigma 1.525 0.037 1.431 1.633 0.001 \n", "\n", " ess_median ess_tail r_hat \n", "Intercept 5919.285 2712.0 1.0 \n", "C(genre, levels=levels)[Action] 6063.053 3142.0 1.0 \n", "rating_sigma 5276.656 3283.0 1.0 " ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import arviz as az\n", "az.summary(idata_eq, stat_focus=\"median\")" ] }, { "cell_type": "markdown", "id": "b5b76140-3133-41a5-8ef9-3c997b69b602", "metadata": {}, "source": [ "### Regression, assuming unequal variances\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "ad794f56-e62e-455d-bbd0-ab0733cd6c5b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Formula: rating ~ 1 + C(genre,levels=levels)\n", " sigma ~ C(genre,levels=levels)\n", " Family: gaussian\n", " Link: mu = identity\n", " sigma = log\n", " Observations: 400\n", " Priors: \n", " target = mu\n", " Common-level effects\n", " Intercept ~ Normal(mu: 0.0, sigma: 5.0)\n", " C(genre, levels=levels) ~ Normal(mu: 0.0, sigma: 1.0)\n", " target = sigma\n", " Common-level effects\n", " sigma_Intercept ~ Normal(mu: 0.0, sigma: 1.0)\n", " sigma_C(genre, levels=levels) ~ Cauchy(alpha: 0.0, beta: 1.0)\n" ] } ], "source": [ "levels = [\"Comedy\", \"Action\"]\n", "formula_uneq = bmb.Formula(\"rating ~ 1 + C(genre,levels=levels)\", \"sigma ~ C(genre,levels=levels)\")\n", "\n", "priors = {\n", " \"Intercept\": bmb.Prior(\"Normal\", mu=0, sigma=5),\n", " \"C(genre, levels=levels)\": bmb.Prior(\"Normal\", mu=0, sigma=1),\n", " \"sigma\": {\"C(genre, levels=levels)\": bmb.Prior(\"Cauchy\", alpha=0, beta=1)},\n", "}\n", "\n", "# Build model\n", "model_uneq = bmb.Model(formula_uneq, priors=priors, data=movies)\n", "\n", "# Get model description\n", "print(model_uneq)\n", "# model_uneq.build()\n", "# model_uneq.graph()" ] }, { "cell_type": "code", "execution_count": 14, "id": "ce4c0ce1-10d1-41c3-84da-0d8e0d04e19f", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [Intercept, C(genre, levels=levels), sigma_Intercept, sigma_C(genre, levels=levels)]\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [8000/8000 00:06<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 17 seconds.\n" ] } ], "source": [ "idata_uneq = model_uneq.fit(draws=1000)" ] }, { "cell_type": "code", "execution_count": 15, "id": "f9261f58-3adc-47ad-8f76-376ec8581694", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
medianmadeti_3%eti_97%mcse_medianess_medianess_tailr_hat
Intercept5.9580.0755.7536.1510.0025501.0803402.01.0
C(genre, levels=levels)[Action]-0.6670.101-0.940-0.3840.0036121.2313284.01.0
sigma_Intercept0.4340.0340.3400.5330.0016062.6343350.01.0
sigma_C(genre, levels=levels)[Action]-0.0230.047-0.1540.1130.0015655.4603044.01.0
sigma[0]1.5070.0501.3771.6600.0016215.9972575.01.0
...........................
sigma[395]1.5440.0531.4061.7030.0016062.6343350.01.0
sigma[396]1.5440.0531.4061.7030.0016062.6343350.01.0
sigma[397]1.5440.0531.4061.7030.0016062.6343350.01.0
sigma[398]1.5440.0531.4061.7030.0016062.6343350.01.0
sigma[399]1.5440.0531.4061.7030.0016062.6343350.01.0
\n", "

404 rows × 8 columns

\n", "
" ], "text/plain": [ " median mad eti_3% eti_97% \\\n", "Intercept 5.958 0.075 5.753 6.151 \n", "C(genre, levels=levels)[Action] -0.667 0.101 -0.940 -0.384 \n", "sigma_Intercept 0.434 0.034 0.340 0.533 \n", "sigma_C(genre, levels=levels)[Action] -0.023 0.047 -0.154 0.113 \n", "sigma[0] 1.507 0.050 1.377 1.660 \n", "... ... ... ... ... \n", "sigma[395] 1.544 0.053 1.406 1.703 \n", "sigma[396] 1.544 0.053 1.406 1.703 \n", "sigma[397] 1.544 0.053 1.406 1.703 \n", "sigma[398] 1.544 0.053 1.406 1.703 \n", "sigma[399] 1.544 0.053 1.406 1.703 \n", "\n", " mcse_median ess_median ess_tail \\\n", "Intercept 0.002 5501.080 3402.0 \n", "C(genre, levels=levels)[Action] 0.003 6121.231 3284.0 \n", "sigma_Intercept 0.001 6062.634 3350.0 \n", "sigma_C(genre, levels=levels)[Action] 0.001 5655.460 3044.0 \n", "sigma[0] 0.001 6215.997 2575.0 \n", "... ... ... ... \n", "sigma[395] 0.001 6062.634 3350.0 \n", "sigma[396] 0.001 6062.634 3350.0 \n", "sigma[397] 0.001 6062.634 3350.0 \n", "sigma[398] 0.001 6062.634 3350.0 \n", "sigma[399] 0.001 6062.634 3350.0 \n", "\n", " r_hat \n", "Intercept 1.0 \n", "C(genre, levels=levels)[Action] 1.0 \n", "sigma_Intercept 1.0 \n", "sigma_C(genre, levels=levels)[Action] 1.0 \n", "sigma[0] 1.0 \n", "... ... \n", "sigma[395] 1.0 \n", "sigma[396] 1.0 \n", "sigma[397] 1.0 \n", "sigma[398] 1.0 \n", "sigma[399] 1.0 \n", "\n", "[404 rows x 8 columns]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(idata_uneq, stat_focus=\"median\")" ] }, { "cell_type": "code", "execution_count": 16, "id": "9b882f87-0a06-473c-9b1f-9bb3f4a98c4e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.9772624837732771" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.exp(az.summary(idata_uneq, stat_focus=\"median\").loc[\"sigma_C(genre, levels=levels)[Action]\",\"median\"])" ] }, { "cell_type": "code", "execution_count": null, "id": "f7e73b05-6479-472c-95f7-4b58a44864f1", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "1a0b138c-5cd3-4733-971f-0c64a146bb12", "metadata": {}, "source": [ "### BEST" ] }, { "cell_type": "code", "execution_count": 17, "id": "0b9d4894-d31b-4bde-870f-0f4641d8a0a6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Formula: rating ~ 1 + C(genre,levels=levels)\n", " sigma ~ C(genre,levels=levels)\n", " Family: t\n", " Link: mu = identity\n", " sigma = log\n", " Observations: 400\n", " Priors: \n", " target = mu\n", " Common-level effects\n", " Intercept ~ Normal(mu: 0.0, sigma: 5.0)\n", " C(genre, levels=levels) ~ Normal(mu: 0.0, sigma: 1.0)\n", " \n", " Auxiliary parameters\n", " nu ~ Exponential(lam: 0.0345)\n", " target = sigma\n", " Common-level effects\n", " sigma_Intercept ~ Normal(mu: 0.0, sigma: 1.0)\n", " sigma_C(genre, levels=levels) ~ Cauchy(alpha: 0.0, beta: 1.0)\n" ] } ], "source": [ "levels = [\"Comedy\", \"Action\"]\n", "formula_best = bmb.Formula(\"rating ~ 1 + C(genre,levels=levels)\", \"sigma ~ C(genre,levels=levels)\")\n", "\n", "priors = {\n", " \"Intercept\": bmb.Prior(\"Normal\", mu=0, sigma=5),\n", " \"C(genre, levels=levels)\": bmb.Prior(\"Normal\", mu=0, sigma=1),\n", " \"sigma\": {\"C(genre, levels=levels)\": bmb.Prior(\"Cauchy\", alpha=0, beta=1)},\n", " \"nu\": bmb.Prior(\"Exponential\", lam=1/29),\n", "}\n", "\n", "# Build model\n", "model_best = bmb.Model(formula_best, priors=priors, family=\"t\", data=movies)\n", "\n", "# Get model description\n", "print(model_best)\n", "# model_best.build()\n", "# model_best.graph()" ] }, { "cell_type": "code", "execution_count": 18, "id": "d953cca3-3f52-4af0-b6ff-d660b36b7752", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [rating_nu, Intercept, C(genre, levels=levels), sigma_Intercept, sigma_C(genre, levels=levels)]\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [8000/8000 00:10<00:00 Sampling 4 chains, 4 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 22 seconds.\n", "There were 4 divergences after tuning. Increase `target_accept` or reparameterize.\n" ] } ], "source": [ "idata_best = model_best.fit(draws=1000)" ] }, { "cell_type": "code", "execution_count": 19, "id": "b94ada55-4676-4bcd-9f14-590c1c14d041", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
medianmadeti_3%eti_97%mcse_medianess_medianess_tailr_hat
Intercept5.9810.0705.7876.1720.0024154.8202915.01.0
C(genre, levels=levels)[Action]-0.6800.100-0.960-0.4020.0033883.9132685.01.0
sigma_Intercept0.3870.0400.2680.4950.0023197.4442080.01.0
sigma_C(genre, levels=levels)[Action]0.0000.051-0.1420.1440.0023651.1212894.01.0
rating_nu32.07913.9559.874105.4160.4553032.6442062.01.0
...........................
sigma[395]1.4720.0591.3071.6400.0023197.4442080.01.0
sigma[396]1.4720.0591.3071.6400.0023197.4442080.01.0
sigma[397]1.4720.0591.3071.6400.0023197.4442080.01.0
sigma[398]1.4720.0591.3071.6400.0023197.4442080.01.0
sigma[399]1.4720.0591.3071.6400.0023197.4442080.01.0
\n", "

405 rows × 8 columns

\n", "
" ], "text/plain": [ " median mad eti_3% eti_97% \\\n", "Intercept 5.981 0.070 5.787 6.172 \n", "C(genre, levels=levels)[Action] -0.680 0.100 -0.960 -0.402 \n", "sigma_Intercept 0.387 0.040 0.268 0.495 \n", "sigma_C(genre, levels=levels)[Action] 0.000 0.051 -0.142 0.144 \n", "rating_nu 32.079 13.955 9.874 105.416 \n", "... ... ... ... ... \n", "sigma[395] 1.472 0.059 1.307 1.640 \n", "sigma[396] 1.472 0.059 1.307 1.640 \n", "sigma[397] 1.472 0.059 1.307 1.640 \n", "sigma[398] 1.472 0.059 1.307 1.640 \n", "sigma[399] 1.472 0.059 1.307 1.640 \n", "\n", " mcse_median ess_median ess_tail \\\n", "Intercept 0.002 4154.820 2915.0 \n", "C(genre, levels=levels)[Action] 0.003 3883.913 2685.0 \n", "sigma_Intercept 0.002 3197.444 2080.0 \n", "sigma_C(genre, levels=levels)[Action] 0.002 3651.121 2894.0 \n", "rating_nu 0.455 3032.644 2062.0 \n", "... ... ... ... \n", "sigma[395] 0.002 3197.444 2080.0 \n", "sigma[396] 0.002 3197.444 2080.0 \n", "sigma[397] 0.002 3197.444 2080.0 \n", "sigma[398] 0.002 3197.444 2080.0 \n", "sigma[399] 0.002 3197.444 2080.0 \n", "\n", " r_hat \n", "Intercept 1.0 \n", "C(genre, levels=levels)[Action] 1.0 \n", "sigma_Intercept 1.0 \n", "sigma_C(genre, levels=levels)[Action] 1.0 \n", "rating_nu 1.0 \n", "... ... \n", "sigma[395] 1.0 \n", "sigma[396] 1.0 \n", "sigma[397] 1.0 \n", "sigma[398] 1.0 \n", "sigma[399] 1.0 \n", "\n", "[405 rows x 8 columns]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(idata_best, stat_focus=\"median\")" ] }, { "cell_type": "markdown", "id": "10f4dbfe-ce60-4055-8224-fc0621c8a1e6", "metadata": {}, "source": [ "### BEST with priors on variables instead of difference\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 42, "id": "6e76947c-5f64-45af-9526-7d6c78456845", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Formula: rating ~ 0 + C(genre,levels=levels)\n", " sigma ~ 0 + C(genre,levels=levels)\n", " Family: t\n", " Link: mu = identity\n", " sigma = log\n", " Observations: 400\n", " Priors: \n", " target = mu\n", " Common-level effects\n", " C(genre, levels=levels) ~ TruncatedNormal(mu: 6.0, sigma: 2.0, lower: 1.0, upper: 10.0)\n", " \n", " Auxiliary parameters\n", " nu ~ Exponential(lam: 0.0345)\n", " target = sigma\n", " Common-level effects\n", " sigma_C(genre, levels=levels) ~ Cauchy(alpha: 0.0, beta: 1.0)\n" ] } ], "source": [ "levels = [\"Comedy\", \"Action\"]\n", "formula_best2 = bmb.Formula(\"rating ~ 0 + C(genre,levels=levels)\",\n", " \"sigma ~ 0 + C(genre,levels=levels)\")\n", "\n", "priors = {\n", " \"C(genre, levels=levels)\": bmb.Prior(\"TruncatedNormal\", mu=6, sigma=2, lower=1, upper=10),\n", " \"sigma\": {\"C(genre, levels=levels)\": bmb.Prior(\"Cauchy\", alpha=0, beta=1)},\n", " \"nu\": bmb.Prior(\"Exponential\", lam=1/29),\n", "}\n", "\n", "# Build model\n", "model_best2 = bmb.Model(formula_best2, priors=priors, family=\"t\", data=movies)\n", "\n", "# Get model description\n", "print(model_best2)\n", "# model_best2.build()\n", "# model_best2.graph()" ] }, { "cell_type": "code", "execution_count": 21, "id": "7f04782d-6a8b-4c25-8f75-0091a92b8af7", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [rating_nu, C(genre, levels=levels), sigma_C(genre, levels=levels)]\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [8000/8000 00:10<00:00 Sampling 4 chains, 7 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 22 seconds.\n", "There were 7 divergences after tuning. Increase `target_accept` or reparameterize.\n" ] } ], "source": [ "idata_best2 = model_best2.fit(draws=1000)" ] }, { "cell_type": "code", "execution_count": 22, "id": "ef71e76d-2176-4b6f-a467-d87c49039806", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
medianmadeti_3%eti_97%mcse_medianess_medianess_tailr_hat
sigma_C(genre, levels=levels)[Comedy]0.3850.0400.2620.4970.0013433.9092076.01.0
sigma_C(genre, levels=levels)[Action]0.3840.0370.2800.4890.0013515.3782336.01.0
rating_nu30.34513.4209.972105.7260.4202940.6612236.01.0
C(genre, levels=levels)[Comedy]5.9900.0705.7766.1900.0023812.9292820.01.0
C(genre, levels=levels)[Action]5.2940.0745.0905.4930.0023737.3402513.01.0
...........................
sigma[395]1.4690.0591.3001.6440.0023433.9092076.01.0
sigma[396]1.4690.0591.3001.6440.0023433.9092076.01.0
sigma[397]1.4690.0591.3001.6440.0023433.9092076.01.0
sigma[398]1.4690.0591.3001.6440.0023433.9092076.01.0
sigma[399]1.4690.0591.3001.6440.0023433.9092076.01.0
\n", "

405 rows × 8 columns

\n", "
" ], "text/plain": [ " median mad eti_3% eti_97% \\\n", "sigma_C(genre, levels=levels)[Comedy] 0.385 0.040 0.262 0.497 \n", "sigma_C(genre, levels=levels)[Action] 0.384 0.037 0.280 0.489 \n", "rating_nu 30.345 13.420 9.972 105.726 \n", "C(genre, levels=levels)[Comedy] 5.990 0.070 5.776 6.190 \n", "C(genre, levels=levels)[Action] 5.294 0.074 5.090 5.493 \n", "... ... ... ... ... \n", "sigma[395] 1.469 0.059 1.300 1.644 \n", "sigma[396] 1.469 0.059 1.300 1.644 \n", "sigma[397] 1.469 0.059 1.300 1.644 \n", "sigma[398] 1.469 0.059 1.300 1.644 \n", "sigma[399] 1.469 0.059 1.300 1.644 \n", "\n", " mcse_median ess_median ess_tail \\\n", "sigma_C(genre, levels=levels)[Comedy] 0.001 3433.909 2076.0 \n", "sigma_C(genre, levels=levels)[Action] 0.001 3515.378 2336.0 \n", "rating_nu 0.420 2940.661 2236.0 \n", "C(genre, levels=levels)[Comedy] 0.002 3812.929 2820.0 \n", "C(genre, levels=levels)[Action] 0.002 3737.340 2513.0 \n", "... ... ... ... \n", "sigma[395] 0.002 3433.909 2076.0 \n", "sigma[396] 0.002 3433.909 2076.0 \n", "sigma[397] 0.002 3433.909 2076.0 \n", "sigma[398] 0.002 3433.909 2076.0 \n", "sigma[399] 0.002 3433.909 2076.0 \n", "\n", " r_hat \n", "sigma_C(genre, levels=levels)[Comedy] 1.0 \n", "sigma_C(genre, levels=levels)[Action] 1.0 \n", "rating_nu 1.0 \n", "C(genre, levels=levels)[Comedy] 1.0 \n", "C(genre, levels=levels)[Action] 1.0 \n", "... ... \n", "sigma[395] 1.0 \n", "sigma[396] 1.0 \n", "sigma[397] 1.0 \n", "sigma[398] 1.0 \n", "sigma[399] 1.0 \n", "\n", "[405 rows x 8 columns]" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(idata_best2, stat_focus=\"median\")" ] }, { "cell_type": "code", "execution_count": 23, "id": "603d0d96-f285-49e3-8d0b-44e4a2aee249", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.4681454416819895" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.exp(az.summary(idata_best2, stat_focus=\"median\").loc[\"sigma_C(genre, levels=levels)[Action]\",\"median\"])" ] }, { "cell_type": "code", "execution_count": null, "id": "b0a0cd8f-7ee8-4973-b3a6-12d10710ba99", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "520ee7a4-c570-4f85-8197-0f7467c00c18", "metadata": {}, "source": [ "## Example from original BEST paper\n", "\n", "Data taken from `BESTexample-original.R` in `BEST.zip`\n", "via https://web.archive.org/web/20170708173718/https://www.indiana.edu/~kruschke/BEST/\n", "\n", "\n", "Steps following Matti Vuorre's [blog post](https://mvuorre.github.io/posts/2017-01-02-how-to-compare-two-groups-with-robust-bayesian-estimation-using-r-stan-and-brms/) see also [src](https://github.com/mvuorre/mvuorre.github.io/tree/main/posts/2017-01-02-how-to-compare-two-groups-with-robust-bayesian-estimation-using-r-stan-and-brms) notebook." ] }, { "cell_type": "markdown", "id": "4fbc66ea-6f70-45f1-ac1b-1643a2d97744", "metadata": {}, "source": [ "### Data" ] }, { "cell_type": "code", "execution_count": 24, "id": "c4bba0e0-c5ab-4421-8030-cc018aec7ca2", "metadata": {}, "outputs": [], "source": [ "treated = np.array([101, 100, 102, 104, 102, 97, 105, 105, 98, 101, 100, 123, 105,\n", " 103, 100, 95, 102, 106, 109, 102, 82, 102, 100, 102, 102, 101,\n", " 102, 102, 103, 103, 97, 97, 103, 101, 97, 104, 96, 103, 124,\n", " 101, 101, 100, 101, 101, 104, 100, 101])\n", "controls = np.array([ 99, 101, 100, 101, 102, 100, 97, 101, 104, 101, 102, 102, 100,\n", " 105, 88, 101, 100, 104, 100, 100, 100, 101, 102, 103, 97, 101,\n", " 101, 100, 101, 99, 101, 100, 100, 101, 100, 99, 101, 100, 102,\n", " 99, 100, 99])\n", "d = pd.DataFrame({\n", " \"group\": [\"treatment\"]*len(treated) + [\"control\"]*len(controls), \n", " \"iq\": np.concatenate([treated, controls])\n", "})" ] }, { "cell_type": "code", "execution_count": 25, "id": "32653699-1c20-4fc0-a535-bbcae3d43cef", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4UAAAJLCAYAAABQRTPQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAB7CAAAewgFu0HU+AABdEklEQVR4nO3dd3wUdf7H8fdusumNAEkQKVJCUwRPBcWCHCfFcoqAigqIohRRUTxsp2LBwiGniIoVVPQQBRuioiACivwERQ5CL2IgCZCEkL7Znd8fIXsJqZtkSzKv5+Phg7jz/e58dnd2dt873/mOxTAMQwAAAAAAU7L6ugAAAAAAgO8QCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAJ9XQCq9+uvv8rpdMpisSgwkJcMAAAAMKuioiIZhiGr1aqePXvWy32SMBoAp9MpSTIMQ3a73cfVAAAAAPC1koxQHwiFDYDFYpFhGK7/DwwMlMVi8WFFQFmGYaioqEgS2yf8D9sn/BXbJvwZ26f/KjlIVJ+vCaGwAQgMDCxzhLBLly4KCgryYUVAWYWFhdq8ebMktk/4H7ZP+Cu2Tfgztk//9fvvv8tut9fraWVMNAMAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGNcpBAAAgF8xDEMOh0OGYfi6FNMquXB9yd9cvL7+WCwWBQQE+NVzSigEAACAX8jOztaxY8eUnZ0tp9Pp63JMrXQg37dvn18FmMbAarUqIiJC0dHRioiI8HU5hEIAAAD4Xnp6ulJTU31dBkoJCQnxdQmNltPpVFZWlrKyshQfH6/Y2Fif1kMoBAAAgE9lZ2e7AmFoaKiaNGmi0NBQjk75kNPpVH5+vqTicGi1MhVJfTEMQ3l5ecrIyFBeXp5SU1MVFBTk0yOGhEIAAAD41LFjxyQVB8LWrVsTQPyA0+l0nVdos9l4TepZUFCQIiMj9ccffygvL0/Hjh3zaSjk1QUAAIDPGIah7OxsSVKTJk0IHzANq9WqJk2aSCo+Wu7LiZV41wEAAMBnHA6Ha1KZ0NBQH1cDeFfJNu90OuVwOHxWB6EQAAAAPlP66AjnEMJsSm/zHCkEAAAAAPgEoRAAAAAATIxQCAAAAAAmxiUpAACNmsVicU3zzflKAACURygEADRqNptNnTp18nUZAAD4LUIhAKBRy96zTEf+TJJU82ugBYQlKOy0AZ4uDQAAv0AoBAA0ao6cFOUe2SlJirDGKSAgwMcVAQDgX5hoBgAAAABMjFAIAAAAACZGKAQAAAAAEyMUAgAAAKjU5s2bdeTIkUqX5+Tk6Oeff/ZiRahvTDQDAAAA+IGVK1fqk08+0bZt25SamqqePXvqoYceUocOHSps/+eff2rVqlVasWKF3nzzTaWkpOjRRx/V+vXrdc011+jhhx92tc3OztaCBQv07bffKjU1VZmZmWrevLl69+6tMWPGqH379q62qampWrVqlbZu3ar/+7//065du7RkyRI1a9bM1WbNmjXaunWrNm3apB9//FHnnXeeevXqVaa+rKwsbdiwQV988YWuvPJKXXzxxfrvf/+rZ599Vlu3bpXD4VDLli01aNAgjR492nVNWXgfoRAAAADwoczMTN15553KycnRc889p/bt22vHjh26++67df311+s///lPmdD26aefasaMGTp8+LAkqWXLliosLNStt96qnTuLZ1vOyclxtf/99991xx13qEuXLpo5c6Zat26tnJwcffzxx3ruuef06aef6tFHH9WwYcMkSRkZGUpJSdHOnTu1a9cuSVJkZGSZmvft26esrCytWbNG+fn5ZQJdRkaGrr32Wh04cEBOp1OS1K9fP61du1Z33HGHIiIiZLPZlJ2drV27dmn27NlatmyZ5s2bp+bNm3vgGUZ1GD4KAAAA+NCkSZO0YcMGzZo1yxX+EhMTNWPGDGVlZenvf/+7evXqpT59+qhv377q1KmTFixYoJtuukmSZLPZNGvWLNlsNg0YMEA2m00xMTGSpAMHDmjs2LGKjY3VnDlz1Lp1a0lSeHi4Ro4cqaeeekp2u12PPPKIVq9eLUnq3Lmz7rjjDt18882uGi0WS5mab7zxRk2ZMkVdunSRpDLXgI2JidFHH32kFStWKCwsTFLxkcXXXntNH3zwgVavXq1169Zp0aJF6tatmyRp165devLJJ+v7qUUNEQoBAAAAH9mwYYPWr1+vU0891RXYSnTr1k3dunWT3W7XyJEjtXbtWn3//ffq3Lmz2rRpo65du0oqHu6Zl5enxYsX68UXX9Tjjz+u4OBgSdKzzz6rzMxM3XDDDQoMLD9I8Morr1Tnzp3ldDr1zDPPlFlms9mqrb+iNhaLRVFRUWrRooWioqIkSdu2bdPcuXPVuXNnV7vu3bvrvffeU6tWrSRJ33zzjY4ePVrtOlH/CIUAAACAj/z3v/+VJNeRvZOdeuqpkqS0tLRK7yMvL0+TJk1yHc0bMmSI7r77bmVkZOi7776TJJ1xxhkV9rVYLBowYICk4qN1u3fvrtXjqEzJEcQRI0YoJCSk3PKwsDDdfvvtkiSn0+l6PuBdhEIAAADAR4qKiiQVTxpTkePHj0sqHk5amZiYGDVt2rTc7Zs3b3ad03fyOYGlnXXWWa6/9+zZU33RbigJqiXDSCty/vnnu/4ufS4kvIdQCAAAAPjImWeeKUk6cuSI1q1bV2bZsWPHtGXLFsXFxenKK6+s9D7Cw8MrvP3YsWOuv7OzsyvtXzpQOhyOGtVdn0pPLhMfH+/19YNQCAAAAPjM2Wef7Rq++cADD+jXX3+VYRjau3ev7rrrLgUHB+uVV16p8khfZUr3SU5OrrRd6WGdp5xyiutvwzAq/Lu+lQTW8PBwnX766R5bDyrHJSkAAAAAH3r++ef10Ucfad68ebrpppsUExOjFi1a6K9//ateeumlWl+/r/T1DX/55Rf169evwnYZGRmSpNjY2DITwZQMbZWkgoKCCvvm5ubWqrbSduzYIal40puSCXLgXYRCAAAAwIcCAwMVHBysuLg4zZ8/v96GUJ566qnq0qWLkpKStHTpUk2ePLnC2UI3bdokSRo+fLiCgoJct5f+e//+/WVCpiStXr1aW7ZsqXOdy5YtU2RkpCZMmFDn+0LtMHwUAAAA8KF3331X999/vxITExUXF1ev933XXXdJklJSUvT222+XW15QUKD33ntPXbt21cSJE8ssa9OmjWuimJdeekkpKSmSio8gLlq0SPPmzat0VtOT7du3r8Lbk5KStHz5cr3yyiv1/thRcxwpBAAAAHxo7dq1korD4ccff6wmTZooICDAtTw0NFRt2rTRqFGjdPbZZ7tuP3TokKTiCWWys7MrHGZ6ySWX6K677tILL7yg559/Xjk5ORo9erSaNGmiffv26bHHHlPLli01Y8aMMkcGJalZs2bq37+/li9frq1bt6pfv36Kj4/XkSNHlJiYqPnz57uO7pXUUpkXX3xRqampuvXWW9W6dWsVFBToq6++0nfffacPP/zQdekN+AahEAAAAPChBx98UCEhIVq2bJlyc3MrPE9v+/btWrlypT7++GO99957+uabb5SZmSmpeKKWc889V+Hh4Zo0aZJGjhxZpu+ECRN01llnaf78+frwww/15ptvqm3btmrdurVuvPFG9evXz3U9wZNNmzZNLVq00LJly5SZmSmn06kbb7xREydOVEREhKvWdevW6ZxzztFFF12kmTNnlrufs846Sxs2bNCnn36q6OhotW7dWsOHD9cLL7zgOhoJ3yEUAgAAAD60Z88e7dixQ48++qiGDBmi4OBgWSwWGYahoqIirV27VnfddZfy8/O1aNEiPfHEE3riiSfcWkfv3r3Vu3dvt2sLCQnRAw88oIceeqjC5R999FGN7mfYsGEaMmSI2+uHd3BOIQAAAOAjs2fP1u23365JkyZpxIgRCgkJcR05s1gsstls6tu3r0aMGCGp8nPzgLogFAIAAAA+YBiG5s+fr6CgIP31r3+tsm3btm0lSTExMZ4vDKZDKAQAAAB84OjRozp+/LhiYmLKTfJysrS0NElSz549vVEaTIZQCAAAAPhASRjMyMhQTk5Ope2cTqeWLVum6OhoXXXVVd4rEKZBKJRUWFjouu4KAAAA4A2BgYG64oorZLfb9dZbb1Xa7oUXXtDevXs1Y8YMhYeHe7FCmIWpZx9NSUnRwoUL9eGHH+rqq6/WlClTyrXZvXu33nrrLa1bt05paWkKDw9X165ddf311+tvf/ubD6oGAABAY/HAAw9o27ZtmjNnjvLy8jRq1CjFx8fL4XBoy5Ytev3117Vx40a9/PLLuvjii31drlsKCgpcl834448/fFsMqmS6UJidna0VK1boiy++0KFDh5Samqpjx45V2Pann37S+PHjlZeXp9jYWDVv3lwpKSlau3at1q5dq7Fjx1YYJAEAAICaiIyM1MKFC7Vw4UJ9+eWXWrx4sQzDUEREhE499VQNGDBATz31lKKionxdqlv69eun9PR05eXlSZJeeeUVvf/++7rmmms0depUH1eHk5kuFH755Zf65ptvNGbMGJ1//vm66aabtH79+grb/vDDD2rVqpWee+45denSRZJ05MgRPfroo/r222/1+uuva+DAgTr99NO9+RAAAADQiNhsNt1444268cYbfV1Kvfn6668VGBjIhekbCNOdUzh8+HC98cYbOv/886ttm5GRoSlTprgCoSQ1a9ZMzz//vE455RRJ0vLlyz1WKwAAANAQ2Ww2AmEDYrpQ6A6r1VpheAwODlbfvn0lSVlZWV6uCgAAAADqj+mGj7pj+vTplS4r+eWjTZs23ioHAAAAAOodobCWtm7dqsDAQA0aNMjr67bb7V5fJ1CV0tsk2yf8hcVikc1mk8PpdN1W+u+qWE+0s9vtMgzDI/UB7DuLFRUVud5nTqdTzhq+T+FZpfd9hmHwuniI0+l0Pdc1/czxxOcSobAWkpKS9Ntvv2nMmDGKj4/3+vq3bdvm9XUCNcX2CX8RERGhTp06KavUDNNHjxypWV9LnKIlHTp0yDVznrvy8vLkcDhq1Rfmw75TCgkJUX5+voqKinxdCk6Sn5/v6xIarZIfRfLz85WUlOSzOgiFbrLb7frnP/+pc845R5MnT/Z1OQAADwiKaC5Jat26da3vY/v27crOzq6vkgAA8BhCoZseeeQROZ1Ovfzyy7LZbD6poXPnzj5bN1ARu93u+pWb7RP+ouTc76joaOWeOEDYtFkzBVirn2MtJDJSkpS95ys5cg65td6A8ARFtBukdu3aMfQUVWLfWayoqEj79u2TVHy00KzPg78pOXolFb8uzCTqGXa7XRaLRaGhoWrbtq0CA6uPZ0lJSfV+RJ1Q6IYZM2Zo165dmjdvniJPfGHwBZvNpqCgIJ+tH6gK2yf8TekQGGC1KiAgoNo+1hNffoy8VDmzD7i1PuuJ9fHFFu4w877TYrG4AofVanW9h+Bbpc8htFgsvC4eYrVaXdu/zWar0WeHJwI6obCGZs+erU2bNuntt99WRESEr8sBAAAAgHpBKKyB2bNn67ffftMbb7yhkJAQX5cDAAAAAPXG9MeBq5sZruQI4SuvvFJhIDQMg1myAAAAADRYpj5SaBiGjpyYnjw9Pb3c8hdeeEEvv/yyYmJidMkll5Trm5eXp7y8PE2cOFGTJk3ySs0AAAAAUJ9MGQrnzJmjlStXKj09XcnJyZKkjz/+WL///rtCQ0M1ffp07d27Vy+//LIkKTMzs8r740ghAACA97z2035lFzS8718RwYG67bw2vi4DKMeUoXDixImaOHFilW06duyo7du3V3tfhYWFTNELAADgRdkFRcpqgKHQW1JSUpSQkODrMirl7/XVxaFDh9SiRQtfl+E2059TWFdBQUFMOw4AAACf++GHH3T11Vdr5syZvi6lQv5eX13s379ft99+u0aMGOHrUmqFUAgAAAA0AuvXr9fWrVsVGxvr61Iq5O/11cX+/fv1/fffKyYmxtel1AqhEAAAAGgEAgOLzwzz10uo1bS+6q4O4A3uzhlSMnLQX5/76pjynEIAAACgsbFai4/3BAQEuG4rLCzUnj17tG3bNq1bt04Wi0VPP/205s6dq/fee09Op1NLly51HeFKSUnRm2++qa1btyozM1NZWVn6y1/+onHjxqlz585l1rd//34tWrRIu3fv1pEjR7R//3517txZ48eP13nnnVej+kru59dff9XOnTv1f//3f7rgggt07bXXavr06fr555+Vm5urDh066Oabb9YVV1yhwsJCzZ49W19//bXr/MTrrrtON998c4VzfWRmZmrevHnasGGDcnJydOTIEXXt2lVjxozRueeeK0natm2bNm/erO3bt2vdunUaPXq0+vXrp+nTp+vHH39UVlaW2rRpo2HDhmnkyJGux1KiZL0nP7aGglAIAAAANGA7d+7Uxx9/rGXLlkmSPvnkE61Zs0aSFB8fr6ZNm2rDhg3asWOHrrzySr388suaPXu2nE6nLBaLIiIiJEnvvfeeXn31Vc2aNUsPPPCA8vLy9Mcff+jOO+/U0KFD9frrr7vCXnJysq644gotWbJE7du3d9Vx++23a8yYMXr99dd1wQUXVFvfWWedpSFDhshms2nVqlXauXOnYmNj9dRTT+maa67RbbfdppUrV2rOnDmaMmWKkpOTtXnzZvXt21ezZs1SUlKSnnvuOT377LM6evSo7rvvvjLPzVdffaVp06bpscce09133y1JSktL06RJkzRy5Eg9/fTTuvrqq+VwOBQZGamvv/5aaWlp+v333/XLL7/oiiuu0JgxY7Rx40bNmjVLTz/9tFJSUnT//fdLko4fP6558+bp+++/lyTt2LFDw4cPl1R81HD+/PkNYlJKQiEAAADQgHXs2FH333+/MjMztWTJEt1xxx0aMmRImTYfffSRHnroIW3YsEGxsbHasGGD0tPTddtttykwMFArVqzQE088oalTp+qcc86R0+mUJLVu3Vp33XWX7r33Xj311FP64osvJBUfEbvmmmtcgbCkjsmTJ2vKlCl66623XKGwJvUlJiZq9erV2rlzp/bt26fPPvtMQUFBkqRu3brp8OHD+s9//qNZs2bpo48+0hlnnOFaFhYWpsmTJ2v+/Pm66667XP22bt2qe++9V9ddd50GDBjgWldcXJwefPBBDR8+XM8884wuu+wydevWTd26ddOCBQuUlpamvXv36p133nEFuq5du6pp06a6++67tXDhQk2ZMkWBgYGKjIzUpEmT1LJlSz3wwAPq16+fnnnmmXp7bb2FcwoBAACARq5kuKPT6dR9992nsLAwnXrqqfryyy8lSW+99ZYk6cwzzyzXt1OnTpKKj/hlZWVJkhISEvToo4+Wa9utWzdJcl0LvDZ69+7tCnYlevXqJUmKjIx0BcKTl9nt9jLXF58/f76KioqqfEyZmZnavXt3ueV//etfyx3hO//88yVJubm5ysjIcPNR+TeOFAIAAAAm0blzZ9eELyUKCwv166+/SpL++c9/KiIiQoZhuI4WOp1OnXLKKYqNjS03CUxhYaE2bdqkLVu2aPv27a6AZbfba11jdHR0udtKznksGepaWlRUlOvv0hPErF+/XpL04osv6r333ivTxzAMtWjRQk2aNCn3fEjF4fNkpdft7kQ0/o5QCAAAAJhEaGhoudsyMzNdIefZZ5/VGWecIafTqby8PFefkydWyc/P14wZM7R48WK1bNlSV1xxhUaMGKGYmBj179+/TjVWFNJKJnA5uQ5JZa4ZXhJkJenw4cOSpClTpmjgwIFu1VDReYClJ5EpvZ7GgFAIAAAAmFjpI21paWk16vPkk09q0aJFuv766/Xoo4+6QlRKSopHaqyNqKgoHT16tMaPyR2GYdT7ffoS5xQCAAAAjYi7gSUkJETt2rWTJP3444/Vti8oKNAnn3wiSbr99tvdnl3TW4Gqa9eukmr2mOpLQw2LhEIAAACgESgJZ7WZBOXGG2+UJC1evFgHDhyosm1OTo7rnMGSc/1KbNy40SP11UbJY1q1apV+++23er3vk4ePevux1TeGjwIAAKBBiQhumF9hPV13fHy8JOnNN99USEiIoqKiFBUVpb59+yo7O1uSXOcJnmzEiBHavHmzlixZolGjRumpp55S9+7dJRVPqrJ8+XKlpqZqzJgxio2NVYcOHbRr1y59/vnnGj58uJxOp5YuXaoNGzZIUrkJaaqrTyqe1bP0v6Xl5+dXuqy00o+vb9++mjBhgl5++WXdfvvteuKJJ9S/f39ZrVYZhqHVq1drw4YNmjx5sqtPVTWUlp6ertatW7v+PyEhQZK0Zs0azZo1S4mJicrNzdWwYcOqvB9/YTEa6jFOE/n999/LzOB0xhlnlJumF/ClwsJCbd68WRLbJ/zPsc1vK3ln8Qx0cXFxZSYKqExwwrkKO22gjm+Zr6Ks/W6tLzCqjSK7japVrTAX9p3F7Ha7du3aJUnq0KFDmUlD4J7s7Gzde++9Wrt2rWw2m3r16qXLLrtM7777rnbt2qWcnBxZLBa1a9dO4eHhmjdvnsLDw8vcx/Lly7Vo0SL997//VVBQkNq3b68mTZrovPPO0xVXXOHaTnfv3q2nnnpKmzZtUseOHXXBBRdo+PDhOn78uAYPHiyLxaKbb75Zl112mU4//fRK63viiSf0zTffaMmSJUpKSlJRUZFCQkLUqVMnjR07Vl27dtXkyZP1559/6ujRo5KKr2nYs2dPPf744666Sy4x0aJFC8XFxemdd95RSEiIJOnnn3/WggULXDOsdurUSVFRUerRo4euueYahYeH69///rd++OEHbdmyRVLx7KPt2rXTgw8+qB49epRbT1xcnFq0aKG3337b9RzOmDFD77//vux2uzp27KhHH320TN+K1Gb7L8kGNpvNFdzrilDYABAK4e/4YgN/RiiEv2LfWYxQ6J+qm30UxUdQK5op1R3+Egp5dQEAAADATXUNhP6EUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAACg0TIMw9cl+L1AXxcAAAAAuCPr9zfktB/3dRlus9oiFdX9Vl+X0ej9+eef2rJli7Zt26b169frvPPO0x133OHrsvwaoRAAAAANitN+XEZhwwuFTl8XYBI7d+7U1q1b9emnn+rQoUPq3bu3r0vyewwfBQAAAOBRqampXhvGeckll2jy5Mk655xzvLK+xoBQCAAAAMAj0tPTdc8996h///6y2+1eXXdAQIBX19eQEQoBAAAAeERGRoaWLl2q0NBQBQUF+bocVIJQCAAAAKBGioqK3Gpvs9kkSSEhIZ4oB/WEiWYAAACARuLnn3/W0qVLdeDAAe3cuVMOh0MtW7bUI488ou7du5dp9/777+uPP/6QYRg6fPiwOnbsqKFDh+ryyy93tfvjjz+UlJSkXbt2af369brooos0duxY/etf/9Ly5cuVnp6uFi1a6IorrtDtt9+u4ODgMvVYLBZJFQ/lPHjwoHbu3KmNGzfqm2++0UcffaTt27dr+vTp2r17tx588EENGzZMkuRwOPTZZ5/ps88+0/Hjx5Wbm6vjx4/rL3/5i0aPHq0ePXp44Nk0D0IhAAAA0MAdOHBATz75pAYNGqRp06bJYrHIMAx99913mjp1apm206dP1+eff64ZM2boggsukCQdPnxYDz30kO699159/fXXmjlzpqTiMBYeHq41a9Zo+/btatOmjR555BH97W9/09ChQ5WUlKRZs2Zpzpw52rlzp2bPnu3q9/rrr+vHH3+UVDyMdPjw4a4a3nzzTS1YsEAHDhzQ119/7XoMY8eOVV5enhwOhytQHj9+XBMmTFBOTo5mzZqlNm3aSJK2bt2qO++8U9dff72mTJmiW265xYPPcONGKAQAAAAasNTUVF177bUaNGiQrrrqKtftFotF/fv3V8+ePRUfHy+pOIzNnz9fM2fOdAVCSWrevLleeuklDRo0SN98841efPFFjR8/Xqeddpq6du2qlStXavv27dqyZYs+/fRT1xHBLl26qEOHDho2bJi++eYbpaamKj4+XgEBARo3bpx69uypn3/+WWeccYbefffdMnXfd999cjqd6tKliyRp1qxZ+uijj9SyZUtNmzZNERERkqQpU6Zo48aN+vbbb9WiRQtX/65du+rVV1/V3//+dz333HPq2LGjLrroIo88x40d5xQCAAAADdjcuXN19OhR/e1vf6tw+c0336ymTZsqPz9fc+bMkdVq1V//+tdy7YKCgjRy5EhJ0nvvvaf8/PxybS688MJyQ0S7d+/uCnCHDh1yq3ar9X9x5KabbtJpp52moKAgPfXUUxo4cKA2bNig77//Xt26dSsTCEt06NBBffv2lSS98cYbbq0b/0MoBAAAABqw77//XpJcRwNP1qdPHwUGBurXX39VTk6OoqKiFBoaWmHbkqOHeXl52rVrV7nlkZGRFfYrCYXuTkRTWrdu3crdtmbNGklSQkJCpf369OkjSdq0aVOt1212hEIAAACgAUtLS5NUHOSqcvTo0WrblQ6WFV1XsOQ8v5OVTCTjdDqrLrYKYWFh5W4rqTk3N7fSfiU1OxyOOq3fzAiFAAAAQANWMpxz//79VbaLjo6WJBUUFOjgwYMVtil9LcFTTz3V7VoMw3C7T1VKat67d2+lbUoef8uWLcsMR0XN8awBAAAADVjXrl0lSStWrKiyXffu3V3XDVy9enWFbUrOCTznnHPUvHnzequxtmHxrLPOkiT9+eeflQbDkppLX0oD7iEUAgAAoEGx2iJlCWp4/1ltFZ+PV1fXX3+9JOnLL7/Utm3bKm0XHR3tCk5vv/22CgsLy7VZtmyZgoKCNGXKlFrVcvLwzZLhphkZGbW6vwsuuEBt27aVJL322msVtlm2bJlat27tmiQH7uOSFAAAAGhQorrf6usS/MrgwYO1ceNGvfvuu7rlllv09NNPl7k0w8aNG9W8eXO1atVKDz30kLZv3+66xt/jjz+uuLg4SdKnn36q999/Xy+88IK6d+9e5tzDknP6qjq3T5LS09PL/H/J+X67du3StGnTdPbZZ+vIkSMaNWqUJCk7O7vMOk6e2dRms+nf//63Ro8ercWLFyshIUG33367QkJCVFhYqOeff15paWmaO3eua6jpyTXn5ORU/ySanMWo74G/qHe///57mRN9zzjjjDLjvQFfKyws1ObNmyWxfcL/HNv8tpJ3rpckxcXFuSZDqEpwwrkKO22gjm+Zr6Ksqs/ROVlgVBtFdhtVq1phLuw7i9ntdtcslx06dHANb4T7vvvuO33wwQfavHmzgoOD1bFjR4WEhOjMM8/UzTff7Hpu8/Ly9NZbb2nZsmVKSUlRhw4dFBISor/85S+69tprFRcXJ6fTqby8PH344YdaunSpkpKSVFRUpJCQEHXq1Eljx44tcwmMfv36KTk5WTExMWrTpo2ee+451xG++fPna86cOcrJyVGrVq00depUnX322brllluUnp6uAwcOSCq+VmJ8fLxuueUWDR48uMxjO3TokF555RWtXr1aubm56tKli2w2m/r166drrrmmzPtnwYIFWrJkSbmax4wZo4EDB3r4VXBPbbb/kmxgs9nUvXv3eqmDUKjinXJ6enqVU936EqEQ/o4vNvBnhEL4K/adxQiF/qkkFEpSaGioWxO4FBUVKSAgoNKZSvE//hIKTT18NCUlRQsXLtSHH36oq6++utKx0wcPHtTs2bO1Y8cOOZ1OHTt2TH379tW4ceNch9sBAAAASIGBpo4YDZLpXrHs7GytWLFCX3zxhQ4dOqTU1FQdO3as0vb79u3TddddpzvvvFPTp0+XxWLR9u3bdeutt2r58uVavHhxvc7MBAAAAADeZLrZR7/88kt99tlnGj16tD7//HN16tSp0rZOp1MTJkxQjx49NGLECNch8E6dOmnq1KlKS0vTtGnTvFU6AAAAANQ70x0pHD58uIYPH16jtt999512796t8ePHl1t26aWXKiQkRCtXrlRGRoaaNGlS36UCAAAAgMeZ7kihO7755htJUmJiYrllQUFB6tKli4qKipSUlOTt0gAAAACgXpjuSKE7duzYIUmVnjPYsmVL/frrr9q/f7/OP/98r9VVeiZSwB+U3ibZPuEvLBaLbDabHKUupOw46aLKlXGemJjb6TTkcDjcWq/1xDrsdruY4BtVYd9ZrKioyPVecTqd5S5+Dt8ovf8yDIPXxUOcTqfrua7p54YnPlsIhVU4fPiwJCkiIqLC5c2aNZNU9qKb3rBt2zavrg9wB9sn/EVERIQ6deqkrFKTiR09cqRGfWMjjytCUnb2cR1LS3NrvWHOaEVL2rNnj9c/H9Bwse+UQkJClJ+fr6KiIl+XgpPk5+f7uoRGq+RHkfz8fJ+OPmT4aBWysrIkVT6tbnBwsCRz/7oHAAAAoGHjSGEVSn6pKioqqvCCsiUX8QwLC/NqXZ07d+bCrvArdrvd9Ss32yf8RcmM0VHR0co9cYCwabNmCqjBBZhDIiMlSRERkQp283q0tujiicfatWvH8FFUiX1nsaKiIu3bt09S8dFCsz4P/qbk6JVU/LpwIXrPsNvtslgsCg0NVdu2bWt0jcekpKR6P6JOKKxCUFCQCgoKdOzYsQrPKywoKJBU+fBST7HZbBWGVMAfsH3C35QOgQFWqwICAqrtYz3x5cdqtdSofZm+J9bHF1u4w8z7TqvV6gocFovF9R6Cb5U+h5DXxXMsFotr+w8KCqpRKPREQOfVrULTpk0l/W8Y6clycnIkSW3atPFaTQAAAI1JQECAK3Dk5eX5uBrAu0q2eWsNf7T0FEJhFdq1aydJSklJqXD5oUOHFBgYWOElKwAAAFA9i8XiGnWVkZHBLJcwDafTqYyMDEnFIw99OUSX4aNVOP/887VmzRpt2rRJffr0Kbd869atuvDCCxUdHe2D6gAAABqH6OhoZWVlKS8vT3/88YeaNGmi0NBQzmPzIafT6TpvzW63M3y0HhmGoby8PGVkZLiOFPo6T5g+FFZ1/amrrrpKL730kpYuXarx48eX2TFt3LhR6enpGjNmjDfKBAAAaLQiIiIUHx+v1NRU5eXlMYzUDxiG4Zosq/R5b6h/8fHxXp+j5GSmjvyGYejIiWtWpaenl1vetGlTPfjgg9q1a5emT5+uwsJCV9snn3xSEydO1LnnnuvVmgEAABqj2NhYtWrVSlFRURyV8hP5+flco9BDrFaroqKi1KpVK8XGxvq6HHMeKZwzZ45Wrlyp9PR0JScnS5I+/vhj/f777woNDdX06dPVsWNHSdKwYcMUFxen119/XQMHDlR8fLyCg4M1duxYDRo0yJcPAwAAoFGJiIhQRESEDMOQw+Hgsi4+ZLfbXRdTb9u2LTMq1yOLpXhma386+mrKUDhx4kRNnDixxu0vvvhiXXzxxR6sCAAAACUsFkuNpuaH55QO5IGBgYTCRo5j8wAAAABgYoRCAAAAADAxQiEAAAAAmBihEAAAAABMjFAIAAAAACZGKAQAAAAAEyMUAgAAAICJEQoBAAAAwMQIhQAAAABgYoRCAAAAADAxQiEAAAAAmBihEAAAAABMjFAIAAAAACZGKAQAAAAAEyMUAgAAAICJEQoBAAAAwMQIhQAAAABgYoRCAAAAADAxQiEAAAAAmBihEAAAAABMjFAIAAAAACZGKAQAAAAAEyMUAgAAAICJBfq6AAAA/I61+OMxICze7a4BES3ruxoAADyKUAgAwEmsQVGSpLDTBvq4EgAAPI9QCABAJb745VcdTP3TrT5ndOis8zp11Ffb0pSSle9W34SoEA3sHOdWHwAA6opQCABAJQ5lZGpPyiG3+pwSf6okKSUrX/sy8jxRFgAA9YqJZgAAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkF+rqAhuSzzz7TokWLlJ2dLYvFouPHj6tPnz4aN26cEhISfF0eAAAAALiNUFgDhmHoH//4h3799VfNmTNHnTp1kiQlJydr/PjxGjZsmBYvXqzmzZv7uFIAAAAAcA/DR2vg888/12effab777/fFQglqWXLlnrssceUlpamhQsX+rBCAAAAAKgdQmENLF++XJLUunXrcsu6desmSTp69KhXawIAAACA+kAorIGioiJJ0k8//VRuWXJysiSpT58+Xq0JAAAAAOoD5xTWwHnnnacVK1Zo5syZCgwM1PXXXy+r1SrDMPTiiy9qyJAh6t+/v9fqsdvtXlsXUBOlt0m2T/gLi8Uim80mh9Ppuq3031UxDMP1r9NpuLXekr5Ow5DD4XCrb0l9drvddT9ovNh3wp+xffovT3w+EApr4LrrrtMPP/yg1atX6/HHH9cHH3ygW2+9VatWrVKPHj00evRor9azbds2r64PcAfbJ/xFRESEOnXqpKxjx1y3HT1ypEZ9ExJyFSqpyF6knJwct9ZbZC8eXZKbm6u0tFS3+sZYYyS11Z49e5Sdne1WXzRs7Dvhz9g+Gz+Gj9ZAUFCQ5s6dqylTpkiSdu7cqalTp2r58uXas2ePUlPd+9AHAAAAAH/BkcIaWrdunZYtW6Y5c+YoISFBb731lpYtW6aFCxdq6dKlmjt3rs4++2yv1NK5c2fZbDavrAuoCbvd7voVke0T/sJisUiSoqKjlXviAGHTZs0UYK3+99DQsDBJUqAtUOHh4W6tN9BW/NEaFhamuLh4t/rGxBSvt127dgwfNQH2nfBnbJ/+KykpyTXnSX0hFNbAqlWrNGnSJM2bN09nnXWWJOn555/Xbbfdpvvvv19JSUmaMmWKvv32WwUGev4ptdlsCgoK8vh6gNpg+4S/KR0CA6xWBQQEVNunJFBaLBZZrRa31lfS12qx1GhdpZXUypcv82HfCX/G9ulfSj5n6hPDR6tRVFSkadOmqXfv3q5AWKJz586aP3++mjRpokOHDmnHjh0+qhIAAAAAaodQWI19+/YpOTlZHTt2rHB5dHS0evXqJUn8ggIAAACgwfF6KPz9999r1C4jI8PDldRMyXjd3NzcStvk5OSoefPmatOmjbfKAgAAAIB64fVQOH78+Bq1u++++yq8WLy3tW/fXnFxcVq1apUKCwvLLU9LS9P69es1YcIEzgEBAAAA0OB4PRTW9MTIZs2a6dlnn/VwNdWz2Wx68sknlZaWpnvvvVfp6emuZXv27NGECRN00003acSIET6sEgAAAABqxy9nH3U4HNqyZYv++OMPX5ciSbr44ou1ZMkSzZ07V0OHDlVMTIwiIyPVvHlzPfjgg+UmoAEAAACAhsJjobCwsFADBw5UeHi4wsLCFBQUJIvFouzsbL300ktV9vvpp5+0c+dORUREeKo8t3Xs2FH/+te/fF0GAAAAANQrj4XCoKAgDRgwQG+//XaZ2y0Wi+bMmVNl35IL9g4dOtRT5QEAAAAA5OHhozfddJOk4oAYEBAgi8WiefPm6dJLL62yn81mU9euXXXttdd6sjwAAAAAMD2PhsJTTjlFU6dOLXPbnj179PTTT3tytQAAAACAGvL67KPXXXddhZd2AAAAAAB4n9dDYa9evRQUFFSjtnfeeafr4vEAAAAAgPrn9VDojp9++kkBAQG+LgMAAAAAGi2vX6fQ6XRq+fLl2rRpk3JycuRwOOR0Ol0zjkrFl6X4888/lZ2dXeOL3QMAAAAA3OfVUJiTk6NRo0Zpy5Yt1bY1DINACAAAAAAe5tVQOHfuXP33v/+VVHyZirZt2yoiIoIhogAAAADgI14NhV9++aUsFovOO+88zZw5U02aNPHm6gEAAAAAJ/HqRDOpqakKCAjQ008/TSAEAAAAAD/g1SOFCQkJcjgcio+P9+ZqAQAAAACV8OqRwssuu0zHjh2Tw+GoUfvt27d7uCIAAAAAMDevhsK77rpLffv21eLFi2vUfvTo0Z4tCAAAAABMzqvDRw8ePKh77rlHDz30kHJyctS7d29FRkZWeOmJQ4cOKScnx5vlAQAAAIDpeDUU3n///frll18kST///LM3Vw0AAAAAqIBXh49eddVVMgyjxv8BAAAAADzLq0cK+/TpI0maN2+eWrVqVWXbtLQ0jRw50htlAQAAAIBpef2SFLGxsTr77LMVGFj1qlu2bKnw8HAvVQYAAAAA5uTV4aOS9O6771YbCEvMmzfPs8UAAAAAgMl5PRS2b9++xm3ffPNND1YCAAAAAPB6KHTHDz/84OsSAAAAAKBR89o5hevXr9c777yjiIiICq9LeLKUlBRlZWV5oTIAAAAAMC+vhcLTTz9dP/74o/Ly8mrU3jCMGoVHAAAAAEDteS0UhoWF6fzzz9eaNWvUpUuXCiebycvL05YtW9S9e3eFhIQoNzfXW+UBAAAAgCl59ZIUZ555ptq1a6d77rmn0jarV6/WO++8oxdeeEFhYWFerA4AAAAAzMerE8107dpV8fHxVba58MILdemll2rUqFEcKQQAAAAAD/NqKOzTp49uuOGGatsNGTJE27dv18svv+yFqgAAAADAvPzykhQBAQFq1qyZli5d6utSAAAAAKBR8+o5hTWVmpqq1NRUWa1+mVkBAAAAoNHwWihMT09XWlqaIiMjK7zUhGEYys/P1x9//KHZs2fL4XCoadOm3ioPAAAAAEzJa6EwOTlZw4YNq/G1By0WiwYPHuzhqgAAAADA3Lw2PvOMM85Qx44dZRhGtf9J0sCBA6u8dAUAAAAAoO68ek5h79691bRpU40dO1Y2m63MMsMwVFRUpKCgIJ122mlq1qyZN0sDAAAAAFPyaijs1KmToqOj1adPH2+uFgAAAABQCa+GwosuukgZGRneXCUAAAAAoApeDYVxcXGKi4tz/f+ePXu0Z88eOZ1OnXLKKercubMCA/3yKhkAAAAA0Cj5JIGtXbtWzzzzjHbt2lXm9tjYWF177bUaP358uXMOAQAAAAD1z+tXh1+0aJFuu+027dq1q9yso0ePHtXLL7+s6667Tunp6d4uDQAAAABMx6tHCnfv3q1p06bJ4XCoc+fOGjBggM466yw1bdpUISEhysjI0G+//aYFCxboH//4h9544w1vlgcAAAAApuPVUPj222/LYrFoxowZuuKKK8otP/XUU3XGGWdo+PDhGjdunFasWKF+/fp5s0QAAAAAMBWvDh9dt26dpk6dWmEgLC04OFiPPvqoFi9e7KXKAAAAAMCcvBoKDx8+rGuuuaZGbdu2bav9+/d7uCIAAAAAMDevhsKYmBiFhobWqG1mZiahEAAAAAA8zKuhMD4+XgcOHKhR23nz5ikyMtLDFQEAAACAuXk1FPbv319PPfWUnE5npW0cDodef/11vfbaa+rRo4f3igMAAAAAE/Lq7KMjRozQggULNHToUN1www06/fTTFRUVpezsbCUnJ+uXX37R0qVLlZKSIkkaNWqUN8sDAAAAANPxaiiMiIjQ7NmzNXbsWD388MMVtjEMQ5J0991369xzz/VmeQAAAABgOl4dPipJ3bt310cffaS//vWvslgsMgyjzH9dunTRq6++qnHjxnm7NAAAAAAwHa8eKSzRqlUrvfTSS8rMzNTmzZt17NgxRUZGqkOHDmrZsqUvSnKbw+HQ559/rs8//1zh4eF64oknFB0d7euyAAAAAMAtHgmFhw4d0g033KALLrhAp512miIiInTNNdfIai17YDImJkYXXnhhuf5r167VF198oaefftoT5dXZjz/+qMcee0xdunTRQw89pHbt2vm6JAAAAACoFY8MH92wYYMOHjyoDRs2KCMjQ6mpqVXOOHqyPn36KDk5Wbt37/ZEeXUyc+ZMjR8/XnfffbdeeOEFAiEAAACABs0jRwp/++03de3aVR988IGCg4NrdR9XX321vv32W7Vv376eq6u9xx57TIsWLdLrr7+u888/39flAAAAAECdeeRIYVJSkm699dZaB0JJOvfcc/XLL7/UY1V1s3jxYn3wwQeaNGkSgRAAAABAo+GRUJicnKxevXrV6T5OOeUU/fHHH/VUUd3k5ubq2WefVfPmzTVmzBhflwMAAAAA9cYjobCgoEBNmzat031YLBbl5+fXU0V18/HHHyszM1MDBgxQUFCQr8sBAAAAgHrjkXMKT55ltDYMw1B6eno9VFN3K1eulCSdeeaZ+vHHH/XBBx8oOTlZ2dnZateuna6//npdfPHFXqvHbrd7bV1ATZTeJtk+4S8sFotsNpscpSY6c9Rw0jPDMFz/Op2GW+st6es0DDkcDrf6ltRnt9td94PGi30n/Bnbp//yxOeDR0JhYGDd7zYlJUUhISH1UE3d/frrr5KkV199VSNHjtTMmTMVFBSkX375Rffdd59uu+02PfTQQxo5cqRX6tm2bZtX1gPUBtsn/EVUVJQ6duyoQkUqrFlHSVKeJNUgFwYZQQqVVFTkVE5OjlvrLbIXSSo+9SAtLdWtvjHWGElttWfPHmVnZ7vVFw0b+074M7bPxs8joTAiIkIZGRlq0qRJre9j3bp1io2Nrceqaic7O1u5ubmSpLvuuksDBgxwLTv77LP173//W8OHD9fzzz+vq666SlFRUb4qFQBQSslkZwlnDq/1fVht/vHjJAAAnuSRUNitWzdt2LBB/fv3r/V9LF68WJ07d67Hqmqn9C+1ffr0Kbf8zDPPVGJionbs2KHffvtNF110kcdr6ty5s2w2m8fXA9SU3W53/YrI9gl/UXIqw+f/t1G79++SJAUHh8hiqb7vmYln6LzE0xQQEKDw8HC31htoK/5oDQsLU1xcvFt9Y2LCJEnt2rVj+KgJsO+EP2P79F9JSUkqKiqq1/v0SCg855xz9P7779c6FH7//ff65Zdf9PDDD9dzZe4LDQ11/V3ZJTZOO+007dixQ1lZWV6pyWazMeEN/BbbJ/zNoYxM7fzzT0lSeHi4rNbqU2GrU4qvkWuRatS+NMuJ1Gm1WBQQEOBW34ATQZYvX+bDvhP+jO3Tv1hq8uummzwy++jAgQO1adMmffzxx2733bdvn6ZOnarAwEBddtllHqjOPdHR0YqOjpYkpaZWfW5IXWdcBQAAAABv80gojIyM1PXXX69HH31U7777bo37/fTTTxoxYoSysrI0YsQIxcTEeKI8t5177rmSpF27dlW4fO/evQoLC1OPHj28WBUAAAAA1J1HQqEk3XnnnWrfvr2mT5+uYcOGadmyZTp+/Hi5doWFhVq7dq3Gjx+vMWPGKD09Xa1bt9akSZM8VZrbbrrpJknF5zmeLCkpSTt27NDNN99cZqgpAAAAADQEHjmnUJKCgoI0Z84cjRkzRps3b9Y999yjgIAAJSQkKCYmRjabTRkZGTp48KDr2ieGYah58+aaO3euIiIiPFWa23r16qVx48bp1Vdf1bx58zR69GhJUnJysv7xj39o0KBBmjBhgm+LBAAAAIBa8FgolKRTTz1V//nPf/Tggw/q+++/V1FRkf78808lJye72pSeXa1Xr1569tlnlZCQ4MmyamXy5Mnq1KmT5s+fr/nz5ys+Pl42m02TJk3SpZde6uvyAAAAAKBWPBoKJSk2Nlavvvqq1q1bpyVLlmj16tVKT08vs/ycc87RsGHDdMEFF3i6nDoZPHiwBg8e7OsyAAAAAKDeeDwUlujdu7d69+4tSSooKFBmZqbCw8P9apgoAAAAAJiN10JhacHBwYqPd++CvgAAAACA+uex2UcBAAAAAP6PUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAOqVpdY9gwJq3xcAgNoiFAIAUI8sVlut+8aGBdVjJQAA1EygrwsAAKA6uXu/liM3xa0+wQnnKqhpFw9VVD1H3hEVHt7lXp/m3STF66ttaUrJynerb0JUiAZ2jnOrDwAAEqEQANAAOHJTVJS1360+tljfBUJJMhwFchZkutfHWShJSsnK176MPA9UBQBAeQwfBQAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKGwlgzD0O7du31dBgAAAADUCaGwFrZv364BAwZo7Nixvi4FAAAAAOqEUOgmh8OhBx54QPv37/d1KQAAAABQZ4RCN7355puKiYnxdRkAAAAAUC8IhW7Yu3evVq9ezbBRAAAAAI0GodANTz75pB577DFZrTxtAAAAABoH0k0NffTRR+rZs6fat2/v61IAAAAAoN4E+rqAhiA9PV1LlizR22+/7etSJEl2u93XJQBllN4m2T5RnywWi2w2m5xOpxwOh1t9DcM48UfZ25zOGvQt9a/TaVTVtH77nqjZaRhuP17HiQdmt9v/99jh19h3wp+xffovT+zjCYU18Oyzz2ry5MkKCgrydSmSpG3btvm6BKBSbJ/+KyAgQKGhoXW6j7y8PLfDSl1ERESoU6dOysjIUO6RNLf6JiTkKlRSUVGR67bc3Nwa9XWc6ONwOJSTk+PWeuvSt8he3Dc3N1dpaalu9Y2xxkhqqz179ig7O9utvvA99p3wZ2yfjR+hsBrr16+XxWLR2Wef7etSAKBOQkND1alTpzrdx/bt2wkcAAA0MoTCKhQVFWnmzJmaM2eOr0spo3PnzrLZbL4uA3Cx2+2uXxHZPv2XxWKRJH25NUWHsvLd6tsiKkSDuyaoXbt2Xh2aWFJzkyZNFGGNc6tvaFiYJCkw8H8fdWFhYa77rErAiT4BAQEKDw93a7116RtoC3TVGRcX71bfmJjix+vt1wi1x74T/ozt038lJSWVGQVTHwiFVVi4cKF27NihG2+8sczt+fnFX6ZSU1M1cOBASdK4ceN01VVXeaUum83mN0NZgZOxffq/tOxCHThW4FafgBOzLvvqS4HValVAQIBbfVzhz1L2Nqu1+lBoKfVvTdrXW98TNVstFrcfr69fI9QN+074M7ZP/1KTHzfdRSiswscff6zc3Fzt3bu3wuVFRUWuZZmZmV6sDAAAAADqB6GwCosXL67w9p9//lkjR45Uy5YttWLFCi9XBQAAAAD1h+sU1kLJbHLenIEPAAAAADyBUFgLhw4dkiRlZGSooMC983IAAAAAwJ8wfLSGdu7cqYceekhFRUXauXOnJKmgoEAXX3yxTjnlFI0cOdJrE80AAAAAQH0hFNZQx44d9eGHH/q6DAAAAACoVwwfBQAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAQCMQFGDxdQkAgAaKUAgAQCMQGxbk6xIAAA1UoK8LAAAA9eerbWlKycp3q0/XhEid27pJrfomRIVoYOc4t/oAAPwLoRAAgEYkJStf+zLy3OqTEBlc674AgIaP4aMAAAAAYGKEQgAAAAAwMUIhAAAAAJgYoRAAAAAATIxQCAAAAAAmRigEAAAAABMjFAIAAACAiREKAQAAAMDECIUAAAAAYGKEQgAAAAAwMUIhAAAAAJgYoRAAAAAATIxQCAAAAAAmRigEAAAAABMjFAIAAACAiREKAQAAAMDECIUAAAAAYGKEQgAAAAAwMUIhAAAAAJgYoRAAAAAATIxQCAAAAAAmRigEAAAAABMjFAIAAACAiREKAQAAAMDEAn1dQEOxdu1avfPOO9q6dasyMjIUGxurc845R7feequ6dOni6/IAAAAAoFY4UlgDc+fO1ZgxY7Rq1SoFBAQoJiZGqamp+uKLLzR06FCtWLHC1yUCAAAAQK0QCmvgm2++0aBBg/TDDz/o+++/15o1a/Tpp5+qXbt2Kioq0kMPPaSCggJflwkAAAAAbiMU1kB2drYef/xxxcXFuW7r3LmzXnjhBVmtVqWnp2vDhg0+rBAAAAAAaodQWI3CwkKdffbZioqKKrcsMTFR7dq1kyQdO3bM26UBAAAAQJ0x0Uw1goKC9NRTT1W63GKxSJLatm3rpYoku93utXUBNVF6m2T79F8Wi0U2m00Op1MOh8Otvg6nU1Lx62sYhifKq1BJzc5a1Oyq0yh724mHUnXfUv86ne493jr1PVGz0zDcf41cfWvx+talr4+2jcaAfSf8Gdun//LEvpZQWAe5ubnat2+fOnbs6NUZSLdt2+a1dQHuYvv0X9HR0erQoYOys48rMzPLrb7ZQcVBYf/+/V4dGREREaFOnTopIyNDuUfS3OqbkJCrUElFRUWu23Jzc2vU13Gij8PhUE5OjlvrrUvfIntx39zcXKWlpbrVN6dZ0Im+ecrMzHCrb15eiCQpOztHaWnuPc8x1hhJbbVnzx5lZ2e71Rf/w74T/ozts/EjFNbBokWLZLfbNWXKFF+XAgDVCgoqDg3RUdFq6ghyq290VGiZ+4D/sQYUnxESER6upk3dOzskIiy8+D6snFUCAGZEKKylP//8U7Nnz9a4cePUt29fr667c+fOstlsXl0nUBW73e76FZHt03+VfOG3ZyerIO0Pt/raba0ltdYpp5yihIQED1RXsZIh+k2aNFGENa6a1mWFhoVJkgID//dRFxYW5rrPqgSc6BMQEKDw8HC31luXvoG2QFedcXHxbvUNCy1+vI68wypI2+VW36Lmp0uKr9V6Y2KK19uuXTuGj7qJfSf8Gdun/0pKSiozCqY+EAprITs7WxMmTNCgQYM0efJkr6/fZrPxaz38FttnA+DIl1Ho5hBQR76ksgHLm6xWqwICAtzq4wp/lrK3Wa3Vh0JLqX9r0r7e+p6o2WqxuP14rSf6Go5C919fZ2Gt1xtw4scGvjDWDftO+DO2T/9Skx833cU4ETcVFBRo/Pjx6tGjhx5//HFflwMAAAAAdcKRQjcUFBRowoQJ6tq1qx544AFflwMAAAAAdUYorKGSQHjGGWfo7rvv9nU5AAAAAFAvCIU1UBIIe/bsqTvuuKPCNg6H48S5KozIBQAAANBwEAqrkZ+fr4kTJ2rNmjXatm2bPvjggzLLHQ6HcnNzVVBQoHfeeUe9evXyUaUAAAAA4D5CYTX+9a9/ac2aNZKkI0eOVNm2vqeGBQAAAABPIxRW4+GHH9bDDz9cZRvDMFRYWMh03AAAAAAaHEJhPbBYLAoODvZ1GQAAAADgNmZFAQAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAEzH/V1/UGCgB+pAeZYGt96gAF/VDACoL4RCADAZS0CQ232aRkZ4oBKczGK1Nbj1xoa5vz0BAPwLP/2i1r7alqaUrHy3+iREhWhg5zgPVdQ45e79Wo7cFLf6BIQlKOy0AR6qCPWlNu+hrgmROrd1kzr1rQve997hyDuiwsO73OsT30NSvPfX27ybpHjTbRtme7yAv/H2Z2hjf/8SClFrKVn52peR5+syGj1HboqKsvb7ugx4QG3eQwmRwXXuWxe8773DcBTIWZDpXien3SfrNZyFksy3bZjt8QL+xtufoY0dw0cBAAAAwMQIhQAAAABgYoRCAAAAADAxQiEAAAAAmBihEAAAAABMjFAIAAAAACZGKAQAAAAAEyMUAgAAAICJEQoBAAAAwMQIhQAAAABgYoRCAAAAADAxQiEAAAAAmBihEAAAAABMjFAIAAAAACZGKAQAAAAAEyMUAgAAAICJEQoBAAAAwMQIhQAAAABgYoRCAAAAADAxQiEAAAAAmBihEAAAAABMjFAIAAAAACZGKAQAAAAAEyMUAgAAAICJBfq6gIZkw4YNeu2115SRkaGsrCwFBgbqyiuv1OjRoxUUFOTr8gAAAADAbYTCGvrqq6/0yCOP6LXXXlOPHj0kSV9++aWmTJmidevW6fXXX1dAQIBviwQAAAAANzF8tAb++OMPTZ06VePGjXMFQkkaPHiwhg4dqrVr12rBggW+KxAAAAAAaolQWAPvvPOO8vPzddlll5VbduWVV0qSlixZ4u2yAAAAAKDOCIU1sHz5csXExCg+Pr7cstNPP10BAQHavn27CgsLfVAdAAAAANQe5xRW4/jx40pJSVGHDh0qXB4SEqKmTZsqLS1NBw4cUPv27T1ek91u9/g6qmKxWGSz2eRwOuVwONzq63A6JRU/BsMwPFFeo1HyPDtr8Txbvfw8l94mfb19NgR1eg+deD2dhvt9nSf6GpKcTve2C8O1XsOr7/u6vA9c6zLK3nainKr7lvrX7efKbH19tG34Sn1+BrLvhD/z1+3TV5+h/rS/8sT6LYavH5Wf2717twYPHqwePXpo4cKFFba5+uqrtXXrVi1atEjdu3ev9xo2btzo842vNIvFosDAQOXZi+Rw88tDgNWiUFugnE6nXz0mf2SxWGS1WuUsKpBh1OBbbJm+VlkDg3me/VTJa1ub95AtwKrgwADl2R2uDyh3++bbi1TkKHKrb1BgoIICA2u13rq87+vyPrAG2GSxBiq/0K4ip3sf/kGBNgXV+rkyW1/fbBu+Upf3b0N8vIC/8dVnaMn7t6ioyG/evxaLRWeddVa93BdHCquRlZUlSQoMrPypCg4OliSPDR/1lw2vRMmvm4GSAt0egGz41a9N/q74Vyyr3B7pbUgOnme/5nA4avceMpyy25116hsguT9bsg/f97V9HzgdhuSwK8BSm8dbl+fKbH3N95lQ6/dvA328gL/xyWeoH75/6zMjEAqrUVRU/ItpVYeYrdbirSosLMwjNVitVjmdTtcROgAAAADmVHK0siSD1AcSRjVKLkqfmZlZaZuCggJJUmRkpEdq6Nmzp0fuFwAAAACYfbQasbGxkv43jLQi2dnZstlsatGihbfKAgAAAIB6QSisRsuWLRUaGqrMzEzXEcGTpaSkKDExkaGdAAAAABocQmE1rFarevXqJYfDoc2bN5dbvnv3buXn5+vyyy/3QXUAAAAAUDeEwhoYMWKEJOnzzz8vt2z58uVq2rSprrrqKi9XBQAAAAB1RyisgYsvvlhXXnmlFi1apC+++MJ1+3//+1+98847mjFjhuvcQwAAAABoSLh4fQ05nU699957WrRokQoLC9WsWTM1adJE48ePV7du3XxdHgAAAADUCqEQAAAAAEyM4aMAAAAAYGKEQgAAAAAwMUIhAAAAAJgYoRAAAAAATIxQCAAAAAAmRigEAAAAABMjFAIAAACAiREKAQAAAMDECIUAAAAAYGKEQgAAAAAwMUIhAAAAAJgYoRAAAAAATIxQCAAAAAAmRigEADRIhYWFSklJ8XUZQIVqun1mZ2fr6NGjXqgIKMa+ExUJ9HUBKOuzzz7TokWLlJ2dLYvFouPHj6tPnz4aN26cEhISyrTdsWOH5syZo+TkZBUUFKigoEADBgzQbbfdpsjISB89AgDwrJSUFC1cuFAffvihrr76ak2ZMqXCdgcPHtTs2bO1Y8cOOZ1OHTt2TH379tW4ceMUFxdXrn1mZqZeeukl/frrrzIMQxkZGTr33HM1btw4nXbaaZ5+WGgkarp97tmzRx988IEWL16sKVOm6Prrr6/0Pvm8R32oyba5e/duvfXWW1q3bp3S0tIUHh6url276vrrr9ff/va3Cu+XfWfjQCj0E4Zh6B//+Id+/fVXzZkzR506dZIkJScna/z48Ro2bJgWL16s5s2bS5I2btyosWPH6tlnn1X//v0lSevXr9e4ceP0/fffa+HChQoLC/PZ40Hj4nA49OGHH+qzzz6TYRjKz89XWFiYbrjhBl122WXl2rv7ZRyoTnZ2tlasWKEvvvhChw4dUmpqqo4dO1Zp+3379um6667TnXfeqenTp8tisWj79u269dZbtXz58jL7U0lKT0/Xtddeq8GDB2vhwoUKDAzUwYMHddttt+maa67RwoUL1bFjR288VDRANd0+09PT9e233+qzzz5TQUGBkpKSZLfbq7xvPu9RF+7sO3/66SeNHz9eeXl5io2NVfPmzZWSkqK1a9dq7dq1Gjt2bLkgyb6zETHgFz799FMjMTHRWL58ebllGzZsMBITE43Zs2cbhmEYx48fN3r37m08/vjj5dq++uqrRmJiojF9+nSP1wxzyM/PN26++Wbj73//u3HgwAHX7d99951x+umnG9OmTSvTfu/evUavXr2MBQsWGE6n0zAMw9i2bZtxwQUXGBdccIGRlpbm1frROCxcuNC45ZZbjLVr1xqGYRg33nijkZiYaMyYMaNcW4fDYQwaNMi4/fbbyy37/PPPjcTERGPixIllbh8zZoxx+eWXl2u/ceNGIzEx0bjqqqvq6ZGgMarp9jlr1izj7rvvNjZt2mQYhmFccsklRmJiovH+++9XeL983qOu3Nl3PvPMM8bll19ubN261XXb4cOHjQkTJhiJiYlGYmKisXnz5jJ92Hc2HpxT6CeWL18uSWrdunW5Zd26dZMk1zkHH330kdLT0ys8QnPllVdKkj799FM5HA5PlQsTmTFjhtatW6fZs2fr1FNPdd3er18/jRo1SgsWLNCyZcskSU6nUxMmTFCPHj00YsQIWSwWSVKnTp00depUpaWladq0aT55HGjYhg8frjfeeEPnn39+tW2/++477d69u8J95KWXXqqQkBCtXLlSGRkZkqStW7dqzZo1Fbbv2bOnWrVqpa1bt2rbtm11fyBolGq6fd59992aNWuWunfvXqP75fMedeXOvjMjI0NTpkxRly5dXLc1a9ZMzz//vE455RRJ//u+KrHvbGwIhX6iqKhIUvGh+5MlJydLkvr06SNJ+uabbyRJiYmJ5dq2aNFCcXFxysjI0J9//umpcmESWVlZ+s9//qPu3burVatW5ZZfe+21kqT58+dLcv/LOOAJVe0jg4KC1KVLFxUVFSkpKana9pJ05plnSpI2b97siXKBSvF5D2+yWq0Vhsfg4GD17dtXUvH3ghLsOxsXQqGfOO+88yRJM2fO1IIFC+R0OiUVn2v44osvasiQIa5zCXbu3KnQ0FBFRERUeF8tW7aUJO3fv98LlaMx++2332S3213b1MlatWqlkJAQbdq0SVlZWW5/GQc8YceOHZJU5pzB0k7eR5a0b9asWYXtS46Qs0+Ft/F5D2+aPn26bDZbhctKRv60adPGdRv7zsaFUOgnrrvuOl144YUqKCjQ448/riuvvFKffPKJ7rnnHvXo0UNPP/20JCk/P19ZWVmVfkBI/3tzZmdne6V2NF6ZmZmSiqevrkxoaKicTqfS0tLc/jIOeMLhw4clqdL95Mn7yOraN23atEx7wBv4vIc/2bp1qwIDAzVo0CDXbew7GxdCoZ8ICgrS3LlzXbM67dy5U1OnTtXy5cu1Z88epaamSvrfYfvAwMonjg0ODpakamc0A6oTHR0tSdqyZUuFy3NyclxDQbOystz+Mg54QnX7yZP3kSXtK/uFnH0qfIHPe/iLpKQk/fbbbxo1apTi4+Ndt7PvbFwIhX5k3bp1WrZsmebMmaOPP/5Yl112mRwOhxYuXKjBgwfrl19+cZ17WPJvRazW4peVKapRV2eddZZCQ0OVnJysVatWlVu+aNEi199NmjRx+8s44AnV7SdP3kfWtH1oaGi91glUhc97+AO73a5//vOfOuecczR58uQyy9h3Ni6EQj+xatUqjR8/Xg8//LD69++v008/Xc8//7yWLFmiLl26KDs7W1OmTFFQUJCksif6nqygoEBS5UdrgJqKjIzUhAkTJEn33XefvvzyS+Xm5uro0aN67733tHLlSlfb5s2bu/1lHPCEkv1kZdfiOnkfWV37kuHTXCQc3sTnPfzBI488IqfTqZdffrncEUH2nY0LodAPFBUVadq0aerdu7fOOuusMss6d+6s+fPnq0mTJjp06JAOHjyogIAAFRQUKD8/v8L7y8nJkVT2ZGCgtm677TY98sgjstlsmjx5snr27KnLLrtMdrvdNdy5Xbt2ioiIcPvLOOAJJeexVPZl+uR9ZEn7yrbbkuHO7FPhTdHR0Xzew6dmzJihXbt2ad68eRUGO/adjQuh0A/s27dPycnJ6tixY4XLo6Oj1atXL0nFR1hKrmWYkpJSYftDhw4pJiZGLVq08EzBMJ0bbrhBa9as0XfffadVq1bpxx9/1M0336xffvlFkjR48GBJ7n8ZBzyhXbt2kqreRwYGBrpmyS1pX3LudkXtJZW5dhfgaTabjc97+Mzs2bO1adMmvf3224qKiqqwDfvOxoVQ6AdKhtrl5uZW2iYnJ0fNmzdXmzZtXNeQ2bRpU7l22dnZ2r9/vy6//HLX9MFAfbBYLDr11FOVkJAgq9Uqp9OpDz/8UCEhIRo6dKgk97+MA55Q1T5SKp5F78ILL3RNpFRd+y1btigxMVGdOnXyQLVA5fi8hy/Mnj1bv/32m954440qR/aw72xcCIV+oH379oqLi9OqVasqnPo/LS1N69ev14QJE2Sz2XTdddfJarXq888/L9d25cqVCgwM1A033OCN0mFi8+bN0549e3TXXXe5fqV298s4UFsOh6PSZVdddZXCwsK0dOlSGYZRZtnGjRuVnp6uMWPGuG7r16+fWrRooW+//VZ5eXll2h88eFBbt27VLbfcUr8PAI1aVdunO+34vEd9q26bKzlC+MorrygkJKTccsMwXAcz2Hc2LoRCP2Cz2fTkk08qLS1N9957r9LT013L9uzZowkTJuimm27SiBEjJBVfGPz222/X6tWr9eabb7oudH/gwAH9+9//1iOPPOI6YgN4wk8//aRZs2ZpyJAhZb5cu/tlHKgNwzB05MgRSSqzvyzRtGlTPfjgg9q1a5emT5/u+rEtPT1dTz75pCZOnKhzzz3X1T4oKEiPP/64srOzdf/997uGOefk5Oixxx7TlVdeqauuusrzDwyNQnXbZ4mCggLXUPvK2vF5j/pU3bb5wgsv6KWXXtLmzZt1ySWXqE+fPq7/zj//fPXs2VNdunTRK6+8Iol9Z2NjMU7+5gaf2blzp+bOnauNGzcqJiZGkZGRat68uUaMGFFuAhpJ+vTTT/XOO+/o2LFjiouLU3h4uMaMGaPzzjvPB9XDDJxOpz744AM988wzGj16tO65555yw5YWLVqkhx9+WCNHjtR9992noKAgpaen69Zbb1W/fv10xx13+Kh6NHRz5szRypUrlZ6eruTkZNftHTt2VGhoqKZPn17m3OxVq1bp9ddf18GDBxUfH6/g4GBde+21ZS6+XNqmTZv00ksvaefOnUpISJDNZtPgwYN13XXXMTwP1arp9vnYY49p8+bNSk1NdV3bNSAgQB06dFBISIhee+01xcTElLlvPu9RFzXZNvfu3atJkybV6P7GjRtX5vIU7DsbB0IhgGplZmZqxYoVeuONNxQbG6vJkyfrL3/5S6Xt3f0yDgAA/F9hYaEsFkulF6xHw0UoBFClQ4cOafr06erQoYMuvPDCCo9aAwAAoOEiFAIAAACAiTHRDAAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAPzE119/rbPPPlv33HOPr0sBAJgIoRAAAD+xePFiHT9+XEuXLlV6erqvywEAmAShEAAAPzFkyBBFRkbqsssuU2xsrK/LAQCYhMUwDMPXRQAAAAAAfIMjhQAAAABgYoRCAAD8DIN4AADeFOjrAgAAgHTo0CElJSVp5cqVWr58udatW1dhu+zsbM2bN0/ffvut9u/fr5CQELVp00YXXXSRmjVrpi+++ELvvPOOl6sHADRkhEIAAHzsk08+0SuvvKLk5GTZ7fYK2zidTv3nP//RV199pTvvvFMTJ06U0+nUunXr9NJLL+mFF16QJPXv39+bpQMAGgEmmgEAwE+88sor+ve//y1J2r59u+v2oqIi3XnnnUpKStInn3yi6OjoMv2cTqcuv/xy7d69WxMnTtSdd97pzbIBAA0c5xQCAOAnmjVrVuHts2fP1nfffadx48aVC4SSZLVaFRISIkkVLgcAoCqEQgAA/ERAQEC52+x2u95//31JUp8+faq9j5JwCABATREKAQDwY0lJScrKypIkxcfHV9u+omAJAEBVCIUAAPixP//80/V3fn5+te2tVj7aAQDu4ZMDAIAGYt++fb4uAQDQCBEKAQDwY61bt3b9vXLlSh9WAgBorAiFAAD4sS5durjOJXzvvfeUnp5erk1WVpaOHDni7dIAAI0EoRAAAD8WEBCgW2+9VZJ07NgxjRo1Sps3b3YtT0pK0uTJk3X8+HFflQgAaOACfV0AAAAodvjwYdffqampriOEN9xwg9auXavvv/9eO3bs0NChQxUbG6uAgACdfvrpmjFjhoYOHarc3FxflQ4AaMA4UggAgI8tWLBAvXr10qxZs1y3XXrpperVq5eys7MVEBCgOXPm6B//+IcSExMVHBwsp9OpG2+8Ua+88opiY2N9WD0AoKGzGIZh+LoIAABQN/369VNycrKefvppDRkyxNflAAAaEI4UAgAAAICJEQoBAAAAwMQIhQAAAABgYoRCAAAAADAxQiEAAI0A88YBAGqLUAgAQAOXnZ2t9PR0SdKhQ4d8XA0AoKHhkhQAADRQ27Zt06hRo5Sfn6/8/HxJksViUXR0tMaNG6ebb77ZxxUCABoCQiEAAAAAmBjDRwEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMDFCIQAAAACYGKEQAAAAAEyMUAgAAAAAJkYoBAAAAAATIxQCAAAAgIkRCgEAAADAxAiFAAAAAGBihEIAAAAAMLH/B9zniPWCQXHQAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 293, "width": 450 } }, "output_type": "display_data" } ], "source": [ "sns.histplot(data=d, x=\"iq\", hue=\"group\");" ] }, { "cell_type": "code", "execution_count": 26, "id": "14ad253d-c7d7-4350-919a-0fc810de02c6", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
iq
group
control100.357143
treatment101.914894
\n", "
" ], "text/plain": [ " iq\n", "group \n", "control 100.357143\n", "treatment 101.914894" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d.groupby(\"group\").mean()" ] }, { "cell_type": "markdown", "id": "2e01bedf-4f0d-4512-b7e2-f5b3bbb2f175", "metadata": {}, "source": [ "### Equal variances t-test" ] }, { "cell_type": "code", "execution_count": 27, "id": "a765681a-981b-4a9c-8311-db747843235b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1.5586953301521096, 0.12269895509665575)" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.stats import ttest_ind\n", "\n", "res_eqvar = ttest_ind(treated, controls, equal_var=True)\n", "res_eqvar.statistic, res_eqvar.pvalue" ] }, { "cell_type": "code", "execution_count": 28, "id": "588bb63c-7154-4ec9-bd34-f2285358f012", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-0.42865302979133335, 3.5441545495481668]" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ci_eqvar = res_eqvar.confidence_interval(confidence_level=0.95)\n", "[ci_eqvar.low, ci_eqvar.high]" ] }, { "cell_type": "code", "execution_count": null, "id": "d9ca14d2-9d26-4423-ad12-72bfe99b4490", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "460c9ad9-81b6-4d90-97ae-295bd445dcea", "metadata": {}, "source": [ "### Unequal variances t-test" ] }, { "cell_type": "code", "execution_count": 29, "id": "a0670103-e5bf-4dea-a6d3-c4998ca63a03", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1.622190457290228, 0.10975381983712831)" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res_uneqvar = ttest_ind(treated, controls, equal_var=False)\n", "res_uneqvar.statistic, res_uneqvar.pvalue" ] }, { "cell_type": "code", "execution_count": 30, "id": "76d367ad-94a2-4353-8115-1d3e9d015d27", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-0.3611847716497789, 3.476686291406612]" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ci_uneqvar = res_uneqvar.confidence_interval(confidence_level=0.95)\n", "[ci_uneqvar.low, ci_uneqvar.high]" ] }, { "cell_type": "code", "execution_count": null, "id": "19a7e676-a01c-4889-b861-36ba87d3e7a6", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "2e5350f9-f505-4452-a779-c9e029b6d60a", "metadata": {}, "source": [ "### Linear model" ] }, { "cell_type": "code", "execution_count": 31, "id": "0df340d4-887e-49ab-bf95-fb959f276884", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: iq R-squared: 0.027\n", "Model: OLS Adj. R-squared: 0.016\n", "Method: Least Squares F-statistic: 2.430\n", "Date: Mon, 19 Aug 2024 Prob (F-statistic): 0.123\n", "Time: 22:59:33 Log-Likelihood: -263.13\n", "No. Observations: 89 AIC: 530.3\n", "Df Residuals: 87 BIC: 535.2\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "=========================================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "-----------------------------------------------------------------------------------------\n", "Intercept 100.3571 0.726 138.184 0.000 98.914 101.801\n", "C(group)[T.treatment] 1.5578 0.999 1.559 0.123 -0.429 3.544\n", "==============================================================================\n", "Omnibus: 46.068 Durbin-Watson: 2.025\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 553.494\n", "Skew: 1.108 Prob(JB): 6.46e-121\n", "Kurtosis: 15.014 Cond. No. 2.69\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "import statsmodels.formula.api as smf\n", "res_ols = smf.ols(\"iq ~ 1 + C(group)\", data=d).fit()\n", "print(res_ols.summary())" ] }, { "cell_type": "markdown", "id": "96f50bd1-8d26-433b-ac3f-1b2c56feea05", "metadata": {}, "source": [ "### Alternative linear model\n", "\n", "Using generalized least squares to reproduce the unequal variance case." ] }, { "cell_type": "code", "execution_count": 32, "id": "2a75cf9c-9ba0-4b1d-a3b8-f9823a120342", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " GLS Regression Results \n", "==============================================================================\n", "Dep. Variable: iq R-squared: 0.029\n", "Model: GLS Adj. R-squared: 0.018\n", "Method: Least Squares F-statistic: 2.629\n", "Date: Mon, 19 Aug 2024 Prob (F-statistic): 0.109\n", "Time: 22:59:33 Log-Likelihood: -248.41\n", "No. Observations: 89 AIC: 500.8\n", "Df Residuals: 87 BIC: 505.8\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "=========================================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "-----------------------------------------------------------------------------------------\n", "Intercept 100.3571 0.388 258.627 0.000 99.586 101.128\n", "C(group)[T.treatment] 1.5578 0.961 1.622 0.109 -0.352 3.467\n", "==============================================================================\n", "Omnibus: 33.582 Durbin-Watson: 2.234\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 346.198\n", "Skew: -0.648 Prob(JB): 6.67e-76\n", "Kurtosis: 12.575 Cond. No. 2.79\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "n_t, var_t = len(treated), treated.var()\n", "n_c, var_c = len(controls), controls.var()\n", "sigma2s = [var_t]*n_t + [var_c]*n_c\n", "\n", "res_gls = smf.gls(\"iq ~ 1 + C(group)\", data=d, sigma=sigma2s).fit()\n", "print(res_gls.summary())" ] }, { "cell_type": "code", "execution_count": null, "id": "62fc825a-650c-4176-9e9f-05c1c43e930a", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "e17037dd-55c0-4d61-b01c-8a30f6aa3501", "metadata": {}, "source": [ "## Bayesian equal variances model" ] }, { "cell_type": "code", "execution_count": 33, "id": "97ca996a-3ff8-4489-b193-060190419e5e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Formula: iq ~ 1 + group\n", " Family: gaussian\n", " Link: mu = identity\n", " Observations: 89\n", " Priors: \n", " target = mu\n", " Common-level effects\n", " Intercept ~ Normal(mu: 101.1798, sigma: 17.17)\n", " group ~ Normal(mu: 0.0, sigma: 23.6275)\n", " \n", " Auxiliary parameters\n", " sigma ~ HalfStudentT(nu: 4.0, sigma: 4.718)\n" ] } ], "source": [ "import bambi as bmb\n", "\n", "mod_eqvar = bmb.Model(\"iq ~ 1 + group\", data=d)\n", "print(mod_eqvar)" ] }, { "cell_type": "code", "execution_count": 34, "id": "7a7b6ec2-d87b-4f47-a6bb-00f1e4c99ae8", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [iq_sigma, Intercept, group]\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [12000/12000 00:08<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 20 seconds.\n" ] } ], "source": [ "idata_eqvar = mod_eqvar.fit(draws=2000)" ] }, { "cell_type": "code", "execution_count": 35, "id": "2c566027-67c4-41d7-8eff-b8a3e09909ec", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_3%hdi_97%
Intercept100.3580.72798.973101.722
group[treatment]1.5541.002-0.3583.365
iq_sigma4.7420.3604.0795.411
\n", "
" ], "text/plain": [ " mean sd hdi_3% hdi_97%\n", "Intercept 100.358 0.727 98.973 101.722\n", "group[treatment] 1.554 1.002 -0.358 3.365\n", "iq_sigma 4.742 0.360 4.079 5.411" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import arviz as az\n", "\n", "az.summary(idata_eqvar, kind=\"stats\")" ] }, { "cell_type": "code", "execution_count": null, "id": "f05fac0c-c19f-46f4-b308-f475e0b2f2c1", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "ed3adf62-60da-4d80-81c4-638c98ecd0c5", "metadata": {}, "source": [ "## Bayesian unequal variances model" ] }, { "cell_type": "code", "execution_count": 36, "id": "65599934-f457-4c75-8450-f60bc7743523", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Formula: iq ~ 1 + group\n", " sigma ~ group\n", " Family: gaussian\n", " Link: mu = identity\n", " sigma = log\n", " Observations: 89\n", " Priors: \n", " target = mu\n", " Common-level effects\n", " Intercept ~ Normal(mu: 101.1798, sigma: 17.17)\n", " group ~ Normal(mu: 0.0, sigma: 23.6275)\n", " target = sigma\n", " Common-level effects\n", " sigma_Intercept ~ Normal(mu: 0.0, sigma: 1.0)\n", " sigma_group ~ Normal(mu: 0.0, sigma: 1.0)\n" ] } ], "source": [ "formula = bmb.Formula(\"iq ~ 1 + group\", \"sigma ~ group\")\n", "mod_uneqvar = bmb.Model(formula, data=d)\n", "print(mod_uneqvar)" ] }, { "cell_type": "code", "execution_count": 37, "id": "c79d67df-82c2-4672-8d0a-3c83af7296cd", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [Intercept, group, sigma_Intercept, sigma_group]\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [12000/12000 00:10<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 22 seconds.\n" ] } ], "source": [ "idata_uneqvar = mod_uneqvar.fit(draws=2000)" ] }, { "cell_type": "code", "execution_count": 38, "id": "d3e409e7-72b5-42bf-bfe1-3b65f7721116", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_3%hdi_97%
Intercept100.3570.39799.617101.105
group[treatment]1.5530.965-0.1773.434
sigma_Intercept0.9390.1090.7411.150
sigma_group[treatment]0.8490.1510.5661.133
sigma[0]6.0130.6374.8817.212
...............
sigma[84]2.5740.2852.0683.121
sigma[85]2.5740.2852.0683.121
sigma[86]2.5740.2852.0683.121
sigma[87]2.5740.2852.0683.121
sigma[88]2.5740.2852.0683.121
\n", "

93 rows × 4 columns

\n", "
" ], "text/plain": [ " mean sd hdi_3% hdi_97%\n", "Intercept 100.357 0.397 99.617 101.105\n", "group[treatment] 1.553 0.965 -0.177 3.434\n", "sigma_Intercept 0.939 0.109 0.741 1.150\n", "sigma_group[treatment] 0.849 0.151 0.566 1.133\n", "sigma[0] 6.013 0.637 4.881 7.212\n", "... ... ... ... ...\n", "sigma[84] 2.574 0.285 2.068 3.121\n", "sigma[85] 2.574 0.285 2.068 3.121\n", "sigma[86] 2.574 0.285 2.068 3.121\n", "sigma[87] 2.574 0.285 2.068 3.121\n", "sigma[88] 2.574 0.285 2.068 3.121\n", "\n", "[93 rows x 4 columns]" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "az.summary(idata_uneqvar, kind=\"stats\")" ] }, { "cell_type": "code", "execution_count": null, "id": "cf380df8-1049-44a8-9298-51bf3cd49258", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "id": "a2e724c9-e853-4d9f-bd02-54473c61d022", "metadata": {}, "source": [ "## Robust Bayesian Estimation" ] }, { "cell_type": "code", "execution_count": 51, "id": "12f66e4a-8e0f-4d2c-aa1e-b0944fbb499c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Formula: iq ~ 1 + group\n", " sigma ~ group\n", " Family: t\n", " Link: mu = identity\n", " sigma = log\n", " Observations: 89\n", " Priors: \n", " target = mu\n", " Common-level effects\n", " Intercept ~ Normal(mu: 101.1798, sigma: 17.17)\n", " group ~ Normal(mu: 0.0, sigma: 23.6275)\n", " \n", " Auxiliary parameters\n", " nu ~ Gamma(alpha: 2.0, beta: 0.1)\n", " target = sigma\n", " Common-level effects\n", " sigma_Intercept ~ Normal(mu: 0.0, sigma: 1.0)\n", " sigma_group ~ Normal(mu: 0.0, sigma: 1.0)\n" ] } ], "source": [ "formula = bmb.Formula(\"iq ~ 1 + group\", \"sigma ~ group\")\n", "mod_robust = bmb.Model(formula, family=\"t\", data=d)\n", "print(mod_robust)" ] }, { "cell_type": "code", "execution_count": 41, "id": "c64e2037-ba8b-4bb5-87c6-1e0bf5dffb37", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [iq_nu, Intercept, group, sigma_Intercept, sigma_group]\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n", "WARNING (pytensor.tensor.blas): Using NumPy C-API based implementation for BLAS functions.\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [8000/8000 00:08<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 20 seconds.\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
meansdhdi_3%hdi_97%
Intercept100.5250.211100.142100.931
group[treatment]1.0260.4140.2481.799
sigma_Intercept0.0090.196-0.3640.361
sigma_group[treatment]0.6240.2520.1561.090
iq_nu1.8270.4871.0902.768
...............
sigma[84]1.0290.2040.6771.416
sigma[85]1.0290.2040.6771.416
sigma[86]1.0290.2040.6771.416
sigma[87]1.0290.2040.6771.416
sigma[88]1.0290.2040.6771.416
\n", "

94 rows × 4 columns

\n", "
" ], "text/plain": [ " mean sd hdi_3% hdi_97%\n", "Intercept 100.525 0.211 100.142 100.931\n", "group[treatment] 1.026 0.414 0.248 1.799\n", "sigma_Intercept 0.009 0.196 -0.364 0.361\n", "sigma_group[treatment] 0.624 0.252 0.156 1.090\n", "iq_nu 1.827 0.487 1.090 2.768\n", "... ... ... ... ...\n", "sigma[84] 1.029 0.204 0.677 1.416\n", "sigma[85] 1.029 0.204 0.677 1.416\n", "sigma[86] 1.029 0.204 0.677 1.416\n", "sigma[87] 1.029 0.204 0.677 1.416\n", "sigma[88] 1.029 0.204 0.677 1.416\n", "\n", "[94 rows x 4 columns]" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "idata_robust = mod_robust.fit(draws=1000)\n", "az.summary(idata_robust, kind=\"stats\")" ] }, { "cell_type": "code", "execution_count": null, "id": "d8f766e1-b53a-42bc-bb1f-f01825179219", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.4" } }, "nbformat": 4, "nbformat_minor": 5 }