diff --git a/examples/howto/decision_analysis.ipynb b/examples/howto/decision_analysis.ipynb new file mode 100644 index 000000000..290c2256f --- /dev/null +++ b/examples/howto/decision_analysis.ipynb @@ -0,0 +1,3039 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "2ef07cb2", + "metadata": {}, + "source": [ + "# Decision analysis with a utility function\n", + "\n", + ":::{post} Dec 15, 2022\n", + ":tags: Bayesian workflow, decision analysis\n", + ":category: beginner, how-to\n", + ":author: Oriol Abril Pla\n", + ":::" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "893d05de", + "metadata": {}, + "outputs": [], + "source": [ + "import arviz as az\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pymc as pm\n", + "\n", + "from scipy import stats\n", + "from xarray_einstats.stats import XrContinuousRV, XrDiscreteRV" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "e2dcb151", + "metadata": {}, + "outputs": [], + "source": [ + "RANDOM_SEED = 4235\n", + "rng = np.random.default_rng(RANDOM_SEED)" + ] + }, + { + "cell_type": "markdown", + "id": "53d11d04", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "In this example we will show how to perform Bayesian decison analysis over a set of discrete decisions. We will use as an example, we will take the role of a shopper, who has done their weekly grocery shopping at a different supermarket every week for the last year, and has logged:\n", + "\n", + "* The supermarket they went to that week\n", + "* The cumulative time it took from exiting home to getting back, so both travel time and unintuitive supermarket layout are penalized (in hours)\n", + "* The money spent (in euros)\n", + "* The number of rubbish bins that were filled that week" + ] + }, + { + "cell_type": "markdown", + "id": "b920026a", + "metadata": {}, + "source": [ + "## Example data" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "6e831402", + "metadata": {}, + "outputs": [], + "source": [ + "obs_supermarket = rng.choice([\"aldi\", \"caprabo\", \"condis\", \"consum\", \"keisy\"], size=52)\n", + "super_id, supermarkets = pd.factorize(obs_supermarket, sort=True)\n", + "time_means = np.array([7, 5, 4.7, 6, 5.2])\n", + "cost_means = np.array([90, 125, 130, 120, 150])\n", + "rubbish_means = np.array([5, 3, 4, 4, 3])\n", + "obs_time = rng.normal(time_means[super_id], 0.2)\n", + "obs_cost = rng.normal(cost_means[super_id], 1)\n", + "obs_rubbish = rng.poisson(rubbish_means[super_id])" + ] + }, + { + "cell_type": "markdown", + "id": "c5a92389", + "metadata": {}, + "source": [ + "## Model definition" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "65694151", + "metadata": {}, + "outputs": [], + "source": [ + "coords = {\n", + " \"supermarket\": supermarkets,\n", + " \"measure\": [\"time\", \"cost\"],\n", + "}\n", + "with pm.Model(coords=coords) as model:\n", + " model.add_coord(\"week\", length=len(super_id), mutable=True)\n", + " supermarket_id = pm.MutableData(\"supermarket_id\", super_id, dims=\"week\")\n", + "\n", + " mu = pm.Normal(\"mu\", 0, 1, dims=(\"measure\", \"supermarket\"))\n", + " sigma = pm.HalfNormal(\"sigma\", 0.5, dims=(\"measure\", \"supermarket\"))\n", + " lam = pm.Exponential(\"lam\", 0.3, dims=\"supermarket\")\n", + "\n", + " pm.Normal(\n", + " \"time\",\n", + " mu=mu[0, supermarket_id],\n", + " sigma=sigma[0, supermarket_id],\n", + " dims=\"week\",\n", + " observed=obs_time,\n", + " )\n", + " pm.Normal(\n", + " \"cost\",\n", + " mu=mu[1, supermarket_id],\n", + " sigma=sigma[1, supermarket_id],\n", + " dims=\"week\",\n", + " observed=obs_cost,\n", + " )\n", + " pm.Poisson(\"rubbish\", mu=lam[supermarket_id], dims=\"week\", observed=obs_rubbish)" + ] + }, + { + "cell_type": "markdown", + "id": "0c08664d", + "metadata": {}, + "source": [ + "## Inference" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "63b7c123", + "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: [mu, sigma, lam]\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 8 seconds.\n" + ] + } + ], + "source": [ + "with model:\n", + " idata = pm.sample()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "77c75650", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "
\n", + "
\n", + "
arviz.InferenceData
\n", + "
\n", + "
    \n", + " \n", + "
  • \n", + " \n", + " \n", + "
    \n", + "
    \n", + "
      \n", + "
      \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
      <xarray.Dataset>\n",
+       "Dimensions:      (chain: 4, draw: 1000, measure: 2, supermarket: 5)\n",
+       "Coordinates:\n",
+       "  * chain        (chain) int64 0 1 2 3\n",
+       "  * draw         (draw) int64 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999\n",
+       "  * measure      (measure) <U4 'time' 'cost'\n",
+       "  * supermarket  (supermarket) <U7 'aldi' 'caprabo' 'condis' 'consum' 'keisy'\n",
+       "Data variables:\n",
+       "    mu           (chain, draw, measure, supermarket) float64 7.076 ... 5.212\n",
+       "    sigma        (chain, draw, measure, supermarket) float64 0.2387 ... 14.65\n",
+       "    lam          (chain, draw, supermarket) float64 4.774 3.181 ... 3.683 2.748\n",
+       "Attributes:\n",
+       "    created_at:                 2022-12-15T03:24:06.590643\n",
+       "    arviz_version:              0.14.0\n",
+       "    inference_library:          pymc\n",
+       "    inference_library_version:  5.0.0\n",
+       "    sampling_time:              8.219048500061035\n",
+       "    tuning_steps:               1000

      \n", + "
    \n", + "
    \n", + "
    • \n", + " \n", + "
  • \n", + " \n", + " \n", + "
    \n", + "
    \n", + "
      \n", + "
      \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
      <xarray.Dataset>\n",
+       "Dimensions:                (chain: 4, draw: 1000)\n",
+       "Coordinates:\n",
+       "  * chain                  (chain) int64 0 1 2 3\n",
+       "  * draw                   (draw) int64 0 1 2 3 4 5 ... 994 995 996 997 998 999\n",
+       "Data variables: (12/17)\n",
+       "    largest_eigval         (chain, draw) float64 nan nan nan nan ... nan nan nan\n",
+       "    smallest_eigval        (chain, draw) float64 nan nan nan nan ... nan nan nan\n",
+       "    max_energy_error       (chain, draw) float64 -0.3631 1.565 ... 0.529 0.5611\n",
+       "    perf_counter_start     (chain, draw) float64 6.252e+04 ... 6.252e+04\n",
+       "    energy                 (chain, draw) float64 4.181e+03 ... 4.184e+03\n",
+       "    perf_counter_diff      (chain, draw) float64 0.001759 0.002019 ... 0.001741\n",
+       "    ...                     ...\n",
+       "    lp                     (chain, draw) float64 -4.172e+03 ... -4.171e+03\n",
+       "    reached_max_treedepth  (chain, draw) bool False False False ... False False\n",
+       "    step_size_bar          (chain, draw) float64 0.5795 0.5795 ... 0.5609 0.5609\n",
+       "    energy_error           (chain, draw) float64 0.1146 -0.1436 ... 0.1856\n",
+       "    step_size              (chain, draw) float64 0.7116 0.7116 ... 0.5416 0.5416\n",
+       "    diverging              (chain, draw) bool False False False ... False False\n",
+       "Attributes:\n",
+       "    created_at:                 2022-12-15T03:24:06.604476\n",
+       "    arviz_version:              0.14.0\n",
+       "    inference_library:          pymc\n",
+       "    inference_library_version:  5.0.0\n",
+       "    sampling_time:              8.219048500061035\n",
+       "    tuning_steps:               1000

      \n", + "
    \n", + "
    \n", + "
    • \n", + " \n", + "
  • \n", + " \n", + " \n", + "
    \n", + "
    \n", + "
      \n", + "
      \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
      <xarray.Dataset>\n",
+       "Dimensions:  (week: 52)\n",
+       "Coordinates:\n",
+       "  * week     (week) int64 0 1 2 3 4 5 6 7 8 9 ... 42 43 44 45 46 47 48 49 50 51\n",
+       "Data variables:\n",
+       "    time     (week) float64 5.855 5.007 5.248 5.042 ... 4.866 4.623 7.34 4.608\n",
+       "    cost     (week) float64 118.5 148.7 127.4 126.0 ... 128.6 132.0 89.79 129.8\n",
+       "    rubbish  (week) int64 8 2 4 4 4 2 6 1 2 3 5 5 4 ... 2 4 3 4 3 9 0 4 3 7 5 4\n",
+       "Attributes:\n",
+       "    created_at:                 2022-12-15T03:24:06.610243\n",
+       "    arviz_version:              0.14.0\n",
+       "    inference_library:          pymc\n",
+       "    inference_library_version:  5.0.0

      \n", + "
    \n", + "
    \n", + "
    • \n", + " \n", + "
  • \n", + " \n", + " \n", + "
    \n", + "
    \n", + "
      \n", + "
      \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
      <xarray.Dataset>\n",
+       "Dimensions:         (week: 52)\n",
+       "Coordinates:\n",
+       "  * week            (week) int64 0 1 2 3 4 5 6 7 8 ... 44 45 46 47 48 49 50 51\n",
+       "Data variables:\n",
+       "    supermarket_id  (week) int32 3 4 1 1 3 3 0 2 1 3 4 ... 2 1 0 3 0 3 1 2 2 0 2\n",
+       "Attributes:\n",
+       "    created_at:                 2022-12-15T03:24:06.611223\n",
+       "    arviz_version:              0.14.0\n",
+       "    inference_library:          pymc\n",
+       "    inference_library_version:  5.0.0

      \n", + "
    \n", + "
    \n", + "
    • \n", + " \n", + "
\n", + "
\n", + " " + ], + "text/plain": [ + "Inference data with groups:\n", + "\t> posterior\n", + "\t> sample_stats\n", + "\t> observed_data\n", + "\t> constant_data" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "idata" + ] + }, + { + "cell_type": "markdown", + "id": "80c6e95a", + "metadata": {}, + "source": [ + "## Posterior predictive sampling" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "717537fa", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Sampling: [cost, rubbish, time]\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "\n", + "
\n", + " \n", + " 100.00% [4000/4000 00:00<00:00]\n", + "
\n", + " " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "with model:\n", + " model.set_dim(\"week\", len(supermarkets))\n", + " model.set_data(\"supermarket_id\", np.arange(len(supermarkets)))\n", + "\n", + " pm.sample_posterior_predictive(idata, extend_inferencedata=True)\n", + "\n", + "# we can't rename dimensions in pymc, so we do it now\n", + "post_pred = idata.posterior_predictive\n", + "post_pred = post_pred.rename(week=\"supermarket\").assign_coords(\n", + " supermarket=idata.posterior[\"supermarket\"]\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "1f42d9c7", + "metadata": {}, + "source": [ + "## Decision analysis" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "93ea711d", + "metadata": {}, + "outputs": [], + "source": [ + "def utility(time, cost, rubbish):\n", + " return -(cost + 20 * time + 2 * (rubbish + 1) ** 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "0bd65d90", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray (supermarket: 5)>\n",
+       "array([-246.72746742, -145.02445363, -162.09484257, -180.34604534,\n",
+       "       -157.80864075])\n",
+       "Coordinates:\n",
+       "  * supermarket  (supermarket) <U7 'aldi' 'caprabo' 'condis' 'consum' 'keisy'
" + ], + "text/plain": [ + "\n", + "array([-246.72746742, -145.02445363, -162.09484257, -180.34604534,\n", + " -157.80864075])\n", + "Coordinates:\n", + " * supermarket (supermarket) " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "_, ax = plt.subplots(1, 2, figsize=(12, 6))\n", + "az.plot_density(\n", + " [util.sel(supermarket=s) for s in supermarkets],\n", + " data_labels=supermarkets,\n", + " labeller=az.labels.NoVarLabeller(),\n", + " ax=ax[0],\n", + ")\n", + "ax[0].set_title(\"Utility function **with** rubbish term\")\n", + "az.plot_density(\n", + " [util2.sel(supermarket=s) for s in supermarkets],\n", + " data_labels=supermarkets,\n", + " labeller=az.labels.NoVarLabeller(),\n", + " ax=ax[1],\n", + ")\n", + "ax[1].set_title(\"Utility function **without** rubbish term\")" + ] + }, + { + "cell_type": "markdown", + "id": "7fa6f2a9", + "metadata": {}, + "source": [ + "## Authors\n", + "* Authored by Oriol Abril Pla in Dec 2022\n", + "\n", + "## References\n", + "\n", + ":::{bibliography}\n", + ":filter: docname in docnames\n", + ":::\n", + "\n", + "## Watermark" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "6ec98064", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Last updated: Thu Dec 15 2022\n", + "\n", + "Python implementation: CPython\n", + "Python version : 3.10.8\n", + "IPython version : 8.7.0\n", + "\n", + "pytensor: 2.8.10\n", + "xarray : 2022.12.0\n", + "\n", + "arviz : 0.14.0\n", + "scipy : 1.9.3\n", + "matplotlib: 3.6.2\n", + "pandas : 1.5.2\n", + "pymc : 5.0.0\n", + "numpy : 1.23.5\n", + "\n", + "Watermark: 2.3.1\n", + "\n" + ] + } + ], + "source": [ + "%load_ext watermark\n", + "%watermark -n -u -v -iv -w -p pytensor,xarray" + ] + }, + { + "cell_type": "markdown", + "id": "1469825e", + "metadata": {}, + "source": [ + ":::{include} ../page_footer.md\n", + ":::" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "59a49c5e", + "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.10.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/myst_nbs/howto/decision_analysis.myst.md b/myst_nbs/howto/decision_analysis.myst.md new file mode 100644 index 000000000..86ea6eee4 --- /dev/null +++ b/myst_nbs/howto/decision_analysis.myst.md @@ -0,0 +1,181 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.13.7 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Decision analysis with a utility function + +:::{post} Dec 15, 2022 +:tags: Bayesian workflow, decision analysis +:category: beginner, how-to +:author: Oriol Abril Pla +::: + +```{code-cell} ipython3 +import arviz as az +import numpy as np +import pandas as pd +import pymc as pm + +from scipy import stats +from xarray_einstats.stats import XrContinuousRV, XrDiscreteRV +``` + +```{code-cell} ipython3 +RANDOM_SEED = 4235 +rng = np.random.default_rng(RANDOM_SEED) +``` + +## Overview + +In this example we will show how to perform Bayesian decison analysis over a set of discrete decisions. We will use as an example, we will take the role of a shopper, who has done their weekly grocery shopping at a different supermarket every week for the last year, and has logged: + +* The supermarket they went to that week +* The cumulative time it took from exiting home to getting back, so both travel time and unintuitive supermarket layout are penalized (in hours) +* The money spent (in euros) +* The number of rubbish bins that were filled that week + ++++ + +## Example data + +```{code-cell} ipython3 +obs_supermarket = rng.choice(["aldi", "caprabo", "condis", "consum", "keisy"], size=52) +super_id, supermarkets = pd.factorize(obs_supermarket, sort=True) +time_means = np.array([7, 5, 4.7, 6, 5.2]) +cost_means = np.array([90, 125, 130, 120, 150]) +rubbish_means = np.array([5, 3, 4, 4, 3]) +obs_time = rng.normal(time_means[super_id], 0.2) +obs_cost = rng.normal(cost_means[super_id], 1) +obs_rubbish = rng.poisson(rubbish_means[super_id]) +``` + +## Model definition + +```{code-cell} ipython3 +coords = { + "supermarket": supermarkets, + "measure": ["time", "cost"], +} +with pm.Model(coords=coords) as model: + model.add_coord("week", length=len(super_id), mutable=True) + supermarket_id = pm.MutableData("supermarket_id", super_id, dims="week") + + mu = pm.Normal("mu", 0, 1, dims=("measure", "supermarket")) + sigma = pm.HalfNormal("sigma", 0.5, dims=("measure", "supermarket")) + lam = pm.Exponential("lam", 0.3, dims="supermarket") + + pm.Normal( + "time", + mu=mu[0, supermarket_id], + sigma=sigma[0, supermarket_id], + dims="week", + observed=obs_time, + ) + pm.Normal( + "cost", + mu=mu[1, supermarket_id], + sigma=sigma[1, supermarket_id], + dims="week", + observed=obs_cost, + ) + pm.Poisson("rubbish", mu=lam[supermarket_id], dims="week", observed=obs_rubbish) +``` + +## Inference + +```{code-cell} ipython3 +with model: + idata = pm.sample() +``` + +```{code-cell} ipython3 +idata +``` + +## Posterior predictive sampling + +```{code-cell} ipython3 +with model: + model.set_dim("week", len(supermarkets)) + model.set_data("supermarket_id", np.arange(len(supermarkets))) + + pm.sample_posterior_predictive(idata, extend_inferencedata=True) + +# we can't rename dimensions in pymc, so we do it now +post_pred = idata.posterior_predictive +post_pred = post_pred.rename(week="supermarket").assign_coords( + supermarket=idata.posterior["supermarket"] +) +``` + +## Decision analysis + +```{code-cell} ipython3 +def utility(time, cost, rubbish): + return -(cost + 20 * time + 2 * (rubbish + 1) ** 2) +``` + +```{code-cell} ipython3 +util = utility(post_pred["time"], post_pred["cost"], post_pred["rubbish"]) +util.mean(dim=("chain", "draw")) +``` + +```{code-cell} ipython3 +def utility_no_rubbish(time, cost, rubbish): + return -(cost + 30 * time) + + +util2 = utility_no_rubbish(post_pred["time"], post_pred["cost"], post_pred["rubbish"]) +``` + +```{code-cell} ipython3 +import matplotlib.pyplot as plt + +_, ax = plt.subplots(1, 2, figsize=(12, 6)) +az.plot_density( + [util.sel(supermarket=s) for s in supermarkets], + data_labels=supermarkets, + labeller=az.labels.NoVarLabeller(), + ax=ax[0], +) +ax[0].set_title("Utility function **with** rubbish term") +az.plot_density( + [util2.sel(supermarket=s) for s in supermarkets], + data_labels=supermarkets, + labeller=az.labels.NoVarLabeller(), + ax=ax[1], +) +ax[1].set_title("Utility function **without** rubbish term") +``` + +## Authors +* Authored by Oriol Abril Pla in Dec 2022 + +## References + +:::{bibliography} +:filter: docname in docnames +::: + +## Watermark + +```{code-cell} ipython3 +%load_ext watermark +%watermark -n -u -v -iv -w -p pytensor,xarray +``` + +:::{include} ../page_footer.md +::: + +```{code-cell} ipython3 + +```