From 31c3b1734d6dc13107b3ed804290bd41223ff407 Mon Sep 17 00:00:00 2001 From: julianlheureux Date: Fri, 29 Mar 2024 21:48:59 -0300 Subject: [PATCH 1/9] Update of the examples section of README.md file and homepage from index.qmd --- .gitignore | 1 + README.md | 49 ++++++++++++++++++++++++++++---------- docs/index.qmd | 64 ++++++++++++++++++++++++++++++++++++++++++-------- 3 files changed, 91 insertions(+), 23 deletions(-) diff --git a/.gitignore b/.gitignore index 7744b9e87..63b0ce404 100644 --- a/.gitignore +++ b/.gitignore @@ -19,3 +19,4 @@ docs/_build .vscode/ pytest.ini /.quarto/ +.Rproj.user diff --git a/README.md b/README.md index 90d07424f..878408fe8 100644 --- a/README.md +++ b/README.md @@ -32,45 +32,68 @@ Bambi requires working versions of ArviZ, formulae, NumPy, pandas and PyMC. Depe In the following two examples we assume the following basic setup ```python +import arviz as az import bambi as bmb import numpy as np import pandas as pd - -data = pd.DataFrame({ - "y": np.random.normal(size=50), - "g": np.random.choice(["Yes", "No"], size=50), - "x1": np.random.normal(size=50), - "x2": np.random.normal(size=50) -}) ``` ### Linear regression +A simple fixed effect model is shown in the example below. + ```python -model = bmb.Model("y ~ x1 + x2", data) -fitted = model.fit() +#### Read in a database from the package content +data = bmb.load_data("sleepstudy") + +# Initialize the fixed effect only model +model = bmb.Model('Reaction ~ Days', data) +print(model) # Get model description + +# Fit the model using 1000 on each chain +results = model.fit(draws=1000) + +# Key summary and diagnostic info on the model parameters +az.summary(results) + +# Use ArviZ to plot the results +az.plot_trace(results) ``` -In the first line we create and build a Bambi `Model`. The second line tells the sampler to start +First, we create and build a Bambi `Model`. Then, the function `model.fit` tells the sampler to start running and it returns an `InferenceData` object, which can be passed to several ArviZ functions such as `az.summary()` to get a summary of the parameters distribution and sample diagnostics or - `az.plot_trace()` to visualize them. - +`az.plot_trace()` to visualize them. ### Logistic regression +In this example we will use a simulated dataset created as shown below. + +```python +data = pd.DataFrame({ + "g": np.random.choice(["Yes", "No"], size=50), + "x1": np.random.normal(size=50), + "x2": np.random.normal(size=50) +}) +``` + Here we just add the `family` argument set to `"bernoulli"` to tell Bambi we are modelling a binary response. By default, it uses a logit link. We can also use some syntax sugar to specify which event we want to model. We just say `g['Yes']` and Bambi will understand we want to model the probability of a `"Yes"` response. But this notation is not mandatory. If we use `"g ~ x1 + x2"`, Bambi will pick one of the events to model and will inform us which one it picked. - ```python model = bmb.Model("g['Yes'] ~ x1 + x2", data, family="bernoulli") fitted = model.fit() ``` +After this, we can evaluate the model as before. + +### More + +There are many additional examples in our [Examples](https://bambinos.github.io/bambi/notebooks/) webpage. + ## Documentation The Bambi documentation can be found in the [official docs](https://bambinos.github.io/bambi/index.html) diff --git a/docs/index.qmd b/docs/index.qmd index 1499cbba8..9eb76054c 100644 --- a/docs/index.qmd +++ b/docs/index.qmd @@ -58,31 +58,75 @@ If you use Conda, you can also install the latest release of Bambi with the foll conda install -c conda-forge bambi ``` -## Usage +## Examples of usage -A simple fixed effects model is shown in the example below. +## Example + +In the following two examples we assume the following basic setup ```python import arviz as az import bambi as bmb +import numpy as np import pandas as pd +``` + +### Linear regression -# Read in a tab-delimited file containing our data -data = pd.read_table('my_data.txt', sep='\t') +A simple fixed effect model is shown in the example below. -# Initialize the fixed effects only model -model = bmb.Model('DV ~ IV1 + IV2', data) +```python +#### Read in a database from the package content +data = bmb.load_data("sleepstudy") -# Fit the model using 1000 on each of 4 chains -results = model.fit(draws=1000, chains=4) +# Initialize the fixed effect only model +model = bmb.Model('Reaction ~ Days', data) +print(model) # Get model description -# Use ArviZ to plot the results -az.plot_trace(results) +# Fit the model using 1000 on each chain +results = model.fit(draws=1000) # Key summary and diagnostic info on the model parameters az.summary(results) + +# Use ArviZ to plot the results +az.plot_trace(results) +``` + +First, we create and build a Bambi `Model`. Then, the function `model.fit` tells the sampler to start +running and it returns an `InferenceData` object, which can be passed to several ArviZ functions +such as `az.summary()` to get a summary of the parameters distribution and sample diagnostics or +`az.plot_trace()` to visualize them. + +### Logistic regression + +In this example we will use a simulated dataset created as shown below. + +```python +data = pd.DataFrame({ + "g": np.random.choice(["Yes", "No"], size=50), + "x1": np.random.normal(size=50), + "x2": np.random.normal(size=50) +}) ``` +Here we just add the `family` argument set to `"bernoulli"` to tell Bambi we are modelling a binary +response. By default, it uses a logit link. We can also use some syntax sugar to specify which event +we want to model. We just say `g['Yes']` and Bambi will understand we want to model the probability +of a `"Yes"` response. But this notation is not mandatory. If we use `"g ~ x1 + x2"`, Bambi will +pick one of the events to model and will inform us which one it picked. + +```python +model = bmb.Model("g['Yes'] ~ x1 + x2", data, family="bernoulli") +fitted = model.fit() +``` + +After this, we can evaluate the model as before. + +### More + +There are many additional examples in our [Examples](https://bambinos.github.io/bambi/notebooks/) webpage. + For a more in-depth introduction to Bambi see our [Quickstart](https://github.com/bambinos/bambi#quickstart) or our set of example notebooks. From c132acd8bb59fe0adddef98765bab733a3d435ef Mon Sep 17 00:00:00 2001 From: julianlheureux <163574230+julianlheureux@users.noreply.github.com> Date: Fri, 29 Mar 2024 21:57:09 -0300 Subject: [PATCH 2/9] Update index.qmd --- docs/index.qmd | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/index.qmd b/docs/index.qmd index 9eb76054c..780523c57 100644 --- a/docs/index.qmd +++ b/docs/index.qmd @@ -60,8 +60,6 @@ conda install -c conda-forge bambi ## Examples of usage -## Example - In the following two examples we assume the following basic setup ```python From f7e4a8fb76fe7d9acac1930b1d95f3bf3c39aa7c Mon Sep 17 00:00:00 2001 From: julianlheureux Date: Fri, 29 Mar 2024 22:46:33 -0300 Subject: [PATCH 3/9] Update index.qmd --- docs/index.qmd | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/index.qmd b/docs/index.qmd index 9eb76054c..780523c57 100644 --- a/docs/index.qmd +++ b/docs/index.qmd @@ -60,8 +60,6 @@ conda install -c conda-forge bambi ## Examples of usage -## Example - In the following two examples we assume the following basic setup ```python From 8087a736124f0a890577eda56ca536bf8cc95480 Mon Sep 17 00:00:00 2001 From: julianlheureux Date: Mon, 1 Apr 2024 13:03:52 -0300 Subject: [PATCH 4/9] Update README.md and index.qmd --- README.md | 51 ++++++++++++++++++++++++++++++++++++++--------- docs/index.qmd | 54 +++++++++++++++++++++++++++++++++++++++----------- 2 files changed, 84 insertions(+), 21 deletions(-) diff --git a/README.md b/README.md index 878408fe8..772262c00 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ Bambi is a high-level Bayesian model-building interface written in Python. It's ## Installation -Bambi requires a working Python interpreter (3.9+). We recommend installing Python and key numerical libraries using the [Anaconda Distribution](https://www.anaconda.com/products/individual#Downloads), which has one-click installers available on all major platforms. +Bambi requires a working Python interpreter (3.10+). We recommend installing Python and key numerical libraries using the [Anaconda Distribution](https://www.anaconda.com/products/individual#Downloads), which has one-click installers available on all major platforms. Assuming a standard Python environment is installed on your machine (including pip), Bambi itself can be installed in one line using pip: @@ -27,7 +27,7 @@ Alternatively, if you want the bleeding edge version of the package you can inst Bambi requires working versions of ArviZ, formulae, NumPy, pandas and PyMC. Dependencies are listed in `pyproject.toml` and should all be installed by the Bambi installer; no further action should be required. -## Example +## Examples In the following two examples we assume the following basic setup @@ -40,15 +40,20 @@ import pandas as pd ### Linear regression -A simple fixed effect model is shown in the example below. +A simple fixed effects model is shown in the example below. ```python -#### Read in a database from the package content +# Read in a dataset from the package content data = bmb.load_data("sleepstudy") -# Initialize the fixed effect only model +# See first rows +data.head() + +# Initialize the fixed effects only model model = bmb.Model('Reaction ~ Days', data) -print(model) # Get model description + +# Get model description +print(model) # Fit the model using 1000 on each chain results = model.fit(draws=1000) @@ -59,8 +64,36 @@ az.summary(results) # Use ArviZ to plot the results az.plot_trace(results) ``` +``` + Reaction Days Subject +0 249.5600 0 308 +1 258.7047 1 308 +2 250.8006 2 308 +3 321.4398 3 308 +4 356.8519 4 308 +``` +``` + Formula: Reaction ~ Days + Family: gaussian + Link: mu = identity + Observations: 180 + Priors: + target = mu + Common-level effects + Intercept ~ Normal(mu: 298.5079, sigma: 261.0092) + Days ~ Normal(mu: 0.0, sigma: 48.8915) + + Auxiliary parameters + sigma ~ HalfStudentT(nu: 4.0, sigma: 56.1721) +``` +``` + mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat +Intercept 251.552 6.658 238.513 263.417 0.083 0.059 6491.0 2933.0 1.0 +Days 10.437 1.243 8.179 12.793 0.015 0.011 6674.0 3242.0 1.0 +Reaction_sigma 47.949 2.550 43.363 52.704 0.035 0.025 5614.0 2974.0 1.0 +``` -First, we create and build a Bambi `Model`. Then, the function `model.fit` tells the sampler to start +First, we create and build a Bambi `Model`. Then, the method `model.fit()` tells the sampler to start running and it returns an `InferenceData` object, which can be passed to several ArviZ functions such as `az.summary()` to get a summary of the parameters distribution and sample diagnostics or `az.plot_trace()` to visualize them. @@ -92,7 +125,7 @@ After this, we can evaluate the model as before. ### More -There are many additional examples in our [Examples](https://bambinos.github.io/bambi/notebooks/) webpage. +For a more in-depth introduction to Bambi see our [Quickstart](https://github.com/bambinos/bambi#quickstart) and check the notebooks in the [Examples](https://bambinos.github.io/bambi/notebooks/) webpage. ## Documentation @@ -102,7 +135,7 @@ The Bambi documentation can be found in the [official docs](https://bambinos.git If you use Bambi and want to cite it please use -``` +```bibtex @article{Capretto2022, title={Bambi: A Simple Interface for Fitting Bayesian Linear Models in Python}, volume={103}, diff --git a/docs/index.qmd b/docs/index.qmd index 780523c57..b0188beae 100644 --- a/docs/index.qmd +++ b/docs/index.qmd @@ -25,7 +25,7 @@ social sciences and other disciplines. ## Dependencies -Bambi is tested on Python 3.9+ and depends on ArviZ, formulae, NumPy, pandas and PyMC +Bambi is tested on Python 3.10+ and depends on ArviZ, formulae, NumPy, pandas and PyMC (see [pyproject.toml](https://github.com/bambinos/bambi/blob/main/pyproject.toml) for version information). @@ -58,7 +58,7 @@ If you use Conda, you can also install the latest release of Bambi with the foll conda install -c conda-forge bambi ``` -## Examples of usage +## Examples In the following two examples we assume the following basic setup @@ -71,15 +71,20 @@ import pandas as pd ### Linear regression -A simple fixed effect model is shown in the example below. +A simple fixed effects model is shown in the example below. ```python -#### Read in a database from the package content +# Read in a dataset from the package content data = bmb.load_data("sleepstudy") -# Initialize the fixed effect only model +# See first rows +data.head() + +# Initialize the fixed effects only model model = bmb.Model('Reaction ~ Days', data) -print(model) # Get model description + +# Get model description +print(model) # Fit the model using 1000 on each chain results = model.fit(draws=1000) @@ -90,8 +95,36 @@ az.summary(results) # Use ArviZ to plot the results az.plot_trace(results) ``` +``` + Reaction Days Subject +0 249.5600 0 308 +1 258.7047 1 308 +2 250.8006 2 308 +3 321.4398 3 308 +4 356.8519 4 308 +``` +``` + Formula: Reaction ~ Days + Family: gaussian + Link: mu = identity + Observations: 180 + Priors: + target = mu + Common-level effects + Intercept ~ Normal(mu: 298.5079, sigma: 261.0092) + Days ~ Normal(mu: 0.0, sigma: 48.8915) + + Auxiliary parameters + sigma ~ HalfStudentT(nu: 4.0, sigma: 56.1721) +``` +``` + mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat +Intercept 251.552 6.658 238.513 263.417 0.083 0.059 6491.0 2933.0 1.0 +Days 10.437 1.243 8.179 12.793 0.015 0.011 6674.0 3242.0 1.0 +Reaction_sigma 47.949 2.550 43.363 52.704 0.035 0.025 5614.0 2974.0 1.0 +``` -First, we create and build a Bambi `Model`. Then, the function `model.fit` tells the sampler to start +First, we create and build a Bambi `Model`. Then, the method `model.fit()` tells the sampler to start running and it returns an `InferenceData` object, which can be passed to several ArviZ functions such as `az.summary()` to get a summary of the parameters distribution and sample diagnostics or `az.plot_trace()` to visualize them. @@ -123,16 +156,13 @@ After this, we can evaluate the model as before. ### More -There are many additional examples in our [Examples](https://bambinos.github.io/bambi/notebooks/) webpage. - -For a more in-depth introduction to Bambi see our -[Quickstart](https://github.com/bambinos/bambi#quickstart) or our set of example notebooks. +For a more in-depth introduction to Bambi see our [Quickstart](https://github.com/bambinos/bambi#quickstart) and check the notebooks in the [Examples](https://bambinos.github.io/bambi/notebooks/) webpage. ## Citation If you use Bambi and want to cite it please use -```bib +```bibtex @article{ Capretto2022, title={Bambi: A Simple Interface for Fitting Bayesian Linear Models in Python}, From 06cbddcc164b91fc9ceaeeb97faf5ec3e75c8eb9 Mon Sep 17 00:00:00 2001 From: julianlheureux <“lheureuxjulian@gmail.com”> Date: Thu, 22 Aug 2024 21:03:34 -0300 Subject: [PATCH 5/9] Adding horseshoe prior and example --- bambi/backend/utils.py | 9 +- docs/notebooks/horseshoe_prior.ipynb | 648 +++++++++++++++++++++++++++ 2 files changed, 656 insertions(+), 1 deletion(-) create mode 100644 docs/notebooks/horseshoe_prior.ipynb diff --git a/bambi/backend/utils.py b/bambi/backend/utils.py index 6827099d9..e7d78535f 100644 --- a/bambi/backend/utils.py +++ b/bambi/backend/utils.py @@ -5,7 +5,14 @@ import pytensor.tensor as pt import pymc as pm -MAPPING = {"Cumulative": pm.Categorical, "StoppingRatio": pm.Categorical} +def Horseshoe(name, tau_nu = 3, lam_nu = 1, dims=None): + tau = pm.HalfStudentT(f"{name}_tau", nu=tau_nu) + lam = pm.HalfStudentT(f"{name}_lam", nu=lam_nu, dims=dims) + beta_raw = pm.Normal(f"{name}_raw", 0, 1, dims=dims) + beta = pm.Deterministic(name, beta_raw * tau ** 2 * lam ** 2, dims=dims) + return beta + +MAPPING = {"Cumulative": pm.Categorical, "StoppingRatio": pm.Categorical, "Horseshoe": Horseshoe} def get_distribution(dist): diff --git a/docs/notebooks/horseshoe_prior.ipynb b/docs/notebooks/horseshoe_prior.ipynb new file mode 100644 index 000000000..a5869e974 --- /dev/null +++ b/docs/notebooks/horseshoe_prior.ipynb @@ -0,0 +1,648 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Horseshoe Prior\n", + "\n", + "In this example, we will use the Horseshoe Prior (Carvalho et al., 2009) to model a large number of variables, with only a few slopes being significantly different from zero. " + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "import arviz as az\n", + "import bambi as bmb\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pymc as pm\n", + "import pytensor.tensor as pt\n", + "\n", + "from matplotlib import pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We are simulating 100 observations with the 50 explanatory variables, also simulated. " + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "D = 50\n", + "D0 = 5\n", + "\n", + "SEED = 123456789 # for reproducibility\n", + "\n", + "rng = np.random.default_rng(SEED)\n", + "\n", + "INTERCEPT = rng.uniform(-3, 3) # simulate an intercept\n", + "\n", + "COEF = np.zeros(D)\n", + "# Simulate the slopes for significant variables\n", + "COEF[:D0] = rng.choice([-1, 1], size=D0) * rng.normal(5, 1, size=D0)\n", + "\n", + "N = 100\n", + "X = rng.normal(size=(N, D))\n", + "SIGMA = 1.\n", + "# Simulate the data\n", + "y = INTERCEPT + X.dot(COEF) + rng.normal(0, SIGMA, size=N)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we create the dataframe and the term name for the set of variables, to define the formula. " + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "df = pd.DataFrame(X)\n", + "df.columns = [f\"x{i}\" for i in range(X.shape[1])]\n", + "df[\"y\"] = y" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'y ~ c(x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, x39, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49)'" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "term_name = \"c(\" + \", \".join([f\"x{i}\" for i in range(X.shape[1])]) + \")\"\n", + "formula = f\"y ~ {term_name}\"\n", + "formula" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we call the Horseshoe prior and create the model" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "clusterpredictors_dim (50)\n", + "\n", + "predictors_dim (50)\n", + "\n", + "\n", + "cluster__obs__ (100)\n", + "\n", + "__obs__ (100)\n", + "\n", + "\n", + "\n", + "predictors_tau\n", + "\n", + "predictors_tau\n", + "~\n", + "HalfStudentT\n", + "\n", + "\n", + "\n", + "predictors\n", + "\n", + "predictors\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "predictors_tau->predictors\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Intercept\n", + "\n", + "Intercept\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "mu\n", + "\n", + "mu\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "Intercept->mu\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "sigma\n", + "\n", + "sigma\n", + "~\n", + "HalfStudentT\n", + "\n", + "\n", + "\n", + "y\n", + "\n", + "y\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "sigma->y\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "predictors->mu\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "predictors_lam\n", + "\n", + "predictors_lam\n", + "~\n", + "HalfStudentT\n", + "\n", + "\n", + "\n", + "predictors_lam->predictors\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "predictors_raw\n", + "\n", + "predictors_raw\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "predictors_raw->predictors\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "mu->y\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "priors = {\n", + " term_name: bmb.Prior(\"Horseshoe\"),\n", + "}\n", + "model = bmb.Model(formula, df, priors=priors)\n", + "model.set_alias({term_name: \"predictors\"})\n", + "\n", + "model.build()\n", + "model.graph()" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (2 chains in 2 jobs)\n", + "NUTS: [sigma, Intercept, predictors_tau, predictors_lam, predictors_raw]\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "1939127c06ab47b69b9e5051e3500949", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
\n"
+      ],
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/html": [
+       "
\n",
+       "
\n" + ], + "text/plain": [ + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 444 seconds.\n", + "There were 64 divergences after tuning. Increase `target_accept` or reparameterize.\n", + "Chain 0 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.\n", + "Chain 1 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.\n", + "We recommend running at least 4 chains for robust computation of convergence diagnostics\n" + ] + } + ], + "source": [ + "idata = model.fit(target_accept = 0.95, chains=2)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "clusterpredictors_dim (50)\n", + "\n", + "predictors_dim (50)\n", + "\n", + "\n", + "cluster__obs__ (100)\n", + "\n", + "__obs__ (100)\n", + "\n", + "\n", + "\n", + "predictors_tau\n", + "\n", + "predictors_tau\n", + "~\n", + "HalfStudentT\n", + "\n", + "\n", + "\n", + "predictors\n", + "\n", + "predictors\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "predictors_tau->predictors\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Intercept\n", + "\n", + "Intercept\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "mu\n", + "\n", + "mu\n", + "~\n", + "Deterministic\n", + "\n", + "\n", + "\n", + "Intercept->mu\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "sigma\n", + "\n", + "sigma\n", + "~\n", + "HalfStudentT\n", + "\n", + "\n", + "\n", + "y\n", + "\n", + "y\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "sigma->y\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "predictors->mu\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "predictors_lam\n", + "\n", + "predictors_lam\n", + "~\n", + "HalfStudentT\n", + "\n", + "\n", + "\n", + "predictors_lam->predictors\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "predictors_raw\n", + "\n", + "predictors_raw\n", + "~\n", + "Normal\n", + "\n", + "\n", + "\n", + "predictors_raw->predictors\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "mu->y\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "priors = {\n", + " term_name: bmb.Prior(\"Horseshoe\", tau_nu = 3, lam_nu = 3),\n", + "}\n", + "model = bmb.Model(formula, df, priors=priors)\n", + "model.set_alias({term_name: \"predictors\"})\n", + "\n", + "model.build()\n", + "model.graph()" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (2 chains in 2 jobs)\n", + "NUTS: [sigma, Intercept, predictors_tau, predictors_lam, predictors_raw]\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "b0b4f67884bd4607a15224df17c185f8", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
\n"
+      ],
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/html": [
+       "
\n",
+       "
\n" + ], + "text/plain": [ + "\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 293 seconds.\n", + "There were 29 divergences after tuning. Increase `target_accept` or reparameterize.\n", + "We recommend running at least 4 chains for robust computation of convergence diagnostics\n" + ] + } + ], + "source": [ + "idata = model.fit(target_accept = 0.97, chains=2)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqcAAAIlCAYAAADonVHXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABawklEQVR4nO3deXhU5eH28fvMTLZJZiYzkz2BsImIyCK4a0XccS9archP0dpqoXV521p3tO7a1lbr1larrVprccENRWUT2SEIYd/CDmFNyJ7M8/4RkxKSkASSzEnm+7muXF5z5sxz7pyZCbdnznPGMsYYAQAAADbgCHcAAAAAoAblFAAAALZBOQUAAIBtUE4BAABgG5RTAAAA2AblFAAAALZBOQUAAIBtUE4BAABgG5RTAAAA2AblFB3aP/7xD1mWVfvjcrmUlZWl0aNHa/Pmza2+veLiYo0bN05Tpkxp9bElacqUKbIsq83Gb46afbp+/fraZUOHDtXQoUNbNM7SpUs1bty4OuM0x8HbWr9+vSzL0jPPPNOicZry2GOP6YMPPqi33A7PQXM899xz6tWrl6Kjo2VZlvbu3RvuSGiB8vJy3XLLLUpPT5fT6dTAgQMlSbt379Y111yjlJQUWZalyy+/XJJkWZbGjRvXom3UvHf+8Y9/tGr2g7311lt69tln23QbiCyucAcAWsNrr72mPn36qKSkRNOmTdPjjz+uqVOnavHixYqPj2+17RQXF+uhhx6SpBaXteY4/vjjNXPmTPXt27fVxz4SL7zwQosfs3TpUj300EMaOnSounXr1qbbOhyPPfaYrrzyytp//GvY9Tk4UE5Ojn75y1/qJz/5ia6//nq5XC55PJ5wx0ILvPjii3r55Zf13HPPafDgwUpISJAk/e53v9P777+vV199VT179lQgEJAkzZw5U1lZWS3aRnp6umbOnKmePXu2ev4DvfXWW1qyZIluv/32Nt0OIgflFJ1Cv379NGTIEEnSWWedpaqqKv3ud7/TBx98oJEjR4Y5XdMqKipkWZa8Xq9OPvnkVhu3uLhYbrf7iMdpj6JWkzXcpbC1n4O2kJubK0m6+eabdeKJJ4Y5TfjUvG9cro73T9mSJUsUFxensWPH1lves2fPen+3Duc1GRMTY/vXMtAQPtZHp1TzBzkvL0+SVFpaqrvvvlvdu3dXdHS0MjMzNWbMmHofhX799dcaOnSogsGg4uLi1LVrV40YMULFxcVav369kpOTJUkPPfRQ7akEN9xwQ+3jV61apWuvvVYpKSmKiYnRMccco7/85S91tlHzsfE///lP/b//9/+UmZmpmJgYrV69utGPlCdMmKBTTjlFbrdbHo9H5557rmbOnFlnnXHjxsmyLC1YsEBXXnml/H5/k0dMZs2apdNOO02xsbHKyMjQ3XffrYqKinrrNfSx/osvvqgBAwYoISFBHo9Hffr00T333COp+tSAq666SlL1/yzU7KuajxeHDh2qfv36adq0aTr11FPldrt14403NrotSQqFQnr00UfVtWtXxcbGasiQIfrqq6/qrHPDDTc0eJS2Zt/UsCxLRUVFev3112uz1WyzNZ6D3Nxc/fjHP5bP51NqaqpuvPFG7du3r16uhrz66qsaMGCAYmNjFQgEdMUVV2jZsmW19w8dOlTXXXedJOmkk06q9xps7HdvTqbmvk+6deumiy++WBMnTtTxxx+vuLg49enTR6+++mq99Q487ebAnwP375G+b5qz3w5l8+bN+ulPf6ouXbooOjpaGRkZuvLKK7V9+/badTZs2KDrrruuTsbf//73CoVCdcYqLy/XI488oj59+igmJkbJyckaPXq08vPza9exLEt/+9vfVFJSUue9YVmWvvzySy1btqzefmroY/2mcjf2sX5L9vfbb7+te++9VxkZGfJ6vTrnnHO0YsWK2vWGDh2qTz75RHl5eXWe3xqH+jsBNMoAHdhrr71mJJm5c+fWWf6nP/3JSDKvvPKKCYVC5vzzzzcul8vcf//95osvvjDPPPOMiY+PN4MGDTKlpaXGGGPWrVtnYmNjzbnnnms++OADM2XKFPPmm2+aUaNGmT179pjS0lIzceJEI8ncdNNNZubMmWbmzJlm9erVxhhjcnNzjc/nM8cdd5x54403zBdffGH+3//7f8bhcJhx48bVZps8ebKRZDIzM82VV15pJkyYYD7++GOza9eu2vsmT55cu/6bb75pJJnzzjvPfPDBB+add94xgwcPNtHR0Wb69Om16z344INGksnOzjZ33XWXmTRpkvnggw8a3Xe5ubnG7Xabvn37mrffftt8+OGH5vzzzzddu3Y1ksy6detq1z3zzDPNmWeeWXv77bffNpLML37xC/PFF1+YL7/80rz00kvml7/8pTHGmB07dpjHHnvMSDJ/+ctfavfVjh07ascLBAKmS5cu5rnnnjOTJ082U6dObXBb69atM5JMly5dzOmnn27Gjx9v3n33XXPCCSeYqKgo8+2339aue/3115vs7Ox6v2vNvqkxc+ZMExcXZ4YPH16bLTc3t87zcyTPwdFHH20eeOABM2nSJPOHP/zBxMTEmNGjRzf6XNSo2Wc//vGPzSeffGLeeOMN06NHD+Pz+czKlStrn7f77rvPSDKvvfZanddgQ5qbqbnvE2OMyc7ONllZWaZv377mjTfeMJ9//rm56qqrjKTa59EYYxYsWFC7f2fOnGlmzJhhjjvuOBMfH2/WrFlT+/sc6fumOfutMZs2bTLp6ekmKSnJ/OEPfzBffvmleeedd8yNN95oli1bZoypfj1nZmaa5ORk89JLL5mJEyeasWPHGknm1ltvrR2rqqrKXHDBBSY+Pt489NBDZtKkSeZvf/ubyczMNH379jXFxcXGmOrX3/Dhw01cXFztvtm2bZuZOXOmGTRokOnRo0ft8n379hljjJFkHnzwwRblrnnvvPbaa7WPa+n+7tatmxk5cqT55JNPzNtvv226du1qjjrqKFNZWVk73mmnnWbS0tLqPNfGNP13AmgM5RQdWk05nTVrlqmoqDCFhYXm448/NsnJycbj8Zht27bVFsqnnnqqzmPfeeed2gJrjDH//e9/jSSTk5PT6Pby8/Pr/SNR4/zzzzdZWVm1/5jUGDt2rImNjTW7d+82xvzvj/4PfvCDemMcXIyqqqpMRkaGOe6440xVVVXteoWFhSYlJcWceuqptctqSsgDDzxw6J32vauvvtrExcWZbdu21S6rrKw0ffr0abKcjh071iQmJh5y/HfffbdeyTtwPEnmq6++avC+hsppRkaGKSkpqV1eUFBgAoGAOeecc2qXNbecGmNMfHy8uf766+ut2xrPwcGvtZ///OcmNjbWhEKhetursWfPntrCfKANGzaYmJgYc+2119Yua+x/yhrS3EzNfZ8YU11OY2NjTV5eXu2ykpISEwgEzM9+9rNGs4wdO9a4XC7z6aef1i470vdNS/ZbQ2688UYTFRVlli5d2ug6v/3tb40kM3v27DrLb731VmNZllmxYoUx5n9lbPz48XXWmzt3rpFkXnjhhdpl119/vYmPj6+3rTPPPNMce+yx9ZYf/HenObkbKqct3d8H79f//Oc/RlJtATXGmIsuuqjB911z/k4ADeFjfXQKJ598sqKiouTxeHTxxRcrLS1Nn332mVJTU/X1119LUr2PPq+66irFx8fXfjQ8cOBARUdH66c//alef/11rV27ttnbLy0t1VdffaUrrrhCbrdblZWVtT/Dhw9XaWmpZs2aVecxI0aMaHLcFStWaMuWLRo1apQcjv+9XRMSEjRixAjNmjVLxcXFLR5XkiZPnqyzzz5bqamptcucTqeuvvrqJh974oknau/evfrxj3+sDz/8UDt37mzWNg/k9/s1bNiwZq//wx/+ULGxsbW3PR6PLrnkEk2bNk1VVVUt3n5zHc5zcOmll9a53b9/f5WWlmrHjh2NbmfmzJkqKSmp9zrt0qWLhg0bVu8UhpZqKlNz3yc1Bg4cqK5du9bejo2NVe/evWtPpTnYE088oeeff14vvfSSLrzwQkmt87450v322Wef6ayzztIxxxzT6Dpff/21+vbtW+/83htuuEHGmNp99/HHHysxMVGXXHJJnd9l4MCBSktLa9UrQDQn98EOZ3839LqR1OjzfKDW+DuByEQ5RafwxhtvaO7cuVq4cKG2bNmi7777TqeddpokadeuXXK5XLXni9awLEtpaWnatWuXJKlnz5768ssvlZKSojFjxqhnz57q2bOn/vSnPzW5/V27dqmyslLPPfecoqKi6vwMHz5ckur9YU5PT2/WuI2tm5GRoVAopD179rR43Jqx09LS6i1vaNnBRo0apVdffVV5eXkaMWKEUlJSdNJJJ2nSpEnN2nZLch4qV1pamsrLy7V///4WjdUSh/McBIPBOrdjYmIkSSUlJYe9nZr7D1dTmZr7PmlsvJoxG/od//Wvf+mee+7RAw88oJtuuql2eWu8b450v+Xn5zc5C37Xrl2Njn9ghu3bt2vv3r2Kjo6u9/ts27atVctZc3If7HD29+G8lmu0xt8JRKaON8URaMAxxxxTO1v/YMFgUJWVlcrPz6/zD68xRtu2bdMJJ5xQu+yMM87QGWecoaqqKs2bN0/PPfecbr/9dqWmpuqaa65pdPt+v19Op1OjRo3SmDFjGlyne/fudW4fOGmgMTX/MGzdurXefVu2bJHD4ZDf72/xuDVjb9u2rd7yhpY1ZPTo0Ro9erSKioo0bdo0Pfjgg7r44ou1cuVKZWdnN/n45uY8VK5t27YpOjq69jI8sbGxKisrq7fekZSCw3kO2mI7SUlJR7yNprbf3PdJS0yaNEk33nijbrjhhtrLsNVojffNke635ORkbdq06ZDrBIPBRseXVLuNpKQkBYNBTZw4scFxWvNyX83JfbDD2d9H6kj/TiAyceQUnd7ZZ58tqfrozYHGjx+voqKi2vsP5HQ6ddJJJ9XOYF2wYIGkxo8auN1unXXWWVq4cKH69++vIUOG1Ptp6EhTU44++mhlZmbqrbfekjGmdnlRUZHGjx9fO3v8cJx11ln66quv6sxIrqqq0jvvvNOiceLj43XhhRfq3nvvVXl5ee1ljlpyhKU53nvvPZWWltbeLiws1EcffaQzzjhDTqdTUvXs8B07dtT5ncrLy/X555/XG6+xo3wHa8vn4ECnnHKK4uLi6r1ON23apK+//rrB12lrOpz3SVNycnI0YsQIDRs2TK+88kq9+1vjfXOk++3CCy/U5MmT68xAP9jZZ5+tpUuX1v4dqPHGG2/IsiydddZZkqSLL75Yu3btUlVVVYO/y9FHH33ILC3RnNwHa6u/U815LzX2dwJoCEdO0emde+65Ov/883XXXXepoKBAp512mr777js9+OCDGjRokEaNGiVJeumll/T111/roosuUteuXVVaWlp7aZxzzjlHUvWRj+zsbH344Yc6++yzFQgElJSUpG7duulPf/qTTj/9dJ1xxhm69dZb1a1bNxUWFmr16tX66KOPas9LawmHw6GnnnpKI0eO1MUXX6yf/exnKisr09NPP629e/fqiSeeOOz9ct9992nChAkaNmyYHnjgAbndbv3lL39RUVFRk4+9+eabFRcXp9NOO03p6enatm2bHn/8cfl8vtojbP369ZMkvfLKK/J4PIqNjVX37t0P6x8/qfp/GM4991zdeeedCoVCevLJJ1VQUFDnaNzVV1+tBx54QNdcc41+/etfq7S0VH/+858bPCf1uOOO05QpU/TRRx8pPT1dHo+nwfLQls/BgRITE3X//ffrnnvu0f/93//pxz/+sXbt2qWHHnpIsbGxevDBB1tlO41p7vukuQoKCjR8+HDFxcXpV7/6lebNm1fn/r59+8rr9R7x++ZI99vDDz+szz77TD/4wQ90zz336LjjjtPevXs1ceJE3XnnnerTp4/uuOMOvfHGG7rooov08MMPKzs7W5988oleeOEF3Xrrrerdu7ck6ZprrtGbb76p4cOH67bbbtOJJ56oqKgobdq0SZMnT9Zll12mK664okX78UhyN6Qt/k4dd9xxeu+99/Tiiy9q8ODBcjgcGjJkSLP+TgANCut0LOAINXfWcklJibnrrrtMdna2iYqKMunp6ebWW281e/bsqV1n5syZ5oorrjDZ2dkmJibGBINBc+aZZ5oJEybUGevLL780gwYNMjExMUZSnRnf69atMzfeeKPJzMw0UVFRJjk52Zx66qnmkUceqV2nZhbsu+++Wy9nQ5cxMsaYDz74wJx00kkmNjbWxMfHm7PPPtvMmDGjzjo1s7Lz8/Ob2Gv/M2PGDHPyySebmJgYk5aWZn7961+bV155pcnZ+q+//ro566yzTGpqqomOjjYZGRnmRz/6kfnuu+/qjP/ss8+a7t27G6fTWWfWcGMzkhvaVs2M4yeffNI89NBDJisry0RHR5tBgwaZzz//vN7jP/30UzNw4EATFxdnevToYZ5//vkGZ+vn5OSY0047zbjdbiOpdptt8RzUvE4P3KeN+dvf/mb69+9voqOjjc/nM5dddlntZa4OHq8ls/Wbk6k57xNjqmfrX3TRRfW2deBzV/O8NfZz4P490vdNc/dbYzZu3GhuvPFGk5aWZqKiompfz9u3b69dJy8vz1x77bUmGAyaqKgoc/TRR5unn366zhUcjDGmoqLCPPPMM2bAgAEmNjbWJCQkmD59+pif/exnZtWqVbXrHels/ebkbmi2fs3yw93fDY25e/duc+WVV5rExERjWVbte625fyeAg1nGHPA5FQAAABBGnHMKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA26CcAgAAwDY6/EX4jTEqLCwMdwwAAAA0wePxNPn11R2+nBYWFsrn84U7BgAAAJqwb98+eb3eQ67T4S/Cz5FTAACAjqE5R047fDkFAABA58GEKAAAANgG5RQAAAC2QTkFAACAbXT42foAOg9jjCorK1VVVRXuKLAJp9Mpl8vV5AQKAJ0H5RSALZSXl2vr1q0qLi4OdxTYjNvtVnp6uqKjo8MdBUA7YLY+gLALhUJatWqVnE6nkpOTFR0dzZEyyBij8vJy5efnq6qqSkcddZQcDs5GAzo7jpwCCLvy8nKFQiF16dJFbrc73HFgI3FxcYqKilJeXp7Ky8sVGxsb7kgA2hj/CwrANjgqhobwugAiC+94AAAA2AblFAAAALZBOQWATsyyLH3wwQdtuo1u3brp2WefbdNtAIgclFMAaAXffvutnE6nLrjgghY/lnIHAP9DOQXQuYSqpHXTpcX/rf5vqH0u6P/qq6/qF7/4hb755htt2LChXbYJAJ0R5RRA57F0gvRsP+n1i6XxN1X/99l+1cvbUFFRkf7zn//o1ltv1cUXX6x//OMf9daZMGGChgwZotjYWCUlJemHP/yhJGno0KHKy8vTHXfcIcuyaq/vOm7cOA0cOLDOGM8++6y6detWe3vu3Lk699xzlZSUJJ/PpzPPPFMLFixodu6XX35ZmZmZCoVCdZZfeumluv766yVJa9as0WWXXabU1FQlJCTohBNO0JdfftnomOvXr5dlWcrJyaldtnfvXlmWpSlTptQuW7p0qYYPH66EhASlpqZq1KhR2rlzZ7OzA+i8KKcAOoelE6T//J9UsKXu8oKt1cvbsKC+8847Ovroo3X00Ufruuuu02uvvaYDv9/kk08+0Q9/+ENddNFFWrhwob766isNGTJEkvTee+8pKytLDz/8sLZu3aqtW7c2e7uFhYW6/vrrNX36dM2aNUtHHXWUhg8frsLCwmY9/qqrrtLOnTs1efLk2mV79uzR559/rpEjR0qS9u/fr+HDh+vLL7/UwoULdf755+uSSy45oqPDW7du1ZlnnqmBAwdq3rx5mjhxorZv364f/ehHhz0mgM6Di/AD6PhCVdLEuyQ19IV3RpIlTfyt1OciyeFs9c3//e9/13XXXSdJuuCCC7R//3599dVXOueccyRJjz76qK655ho99NBDtY8ZMGCAJCkQCMjpdMrj8SgtLa1F2x02bFid2y+//LL8fr+mTp2qiy++uMnHBwIBXXDBBXrrrbd09tlnS5LeffddBQKB2tsDBgyozSpJjzzyiN5//31NmDBBY8eObVHeGi+++KKOP/54PfbYY7XLXn31VXXp0kUrV65U7969D2tcAJ0DR04BdHx539Y/YlqHkQo2V6/XylasWKE5c+bommuukSS5XC5dffXVevXVV2vXycnJqS17rWnHjh265ZZb1Lt3b/l8Pvl8Pu3fv79FRzVHjhyp8ePHq6ysTJL05ptv6pprrpHTWV3ii4qK9Jvf/EZ9+/ZVYmKiEhIStHz58iM6cjp//nxNnjxZCQkJtT99+vSRVH0aAYDIxpFTAB3f/u2tu14L/P3vf1dlZaUyMzNrlxljFBUVpT179sjv9ysuLq7F4zocjjqnBkhSRUVFnds33HCD8vPz9eyzzyo7O1sxMTE65ZRTVF5e3uztXHLJJQqFQvrkk090wgknaPr06frDH/5Qe/+vf/1rff7553rmmWfUq1cvxcXF6corr2x0GzXf5nRg9oNzh0IhXXLJJXryySfrPT49Pb3Z2QF0TpRTAB1fQmrrrtdMlZWVeuONN/T73/9e5513Xp37RowYoTfffFNjx45V//799dVXX2n06NENjhMdHa2qqrpXFUhOTta2bdtkjKmdJHXgJCNJmj59ul544QUNHz5ckrRx48YWTyqKi4vTD3/4Q7355ptavXq1evfurcGDB9fZxg033KArrrhCUvU5qOvXr290vOTkZEnV55UOGjSowdzHH3+8xo8fr27dusnl4p8hAHXxsT6Aji/7VMmbIclqZAVL8mZWr9eKPv74Y+3Zs0c33XST+vXrV+fnyiuv1N///ndJ0oMPPqi3335bDz74oJYtW6bFixfrqaeeqh2nW7dumjZtmjZv3lxbLocOHar8/Hw99dRTWrNmjf7yl7/os88+q7P9Xr166Z///KeWLVum2bNna+TIkYd1lHbkyJH65JNP9Oqrr9aeO3vgNt577z3l5ORo0aJFuvbaa+vN7j9QXFycTj75ZD3xxBNaunSppk2bpvvuu6/OOmPGjNHu3bv14x//WHPmzNHatWv1xRdf6MYbb6xX0gFEHsopgI7P4ZQuqPmI+OCC+v3tC55o9clQf//733XOOefI5/PVu2/EiBHKycnRggULNHToUL377ruaMGGCBg4cqGHDhmn27Nm16z788MNav369evbsWXvk8ZhjjtELL7ygv/zlLxowYIDmzJmjX/3qV3W28eqrr2rPnj0aNGiQRo0apV/+8pdKSUlp8e8xbNgwBQIBrVixQtdee22d+/74xz/K7/fr1FNP1SWXXKLzzz9fxx9//CHHe/XVV1VRUaEhQ4botttu0yOPPFLn/oyMDM2YMUNVVVU6//zz1a9fP912223y+Xy1pwUAiFyWOfikJgBoZ6WlpVq3bp26d++u2NjYwx9o6YTqWfsHTo7yZlYX076XHnlQhEWrvT4AdAic7AOg8+h7afXlovK+rZ78lJBa/VF+G1w+CgDQNiinADoXh1Pqfka4UwAADhMn9wAAAMA2KKcAAACwDcopANtgfiYawusCiCyUUwBhFxUVJUkqLi4OcxLYUc3rouZ1AqBzY0IUgLBzOp1KTEzUjh07JElut7v2W5EQuYwxKi4u1o4dO5SYmCink6suAJGA65wCsAVjjLZt26a9e/eGOwpsJjExUWlpafwPCxAhKKcAbKWqqkoVFRXhjgGbiIqK4ogpEGEopwAAALANJkQBAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAANCIxx57TOnpGQqmd9FDD/8u3HEigmWMMeEOAQAAYDeTJ0/WOcMvlWfQxXL6klVVsFPvPTFWw844NdzROjWOnAIAADTg7bffVkxWX3Uf8gN165KpqJTueuE/n4c7VqdHOQUAAGjA1GnTFJ3aUz2zu+iE/n1lKsuUs2G3+NC5bVFOAQAADlJWVqZV6zbImRBQt4wUZWd3lSnbrz1l0pqNW8Mdr1OjnAIAABxkw4YNcroTFRWXoKAvQbGxcfInuOWIcWvyrIXhjtepUU4BAAAOkpeXJ0ecR/Een2JcTklSRnqqZDk1L3dVmNN1bpRTAACAg9SUU7fbLcuyJEnpaemSjJat3xLecJ0c5RQAAOAgeXl5csQkyO121y5LT0+XqShV3q6iMCbr/CinAAAAB8nLy5MzPlEJ7rjaZWlpaQqVl2p/paUNW7aHMV3nRjkFAAA4SF5enpxun3wJ/ztyGhMTI198nKyoWE2fuyiM6To3yikAAMBB1m/YKEdMvBJ9CXWWp6cmyXJFa+6SFWFK1vlRTgEAAA5QVVWlrfl7ZLmiFPD56tyXkZ4uSVqyemM4okUEyikAAMABtmzZopAzWg5XtPwHHzlNT5OpKNfa/MIwpev8KKcAAAAHyMvLkxUdJ7fHqxiXq859aWnpCpUXa1+FQ9t27glTws6NcgoAAHCAvLw8OaLdcrvj5XRYde5zu93yxLjkiInXlDk54QnYyVFOAQAADpCXlycrxl3nGqcHSk9NlhwuzV7MN0W1BcopAADAAaovwH+IcpqeLpmQvlu3tZ2TRQbKKQAAwAHWrVsnp9snb0J8g/enp6crVFasDbtLVFkVaud0nR/lFAAA4ABr166VM96vRG9Cg/enp6cpVF6svUVl2rRzXzun6/wopwAAAN+rrKxU3qbNcsTEK5jobXCdhASP4qNdUlSsvl2Y284JOz/KKQAAwPc2bdqkkCtWzpg4BXwNl1NJSktPlWRpfu7K9gsXISinAAAA31u6dKmccV4lJAYUF+1sdL30tHRJRkvXbmq/cBGCcgoAAPC9qVOnyhHnUzApWVHOxmtSenq6TEWp1u0oaMd0kYFyCgAA8L0pU6bIGZ+o5KSkQ66Xnp6mUEWZdu0v0/7i0nZKFxkopwAAAJIKCgo0f/58uQJZykpPPeS6iYmJcllGxhGlBZx32qoopwAAAJK++eYbGVesApndlRbwNbG2paREr6yoGC1YSjltTZRTAAAAVX+k7/KnKzmjqzyxribXT01JliQtXb2hraNFFMopAACApMmTJ8sVyFJScopch5gMVSM5OVmyLK3ayNeYtibKKQAAiHj79u3TgpxFikntqW6Zhz7ftEZKSrJMRbk27i5q43SRhXIKAAAi3jfffCOHJ1m+tC7KTE5s1mNSUlJkKkq0p0wqLCpp24ARhHIKAAAi3rRp0+RKTFNSarriohq/+P6BPB6PoiwjKypWMxcuaeOEkYNyCgAAIt6cOXMUFeyi5GBAlmU181GWMtKSZLliNGV2TlvGiyiUUwAAENGqqqo0b+EiufwZ6pqW3KLHZnfpKlmW5ixZ1UbpIg/lFAAARLTly5er1BGn6HivuqYf+puhDta1a1eZ8lKt2FagUCjURgkjC+UUAABEtDlz5sjpCSqQlKLYqKavb3qgzMxMqbJExYrmm6JaCeUUAABEtDlz5sjlSZLf72/B+abVoqKilJEckCMmXh9++U0bJYwslFMAABDR5syZI5c/Q2lJgcN6fK8ePSRZmrYgt3WDRSjKKQAAiFglJSVavGylnPGJ6pKeclhj9OzZQ6aiVMu3Fqi8vLyVE0YeyikAAIhYOTk5MrFeuT1+pQZ9hzVGenq6olSlqugEffXN7FZOGHkopwAAIGLNnj1bzni//MkpinE17+L7B7Msh7plpcqKduvjKbNaOWHkoZwCAICINWPGjOqZ+oGWXHy/vl49ekiWpRmLlrdiushEOQUAABHJGKPp06crKthFmaktu77pwXr06CFTXqINBSHt2bu3dQJGKMopAACISDk5OcovKFGML0XdMtOOaKzExER5Yp1yxPs1/rPJrZQwMlFOAQBAxAmFQnrggQcUFeyijG495Y+PPeIxe3XrIisqTv/5fHorJIxclFMAABBxnnjiCX382UQl9BysY4/pI6fj8M83rTF40CCZilLNW5uvLVu2tELKyEQ5BQAAEWXu3Lm6//775Yzz6ZRzL1Hvrkf2kX6N1NQ0pSQHZZxRevTRR1tlzEhEOQUAABGjsrJSP/3pTxUKhXTZ5ZfpqF695LAsWaZKWfvm6+j8z5W1b74sU3VY4w8ZfLwk6eWXX1ZOTk4rJo8crnAHAAAAaC8vvfSScnJy5Pf79dC4h/TqvJ06es8UXbDxj/KU76hdrzA6RVN6/D+tDg5r0fhpaenq3bu35kyp0jXXXKN58+YpISGhtX+NTs0WR05feOEFde/eXbGxsRo8eLCmT+dEYgAA0LpycnJ01113SZIeffRRJSUl6biCabpyzd1KOKCYSlJC+Q5dvPwu9dr1dYu3c+655ykzM1MrVqzQzTffLGNMq+SPFGEvp++8845uv/123XvvvVq4cKHOOOMMXXjhhdqwYUO4owEAgA4oFAqpqKiozs+ECRN01llnqbi4WMOGDdPIkSNVvL9Ql2//syTp4OlQNbfPXPt7VZaXqKKivFk/lZWVcrmceu211+R0OvXvf/9bjz76aL08FNbGWSbMe+ekk07S8ccfrxdffLF22THHHKPLL79cjz/+eBiTAQCAjujtt9/Wtdde2+R6w/r49dXVTZ9bevH8k/XNnuZdpN/pCSpUul97p7x2yPU+/PBDXXrppc0aM9KE9ZzT8vJyzZ8/X7/97W/rLD/vvPP07bffhikVAADoyAoKCpq1XlozTwXNCHoU5ezarHUdMW6VlxU3ud7WrVubt/EIFNZyunPnTlVVVSk1NbXO8tTUVG3bti1MqQAAQEc2atQonX766dq3b58KCgrkdDrVp08fBQKBOuuZ9d9I7/6oyfFGj7xGlyad0OztB+Oj1T/zT7W3t2/frvnz56u4uFiJiYlKS0tT3759m/8LRRhbzNa3rLpnehhj6i0DAABoDrfbrWOPPbbpFY85R/JmSAVbJTV0lqMleTN0ytlXSA7nYefp0aOHevTocdiPjzRhnRCVlJQkp9NZ7yjpjh076h1NBQAAaFUOp3TBk9/faGRK1AVPHFExRcuFtZxGR0dr8ODBmjRpUp3lkyZN0qmnnhqmVAAAIGL0vVT60RuSN73ucm9G9fK+TFpqb2Gfrf/OO+9o1KhReumll3TKKafolVde0V//+lfl5uYqOzs7nNEAAECkCFVJed9K+7dLCalS9qkcMQ2TsJ9zevXVV2vXrl16+OGHtXXrVvXr10+ffvopxRQAALQfh1Pqfka4U0A2OHIKAAAA1Aj7N0QBAAAANSinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAHZAxRp3xokuUUwAAgA7mv//9r+K9Pg0+daj27dsX7jitiuucAgAAdCDGGHXr1k27g/0UnZStq45P00tPPBDuWK2GI6cAAAAdyLJly7R5937FpB0lV2Ka3p+9RsXFJeGO1WoopwAAAB3I1KlT5fKlKiWzq2LL9qgqIVmv/vfjcMdqNZRTAACADiQnJ0dOT5KCSUnqd8xRslzR+mja/HDHajWUUwAAgA5kxYoVcvlSFfB61Puo3jLlJfpu416VlVeEO1qroJwCAAB0ICtXrZbTE1Sy36usrCxFmQpVRLn1yeRvwx2tVVBOAQAAOojCwkJt31skR1SsUpL8cjgc6pmdJcsVq4+nzg53vFZBOQUAAOggVq5cKWecR7EJXiUmuCVJPXv0kEyV5q7YGOZ0rYNyCgAA0EGsXLlSjliPPF6fXA5LktS9e3dVlezX1iKjrfm7w5zwyFFOAQAAOoiVK1fKEZegBI9HllVdTn0+n3zuKFkx8Xp/0vQwJzxylFMAAIAOYsWKFXK4E+VJSKizvGf3bMnh1NdzF4cpWeuhnAIAAHQQK1eulNMTVMBXt5z26N5DprJc363bHqZkrYdyCgAA0AEYY7Ry9Rq53D4l+X117svKylSorFi7S6Vd+/aHKWHroJwCAAB0ANu2bVNRhZEVFaOUYGKd+zwej2KdkqJi9e3CJWHJ11oopwAAAB3AihUr5Ij1KMEXUHxMzEH3WkpNDspyRWnudyvCkq+1UE4BAAA6gBUrVsgZlyCP16cop1Xv/rS0VMlIS9dtCkO61kM5BQAA6ACqj5x65TngMlIHCgaCMqFKbdxZGIZ0rYdyCgAA0AGsWLFCjgS/PJ6EBu8PBoMyFWXK318hY0w7p2s9lFMAAIAOYOmyZXJ5U5SU6G3w/mAwIFNZpoLSCu0pLGnndK2HcgoAAGBzu3fvVt6mrXK6vcpMTWpwHY/HI6dCkjNKi1eubeeErYdyCgAAYHPz5s2TM94vTyBZfk98I2tZCvi8slzRyl21rl3ztSbKKQAAgM3Nmzev+puhgimKcTVe34JBvyRpVd7m9orW6iinAAAANjdnzhy5vMny+xMbnKlfI+APSJal9Vvz2zFd66KcAgAA2FgoFNI338xQVFJ2o+eb1ggEAjJVldq8s6Cd0rU+yikAAICNrVixQntLqxTtDapbZuoh1w0EAjKV5crfX95O6Vof5RQAAMDGpk+fLqcvRcHUDPniDv7a0rqqLydVrsKyKu0vLm2nhK2LcgoAAGBjs2fPlsuTpGBSUA5H4+ebSlJCQoJclpEcLi3poJeTopwCAADYWG5urlz+TCUl+pqxtqVET3z15aRWd8zLSVFOAQAAbMoYo6UrVsnlCSotyd+sxwT8iZLDqVXrN7VtuDZCOQUAALCpjRs3qjjklDPGrfTkQLMeEwxWr7du0/a2jNZmKKcAAAA2lZubK6fbJ48/qPjY6GY9JuCvLqcbtu9uy2hthnIKAABgU7m5uXK4ffJ6vXIc4uL7B6q+nFSFtheWtXG6tkE5BQAAsKnc3Fw546rLaXNVl9MyFZQblZV1vOudUk4BAABsKjc3t/oap4nNL6cej0cOUyXjjNHS1evbLlwboZwCAADYkDFGS5evkCs+UalJic1+nGVVX07K4YrW4g54rVPKKQAAgA1t2LBBJSGHnNFxSg02b6Z+jUCiV3JGafnajW2Uru1QTgEAAGxoyZIl1eebBpIVHxPVoscGAwHJktZu7niXk6KcAgAA2NCSJUvkcHvl9fnkcrassgUCASkU0ob8fW2Uru24wh0AAAAA9S1ZskROt09eb3O+trSumhn7WwtCbZCsbXHkFAAAwIZyc3Pl8qW2aKZ+jbS0NIXKS1VQVqUtO3a1Qbq2QzkFAACwmfLyci1dvlJOT7JSk/wtfnxcXJy88TGyomI1dc7CNkjYdiinAAAANrNo0SJVuuIU6/Ep/TDKqSSlpyTLckVrzncrWjld26KcAgAA2MzMmTPl9AQVTM2QO+bwpghlpKdLMlq0Kq91w7UxyikAAIDNzJgxQy5PsoKBgByWdVhjZGd3lSkv07Ite1VZWdnKCdsO5RQAAMBGKioq9MUXXygqpZsy01IOe5zMzCy5TIUqnG59/e28VkzYtiinAAAANjJjxgwVlEvxSRnq1TX9sMdxOBzqlpkqK8atf3/6dSsmbFuUUwAAABv517/+JcsZpawu2YpxHdkl6fsfd5wUqtIXc5Z2mI/2KacAAAA2sXfvXr311luSpF69eh7xeL1791a006H9ZSF9+umnRzxee6CcAgAA2MQf/vAHlZSUqE+fo5WcfPjnm9ZwOp3qdVQvSdKjjz4qY8wRj9nWKKcAAAA2sHbtWv3+97+XJN1xxx2yLMkyVcraN19H53+urH3zZZmqFo97bN++ckW5NGfOHH300UetHbvVhb2cTps2TZdccokyMjJkWZY++OCDcEcCAABoVxUVFfq///s/FRcX68wzz9T5F1yg4wqmacyiK3TVkls0fOV9umrJLbpp3qXqtatlk5tiY+N08sknS5J+/vOfa8+ePW3xK7SasJfToqIiDRgwQM8//3y4owAAALSLjRs3atq0aVq9erXy8vI0evRozZgxQx6PR//4xz8Ut+ZTjd78gDwVO+o8LqF8hy5efleLC+qZZ56po446Sps3b9aIESO0detW237EbxkbJbMsS++//74uv/zycEcBAABoM3/+859122231VlmWZb+85//6MLzz1P0SyfIVbRNDV1+30gqjE7RywP+K2M5m9zW5r2lykqM1am+fTr//PNVVFQkqfpSUwUFBYqPj2+F36j1HNn1CQAAANAqjDG66qqrdGa2U1NuaLwwWpK85Ts07dVxmlGQ2uS4Ll+qyrev0b5v3qyzPBQK2fLoadg/1gcAAIg0N910U6P3pXua93WlabHlspyuJn+q9u9SVdG+BsewDvOrUdsSR04BAADamdvt1v79+xu8z7HhW+mdK5sc4+nf/lKPZZ3SrO1Fuyy5HPWPSbrd7mY9vj1RTgEAANqZZVmNn+t59DDJmyEVbFX1Gab1Hi15M5R47DmSo+lzTjsaPtYHAACwE4dTuuDJ728c/LH797cveKJTFlPJBuV0//79ysnJUU5OjiRp3bp1ysnJ0YYNG8IbDAAAIFz6Xir96A3Jm153uTejennfS8OTqx2E/VJSU6ZM0VlnnVVv+fXXX69//OMf7R8IAADALkJVUt630v7tUkKqlH1qpz1iWiPs5RQAAACoEfaP9QEAAIAalFMAAADYBuUUAAAAtkE5BQAAgG00u5xed911KikpkSRt3LixzQIBAAAgcjX7G6ISEhJUVlamuLg4ZWdny+/3a8CAARowYIAGDhyoAQMG6Nhjj1VUVFRb5gUAAEAndliXksrLy1NOTo4WLVpU+9/169fL5XKpT58+WrRoUVtkBQAAQCfXatc5LSwsVE5Ojr777juNGTOmNYZEGI0fP1579+7V6NGj5XBwajIAAGgfXIQf9Xz++ee68LIRiu3aX/f//DrdffvPwx0JAABECA6JoZ6XXnpJcd0GKaH/OfrrZ/MUCoXCHQkAAEQIyinqqKqq0pSp0xSd0VsyRvscHk2ZNT/csQAAQISgnKKOxYsXq7DKpVhfsrr44+SIjdd/v/gm3LEAAECEoJyijqlTp8rlTVZSarqO6dVDClVp+ndrwh0LAABECMop6pgxY4ac3mQlJSWpd+9eqirZry3F0tqNW8MdDQAARADKKeqYOXOmooJdlJESlMfjVXKiW45Yj9755KtwRwMAABGAcopamzZt0pb8PXLFJyo7I1WSdHSvXpJladLsxWFOBwAAIgHlFLVmz54tpydJ/pQMJSbESZJ69uwhU16ilTv2c0kpAADQ5iinqDVr1iw5vUkKBIOKcla/NNLTM2TKi1VqorR41frwBgQAAJ0e5RS1Zs2aJVdiupKCwdplUVFRSvJ5ZEXHafKshWFMBwAAIgHlFJKkiooKzVuwUFH+DHVJT6pzX1ZmhmRZWrh8bZjSAQCASEE5hSRp0aJFqnC5FesNKCslWOe+lJQUKRTS6i27wpQOAABECle4A8AeZs2aJZcnWcGUNLmj674skpOTFaoo1ZYCS8YYWZYVppQAAKCz48gpJNVMhkpWIBisVz5TUlJkKkpVWFqlHXsKw5QQAABEAsopJEkzZ81SVDBL6UmBevfFx8cr1uWQXNGat3hZGNIBAIBIQTmFCgsLtW7jVjnj/eqantzgOknBRFnOaOUsW9PO6QAAQCShnELLli2TMyEgt8+vgC+hwXVSk1MkS1q+bmM7pwMAAJGEcgotXbpUzgS/vD6/op0NvySSk5OlUEhrtuS3czoAABBJmK2P6nIa75fH62l0Jn5KSrJMZZm27K1o53QAACCScOQUys3NlcuXqoDX2+g6yckpClWUqbDcKH/PvnZMBwAAIgnlFMpdtlzOhIBSAr5G13G73XJHO2VFxWrOoqXtmA4AAEQSymmEKyoq0sYdu+WIilX6Qd8MdbAkv0+WK1o5y1a3UzoAABBpKKcRbvny5XLG+RSb4FPA2/BM/RopKcmSZWnpWmbsAwCAtkE5jXC5ublyxHnl9fnkamSmfo2U5JTvZ+zvbKd0AAAg0lBOI9zSpUvldHvlPcRkqBrJyckKVZRpy75yGWPaIR0AAIg0lNMIV32N04ASvZ4m101JSZapKFVRpbR15552SAcAACIN5TTC5S5dKqc3WcmHmKlfIzY2TgmxLjmi4vTt/MXtkA4AAEQaymkEKy4u1vpN2+SIdistKdCsx6QmJ0muKM1dsqKN0wEAgEhEOY1gK1askCPOo9h4r/y+Q8/Ur5GeliZJWrKGGfsAAKD1UU4j2NKlS6tn6vsDinU5m/WY9PR0mcoKrdm2t23DAQCAiEQ5jWC15dTrk2VZzXpMWlqaTHmJdhZXaX9xSRsnBAAAkYZyGsFyc3PljE+U19P0TP0aiYk+RTuN5IrRzAVL2jAdAACIRJTTCLZ06VK5vClKCjR9jdP/sZQaDMiKitXsRUvbLBsAAIhMlNMIVVpaqjXrN8gZn6i0JH+LHpuenipZ0qKV69smHAAAiFiU0wi1cuVKWbEexcR7FPQ1/2N9SUpPS5dCRqs27WijdAAAIFJRTiNU9deW+uTzJysuytWix6alpSlUUaotBZWqrKxso4QAACASUU4jVG5urhzxifJ4vXI4mjdTv0ZSUlDOUIWqXDFauHRlGyUEAACRiHIaob777js54wPyelv2kb4kWZZDyX6vHNFxmsHXmAIAgFZEOY1Q8+fPlysxrcWToWqkp6VIllMLl61p5WQAACCSUU4j0LZt27Qlf7dc8T5lpaUc1hjpaemSjJblbW/dcAAAIKJRTiPQ/Pnz5fQE5Q2mKjEh9rDGSEtLU6i8RJv2lSsUCrVyQgAAEKkopxFo+vTpcnmSFAgmK8blPKwxUlJS5KgqU7kjWnO+W9bKCQEAQKSinEagSZMmyeXPVGpq8mGP4XK5lJWaLCvarQlfftOK6QAAQCSjnEaYzZs3a+HChXIlpqp7VsYRjdWrZw/JsjRlfm4rpQMAAJGOchph/vnPf8oYo8zMLMW73Uc0Vu/evRUqK9aKHSVav2FjKyUEAACRjHIaQUpLS/X8889Lko499tgjHi8pKUmZAY+shICefPlfRzweAAAA5TSCvPDCC9q8ebOysrJ0zDF9WmXMk08cIsuy9N+pi7Rz585WGRMAAEQuymmEWLt2re6//35J0gMPPCCns/FZ+papUta++To6/3Nl7Zsvy1Q1um6fPn0USHCrKtqjJ554otVzAwCAyBL2cvr444/rhBNOkMfjUUpKii6//HKtWLEi3LE6lfz8fF122WUqLi7WmWeeqZtuuqnRdXvt+lo3zbtUVy25RcNX3qerltyim+Zdql67vm5wfcuyNOSEwZKk559/XmvW8I1RAADg8IW9nE6dOlVjxozRrFmzNGnSJFVWVuq8885TUVFRuKN1eJMnT9aYMWM0YMAALVmyROnp6frnP/8ph6Php73Xrq918fK7lFC+o87yhPIdunj5XY0W1KzMLHXJ7qqysjLddNNNKisra/XfBQAARAbLGGPCHeJA+fn5SklJ0dSpU/WDH/wg3HE6tKeeekp33XWXJKlnz55655131KdP9bmmT36xWkVllUr1xkiq/ij/Z4tGyFOeL6uBsYykwugUvTzgvzJW3VMCtheUqaJ0v/728wtVUlKirl27avDgwQoGg3rllVdkWQ2NCAAAUJ/tyunq1at11FFHafHixerXr1+443RoM2bM0Omnn97gff6zbpIV7VbV/l2SpNP9O/Xx4FlNjnnx/JP1zZ6kOsucCUGZ8hLtmfy3euvv379f8fHxh5EeAABEIle4AxzIGKM777xTp59+OsW0FQwcOLDxOy1LLm+SHLHVxTEjsXkfxWcEPYpydq2zzBEdp4pdmw43JgAAQC1bHTkdM2aMPvnkE33zzTfKysoKd5wOzxij4uLiBu/7bnOB9hZX1N4O5M/VSTNubHLM2ae9qt3JJ9RbnuiOUv9Mb73lbrebj/UBAECz2aac/uIXv9AHH3ygadOmqXv37uGOE3lCVdKz/aSCrao+w/RgluTNkG5fLDkavwwVAADAkQj7bH1jjMaOHav33ntPX3/9NcU0XBxO6YInv79x8JHO729f8ATFFAAAtKmwHzn9+c9/rrfeeksffvihjj766NrlPp9PcXFxYUwWoZZOkCbeJRVs+d8yb2Z1Me17afhyAQCAiBD2ctrY+YivvfaabrjhhvYNg2qhKinvW2n/dikhVco+lSOmAACgXYS9nAIAAAA1wn7OKQAAAFCDcgoAAADboJwCAADANiinAAAAsA3KKQAAAGyDcgoAAADboJwCAADANiinAAAAsA3KKQAAwGHKz8/X+vXrwx2jU6GcAgAAHIbc3Fwd1bu3ep94lh56/Klwx+k0KKcAAACHYdy4cSoL9JT3pBF69pOFWrF6TbgjdQqWMcaEOwQAAEBHsnfvXqVmZCr+5GsVn5Su0vIK/SAtpPEvPRnuaB0eR04BAABa6Ntvv5WJT5IvvasuP/cMmaoKTVmxQ5u2bA13tA6PcgoAANBCc+bMkcuXpkBSsnp066qUhCg5vKl67MV/hjtah0c5BQAAaKHZs2fLFchUcsAvydIpJ54gyWji/FUKhULhjtehUU4BAABawBijOfMWyJWYqi5pSZKk3r17S2WFKpBb3y5YEuaEHRvlFAAAoAXWrl2rfeWSKzZBXdJTJUnR0dHKTAnIinbri28XhDlhx0Y5BQAAaIHZs2fLmeCXPylVCbFRtcu7ZGZJlqV5S7mk1JGgnAIAALTAnDlz5IwPKBAIyLKs2uWZmZkyFWVam78/jOk6PsopAABAC8yePVsuX6qSg/46y5OTkxUqL9XukioVl5aHKV3HRzkFAABoprKyMi3IWSRXYqqy0pLr3Of3J8ppKhRyuLRoxdowJez4KKcAAADNtGDBAlVFJcjtDSgjue6RU8tyyO/1yHLFaNGy1WFK2PFRTgEAAJpp5syZcnmDCqamyx3tqnd/8PuP+ldv5JuiDhflFAAAoJlmzpwppzdVgUCwzmSoGv5Ev2RZWr81PwzpOgfKKQAAQDMYYzR9+nRFJ3dVZmqwwXX8fr9MVaW27GLG/uGinAIAADTD6tWrlb+/TFGeJHXPSmtwncTERJnKcu0sYrb+4aKcAgAANMP06dPl8qYokJapRHdsg+v4/dXltKC0UqXlle2csHOgnAIAADTDN998I5c3WUlJSXI46p9vKkk+X3U5rQxJ67fuaOeEnQPlFAAAoBlycnLk8mcq9aCL7x/I5XLJ446R5YrWstXr2y9cJ0I5BQAAaEJVVZWWrVwjlzdJ6cmBQ67r83olh0ur8ja3U7rOhXIKAADQhLVr16rCGauouHilJTV+5FSqPu9UllHe5u3tE66ToZwCAAA0YcmSJXLGJ8rrT2rw4vsH8if6JVnauH1X+4TrZCinAAAATViyZImcbp+8Xm+DF98/UO21TndzrdPDcejqDwAAAC1ZskSOhIB8Pm+T69Zc63RXkWmHZJ0PR04BAACasHjJErm8KUr2+5pct/papxUqLKtSSVlFO6TrXCinAAAAh1BWVqbVa/PkjPMqPbnhry09kMfjkcNUSg6XVqzd0A4JOxfKKQAAwCGsXLlSJiZB0e4EJfmb/ljfshzyJLhluaK0fG1eOyTsXCinAAAAh5CbmyuH2ydfIEmxUc5mPcbv9UiOKK3mWqctRjkFAAA4hNqZ+j6fHE3M1K/h9ydKktZv3taGyTonyikAAMAh1MzUT/Q2/ZF+DX+iX7KkjTt2t2GyzolyCgAAcAiLlyxRlC9NyYGmZ+rX8Pv9UlWVtnKt0xajnAIAADSiqKhI6zdtlSPOo4yUpmfq10hMTFSosly7iyvbMF3nRDkFAABoxLJly+SI8ynO41PAm9Dsx1Vf67RcJVWW9hZw9LQlKKcAAACNWLJkiZzxPnn9QcW4ml+b4uLiFOWQrKgYLVm1tg0Tdj6UUwAAgEYsWbJEDneifF6frGbO1K9myeeJk+WK1rI1G9ssX2dEOQUAAGhEbm6unJ4kJfqaP1O/RjAxUbIcWrqGb4lqCcopAABAI5YsWSKXN1kpLZipXyM9PV2ypMWr+ZaolqCcAgAANGDv3r3anL9HjtgEpacEWvz4zMwMmYoKrdm+rw3SdV6UUwAAgAbUfDNUgi+gxHh3ix+fkZEhU16sggqnVq3nvNPmopwCAAA0YPHixdVfW5roV3QLZurXiI2NU1rQJ0dsgt6a8EUbJOycKKcAAAANWLx4sRxun3yJiS2cqf8/fXofJTkcmvhtTuuG68QopwAAAA1YvHixXN5k+Q9jpn6NY445RqHSIq3cWaYNm7a0YrrOi3IKAABwkFAopO8WL5bLm6q0YOJhjxMMBpWWGC9HvF8vvfV+6wXsxCinAAAAB1m6dKkKi0oV5Y5XSlLwiMYa2L+f5HDq/a++baV0nRvlFAAA4CAzZsyQJKWkpMjhOLK6dOyxx8phObRh42YtWrSoNeJ1apRTAACAg3z88ceSpIz0jCMeKy4uTukZ6ZKk119//YjH6+wopwAAAAfYvHmzvvii+tJP3bt3b5Uxs7tmS5L++te/assWJkYdCuUUAADge8YY3XHHHSovL9cpp54if8Df7MdapkpZ++br6PzPlbVvvixTVXtfWlqaevbqpf3792vs2LEyxrRF/E4h7OX0xRdfVP/+/eX1euX1enXKKafos88+C3csAAAQgR599FG9++67crlcevyxx5v9uF67vtZN8y7VVUtu0fCV9+mqJbfopnmXqteuryVJlmXpxtGj5XK59P777+uvf/1rW/0KHV7Yy2lWVpaeeOIJzZs3T/PmzdOwYcN02WWXKTc3N9zRAABABFi5cqXGjh2roUOH6v7775ck/elPf9IJJ5zQrMf32vW1Ll5+lxLKd9RZnlC+Qxcvv6u2oHbrlq3HH68uvL/85S/1yiuvKC8vT/v372/F36bjs4wNjysHAgE9/fTTuummm8IdBQAAdHLz5s2rU0THjRunX/3qVyqpqNKTn6+W02kp4I5q8LGWqdLPFo2QpzxfDX2HlJFUGJ2iX6W/rrOPSdXwY1M0cuRIffTRR7Xr/OY3v9GTTz7Zyr9Vx+UKd4ADVVVV6d1331VRUZFOOeWUcMcBAAARIDMzs87tcePGady4cbJcMfIP+4kcsfGqKi5o8LGnebfr18flNzq2JclbvkP509/QjU8vV/GyqfXWcblsVcfCzhZ7Y/HixTrllFNUWlqqhIQEvf/+++rbt2+4YwEAgAjg9Tb89aQmVKmq4n2SjKxGrnWaFlParG0klW6QKS9p8L677767WWNEClt8rF9eXq4NGzZo7969Gj9+vP72t79p6tSpFFQAANDmjDEqLi5u8L7KUEiVVY1XJdfGmfKN/1GT29g34j9ydj9NllX/w3+3293g8khli3J6sHPOOUc9e/bUyy+/HO4oAAAAjQtVSc/2kwq2qvoM04NZkjdDun2x5HC2d7oOKeyz9RtijFFZWVm4YwAAAByawyldUDOZ6eCjn9/fvuAJimkLhP2c03vuuUcXXnihunTposLCQv373//WlClTNHHixHBHAwAAaFrfS6UfvSFNvEsqOODbn7wZ1cW076Xhy9YBhb2cbt++XaNGjdLWrVvl8/nUv39/TZw4Ueeee264owEAADRP30ulPhdJed9K+7dLCalS9qkcMT0MtjznFAAAAJHJluecAgAAIDJRTgEAAGAblFMAAADYBuUUAAAAtkE5BQAAgG1QTsNg586dOuGEE9T96GO1YOHCcMcBAACwDcppGNx///36bnOBCrJ/oNuf+lu44wAAANgG5bSdlZWV6c233lJczxPlCnTRst0hFewvDncsAAAAW6CctrPFixerxOFWfEoXRZXskomO1xcz5oY7FgAAgC1QTtvZggUL5PKlyp+UovSgV5YrWvOWrAp3LAAAAFugnLazBQsWyOlNkd+fqLTUVElGS9dtDncsAAAAW6CctrMFCxcqyp+hlECiUlNTZSoqtHb7vnDHAgAAsAXKaTuqqKjQ4mWr5IjzKDM1SSkpKTIVpcovqlBFZVW44wEAAIQd5bQdLV++XJWuOEW7E5SWFFAg4JepLFeVnFqxno/2AQAAKKftaOHChXK6vUoMJis2yimn0yVvQqysqBjlLFsd7ngAAABhRzltRwsWLJDDnahEX6Isy5IkJQcCkuXU0jUbwpwOAAAg/Cin7Wju3LlyeZIVDCTWLgsmBSVLWr1xW/iCAQAA2ATltJ1UVFRowcIcuRLTlJWaVLs8KZgkU1mhvPyCMKYDAACwB8ppO5kzZ44qXHGK9fqVkRKsXR4MBmUqy5W/v1zGmDAmBAAACD/KaTt58cUXFeXPVHqXbkqIcdUury6nZSooKdeewpIwJgQAAAg/ymk7mDVrlt76z38Vk91ffXofVTsZSpI8ngRFOYzkjNLilWvCmBIAACD8KKft4Omnn1Zs9kD1Ov4M9e+ZddC9lrwJ8bJc0Vq5blNY8gEAANgF5bSNVVZW6suvvlZMl2PV7+geinbV3+V+n0+StG7z9vaOBwAAYCuU0zaWl5enopBL0Z6AemSlNbiOL9EnWdLG7TvbOR0AAIC9UE7b2Pr16+WM9yve61d8TFSD6yT6EqWQ0ZadXE4KAABENlfTq+BI5OXlyRHnUXx8vBwHTIQ6kM/nk6mqUH5hqJ3TAQAA2AtHTtvY+vXr5Yj1yu12N7qOz+eTqazQ3tJKrnUKAAAiGuW0jeXl5cmZ4JfHHdfoOomJPpmqchWXVWpvUWk7pgMAALAXymkbW79+vZxun3ye+EbXiY9PkMOEJKdLq9ZzOSkAABC5KKdtLG/jZjli3Er0JjS6jmVZ8sTHyXJGadX6je2YDgAAwF4op22osrJSW3bukeWMUjDRd8h1E31eyeHUmg1b2ikdAACA/VBO29DmzZtlnDFyRsXIf4gjp5IUCPglSWs3bWuPaAAAALZEOW1D69evlxUdJ7fHq+go5yHXDQaCkjHasG1XO6UDAACwH65z2oby8vLkiHHL7XY3eo3TGoFAQKayQlv2lrdTOgAAAPvhyGkbysvLkyParfj4xmfq1wgGAzKVZdpdElJVFRfjBwAAkYly2oaqL8AfL/chrnFaw+/3y1SWK2Q5tTJvczukAwAAsB/KaRtat26dnO5EeZtx5NTpdMmXECfLFaNZC5e0QzoAAAD7oZy2obVr18rh9snv8zRr/cyMNMnh0uzvlrdxMgAAAHuinLaR8vJybdq6XY4Yt4KJ3mY9JisjU5LR4tVciB8AAEQmymkbycvLk6Ldiop1N3mN0xpZWZkKlZdo7e4SVVZWtXFCAAAA+6GctpE1a9bIEZugeJ9fMU1c47RGenq6XKEyVThiNX3+d22cEAAAwH4op21k7dq1csR6FJ/gkctx6Guc1nA4nMpITZIVHafPv5nfxgkBAADsh3LaRqrLaYIS4uNlNXEB/gN165otGaNZuWvbMB0AAIA9UU7byMqVK+WI98uT0LzzTWt065atUFmxVu0oUgXnnQIAgAhDOW0jy5Ytk8ub1OyZ+jUyMzPlqCxVmXFo1qJlbZQOAADAniinbaCsrEzrNmySM86r1GBiix7rdLqUnhqUFRWrL2fMbZuAAAAANkU5bQOrVq2SYhIU7U5QMLF5F+A/UHaXrpJlaU7uqjZIBwAAYF+U0zawfPlyOeK88vqDio1ytfjxXbt2kaks1/JNu9sgHQAAgH1RTtvAsmXL5HR75fF65WzmZaQOlJXVRaa8RAWVTuVt2tIGCQEAAOyJctoGli1bJoc7UR5PyyZD1YiNjVXA45Yj2q3Pps5u5XQAAAD2RTltA0uWLJHLm6Kgr+Xnm9bokpUuuaL0zcKlrZgMAADA3iinraykpERLV6yUyxNUZmrSYY/TrWu2FKrSvJUbWzEdAACAvVFOW9miRYukmAS5vYlK9h/ex/qS1LNnT4VKi5Rf5tSy1XxbFAAAiAyU01Y2ffp0yXIokJQsh+Pwd6/b7VZGkk+O2AS989GXrZgQAADAviinrezTTz+VJGWkpx/xWH16HyU5nJo8ddoRjwUAANARUE5b0fLlyzV16lRJ1V9DeqSOPvpoSdKcOXO0devWIx4PAADA7iinrWTbtm26/LJL9IOuDv1uxDEaaK2SZaqOaMxgMKhgMKhQKKTXXnutlZICAADYl63K6eOPPy7LsnT77beHO0qL7Nq1S09cf5omDd+mKTfE6+6uORq78Q6NWXSFeu36+ojG7t69uyTpmWee0Y4dO1ojLgAAgG3ZppzOnTtXr7zyivr37x/uKC2ye/duPXnD6frDyfnK9NbdnZ7yHbp4+V1HVFC7dOmio3r31p49e3T55ZerqKjoSCMDAADYli3K6f79+zVy5Ej99a9/ld/vD3ecRhljVFRUpH379mn9+vV67bXXdNKJQ/SLXpskWTr4m0prbp659veqLC9RRUV5i39CVSE9cP/98vl8mjlzpk477TRNmTJFs2bN0nvvvadly5a1924AAABoM5YxxoQ7xPXXX69AIKA//vGPGjp0qAYOHKhnn3023LHqKSoqUkJCQp1lZ2Y7NeWG+CYfe/H8k/XNnpZflD8qOVsFM/+j8u1rGrz/t7/9rR5//PEWjwsAAGBHrnAH+Pe//60FCxZo7ty54Y5yWNI9VtMrScoIehTl7NryDVRVHvLuQCDQ8jEBAABsKqzldOPGjbrtttv0xRdfKDY2NpxRmsXtdmvt2rVyOBzy+XyKioqSY8O30jtXNvnYn1z3Y/0w+YQWb9NhWRqYda/iYxp+qtxud4vHBAAAsKuwfqz/wQcf6IorrpDT6axdVlVVJcuy5HA4VFZWVuc+WwpVSc/2kwq2SmpoV1qSN0O6fbHksPnvAgAAEGZhLaeFhYXKy8urs2z06NHq06eP7rrrLvXr1y9MyVpo6QTpP//3/Y0Dd+f3H/n/6A2p76XtnQoAAKDDCevH+h6Pp14BjY+PVzAY7DjFVKounj96Q5p4l1Sw5X/LvRnSBU9QTAEAAJop7BOiOo2+l0p9LpLyvpX2b5cSUqXsU/koHwAAoAVscSkpAAAAQLLJRfgBAAAAiXIKAAAAG6GcAgAAwDYopwAAALANyikAAABsg3IKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA23CFO8CRMsaosLAw3DEAAADQBI/HI8uyDrlOhy+nhYWF8vl84Y4BAACAJuzbt09er/eQ61jGGNNOedoER06rFRQUqEuXLtq4cWOTTzrqY/8dPvbdkWH/HRn235Fh/x0+9t3hiYgjp5Zl8aI4gNfrZX8cAfbf4WPfHRn235Fh/x0Z9t/hY9+1PiZEAQAAwDYopwAAALANymknERMTowcffFAxMTHhjtIhsf8OH/vuyLD/jgz778iw/w4f+67tdPgJUQAAAOg8OHIKAAAA26CcAgAAwDYopwAAALANyikAAABsg3LaSX3yySc66aSTFBcXp6SkJP3whz8Md6QOp6ysTAMHDpRlWcrJyQl3nA5h/fr1uummm9S9e3fFxcWpZ8+eevDBB1VeXh7uaLb1wgsvqHv37oqNjdXgwYM1ffr0cEeyvccff1wnnHCCPB6PUlJSdPnll2vFihXhjtVhPf7447IsS7fffnu4o3QYmzdv1nXXXadgMCi3262BAwdq/vz54Y7VaVBOO6Hx48dr1KhRGj16tBYtWqQZM2bo2muvDXesDuc3v/mNMjIywh2jQ1m+fLlCoZBefvll5ebm6o9//KNeeukl3XPPPeGOZkvvvPOObr/9dt17771auHChzjjjDF144YXasGFDuKPZ2tSpUzVmzBjNmjVLkyZNUmVlpc477zwVFRWFO1qHM3fuXL3yyivq379/uKN0GHv27NFpp52mqKgoffbZZ1q6dKl+//vfKzExMdzROg0uJdXJVFZWqlu3bnrooYd00003hTtOh/XZZ5/pzjvv1Pjx43Xsscdq4cKFGjhwYLhjdUhPP/20XnzxRa1duzbcUWznpJNO0vHHH68XX3yxdtkxxxyjyy+/XI8//ngYk3Us+fn5SklJ0dSpU/WDH/wg3HE6jP379+v444/XCy+8oEceeUQDBw7Us88+G+5Ytvfb3/5WM2bM4FOONsSR005mwYIF2rx5sxwOhwYNGqT09HRdeOGFys3NDXe0DmP79u26+eab9c9//lNutzvccTq8ffv2KRAIhDuG7ZSXl2v+/Pk677zz6iw/77zz9O2334YpVce0b98+SeJ11kJjxozRRRddpHPOOSfcUTqUCRMmaMiQIbrqqquUkpKiQYMG6a9//Wu4Y3UqlNNOpubo1Lhx43Tffffp448/lt/v15lnnqndu3eHOZ39GWN0ww036JZbbtGQIUPCHafDW7NmjZ577jndcsst4Y5iOzt37lRVVZVSU1PrLE9NTdW2bdvClKrjMcbozjvv1Omnn65+/fqFO06H8e9//1sLFizgCP1hWLt2rV588UUdddRR+vzzz3XLLbfol7/8pd54441wR+s0KKcdxLhx42RZ1iF/5s2bp1AoJEm69957NWLECA0ePFivvfaaLMvSu+++G+bfInyau/+ee+45FRQU6O677w53ZFtp7v470JYtW3TBBRfoqquu0k9+8pMwJbc/y7Lq3DbG1FuGxo0dO1bfffed3n777XBH6TA2btyo2267Tf/6178UGxsb7jgdTigU0vHHH6/HHntMgwYN0s9+9jPdfPPNdU7PwZFxhTsAmmfs2LG65pprDrlOt27dVFhYKEnq27dv7fKYmBj16NEjoidZNHf/PfLII5o1a1a970oeMmSIRo4cqddff70tY9pWc/dfjS1btuiss87SKaecoldeeaWN03VMSUlJcjqd9Y6S7tixo97RVDTsF7/4hSZMmKBp06YpKysr3HE6jPnz52vHjh0aPHhw7bKqqipNmzZNzz//vMrKyuR0OsOY0N7S09Pr/BsrVZ8rPn78+DAl6nwopx1EUlKSkpKSmlxv8ODBiomJ0YoVK3T66adLkioqKrR+/XplZ2e3dUzbau7++/Of/6xHHnmk9vaWLVt0/vnn65133tFJJ53UlhFtrbn7T6q+xMpZZ51Ve9Te4eADmoZER0dr8ODBmjRpkq644ora5ZMmTdJll10WxmT2Z4zRL37xC73//vuaMmWKunfvHu5IHcrZZ5+txYsX11k2evRo9enTR3fddRfFtAmnnXZavUuXrVy5MqL/jW1tlNNOxuv16pZbbtGDDz6oLl26KDs7W08//bQk6aqrrgpzOvvr2rVrndsJCQmSpJ49e3Jkphm2bNmioUOHqmvXrnrmmWeUn59fe19aWloYk9nTnXfeqVGjRmnIkCG1R5k3bNjAObpNGDNmjN566y19+OGH8ng8tUeffT6f4uLiwpzO/jweT73zc+Pj4xUMBjlvtxnuuOMOnXrqqXrsscf0ox/9SHPmzNErr7zCp0StiHLaCT399NNyuVwaNWqUSkpKdNJJJ+nrr7+W3+8PdzR0cl988YVWr16t1atX1yvzXLWuvquvvlq7du3Sww8/rK1bt6pfv3769NNPOQLThJpz+4YOHVpn+WuvvaYbbrih/QMhopxwwgl6//33dffdd+vhhx9W9+7d9eyzz2rkyJHhjtZpcJ1TAAAA2AYngwEAAMA2KKcAAACwDcopAAAAbINyCgAAANugnAIAAMA2KKcAAACwDcopAAAAbINyCgAAANugnAIAAMA2KKcAAACwDcopANjE66+/rr59+8rtdqtPnz76+OOPwx0JANod5RQAbOD999/XmDFjdN9992nJkiW68MILdcstt4Q7FgC0O8sYY8IdAgAi3emnn65hw4bp4YcfliRNmjRJV111lfbu3RveYADQzjhyCgBhVlhYqJkzZ+qiiy6qXTZx4kQNHDgwfKEAIExc4Q4AAJFu0aJFsixL/fv3V3Fxsd58800999xzGj9+fLijAUC7o5wCQJjl5OSoT58+ysnJ0amnnipJuuKKK2qPpH755ZdavHix7rjjjnDGBIB2wTmnABBmP/nJT1RWVqYXXnhBy5Yt08yZM3Xvvffqjjvu0O9+97twxwOAdsU5pwAQZjk5ORo0aJA8Ho9OPPFE3XbbbRo1apRmzZolSbrwwgu1bNmyMKcEgPZBOQWAMKqsrFRubq769OlTZ/miRYt0xhlnSJJWrVqlo446KhzxAKDdcc4pAITR8uXLVVpaqkceeUTp6elyu9168cUXtW7dOt18883at2+fEhIS5HLx5xpAZOCvHQCEUU5OjtLT0xUfH68zzjhD8fHxOv300zV58mSlp6drxowZOvbYY8MdEwDaDeUUAMIoJydHJ510kt5///0G71+yZIn69evXzqkAIHw45xQAwignJ0f9+/dv9P7c3FzKKYCIwqWkACCMkpOT9dJLL2nEiBHhjgIAtkA5BQAAgG3wsT4AAABsg3IKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA26CcAgAAwDYopwAAALANyikAAABsg3IKAAAA2/j/hLrkmo+6CUcAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ax, = az.plot_forest(\n", + " idata, \n", + " var_names=[\"predictors\"], \n", + " coords={\"predictors_dim\": range(D0)},\n", + " kind='ridgeplot',\n", + " ridgeplot_truncate=False, \n", + " ridgeplot_alpha=0.5,\n", + " hdi_prob=0.95, \n", + " combined=True,\n", + " figsize=(8, 6)\n", + ")\n", + "ax.scatter(COEF[:D0][::-1], ax.get_yticks(), c='C1', label=\"Actual value\");\n", + "ax.set_xlabel(r\"$\\beta_i$\");\n", + "ax.set_ylim(bottom=None, top=1.55 * ax.get_yticks().max())\n", + "ax.set_yticklabels(range(D0)[::-1]);\n", + "ax.set_ylabel(r\"$i$\");\n", + "ax.legend(loc='upper center');\n", + "ax.set_title(\"Posterior distribution of nonzero coefficients\");" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Last updated: Thu Aug 22 2024\n", + "\n", + "Python implementation: CPython\n", + "Python version : 3.11.9\n", + "IPython version : 8.24.0\n", + "\n", + "pandas : 2.2.2\n", + "numpy : 1.26.4\n", + "bambi : 0.14.1.dev12+g64e57423.d20240730\n", + "arviz : 0.18.0\n", + "matplotlib: 3.8.4\n", + "pymc : 5.16.1\n", + "pytensor : 2.23.0\n", + "\n", + "Watermark: 2.4.3\n", + "\n" + ] + } + ], + "source": [ + "%load_ext watermark\n", + "%watermark -n -u -v -iv -w" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "bambi-dev", + "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.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From b5c8d4f46550f0fb234047a6c40835948e887a45 Mon Sep 17 00:00:00 2001 From: julianlheureux <“lheureuxjulian@gmail.com”> Date: Thu, 22 Aug 2024 21:32:18 -0300 Subject: [PATCH 6/9] Adding horseshoe prior and example2 --- bambi/backend/utils.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/bambi/backend/utils.py b/bambi/backend/utils.py index e7d78535f..94915e067 100644 --- a/bambi/backend/utils.py +++ b/bambi/backend/utils.py @@ -5,14 +5,16 @@ import pytensor.tensor as pt import pymc as pm -def Horseshoe(name, tau_nu = 3, lam_nu = 1, dims=None): + +def horseshoe(name, tau_nu=3, lam_nu=1, dims=None): tau = pm.HalfStudentT(f"{name}_tau", nu=tau_nu) - lam = pm.HalfStudentT(f"{name}_lam", nu=lam_nu, dims=dims) + lam = pm.HalfStudentT(f"{name}_lam", nu=lam_nu, dims=dims) beta_raw = pm.Normal(f"{name}_raw", 0, 1, dims=dims) - beta = pm.Deterministic(name, beta_raw * tau ** 2 * lam ** 2, dims=dims) + beta = pm.Deterministic(name, beta_raw * tau**2 * lam**2, dims=dims) return beta -MAPPING = {"Cumulative": pm.Categorical, "StoppingRatio": pm.Categorical, "Horseshoe": Horseshoe} + +MAPPING = {"Cumulative": pm.Categorical, "StoppingRatio": pm.Categorical, "Horseshoe": horseshoe} def get_distribution(dist): From 25678d4600a7420a05f5580c9f8d9b988c2fa822 Mon Sep 17 00:00:00 2001 From: julianlheureux <“lheureuxjulian@gmail.com”> Date: Fri, 23 Aug 2024 00:16:13 -0300 Subject: [PATCH 7/9] Update horseshoe notebook --- docs/notebooks/horseshoe_prior.ipynb | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/docs/notebooks/horseshoe_prior.ipynb b/docs/notebooks/horseshoe_prior.ipynb index a5869e974..3d5e39d3a 100644 --- a/docs/notebooks/horseshoe_prior.ipynb +++ b/docs/notebooks/horseshoe_prior.ipynb @@ -29,7 +29,29 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We are simulating 100 observations with the 50 explanatory variables, also simulated. " + "Here is what we did:\n", + "\n", + "* We defined an intercept.\n", + "* We defined a vector of 50 betas, 5 of which were drawn from a normal(5,1) distribution, and then assigned a random sign.\n", + "* We created the design matrix with normal(0,1) entries and set $\\sigma$ to 1.\n", + "* We calculated the deterministic means $\\mu$ using the intercept and the design matrix multiplied by the betas.\n", + "* We simulated 100 response variables (observations) from a normal distribution with mean $\\mu$ and standard deviation $\\sigma$.\n", + "\n", + "Next, we proceeded with the Bayesian estimation of the model. We proposed the horseshoe prior, for which the following parameters were calculated:\n", + "\n", + "$$\\mu_i = \\alpha + \\beta_1 x_{1i} + \\beta_2 x_{2i} + ... + \\beta_p x_{pi}$$\n", + "\n", + "$$y_i \\sim N(\\mu_i, \\sigma^2)$$\n", + "\n", + "$$\\alpha \\sim N(0,1)$$\n", + "\n", + "$$\\beta_j \\sim N(0,\\lambda_j^2 \\tau^2)$$\n", + "\n", + "$$\\lambda_j \\sim C^+(0,1)$$\n", + "\n", + "$$\\tau \\sim T^+(df=3)$$\n", + "\n", + "$$\\sigma^2 \\sim N^+(0,1)$$" ] }, { @@ -62,7 +84,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Then we create the dataframe and the term name for the set of variables, to define the formula. " + "Here we create the dataframe and the term name for the set of variables, to define the formula. " ] }, { From 025ae0ee907b4fdf6da6502766e328046fb1fcaa Mon Sep 17 00:00:00 2001 From: julianlheureux <“lheureuxjulian@gmail.com”> Date: Wed, 28 Aug 2024 15:00:47 -0300 Subject: [PATCH 8/9] Add docstring --- bambi/backend/utils.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/bambi/backend/utils.py b/bambi/backend/utils.py index 94915e067..625afad3f 100644 --- a/bambi/backend/utils.py +++ b/bambi/backend/utils.py @@ -7,6 +7,26 @@ def horseshoe(name, tau_nu=3, lam_nu=1, dims=None): + """Simulate a beta coefficient value with a horseshoe prior. + This is an internal function which is not supposed to be used by users. + This will be used only when a horseshoe prior is called for beta coefficients. + + Parameters + ---------- + name: str + is the name of the parameters as registered in the PyMC model + tau_nu: int, float + Degrees of freedom of tau. Default: 3 + lam_nu: int, float + Degrees of freedom of lam. Default: 1 (equivalent to a HalfCauchy) + dims: str + dimensions passed to PyMC. Default: None + + Returns + ------ + np.ndarray + Array with the beta coefficient simulated. + """ tau = pm.HalfStudentT(f"{name}_tau", nu=tau_nu) lam = pm.HalfStudentT(f"{name}_lam", nu=lam_nu, dims=dims) beta_raw = pm.Normal(f"{name}_raw", 0, 1, dims=dims) From 9d547dedb29fc276de64085ec7d63c86d64d272b Mon Sep 17 00:00:00 2001 From: julianlheureux <“lheureuxjulian@gmail.com”> Date: Wed, 28 Aug 2024 15:05:16 -0300 Subject: [PATCH 9/9] Add docstring2 --- bambi/backend/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/bambi/backend/utils.py b/bambi/backend/utils.py index 625afad3f..50f772a00 100644 --- a/bambi/backend/utils.py +++ b/bambi/backend/utils.py @@ -7,7 +7,7 @@ def horseshoe(name, tau_nu=3, lam_nu=1, dims=None): - """Simulate a beta coefficient value with a horseshoe prior. + """Simulate a beta coefficient value with a horseshoe prior. This is an internal function which is not supposed to be used by users. This will be used only when a horseshoe prior is called for beta coefficients. @@ -25,7 +25,7 @@ def horseshoe(name, tau_nu=3, lam_nu=1, dims=None): Returns ------ np.ndarray - Array with the beta coefficient simulated. + Array with the beta coefficient simulated. """ tau = pm.HalfStudentT(f"{name}_tau", nu=tau_nu) lam = pm.HalfStudentT(f"{name}_lam", nu=lam_nu, dims=dims)