diff --git a/examples/20231102_multinom.ipynb b/examples/20231102_multinom.ipynb deleted file mode 100644 index dbcc329..0000000 --- a/examples/20231102_multinom.ipynb +++ /dev/null @@ -1,280 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "5053db7f-e3d3-4e5b-923f-ec6554fe6466", - "metadata": {}, - "outputs": [], - "source": [ - "import pandas as pd\n", - "import pymc as pm\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import seaborn as sns\n", - "import numpy as np\n", - "import arviz as az\n", - "import statsmodels.api as sm\n", - "\n", - "import matplotlib.dates as mdates\n", - "import matplotlib.ticker as ticker\n", - "import matplotlib.cm as cm\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "0465678d-155c-4202-8578-382e8773a51f", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Current working directory: /Users/gneiss/Documents/School/CBG2/brexit/TSmodelling/covvfit/examples\n" - ] - } - ], - "source": [ - "import os\n", - "cwd = os.getcwd()\n", - "print(\"Current working directory:\", cwd)\n" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "05b32867-adff-441a-8dd3-34cd1b8ba3b3", - "metadata": {}, - "outputs": [], - "source": [ - "import sys\n", - "\n", - "# Modify this path to match your directory structure\n", - "path_to_covvfit = os.path.abspath(os.path.join(cwd, '../src'))\n", - "sys.path.append(path_to_covvfit)\n", - "\n", - "# Now you can import your modules\n", - "import covvfit as cv # Replace with the actual name of your Python file without .py\n" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "id": "82522dae-866b-445e-a2a2-982c182bd78b", - "metadata": {}, - "outputs": [], - "source": [ - "# import importlib\n", - "# importlib.reload(cv)\n", - "# this does not work" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "0a168eb1-2a4d-467d-b0d2-975c866444c0", - "metadata": {}, - "outputs": [], - "source": [ - "variants = [\n", - "# 'B.1.1.7', 'B.1.351', 'P.1', 'undetermined',\n", - " 'B.1.617.2', 'BA.1', 'BA.2', 'BA.4', 'BA.5', 'BA.2.75',\n", - " 'BQ.1.1', 'XBB.1.5', 'XBB.1.9', 'XBB.1.16', 'XBB.2.3', 'EG.5', \"BA.2.86\"\n", - "]\n", - "\n", - "cities = ['Lugano (TI)', 'Zürich (ZH)', 'Chur (GR)', 'Altenrhein (SG)',\n", - " 'Laupen (BE)', 'Genève (GE)', 'Basel (BS)', 'Porrentruy (JU)',\n", - " 'Lausanne (VD)', 'Bern (BE)', 'Luzern (LU)', 'Solothurn (SO)',\n", - " 'Neuchâtel (NE)', 'Schwyz (SZ)']" - ] - }, - { - "cell_type": "code", - "execution_count": 17, - "id": "16c0fdb0-d3f4-430b-bfa5-19be50154d53", - "metadata": {}, - "outputs": [], - "source": [ - "data = cv.load_data('../data/robust_deconv2_noisy13.csv')\n", - "data2 = cv.preprocess_df(data, cities, variants, date_min='2021-11-01')\n", - "\n", - "ts_lst, ys_lst = cv.make_data_list(data2, cities, variants)" - ] - }, - { - "cell_type": "code", - "execution_count": 30, - "id": "d08021ba-e4e9-4a98-b6e2-393786493894", - "metadata": {}, - "outputs": [], - "source": [ - "from pymc.distributions.dist_math import factln\n", - "\n", - "# model for just one city\n", - "def create_model5(\n", - " ts_lst,\n", - " ys_lst,\n", - " n=1.0,\n", - " coords={\n", - "# \"cities\":cities,\n", - " \"variants\":variants,\n", - " },\n", - " n_pred=60\n", - "):\n", - " ts_pred = np.arange(n_pred) + ts_lst.max()\n", - " with pm.Model(coords=coords) as model:\n", - "# sigma_var = pm.InverseGamma(\"sigma\", alpha=2.1, beta=0.015, dims=[\"cities\",\"variants\"])\n", - " midpoint_var = pm.Normal(\"midpoint\", mu=0.0, sigma=500.0, dims=\"variants\")\n", - "# midpoint_sig = pm.InverseGamma(\"midpoint_sig\", alpha=7.0, beta=60.0)\n", - " rate_var = pm.Gamma(\"rate\", mu=0.15, sigma=0.1, dims=\"variants\")\n", - "# rate_sig = pm.InverseGamma(\"rate_sigma\", alpha=2.0005, beta=0.05)\n", - " n_eff_inv = pm.InverseGamma(\"n_eff_inv\", alpha=20.0, beta=2.0) \n", - " n_eff = pm.Deterministic(\"n_eff\", 1/n_eff_inv)\n", - "# n_eff = pm.TruncatedNormal(\"n_eff\", mu=10, sigma=10, lower=1.0)\n", - "# n_eff = pm.Gamma(\"n_eff\", alpha=1000, beta=100)\n", - " \n", - " # Kaan's trick to avoid overflows\n", - " def softmax(x, rates, midpoints):\n", - " E = rates[:, None] * (x - midpoints[:, None])\n", - " E_max = E.max(axis=0)\n", - " un_norm = pm.math.exp(E - E_max)\n", - " return un_norm / (pm.math.sum(un_norm, axis=0))\n", - " \n", - " ys_smooth = pm.Deterministic(f\"ys_ideal\",softmax(ts_lst, rate_var, midpoint_var), dims=\"variants\")\n", - " ys_pred = pm.Deterministic(f\"ys_pred\",softmax(ts_pred, rate_var, midpoint_var), dims=\"variants\")\n", - "# ys_wiggly = pm.Beta(f\"ys_wiggly\", mu=ys_smooth, nu=n_eff)\n", - " \n", - " # make Multinom/n likelihood\n", - " def log_likelihood(y, p, n):\n", - " return n*pm.math.sum(y * pm.math.log(p) - factln(n*y), axis=0) + pm.math.log(n) + factln(n)\n", - "\n", - " ys_noisy = pm.DensityDist(\n", - " f\"ys_noisy\",\n", - " ys_smooth,\n", - " n_eff,\n", - " logp=log_likelihood,\n", - " observed=ys_lst,\n", - " )\n", - " \n", - " return model" - ] - }, - { - "cell_type": "code", - "execution_count": 18, - "id": "d5047fa2-701e-4f7a-811c-43f865df6e1b", - "metadata": {}, - "outputs": [], - "source": [ - "ys = ys_lst[1]\n", - "ys = ys / ys.sum(0)\n", - "ts = ts_lst[1]" - ] - }, - { - "cell_type": "code", - "execution_count": 28, - "id": "c8f5c012-84fc-43e0-b230-95dcbe3e8167", - "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 4 jobs)\n", - "NUTS: [sigma, midpoint, rate]\n" - ] - }, - { - "data": { - "text/html": [ - "\n", - "\n" - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/html": [ - "\n", - "
\n", - " \n", - " 24.15% [483/2000 03:05<09:42 Sampling 2 chains, 0 divergences]\n", - "
\n", - " " - ], - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "ename": "ValueError", - "evalue": "Not enough samples to build a trace.", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[28], line 4\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m create_model(ts, ys, coords\u001b[38;5;241m=\u001b[39m{\n\u001b[1;32m 2\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mvariants\u001b[39m\u001b[38;5;124m\"\u001b[39m:variants,\n\u001b[1;32m 3\u001b[0m }):\n\u001b[0;32m----> 4\u001b[0m idata_posterior \u001b[38;5;241m=\u001b[39m pm\u001b[38;5;241m.\u001b[39msample(random_seed\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m66\u001b[39m, chains\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m2\u001b[39m, tune\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m500\u001b[39m, draws\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m500\u001b[39m)\n", - "File \u001b[0;32m~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:702\u001b[0m, in \u001b[0;36msample\u001b[0;34m(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs)\u001b[0m\n\u001b[1;32m 698\u001b[0m t_sampling \u001b[38;5;241m=\u001b[39m time\u001b[38;5;241m.\u001b[39mtime() \u001b[38;5;241m-\u001b[39m t_start\n\u001b[1;32m 700\u001b[0m \u001b[38;5;66;03m# Packaging, validating and returning the result was extracted\u001b[39;00m\n\u001b[1;32m 701\u001b[0m \u001b[38;5;66;03m# into a function to make it easier to test and refactor.\u001b[39;00m\n\u001b[0;32m--> 702\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m _sample_return(\n\u001b[1;32m 703\u001b[0m run\u001b[38;5;241m=\u001b[39mrun,\n\u001b[1;32m 704\u001b[0m traces\u001b[38;5;241m=\u001b[39mtraces,\n\u001b[1;32m 705\u001b[0m tune\u001b[38;5;241m=\u001b[39mtune,\n\u001b[1;32m 706\u001b[0m t_sampling\u001b[38;5;241m=\u001b[39mt_sampling,\n\u001b[1;32m 707\u001b[0m discard_tuned_samples\u001b[38;5;241m=\u001b[39mdiscard_tuned_samples,\n\u001b[1;32m 708\u001b[0m compute_convergence_checks\u001b[38;5;241m=\u001b[39mcompute_convergence_checks,\n\u001b[1;32m 709\u001b[0m return_inferencedata\u001b[38;5;241m=\u001b[39mreturn_inferencedata,\n\u001b[1;32m 710\u001b[0m keep_warning_stat\u001b[38;5;241m=\u001b[39mkeep_warning_stat,\n\u001b[1;32m 711\u001b[0m idata_kwargs\u001b[38;5;241m=\u001b[39midata_kwargs \u001b[38;5;129;01mor\u001b[39;00m {},\n\u001b[1;32m 712\u001b[0m model\u001b[38;5;241m=\u001b[39mmodel,\n\u001b[1;32m 713\u001b[0m )\n", - "File \u001b[0;32m~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:733\u001b[0m, in \u001b[0;36m_sample_return\u001b[0;34m(run, traces, tune, t_sampling, discard_tuned_samples, compute_convergence_checks, return_inferencedata, keep_warning_stat, idata_kwargs, model)\u001b[0m\n\u001b[1;32m 731\u001b[0m \u001b[38;5;66;03m# Pick and slice chains to keep the maximum number of samples\u001b[39;00m\n\u001b[1;32m 732\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m discard_tuned_samples:\n\u001b[0;32m--> 733\u001b[0m traces, length \u001b[38;5;241m=\u001b[39m _choose_chains(traces, tune)\n\u001b[1;32m 734\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m 735\u001b[0m traces, length \u001b[38;5;241m=\u001b[39m _choose_chains(traces, \u001b[38;5;241m0\u001b[39m)\n", - "File \u001b[0;32m~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/backends/base.py:601\u001b[0m, in \u001b[0;36m_choose_chains\u001b[0;34m(traces, tune)\u001b[0m\n\u001b[1;32m 599\u001b[0m lengths \u001b[38;5;241m=\u001b[39m [\u001b[38;5;28mmax\u001b[39m(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;28mlen\u001b[39m(trace) \u001b[38;5;241m-\u001b[39m tune) \u001b[38;5;28;01mfor\u001b[39;00m trace \u001b[38;5;129;01min\u001b[39;00m traces]\n\u001b[1;32m 600\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28msum\u001b[39m(lengths):\n\u001b[0;32m--> 601\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNot enough samples to build a trace.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 603\u001b[0m idxs \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39margsort(lengths)\n\u001b[1;32m 604\u001b[0m l_sort \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray(lengths)[idxs]\n", - "\u001b[0;31mValueError\u001b[0m: Not enough samples to build a trace." - ] - } - ], - "source": [ - "with create_model(ts, ys, coords={\n", - " \"variants\":variants,\n", - " }):\n", - " idata_posterior = pm.sample(random_seed=66, chains=2, tune=500, draws=500)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "pymc_env", - "language": "python", - "name": "pymc_env" - }, - "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.3" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/examples/20231102_multinom.qmd b/examples/20231102_multinom.qmd new file mode 100644 index 0000000..fa1d0c8 --- /dev/null +++ b/examples/20231102_multinom.qmd @@ -0,0 +1,109 @@ +## Experimentation with multinomial model + +```{python} +import pandas as pd +import pymc as pm + +import matplotlib.pyplot as plt +import seaborn as sns +import numpy as np +import arviz as az +import statsmodels.api as sm + +import matplotlib.dates as mdates +import matplotlib.ticker as ticker +import matplotlib.cm as cm + +import covvfit as cv +``` + +Load the data: +```{python} +data_path = '../private/data/robust_deconv2_noisy13.csv' + +variants = [ +# 'B.1.1.7', 'B.1.351', 'P.1', 'undetermined', + 'B.1.617.2', 'BA.1', 'BA.2', 'BA.4', 'BA.5', 'BA.2.75', + 'BQ.1.1', 'XBB.1.5', 'XBB.1.9', 'XBB.1.16', 'XBB.2.3', 'EG.5', "BA.2.86" +] + +cities = ['Lugano (TI)', 'Zürich (ZH)', 'Chur (GR)', 'Altenrhein (SG)', + 'Laupen (BE)', 'Genève (GE)', 'Basel (BS)', 'Porrentruy (JU)', + 'Lausanne (VD)', 'Bern (BE)', 'Luzern (LU)', 'Solothurn (SO)', + 'Neuchâtel (NE)', 'Schwyz (SZ)'] + + +data = cv.load_data(data_path) +data2 = cv.preprocess_df(data, cities, variants, date_min='2021-11-01') + +ts_lst, ys_lst = cv.make_data_list(data2, cities, variants) +``` + + +Let's load one city only: +```{python} +ys = ys_lst[1] +ys = ys / ys.sum(0) +ts = ts_lst[1] +``` + +Now we can create model for this one city: +```{python} +from pymc.distributions.dist_math import factln + +# model for just one city +def create_model5( + ts_lst, + ys_lst, + n=1.0, + coords={ +# "cities":cities, + "variants":variants, + }, + n_pred=60 +): + ts_pred = np.arange(n_pred) + ts_lst.max() + with pm.Model(coords=coords) as model: +# sigma_var = pm.InverseGamma("sigma", alpha=2.1, beta=0.015, dims=["cities","variants"]) + midpoint_var = pm.Normal("midpoint", mu=0.0, sigma=500.0, dims="variants") +# midpoint_sig = pm.InverseGamma("midpoint_sig", alpha=7.0, beta=60.0) + rate_var = pm.Gamma("rate", mu=0.15, sigma=0.1, dims="variants") +# rate_sig = pm.InverseGamma("rate_sigma", alpha=2.0005, beta=0.05) + n_eff_inv = pm.InverseGamma("n_eff_inv", alpha=20.0, beta=2.0) + n_eff = pm.Deterministic("n_eff", 1/n_eff_inv) +# n_eff = pm.TruncatedNormal("n_eff", mu=10, sigma=10, lower=1.0) +# n_eff = pm.Gamma("n_eff", alpha=1000, beta=100) + + # Kaan's trick to avoid overflows + def softmax(x, rates, midpoints): + E = rates[:, None] * (x - midpoints[:, None]) + E_max = E.max(axis=0) + un_norm = pm.math.exp(E - E_max) + return un_norm / (pm.math.sum(un_norm, axis=0)) + + ys_smooth = pm.Deterministic(f"ys_ideal",softmax(ts_lst, rate_var, midpoint_var), dims="variants") + ys_pred = pm.Deterministic(f"ys_pred",softmax(ts_pred, rate_var, midpoint_var), dims="variants") +# ys_wiggly = pm.Beta(f"ys_wiggly", mu=ys_smooth, nu=n_eff) + + # make Multinom/n likelihood + def log_likelihood(y, p, n): + return n*pm.math.sum(y * pm.math.log(p) - factln(n*y), axis=0) + pm.math.log(n) + factln(n) + + ys_noisy = pm.DensityDist( + f"ys_noisy", + ys_smooth, + n_eff, + logp=log_likelihood, + observed=ys_lst, + ) + + return model + + + +with create_model(ts, ys, coords={ + "variants":variants, + }): + idata_posterior = pm.sample(random_seed=65, chains=2, tune=500, draws=500) +``` +