From 1e69861060115073aff6a3d8930d1c1b96d28597 Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Sun, 18 Nov 2018 18:31:06 -0500 Subject: [PATCH 01/11] Transmitter up --- .idea/causal_demon.iml | 14 ++ .idea/inspectionProfiles/Project_Default.xml | 14 ++ .idea/libraries/R_User_Library.xml | 6 + .idea/misc.xml | 4 + .idea/modules.xml | 9 + .idea/vcs.xml | 6 + .idea/workspace.xml | 230 +++++++++++++++++++ mapk_demon/__init__.py | 0 mapk_demon/inference.py | 0 mapk_demon/receivers.py | 0 mapk_demon/transmitters.py | 92 ++++++++ 11 files changed, 375 insertions(+) create mode 100644 .idea/causal_demon.iml create mode 100644 .idea/inspectionProfiles/Project_Default.xml create mode 100644 .idea/libraries/R_User_Library.xml create mode 100644 .idea/misc.xml create mode 100644 .idea/modules.xml create mode 100644 .idea/vcs.xml create mode 100644 .idea/workspace.xml create mode 100644 mapk_demon/__init__.py create mode 100644 mapk_demon/inference.py create mode 100644 mapk_demon/receivers.py create mode 100644 mapk_demon/transmitters.py diff --git a/.idea/causal_demon.iml b/.idea/causal_demon.iml new file mode 100644 index 0000000..c1866a5 --- /dev/null +++ b/.idea/causal_demon.iml @@ -0,0 +1,14 @@ + + + + + + + + + + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/Project_Default.xml b/.idea/inspectionProfiles/Project_Default.xml new file mode 100644 index 0000000..f424b52 --- /dev/null +++ b/.idea/inspectionProfiles/Project_Default.xml @@ -0,0 +1,14 @@ + + + + \ No newline at end of file diff --git a/.idea/libraries/R_User_Library.xml b/.idea/libraries/R_User_Library.xml new file mode 100644 index 0000000..71f5ff7 --- /dev/null +++ b/.idea/libraries/R_User_Library.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..b54da92 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,4 @@ + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..8fd0296 --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,9 @@ + + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 0000000..35eb1dd --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/.idea/workspace.xml b/.idea/workspace.xml new file mode 100644 index 0000000..b052bbc --- /dev/null +++ b/.idea/workspace.xml @@ -0,0 +1,230 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Python + + + + + PyPep8Inspection + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 1542574426773 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - \ No newline at end of file From b05f643fc52c16b844ae80c0e491fce96ec3d5ee Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Sun, 18 Nov 2018 18:41:44 -0500 Subject: [PATCH 03/11] gitignored --- .idea/causal_demon.iml | 14 -------------- 1 file changed, 14 deletions(-) delete mode 100644 .idea/causal_demon.iml diff --git a/.idea/causal_demon.iml b/.idea/causal_demon.iml deleted file mode 100644 index c1866a5..0000000 --- a/.idea/causal_demon.iml +++ /dev/null @@ -1,14 +0,0 @@ - - - - - - - - - - - - - \ No newline at end of file From aea22184b835c52fad2d037ab754617bee65ea93 Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Sun, 18 Nov 2018 18:42:25 -0500 Subject: [PATCH 04/11] gitignored --- .gitignore | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index f8ab63a..a67167b 100644 --- a/.gitignore +++ b/.gitignore @@ -23,7 +23,6 @@ wheels/ *.egg-info/ .installed.cfg *.egg -*.xml MANIFEST # PyInstaller @@ -98,4 +97,8 @@ venv.bak/ *.pdf # avoid uploading data -notebooks/data/* \ No newline at end of file +notebooks/data/* + +# avoid PyCharm files +*.iml +*.xml \ No newline at end of file From 04c7351f4dd66c40405f847065be14b45faafa56 Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Sun, 18 Nov 2018 19:20:18 -0500 Subject: [PATCH 05/11] gitignore updated --- mapk_demon/transmitters.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mapk_demon/transmitters.py b/mapk_demon/transmitters.py index 24a3007..9048999 100644 --- a/mapk_demon/transmitters.py +++ b/mapk_demon/transmitters.py @@ -89,4 +89,4 @@ def mapk_signaling(noise_dists): noise_vars = ['N_3k', 'N_2k', 'N_k'] noise_prior = {N: LogNormal(0, 10) for N in noise_vars} - print(mapk_signaling(noise_dists=noise_prior)) \ No newline at end of file + print(mapk_signaling(noise_dists=noise_prior)) From 6084917ca6350f76a72f9e8e624b510a441aa2c9 Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Thu, 13 Dec 2018 14:00:50 -0500 Subject: [PATCH 06/11] Inference Done --- .gitignore | 66 ++++++++++++++- mapk_demon/receivers.py | 160 +++++++++++++++++++++++++++++++++++++ mapk_demon/transmitters.py | 27 ++++--- 3 files changed, 239 insertions(+), 14 deletions(-) diff --git a/.gitignore b/.gitignore index a67167b..49c06e6 100644 --- a/.gitignore +++ b/.gitignore @@ -101,4 +101,68 @@ notebooks/data/* # avoid PyCharm files *.iml -*.xml \ No newline at end of file +*.xml + +# User-specific stuff +.idea/**/workspace.xml +.idea/**/tasks.xml +.idea/**/usage.statistics.xml +.idea/**/dictionaries +.idea/**/shelf + +# Generated files +.idea/**/contentModel.xml + +# Sensitive or high-churn files +.idea/**/dataSources/ +.idea/**/dataSources.ids +.idea/**/dataSources.local.xml +.idea/**/sqlDataSources.xml +.idea/**/dynamic.xml +.idea/**/uiDesigner.xml +.idea/**/dbnavigator.xml + +# Gradle +.idea/**/gradle.xml +.idea/**/libraries + +# Gradle and Maven with auto-import +# When using Gradle or Maven with auto-import, you should exclude module files, +# since they will be recreated, and may cause churn. Uncomment if using +# auto-import. +# .idea/modules.xml +# .idea/*.iml +# .idea/modules + +# CMake +cmake-build-*/ + +# Mongo Explorer plugin +.idea/**/mongoSettings.xml + +# File-based project format +*.iws + +# IntelliJ +out/ + +# mpeltonen/sbt-idea plugin +.idea_modules/ + +# JIRA plugin +atlassian-ide-plugin.xml + +# Cursive Clojure plugin +.idea/replstate.xml + +# Crashlytics plugin (for Android Studio and IntelliJ) +com_crashlytics_export_strings.xml +crashlytics.properties +crashlytics-build.properties +fabric.properties + +# Editor-based Rest Client +.idea/httpRequests + +# Android studio 3.1+ serialized cache file +.idea/caches/build_file_checksums.ser \ No newline at end of file diff --git a/mapk_demon/receivers.py b/mapk_demon/receivers.py index e69de29..86cf0b6 100644 --- a/mapk_demon/receivers.py +++ b/mapk_demon/receivers.py @@ -0,0 +1,160 @@ +from __future__ import division + +from functools import partial + +import numpy as np +import pyro +from pyro import sample +from pyro.distributions import Normal, Uniform +from torch import tensor +from mapk_demon.transmitters import mapk_signaling +from pyro.infer.mcmc import MCMC +from pyro.infer.mcmc.nuts import NUTS +from pyro.infer import EmpiricalMarginal +from causal_demon.inference import infer_dist +import torch +import matplotlib.pyplot as plt + +def infer_mcmc(model): + nuts_kernel = NUTS(model, adapt_step_size=True) + return MCMC(nuts_kernel, num_samples=500, warmup_steps=300) + + +def g1(a): + return a / (a + 1) + + +def g2(a): + return a ** 2 / (a ** 2 + a + 1) + + +def f(mu, N): + return N * mu ** (0.5) + mu + + +# Parameters + +params = { + "alpha_3k": Normal(500, 1), + "alpha_2k": Normal(500, 1), + "alpha_k": Normal(500, 1), + "nu_3k": Normal(.15, .05), + "nu_2k": Normal(.15, .05), + "nu_k": Normal(.15, .05) +} +''' +params = { + "alpha_3k": 500, + "alpha_2k": 500, + "alpha_k": 500, + "nu_3k": .15, + "nu_2k": .15, + "nu_k": .15 +} +''' + +hyperparameters = { + "t_3k": 1.2, + "t_2k": 1.2, + "t_k": 0.003 +} + + +def f_map3k(E1, N_3k, params): + total_3k = hyperparameters['t_3k'] + alpha_3k = sample("alpha_3k",params['alpha_3k']) + nu_3k = sample("nu_3k", params['nu_3k']) + + map3k_mu = total_3k * g1(E1 * (alpha_3k/nu_3k)) + + map3k = sample("map3k", Normal(f(map3k_mu, N_3k), 1e-6)) + return map3k + + +def f_map2k(map3k, N_2k, params): + total_2k = hyperparameters['t_2k'] + alpha_2k = sample("alpha_2k", params['alpha_2k']) + nu_2k = sample("nu_2k", params['nu_2k']) + + map2k_mu = total_2k * g2(map3k * (alpha_2k / nu_2k)) + + map2k = sample("map2k", Normal(f(map2k_mu, N_2k), 1e-6)) + return map2k + + +def f_mapk(map2k, N_k, params): + total_k = hyperparameters['t_k'] + alpha_k = sample("alpha_k",params['alpha_k']) + nu_k = sample("nu_k", params['nu_k']) + + mapk_mu = total_k * g2(map2k * (alpha_k / nu_k)) + mapk = sample("mapk", Normal(f(mapk_mu, N_k), 1e-6)) + return mapk + + +def mapk_signaling_receiver(noise_dists): + + with pyro.iarange("model"): + E1 = sample('E1', Uniform(1.5e-5, 6e-5)) + N_3k = sample('N_3k', noise_dists['N_3k']) + N_2k = sample('N_2k', noise_dists['N_2k']) + N_k = sample('N_k', noise_dists['N_k']) + + map3k = f_map3k(E1, N_3k, params) + map2k = f_map2k(map3k,N_2k, params) + mapk = f_mapk(map2k, N_k, params) + + return { + 'map3k': map3k, + 'map2k': map2k, + 'mapk': mapk + } + +def plot(data): + plt.boxplot(data['erk']) + plt.title("Empirical Distribution of Erk") + plt.xlabel("Condition") + plt.ylabel("Concentration") + plt.show() + + +if __name__ == "__main__": + + # Noise Distributions + noise_vars = ['N_3k', 'N_2k', 'N_k'] + noise_priors = {N: Normal(0, 1) for N in noise_vars} + + parameters = ["alpha_3k","alpha_2k","alpha_k","nu_3k","nu_2k","nu_k"] + + receiver_model = mapk_signaling_receiver(noise_dists=noise_priors) + + # Simulate data from transmitter model + evidence = [] + + for i in range(10): + evidence.append(mapk_signaling(noise_dists=noise_priors)) + + # condition the receiver model with data + conditioned = pyro.condition(mapk_signaling_receiver, data = evidence) + + print(conditioned(noise_priors)) + + # mcmc inference + #nuts_kernel = NUTS(conditioned, adapt_step_size=True) + #mcmc_run = MCMC(nuts_kernel, num_samples=50, warmup_steps=10).run(noise_priors) + + # Importance Sampling + posterior = pyro.infer.Importance(conditioned, num_samples=100).run(noise_priors) + + #print(posterior) + posteriors = { + n: EmpiricalMarginal(posterior, sites=n) + for n in parameters + } + + print("Alpha_3k inference") + alpha_3k_marginal = posteriors['alpha_3k'] + alpha_3k_samples = [alpha_3k_marginal().item() for _ in range(500)] + + print("mean:",np.mean(alpha_3k_samples)) + print("std:",np.std(alpha_3k_samples)) diff --git a/mapk_demon/transmitters.py b/mapk_demon/transmitters.py index 9048999..24b2149 100644 --- a/mapk_demon/transmitters.py +++ b/mapk_demon/transmitters.py @@ -2,8 +2,9 @@ from functools import partial +import pyro from pyro import sample -from pyro.distributions import LogNormal +from pyro.distributions import Normal, Uniform from torch import tensor @@ -12,14 +13,14 @@ def g1(a): def g2(a): - return a**2 / (a**2 + a + 1) + return a ** 2 / (a ** 2 + a + 1) def f(mu, N): - return N ** (0.5) * mu + mu + return N * mu ** (0.5) + mu -constants = { +hyperparameters = { "t_3k": 1.2, "t_2k": 1.2, "t_k": 0.003 @@ -36,7 +37,7 @@ def f(mu, N): def f_map3k(E1, N_3k): - total_3k = constants['t_3k'] + total_3k = hyperparameters['t_3k'] alpha_3k = params['alpha_3k'] nu_3k = params['nu_3k'] @@ -46,7 +47,7 @@ def f_map3k(E1, N_3k): def f_map2k(map3k, N_2k): - total_2k = constants['t_2k'] + total_2k = hyperparameters['t_2k'] alpha_2k = params['alpha_2k'] nu_2k = params['nu_2k'] @@ -56,7 +57,7 @@ def f_map2k(map3k, N_2k): def f_mapk(map2k, N_k): - total_k = constants['t_k'] + total_k = hyperparameters['t_k'] alpha_k = params['alpha_k'] nu_k = params['nu_k'] @@ -67,10 +68,11 @@ def f_mapk(map2k, N_k): def mapk_signaling(noise_dists): - E1 = 1. - N_3k = sample('N_3k', noise_dists['N_3k']) - N_2k = sample('N_2k', noise_dists['N_2k']) - N_k = sample('N_k', noise_dists['N_k']) + with pyro.iarange("model"): + E1 = sample('E1', Uniform(1.5e-5, 6e-5)) + N_3k = sample('N_3k', noise_dists['N_3k']) + N_2k = sample('N_2k', noise_dists['N_2k']) + N_k = sample('N_k', noise_dists['N_k']) map3k = f_map3k(E1, N_3k) map2k = f_map2k(map3k,N_2k) @@ -82,11 +84,10 @@ def mapk_signaling(noise_dists): 'mapk': mapk } - if __name__ == "__main__": # Noise Distributions noise_vars = ['N_3k', 'N_2k', 'N_k'] - noise_prior = {N: LogNormal(0, 10) for N in noise_vars} + noise_prior = {N: Normal(0, 1) for N in noise_vars} print(mapk_signaling(noise_dists=noise_prior)) From 2d9410f1896e7df19ba67e5a2ef7a90270132a84 Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Fri, 15 Feb 2019 13:48:35 -0500 Subject: [PATCH 07/11] mapk inference working --- mapk_demon/svi_test.py | 11 + notebooks/Mapk Counterfactual Inference.ipynb | 214 ++++++++++++++++++ run_mapk.py | 11 + 3 files changed, 236 insertions(+) create mode 100644 mapk_demon/svi_test.py create mode 100644 notebooks/Mapk Counterfactual Inference.ipynb create mode 100644 run_mapk.py diff --git a/mapk_demon/svi_test.py b/mapk_demon/svi_test.py new file mode 100644 index 0000000..edc1699 --- /dev/null +++ b/mapk_demon/svi_test.py @@ -0,0 +1,11 @@ +''' +svi_test +Author: Kaushal Paneri +Project: causal_demon +Date of Creation: 1/2/19 +''' + +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals \ No newline at end of file diff --git a/notebooks/Mapk Counterfactual Inference.ipynb b/notebooks/Mapk Counterfactual Inference.ipynb new file mode 100644 index 0000000..d9c7ea1 --- /dev/null +++ b/notebooks/Mapk Counterfactual Inference.ipynb @@ -0,0 +1,214 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Mapk Counterfactual Inference" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.path.append('../')\n", + "\n", + "import pyro\n", + "from pyro import sample\n", + "from pyro.distributions import Normal, Uniform\n", + "import torch\n", + "from torch.distributions import constraints\n", + "from tqdm import tqdm_notebook as tqdm\n", + "import pyro.optim\n", + "\n", + "from pyro.infer.mcmc import MCMC\n", + "from pyro.infer.mcmc.nuts import NUTS, HMC\n", + "from pyro.infer import EmpiricalMarginal\n", + "\n", + "from mapk_demon.inference import infer_dist\n", + "from mapk_demon.receivers import mapk_receiver\n", + "\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "\n", + "pyro.set_rng_seed(101)\n", + "sample_size = 1000\n", + "\n", + "def hist(marginal, name):\n", + " plt.hist([marginal() for _ in range(5000)])\n", + " plt.title(\"Marginal Histogram of {}\".format(name))\n", + " plt.xlabel(\"concentration\")\n", + " plt.ylabel(\"#\")" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'map2k': tensor(12.1773), 'map3k': tensor(11.9624), 'mapk': tensor(1.3070)}" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "noise_vars = ['N_3k', 'N_2k', 'N_k']\n", + "noise_dists = {N: Normal(10., 1.) for N in noise_vars}\n", + "\n", + "mapk_receiver(noise_dists)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Condition the model on observed map3k, map2k, mapk" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "conditioned_mapk = pyro.condition(mapk_receiver, data={'map2k': 13.8113, 'map3k': 11.0504, 'mapk': 0.2803})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Infer the noise distributions N_3k, N_2k, N_k" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "noise_priors = {N: Normal(0., 1.) for N in noise_vars}\n", + "\n", + "posterior = infer_dist(conditioned_mapk, noise_priors)\n", + "\n", + "noise_marginals = {\n", + " n: EmpiricalMarginal(posterior, sites=n)\n", + " for n in noise_vars\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGaZJREFUeJzt3Xu4HXV97/H3RyIoKIRL9EhCDWK8\ntfVCc8D78Yg3RIXnHLFYWgKiWKVqtRaDPRW1j4rVR3uw1XOooPGGF4oaxRsFrVoLJYggiBwicolB\nCHIHFZHv+WN+kcVmXzJh7z07yfv1POtZM7/5zcx3LcL67PnNWjOpKiRJ2lD3GboASdKmxeCQJPVi\ncEiSejE4JEm9GBySpF4MDklSLwaH5rQkb07y4WnaViV5+L1Y/6tJlk1HLZuqJE9JckmSW5IcMIv7\nPTTJd2drf5qcwaENkuSyJLcn2WVM+w/aB/LimdhvVb2zql4+E9seleRbSV4+pu0ZSdaM1LJvVa3Y\ngG3dq4Ca494O/GNVPaCqvjB2Yft3cnWS7UbaXp7kW5NtNMmDkpyUZG2SG5P8e5K9p798TQeDQ338\nFHjp+pkkfwjcf2M3lmTedBS1JZkD79lDgQun6DMPeF3P7T4AOBv4I2AnYAVwapIH9K5QM87gUB8f\nBw4ZmV8GfGy0Q5L9kpyb5KYkVyZ568iyxe2v8cOTXAGc0doPSXJ5kl8k+dv2V+uz2rK3JvnEmPWX\nJbkiybVJ/mZk+3sl+Y8kNyS5Ksk/Jtl6ul786FFJkocn+bf21/G1ST7T2r/dup/XhnP+uLW/Isnq\nJNclWZlk15HtPifJxW1bH2zbXb+fQ9tf3+9Pch3w1iR7JDmjvV/XJvlkkvkj27ssyV8nOT/JrUlO\nSPLgNtR2c5J/TbLjJK9z3FqT/AR4GPCl9tq2mWAT7wHeOFrTVKrq0qp6X1VdVVW/rarjga2BR05Q\n43uSfDfJDhu6D00fg0N9nAlsn+TRSbYC/hj4xJg+t9KFy3xgP+BV44yF/zfg0cBzkzwG+CBwMPAQ\nYAdg4RR1PJXuA2Uf4C1JHt3afwu8HtgFeFJb/uq+L3ID/R3wDWBHYBHwAYCqenpb/rg2nPOZJM8E\n3gW8hO41Xg58GqAN/Z0MHA3sDFwMPHnMvvYGLgUeBLwDSNvernTv427AW8es8z+BZwOPAF4IfBV4\nM917cx/gteO9qMlqrao9gCuAF7bX9usJ3ptVwLeAN06wfEpJHk8XHKvHtN8nyT8DjwWeU1U3buw+\ntPEMDvW1/qjj2cCPgZ+NLqyqb1XVD6vqzqo6HziJLihGvbWqbq2qXwIvBr5UVd+tqtuBtwBTXUDt\nbVX1y6o6DzgPeFzb9zlVdWZV3VFVlwH/d5x9T+a4drRyQ5IbgC9P0vc3dMM2u1bVr6pqshO3BwMn\nVtX324ft0cCT2nmh5wMXVtUpVXUHcBzw8zHrr62qD7TX9cuqWl1Vp1XVr6tqHfC+cV7nB6rq6qr6\nGfAd4KyqOrft//PAEzai1j7eArwmyYKe65Fke7p/Z28bEwz3pfv3tBNdeN3Wd9uaHgaH+vo48CfA\noYwZpgJIsneSbyZZl+RG4M/p/soddeXI9K6j8+3D4BdT1DD6wXob3fg4SR6R5MtJfp7kJuCd4+x7\nMq+tqvnrH8ALJul7FN1f/v+Z5MIkL5uk7650f7kDUFW30L3Ghdzz9RewZsz6o+/X+hPJn07ys/Y6\nP8E9X+fVI9O/HGd+onMHk9W6warqArrgXd5nvST3B74EnFlV7xqz+OHA/nSBcnuf7Wp6GRzqpaou\npztJ/nzglHG6fApYCexWVTsA/4fuA/ZumxmZvopuqAf43QfHzhtZ3ofojoKWVNX2dEMzY/c9Larq\n51X1iqraFXgl8MFM/E2qtXRHJwCk+8bRznRHa2Nff0bn1+9uzPy7Wttj2+v8U6bvdU5Wa1/HAK9g\nA0OnnTP5QtvXK8fpchFwGPDVJOOe+9DsMDi0MQ4HnllVt46z7IHAdVX1qyR70R2dTOZk4IVJntxO\nZL+Njf8QfCBwE3BLkkcBr9rI7UwpyYFJ1n/AX0/3Qf7bNn813Unk9T4FHJbk8e3D8Z10Q0eXAacC\nf5jkgHTfmDoS+C9T7P6BwC3ADUkWAn89Ha9pA2rtpapWA59hgvMpo5Lcl+7fwi+BQ6rqzgm2eRLd\nHwT/mmSPvjVpehgc6q2qflJVqyZY/Grg7Uluphvn/uwU27oQeA3dCdirgJuBa4CJTrxO5o10QXUz\n8M90H1oz5b8CZyW5he4I63VV9dO27K3Ainau5CVVdTrwt8C/0L3GPYCDAKrqWuBA4O/phoQeQ3dy\nebLX/zZgT+BGuuAZ78hvo0xW60Z6O7DdlL26LwS8AHgOXSDe0h5PG6fGFW27Z2zEuRdNg3gjJ80l\n6b63fwPdcNNPp+q/uUlyH7pzHAdX1TeHrkcaj0ccGlySFybZto2nvxf4IXDZsFXNniTPTTK/DQ2t\nPy9z5sBlSRMyODQX7E93UnYtsAQ4qLasQ+EnAT8BrqX7zcUB7avKm5UkTxsZgrrbY+ja1I9DVZKk\nXjzikCT1MvQF02bELrvsUosXLx66DEnapJxzzjnXVtWUv/bfLINj8eLFrFo10bdFJUnjSXL51L0c\nqpIk9WRwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9bJZ/nJcmsri5acOtu/L\njt1vsH1L08EjDklSLzMWHElOTHJNkgtG2t6T5MdJzk/y+STzR5YdnWR1kouTPHek/XmtbXWS5TNV\nryRpw8zkEcdHgeeNaTsN+IOqeizw/4CjAZI8hu6+xr/f1vlgkq2SbAX8E7Av3b2YX9r6SpIGMmPB\nUVXfBq4b0/aNqrqjzZ4JLGrT+wOfrqpft/tMrwb2ao/VVXVpVd0OfLr1lSQNZMhzHC8DvtqmFwJX\njixb09omar+HJEckWZVk1bp162agXEkSDBQcSf4GuAP45PqmcbrVJO33bKw6vqqWVtXSBQumvA+J\nJGkjzfrXcZMsA14A7FN33fB8DbDbSLdFwNo2PVG7JGkAs3rEkeR5wJuAF1XVbSOLVgIHJdkmye7A\nEuA/gbOBJUl2T7I13Qn0lbNZsyTp7mbsiCPJScAzgF2SrAGOofsW1TbAaUkAzqyqP6+qC5N8FvgR\n3RDWkVX127advwC+DmwFnFhVF85UzZKkqc1YcFTVS8dpPmGS/u8A3jFO+1eAr0xjaZKke8FfjkuS\nejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaH\nJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSL/OGLkDa0ixefuog+73s2P0G2a82\nPx5xSJJ6MTgkSb0YHJKkXgwOSVIvMxYcSU5Mck2SC0badkpyWpJL2vOOrT1JjkuyOsn5SfYcWWdZ\n639JkmUzVa8kacPM5BHHR4HnjWlbDpxeVUuA09s8wL7AkvY4AvgQdEEDHAPsDewFHLM+bCRJw5ix\n4KiqbwPXjWneH1jRplcAB4y0f6w6ZwLzkzwEeC5wWlVdV1XXA6dxzzCSJM2i2T7H8eCqugqgPT+o\ntS8Erhzpt6a1TdR+D0mOSLIqyap169ZNe+GSpM5cOTmecdpqkvZ7NlYdX1VLq2rpggULprU4SdJd\nZjs4rm5DULTna1r7GmC3kX6LgLWTtEuSBjLbwbESWP/NqGXAF0faD2nfrnoicGMbyvo68JwkO7aT\n4s9pbZKkgczYtaqSnAQ8A9glyRq6b0cdC3w2yeHAFcCBrftXgOcDq4HbgMMAquq6JH8HnN36vb2q\nxp5wlyTNohkLjqp66QSL9hmnbwFHTrCdE4ETp7E0SdK9MFdOjkuSNhEGhySpF4NDktSLwSFJ6sXg\nkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknqZscuqSxti8fJT\nhy5BUk8ecUiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSehkkOJK8PsmF\nSS5IclKS+yXZPclZSS5J8pkkW7e+27T51W354iFqliR1Zj04kiwEXgssrao/ALYCDgLeDby/qpYA\n1wOHt1UOB66vqocD72/9JEkDGWqoah5w/yTzgG2Bq4BnAie35SuAA9r0/m2etnyfJJnFWiVJI2Y9\nOKrqZ8B7gSvoAuNG4Bzghqq6o3VbAyxs0wuBK9u6d7T+O89mzZKkuwwxVLUj3VHE7sCuwHbAvuN0\nrfWrTLJsdLtHJFmVZNW6deumq1xJ0hhDDFU9C/hpVa2rqt8ApwBPBua3oSuARcDaNr0G2A2gLd8B\nuG7sRqvq+KpaWlVLFyxYMNOvQZK2WEMExxXAE5Ns285V7AP8CPgm8OLWZxnwxTa9ss3Tlp9RVfc4\n4pAkzY4hznGcRXeS+/vAD1sNxwNvAt6QZDXdOYwT2ionADu39jcAy2e7ZknSXQa5A2BVHQMcM6b5\nUmCvcfr+CjhwNuqSJE3NX45LknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBI\nknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoZ5A6Akmbf4uWnDrbvy47d\nb7B9a/pt0BFHkv81Mr3NzJUjSZrrJg2OJEcleRLw4pHm/5jZkiRJc9lUQ1UXAwcCD0vyHeAiYOck\nj6yqi2e8OknSnDPVUNX1wJuB1cAzgONa+/Ik35vBuiRJc9RURxzPA44B9gDeB5wH3FpVh810YZKk\nuWnSI46qenNV7QNcBnyCLmgWJPluki/NQn2SpDlmQ7+O+/WqOhs4O8mrquqpSXaZycIkSXPTBn0d\nt6qOGpk9tLVdu7E7TTI/yclJfpzkoiRPSrJTktOSXNKed2x9k+S4JKuTnJ9kz43dryTp3uv9y/Gq\nOm8a9vu/ga9V1aOAx9F9W2s5cHpVLQFOb/MA+wJL2uMI4EPTsH9J0kaa9UuOJNkeeDpwAkBV3V5V\nNwD7AytatxXAAW16f+Bj1TkTmJ/kIbNctiSpGeJaVQ8D1gEfSXJukg8n2Q54cFVdBdCeH9T6LwSu\nHFl/TWu7myRHJFmVZNW6detm9hVI0hZsiOCYB+wJfKiqngDcyl3DUuPJOG11j4aq46tqaVUtXbBg\nwfRUKkm6hyGCYw2wpqrOavMn0wXJ1euHoNrzNSP9dxtZfxGwdpZqlSSNMevBUVU/B65M8sjWtA/w\nI2AlsKy1LQO+2KZXAoe0b1c9Ebhx/ZCWJGn2DXVZ9dcAn0yyNXApcBhdiH02yeHAFXTXyAL4CvB8\nusue3Nb6SpIGMkhwVNUPgKXjLNpnnL4FHDnjRUmaMUPdC8T7gMwM7wAoSerF4JAk9WJwSJJ6MTgk\nSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReD\nQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSeplsOBIslWSc5N8\nuc3vnuSsJJck+UySrVv7Nm1+dVu+eKiaJUnDHnG8DrhoZP7dwPuraglwPXB4az8cuL6qHg68v/WT\nJA1kkOBIsgjYD/hwmw/wTODk1mUFcECb3r/N05bv0/pLkgYw1BHHPwBHAXe2+Z2BG6rqjja/BljY\nphcCVwK05Te2/neT5Igkq5KsWrdu3UzWLklbtFkPjiQvAK6pqnNGm8fpWhuw7K6GquOramlVLV2w\nYME0VCpJGs+8Afb5FOBFSZ4P3A/Ynu4IZH6See2oYhGwtvVfA+wGrEkyD9gBuG72y5YkwQBHHFV1\ndFUtqqrFwEHAGVV1MPBN4MWt2zLgi216ZZunLT+jqu5xxCFJmh1z6XccbwLekGQ13TmME1r7CcDO\nrf0NwPKB6pMkMcxQ1e9U1beAb7XpS4G9xunzK+DAWS1MkjShuXTEIUnaBBgckqReBh2q0tyxePmp\nQ5cgaRPhEYckqReDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8Eh\nSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJwSJJ6mfXgSLJb\nkm8muSjJhUle19p3SnJakkva846tPUmOS7I6yflJ9pztmiVJdxniiOMO4K+q6tHAE4EjkzwGWA6c\nXlVLgNPbPMC+wJL2OAL40OyXLElab9aDo6quqqrvt+mbgYuAhcD+wIrWbQVwQJveH/hYdc4E5id5\nyCyXLUlqBj3HkWQx8ATgLODBVXUVdOECPKh1WwhcObLamtY2dltHJFmVZNW6detmsmxJ2qINFhxJ\nHgD8C/CXVXXTZF3Haat7NFQdX1VLq2rpggULpqtMSdIYgwRHkvvShcYnq+qU1nz1+iGo9nxNa18D\n7Day+iJg7WzVKkm6uyG+VRXgBOCiqnrfyKKVwLI2vQz44kj7Ie3bVU8Eblw/pCVJmn3zBtjnU4A/\nA36Y5Aet7c3AscBnkxwOXAEc2JZ9BXg+sBq4DThsdsuVJI2a9eCoqu8y/nkLgH3G6V/AkTNalCRp\ng/nLcUlSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSL0NcckSSZsXi5acO\ntu/Ljt1vsH3PNI84JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxa/jziFDfnVQkjaURxySpF48\n4pCkGTDUCMJs/PDQIw5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknrZZIIjyfOSXJxkdZLlQ9cj\nSVuqTeJ3HEm2Av4JeDawBjg7ycqq+tFM7M9fcEvSxDaVI469gNVVdWlV3Q58Gth/4JokaYu0SRxx\nAAuBK0fm1wB7j3ZIcgRwRJu9JcnFwC7AtbNS4abJ92divjcT872Z2ODvTd59r1Z/6IZ02lSCI+O0\n1d1mqo4Hjr/bSsmqqlo6k4Vtynx/JuZ7MzHfm4ltKe/NpjJUtQbYbWR+EbB2oFokaYu2qQTH2cCS\nJLsn2Ro4CFg5cE2StEXaJIaqquqOJH8BfB3YCjixqi7cgFWPn7rLFs33Z2K+NxPzvZnYFvHepKqm\n7iVJUrOpDFVJkuYIg0OS1MtmHRxJtkpybpIvD13LXJLksiQ/TPKDJKuGrmcuSTI/yclJfpzkoiRP\nGrqmuSLJI9u/mfWPm5L85dB1zRVJXp/kwiQXJDkpyf2GrmmmbNbnOJK8AVgKbF9VLxi6nrkiyWXA\n0qryR1xjJFkBfKeqPty+wbdtVd0wdF1zTbsM0M+Avavq8qHrGVqShcB3gcdU1S+TfBb4SlV9dNjK\nZsZme8SRZBGwH/DhoWvRpiHJ9sDTgRMAqup2Q2NC+wA/MTTuZh5w/yTzgG3ZjH9rttkGB/APwFHA\nnUMXMgcV8I0k57RLtajzMGAd8JE2xPnhJNsNXdQcdRBw0tBFzBVV9TPgvcAVwFXAjVX1jWGrmjmb\nZXAkeQFwTVWdM3Qtc9RTqmpPYF/gyCRPH7qgOWIesCfwoap6AnAr4CX8x2hDeC8CPjd0LXNFkh3p\nLry6O7ArsF2SPx22qpmzWQYH8BTgRW0s/9PAM5N8YtiS5o6qWtuerwE+T3f1YXWXtllTVWe1+ZPp\ngkR3ty/w/aq6euhC5pBnAT+tqnVV9RvgFODJA9c0YzbL4Kiqo6tqUVUtpjukPqOqNtv07yPJdkke\nuH4aeA5wwbBVzQ1V9XPgyiSPbE37ADNyz5dN3EtxmGqsK4AnJtk2Sej+7Vw0cE0zZpO45Iim1YOB\nz3f/tpkHfKqqvjZsSXPKa4BPtuGYS4HDBq5nTkmyLd0N1V45dC1zSVWdleRk4PvAHcC5bMaXH9ms\nv44rSZp+m+VQlSRp5hgckqReDA5JUi8GhySpF4NDktSLwSENIMniJH+yEevNT/Lqkfld29dApVlj\ncEjDWAyMGxztInkTmQ/8Ljiqam1VvXh6S5MmZ3Boi5LkkCTnJzkvyceTPDTJ6a3t9CS/1/p9NMlx\nSb6X5NIkLx7ZxlHtfibnJTm2te2R5GvtwpHfSfKoKbZzLPC0dl+L1yc5NMnnknyJ7gKUD2j1fL/t\na/+R9fZo672nHblc0PZ1vyQfaf3PTfLfW/uhSU5p9V2S5O9n5c3W5quqfPjYIh7A7wMXA7u0+Z2A\nLwHL2vzLgC+06Y/SXcTvPsBjgNWtfV/ge3T36QDYqT2fDixp03vTXeZmsu08A/jySG2H0l0ra/32\n5tHdRwZgF2A1ELojlQtG1vvdPPBXwEfa9KPoLoNxv7btS4Ed2vzlwG5D//fwsek+vOSItiTPBE6u\ndgOrqrqu3eHvf7TlHwdG/xr/QlXdCfwoyYNb27PoPpxvG9nGA+guaPe5dikXgG2m2M54Tquq69p0\ngHe2KxffCSyku1zMZJ4KfKDV9eMklwOPaMtOr6obAZL8CHgocOUU25PGZXBoSxK6e5FMZnT5r8es\nO9E27gPcUFWPn2Cb421nPLeOTB8MLAD+qKp+0670PNWtSCfb9mgNv8X/93UveI5DW5LTgZck2Rkg\nyU50w04HteUH093+czLfAF7WLvZHkp2q6ibgp0kObG1J8rgptnMz8MBJlu9Ad0+Z37RzFQ/dgPW+\n3V4DSR4B/B7d0Jw0rQwObTGq6kLgHcC/JTkPeB/wWuCwJOcDfwa8boptfA1YCaxK8gPgjW3RwcDh\nbbsX0t3UZzLnA3e0E+yvH2f5J4GlSVa1bf+47f8XwL8nuSDJe8as80FgqyQ/BD4DHFpVv0aaZl4d\nV5LUi0cckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknr5/1TqS69Z5sPcAAAAAElFTkSu\nQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "hist(noise_marginals['N_2k'], 'N_2k')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. do Map3k = 5.678" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "do_map3k = pyro.do(mapk_receiver, data={'map3k': 5.678})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Infer Mapk" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "ename": "RuntimeError", + "evalue": "One of the differentiated Tensors appears to not have been used in the graph. Set allow_unused=True if this is the desired behavior.", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmapk_do_dist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minfer_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdo_map3k\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnoise_marginals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mmapk_marginal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mEmpiricalMarginal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmapk_do_dist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msites\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'mapk'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Desktop/Git/causal_demon/mapk_demon/inference.py\u001b[0m in \u001b[0;36minfer_dist\u001b[0;34m(prog, n_dist, type)\u001b[0m\n\u001b[1;32m 19\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mtype\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'mcmc'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 20\u001b[0m \u001b[0mhmc_kernel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mHMC\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprog\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstep_size\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.9\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m \u001b[0mposterior\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMCMC\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhmc_kernel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_samples\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwarmup_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m50\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_dist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 22\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mposterior\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/infer/abstract_infer.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 81\u001b[0m \"\"\"\n\u001b[1;32m 82\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_init\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 83\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mtr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlogit\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpoutine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mblock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_traces\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 84\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexec_traces\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 85\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_weights\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlogit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/infer/mcmc/mcmc.py\u001b[0m in \u001b[0;36m_traces\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 35\u001b[0m \u001b[0mlogging_interval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mceil\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_samples\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_samples\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 37\u001b[0;31m \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkernel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 38\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mlogging_interval\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 39\u001b[0m \u001b[0mstage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"WARMUP\"\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m\"SAMPLE\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/infer/mcmc/hmc.py\u001b[0m in \u001b[0;36msample\u001b[0;34m(self, trace)\u001b[0m\n\u001b[1;32m 232\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_potential_energy\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 233\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep_size\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 234\u001b[0;31m self.num_steps)\n\u001b[0m\u001b[1;32m 235\u001b[0m \u001b[0;31m# apply Metropolis correction.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 236\u001b[0m \u001b[0menergy_proposal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_energy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz_new\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr_new\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/ops/integrator.py\u001b[0m in \u001b[0;36mvelocity_verlet\u001b[0;34m(z, r, potential_fn, step_size, num_steps)\u001b[0m\n\u001b[1;32m 22\u001b[0m \u001b[0mz_next\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0mr_next\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 24\u001b[0;31m \u001b[0mgrads\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpotential_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz_next\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 25\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 26\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum_steps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/ops/integrator.py\u001b[0m in \u001b[0;36m_grad\u001b[0;34m(potential_fn, z)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[0mnode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrequires_grad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 65\u001b[0m \u001b[0mpotential_energy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpotential_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 66\u001b[0;31m \u001b[0mgrads\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgrad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpotential_energy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz_nodes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 67\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mnode\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mz_nodes\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0mnode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrequires_grad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/torch/autograd/__init__.py\u001b[0m in \u001b[0;36mgrad\u001b[0;34m(outputs, inputs, grad_outputs, retain_graph, create_graph, only_inputs, allow_unused)\u001b[0m\n\u001b[1;32m 143\u001b[0m return Variable._execution_engine.run_backward(\n\u001b[1;32m 144\u001b[0m \u001b[0moutputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad_outputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mretain_graph\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcreate_graph\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 145\u001b[0;31m inputs, allow_unused)\n\u001b[0m\u001b[1;32m 146\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 147\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mRuntimeError\u001b[0m: One of the differentiated Tensors appears to not have been used in the graph. Set allow_unused=True if this is the desired behavior." + ] + } + ], + "source": [ + "mapk_do_dist = infer_dist(do_map3k, noise_marginals)\n", + "mapk_marginal = EmpiricalMarginal(mapk_do_dist, sites = 'mapk')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/run_mapk.py b/run_mapk.py new file mode 100644 index 0000000..ef03c0e --- /dev/null +++ b/run_mapk.py @@ -0,0 +1,11 @@ +''' +run_mapk +Author: Kaushal Paneri +Project: causal_demon +Date of Creation: 12/29/18 +''' + +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function +from __future__ import unicode_literals \ No newline at end of file From 37562779273490b94c203139f66ee8aa9e912b3d Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Fri, 15 Feb 2019 13:48:54 -0500 Subject: [PATCH 08/11] mapk inference working --- mapk_demon/inference.py | 25 ++++++++ mapk_demon/receivers.py | 121 +++++++++--------------------------- mapk_demon/svi_test.py | 133 +++++++++++++++++++++++++++++++++++++++- 3 files changed, 185 insertions(+), 94 deletions(-) diff --git a/mapk_demon/inference.py b/mapk_demon/inference.py index e69de29..aa3575a 100644 --- a/mapk_demon/inference.py +++ b/mapk_demon/inference.py @@ -0,0 +1,25 @@ +import pyro +from pyro.infer.mcmc import MCMC +from pyro.infer.mcmc.nuts import HMC + + +def infer_dist(prog, n_dist, type='mcmc'): + """Obtain the unique distribution entailed by a SCM program. + + Do inference on a SCM program and obtain a object representing the + probability distribution entailed by the SCM. + + This implementation depends on simple importance sampling with 5000 + samples. + + `prog`: the subroutine encoding the SCM. + `n_dist`: a dictionary containing distributions for each + noise object. + """ + if type == 'mcmc': + hmc_kernel = HMC(prog, step_size=0.9, num_steps=4) + posterior = MCMC(hmc_kernel, num_samples=1000, warmup_steps=50).run(n_dist) + return posterior + + else: + return 0 diff --git a/mapk_demon/receivers.py b/mapk_demon/receivers.py index 86cf0b6..05eabfe 100644 --- a/mapk_demon/receivers.py +++ b/mapk_demon/receivers.py @@ -2,22 +2,12 @@ from functools import partial -import numpy as np import pyro from pyro import sample from pyro.distributions import Normal, Uniform -from torch import tensor -from mapk_demon.transmitters import mapk_signaling -from pyro.infer.mcmc import MCMC -from pyro.infer.mcmc.nuts import NUTS -from pyro.infer import EmpiricalMarginal -from causal_demon.inference import infer_dist -import torch -import matplotlib.pyplot as plt -def infer_mcmc(model): - nuts_kernel = NUTS(model, adapt_step_size=True) - return MCMC(nuts_kernel, num_samples=500, warmup_steps=300) +import pyro.optim +import torch def g1(a): @@ -32,27 +22,6 @@ def f(mu, N): return N * mu ** (0.5) + mu -# Parameters - -params = { - "alpha_3k": Normal(500, 1), - "alpha_2k": Normal(500, 1), - "alpha_k": Normal(500, 1), - "nu_3k": Normal(.15, .05), - "nu_2k": Normal(.15, .05), - "nu_k": Normal(.15, .05) -} -''' -params = { - "alpha_3k": 500, - "alpha_2k": 500, - "alpha_k": 500, - "nu_3k": .15, - "nu_2k": .15, - "nu_k": .15 -} -''' - hyperparameters = { "t_3k": 1.2, "t_2k": 1.2, @@ -62,40 +31,55 @@ def f(mu, N): def f_map3k(E1, N_3k, params): total_3k = hyperparameters['t_3k'] - alpha_3k = sample("alpha_3k",params['alpha_3k']) - nu_3k = sample("nu_3k", params['nu_3k']) + alpha_3k = params['alpha_3k'].rsample() + nu_3k = params['nu_3k'].rsample() map3k_mu = total_3k * g1(E1 * (alpha_3k/nu_3k)) - map3k = sample("map3k", Normal(f(map3k_mu, N_3k), 1e-6)) + map3k = sample("map3k", Normal(f(map3k_mu, N_3k), 1.)) return map3k def f_map2k(map3k, N_2k, params): total_2k = hyperparameters['t_2k'] - alpha_2k = sample("alpha_2k", params['alpha_2k']) - nu_2k = sample("nu_2k", params['nu_2k']) + alpha_2k = params['alpha_2k'].rsample() + nu_2k = params['nu_2k'].rsample() map2k_mu = total_2k * g2(map3k * (alpha_2k / nu_2k)) - map2k = sample("map2k", Normal(f(map2k_mu, N_2k), 1e-6)) + map2k = sample("map2k", Normal(f(map2k_mu, N_2k), 1.)) return map2k def f_mapk(map2k, N_k, params): total_k = hyperparameters['t_k'] - alpha_k = sample("alpha_k",params['alpha_k']) - nu_k = sample("nu_k", params['nu_k']) + alpha_k = params['alpha_k'].rsample() + nu_k = params['nu_k'].rsample() mapk_mu = total_k * g2(map2k * (alpha_k / nu_k)) - mapk = sample("mapk", Normal(f(mapk_mu, N_k), 1e-6)) + mapk = sample("mapk", Normal(f(mapk_mu, N_k), 1.)) return mapk -def mapk_signaling_receiver(noise_dists): +def mapk_receiver(noise_dists): + # Parameters + al_m = torch.tensor(700.) + al_s = torch.tensor(1.) + + nu_m = torch.tensor(.15) + nu_s = torch.tensor(.05) + + params = { + "alpha_3k": Normal(al_m, al_s), + "alpha_2k": Normal(al_m, al_s), + "alpha_k": Normal(al_m, al_s), + "nu_3k": Normal(nu_m, nu_s), + "nu_2k": Normal(nu_m, nu_s), + "nu_k": Normal(nu_m, nu_s) + } with pyro.iarange("model"): - E1 = sample('E1', Uniform(1.5e-5, 6e-5)) + E1 = Uniform(1.5e-5, 10.).rsample() N_3k = sample('N_3k', noise_dists['N_3k']) N_2k = sample('N_2k', noise_dists['N_2k']) N_k = sample('N_k', noise_dists['N_k']) @@ -109,52 +93,3 @@ def mapk_signaling_receiver(noise_dists): 'map2k': map2k, 'mapk': mapk } - -def plot(data): - plt.boxplot(data['erk']) - plt.title("Empirical Distribution of Erk") - plt.xlabel("Condition") - plt.ylabel("Concentration") - plt.show() - - -if __name__ == "__main__": - - # Noise Distributions - noise_vars = ['N_3k', 'N_2k', 'N_k'] - noise_priors = {N: Normal(0, 1) for N in noise_vars} - - parameters = ["alpha_3k","alpha_2k","alpha_k","nu_3k","nu_2k","nu_k"] - - receiver_model = mapk_signaling_receiver(noise_dists=noise_priors) - - # Simulate data from transmitter model - evidence = [] - - for i in range(10): - evidence.append(mapk_signaling(noise_dists=noise_priors)) - - # condition the receiver model with data - conditioned = pyro.condition(mapk_signaling_receiver, data = evidence) - - print(conditioned(noise_priors)) - - # mcmc inference - #nuts_kernel = NUTS(conditioned, adapt_step_size=True) - #mcmc_run = MCMC(nuts_kernel, num_samples=50, warmup_steps=10).run(noise_priors) - - # Importance Sampling - posterior = pyro.infer.Importance(conditioned, num_samples=100).run(noise_priors) - - #print(posterior) - posteriors = { - n: EmpiricalMarginal(posterior, sites=n) - for n in parameters - } - - print("Alpha_3k inference") - alpha_3k_marginal = posteriors['alpha_3k'] - alpha_3k_samples = [alpha_3k_marginal().item() for _ in range(500)] - - print("mean:",np.mean(alpha_3k_samples)) - print("std:",np.std(alpha_3k_samples)) diff --git a/mapk_demon/svi_test.py b/mapk_demon/svi_test.py index edc1699..1fe3d84 100644 --- a/mapk_demon/svi_test.py +++ b/mapk_demon/svi_test.py @@ -8,4 +8,135 @@ from __future__ import absolute_import from __future__ import division from __future__ import print_function -from __future__ import unicode_literals \ No newline at end of file +from __future__ import unicode_literals + +import numpy as np +import pandas as pd +import pyro +from pyro import sample +from pyro.distributions import Normal, Uniform +from torch import tensor +from mapk_demon.transmitters import mapk_signaling +from pyro.infer.mcmc import MCMC +from pyro.infer.mcmc.nuts import NUTS, HMC +from pyro.infer import EmpiricalMarginal +#from causal_demon.inference import infer_dist +import torch +import matplotlib.pyplot as plt +import json +from torch.distributions import constraints +from tqdm import tqdm + +import pyro.optim + + +def test_transmitter(inp): + al_m = torch.tensor(inp) + al_s = torch.tensor(1.) + + nu_m = torch.tensor(.15) + nu_s = torch.tensor(.05) + + params = { + "alpha_3k": Normal(al_m, al_s), + "alpha_2k": Normal(al_m, al_s), + "alpha_k": Normal(al_m, al_s), + "nu_3k": Normal(nu_m, nu_s), + "nu_2k": Normal(nu_m, nu_s), + "nu_k": Normal(nu_m, nu_s) + } + return { + 'alpha_3k': sample("alpha_3k", params['alpha_3k']), + 'alpha_2k': sample("alpha_2k", params['alpha_2k']), + 'alpha_k': sample("alpha_k", params['alpha_k']) + } + + +def test_receiver(inp): + al_m = sample("al_m", Normal(inp, 1)) + al_s = torch.tensor(1.) + + nu_m = torch.tensor(.15) + nu_s = torch.tensor(.05) + + params = { + "alpha_3k": Normal(al_m, al_s), + "alpha_2k": Normal(al_m, al_s), + "alpha_k": Normal(al_m, al_s), + "nu_3k": Normal(nu_m, nu_s), + "nu_2k": Normal(nu_m, nu_s), + "nu_k": Normal(nu_m, nu_s) + } + return { + 'alpha_3k': sample("alpha_3k", params['alpha_3k']), + 'alpha_2k': sample("alpha_2k", params['alpha_2k']), + 'alpha_k': sample("alpha_k", params['alpha_k']) + } + + +def receiver_guide(inp): + al_m = pyro.param("al_m", Normal(inp, 1)) + al_s = pyro.param("al_s", torch.tensor(1.), constraint=constraints.positive) + + nu_m = pyro.param("nu_m", torch.tensor(.15)) + nu_s = pyro.param('nu_m', torch.tensor(.05), constraint=constraints.positive) + + params = { + "alpha_3k": Normal(al_m, al_s), + "alpha_2k": Normal(al_m, al_s), + "alpha_k": Normal(al_m, al_s), + "nu_3k": Normal(nu_m, nu_s), + "nu_2k": Normal(nu_m, nu_s), + "nu_k": Normal(nu_m, nu_s) + } + + return { + 'alpha_3k': sample("alpha_3k", params['alpha_3k']), + 'alpha_2k': sample("alpha_2k", params['alpha_2k']), + 'alpha_k': sample("alpha_k", params['alpha_k']) + } + + +if __name__ == "__main__": + sample_size = 10000 + + ground = 500. + inp = 700. + output = test_transmitter(ground) + + # Noise Distributions + noise_vars = ['N_3k', 'N_2k', 'N_k'] + noise_priors = {N: Normal(0, 1) for N in noise_vars} + parameters = ["alpha_3k", "alpha_2k", "alpha_k", "nu_3k", "nu_2k", "nu_k"] + + # Simulate data from transmitter model + evidence = [] + + for i in range(sample_size): + evidence.append(test_transmitter(ground)) + + + # parameter settings + print(evidence) + + #receiver_model = mapk_signaling_receiver(noise_dists=noise_priors) + + # condition the receiver model with data + conditioned = pyro.condition(test_receiver, data = output) + + print("conditioned:: ", conditioned(inp)) + + svi = pyro.infer.SVI(model=conditioned, + guide=receiver_guide, + optim=pyro.optim.SGD({"lr": 0.0001, "momentum": 0.1}), + loss=pyro.infer.Trace_ELBO()) + + losses, a = [], [] + num_steps = 2500 + for t in tqdm(range(num_steps)): + losses.append(svi.step(inp)) + a.append(pyro.param("al_m").item()) + print(pyro.get_param_store().get_param('al_m')) + + print(a) + print(losses) From bfb04155414526d273431f2cb306eed314c8495d Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Fri, 15 Feb 2019 14:08:56 -0500 Subject: [PATCH 09/11] mapk inference working --- notebooks/Mapk Counterfactual Inference.ipynb | 49 +++++++++++++------ 1 file changed, 33 insertions(+), 16 deletions(-) diff --git a/notebooks/Mapk Counterfactual Inference.ipynb b/notebooks/Mapk Counterfactual Inference.ipynb index d9c7ea1..de5395b 100644 --- a/notebooks/Mapk Counterfactual Inference.ipynb +++ b/notebooks/Mapk Counterfactual Inference.ipynb @@ -9,8 +9,10 @@ }, { "cell_type": "code", - "execution_count": 1, - "metadata": {}, + "execution_count": 2, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "import sys\n", @@ -44,18 +46,25 @@ " plt.ylabel(\"#\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let $$N_{3k}, N_{2k}, N_k \\sim Normal\\ (10, 1)$$." + ] + }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "{'map2k': tensor(12.1773), 'map3k': tensor(11.9624), 'mapk': tensor(1.3070)}" + "{'map2k': tensor(12.9241), 'map3k': tensor(12.4646), 'mapk': tensor(-0.2236)}" ] }, - "execution_count": 6, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -76,7 +85,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 17, "metadata": { "collapsed": true }, @@ -89,13 +98,19 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 2. Infer the noise distributions N_3k, N_2k, N_k" + "## 2. Infer the noise distributions N_3k, N_2k, N_k\n", + "\n", + "Let there be a prior over Noise distributions. $$N_{3k}, N_{2k}, N_k \\sim Normal\\ (0, 1)$$.\n", + "\n", + "I have used HMC and MCMC for inference." ] }, { "cell_type": "code", - "execution_count": 8, - "metadata": {}, + "execution_count": 18, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "noise_priors = {N: Normal(0., 1.) for N in noise_vars}\n", @@ -110,14 +125,14 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 20, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGaZJREFUeJzt3Xu4HXV97/H3RyIoKIRL9EhCDWK8\ntfVCc8D78Yg3RIXnHLFYWgKiWKVqtRaDPRW1j4rVR3uw1XOooPGGF4oaxRsFrVoLJYggiBwicolB\nCHIHFZHv+WN+kcVmXzJh7z07yfv1POtZM7/5zcx3LcL67PnNWjOpKiRJ2lD3GboASdKmxeCQJPVi\ncEiSejE4JEm9GBySpF4MDklSLwaH5rQkb07y4WnaViV5+L1Y/6tJlk1HLZuqJE9JckmSW5IcMIv7\nPTTJd2drf5qcwaENkuSyJLcn2WVM+w/aB/LimdhvVb2zql4+E9seleRbSV4+pu0ZSdaM1LJvVa3Y\ngG3dq4Ca494O/GNVPaCqvjB2Yft3cnWS7UbaXp7kW5NtNMmDkpyUZG2SG5P8e5K9p798TQeDQ338\nFHjp+pkkfwjcf2M3lmTedBS1JZkD79lDgQun6DMPeF3P7T4AOBv4I2AnYAVwapIH9K5QM87gUB8f\nBw4ZmV8GfGy0Q5L9kpyb5KYkVyZ568iyxe2v8cOTXAGc0doPSXJ5kl8k+dv2V+uz2rK3JvnEmPWX\nJbkiybVJ/mZk+3sl+Y8kNyS5Ksk/Jtl6ul786FFJkocn+bf21/G1ST7T2r/dup/XhnP+uLW/Isnq\nJNclWZlk15HtPifJxW1bH2zbXb+fQ9tf3+9Pch3w1iR7JDmjvV/XJvlkkvkj27ssyV8nOT/JrUlO\nSPLgNtR2c5J/TbLjJK9z3FqT/AR4GPCl9tq2mWAT7wHeOFrTVKrq0qp6X1VdVVW/rarjga2BR05Q\n43uSfDfJDhu6D00fg0N9nAlsn+TRSbYC/hj4xJg+t9KFy3xgP+BV44yF/zfg0cBzkzwG+CBwMPAQ\nYAdg4RR1PJXuA2Uf4C1JHt3afwu8HtgFeFJb/uq+L3ID/R3wDWBHYBHwAYCqenpb/rg2nPOZJM8E\n3gW8hO41Xg58GqAN/Z0MHA3sDFwMPHnMvvYGLgUeBLwDSNvernTv427AW8es8z+BZwOPAF4IfBV4\nM917cx/gteO9qMlqrao9gCuAF7bX9usJ3ptVwLeAN06wfEpJHk8XHKvHtN8nyT8DjwWeU1U3buw+\ntPEMDvW1/qjj2cCPgZ+NLqyqb1XVD6vqzqo6HziJLihGvbWqbq2qXwIvBr5UVd+tqtuBtwBTXUDt\nbVX1y6o6DzgPeFzb9zlVdWZV3VFVlwH/d5x9T+a4drRyQ5IbgC9P0vc3dMM2u1bVr6pqshO3BwMn\nVtX324ft0cCT2nmh5wMXVtUpVXUHcBzw8zHrr62qD7TX9cuqWl1Vp1XVr6tqHfC+cV7nB6rq6qr6\nGfAd4KyqOrft//PAEzai1j7eArwmyYKe65Fke7p/Z28bEwz3pfv3tBNdeN3Wd9uaHgaH+vo48CfA\noYwZpgJIsneSbyZZl+RG4M/p/soddeXI9K6j8+3D4BdT1DD6wXob3fg4SR6R5MtJfp7kJuCd4+x7\nMq+tqvnrH8ALJul7FN1f/v+Z5MIkL5uk7650f7kDUFW30L3Ghdzz9RewZsz6o+/X+hPJn07ys/Y6\nP8E9X+fVI9O/HGd+onMHk9W6warqArrgXd5nvST3B74EnFlV7xqz+OHA/nSBcnuf7Wp6GRzqpaou\npztJ/nzglHG6fApYCexWVTsA/4fuA/ZumxmZvopuqAf43QfHzhtZ3ofojoKWVNX2dEMzY/c9Larq\n51X1iqraFXgl8MFM/E2qtXRHJwCk+8bRznRHa2Nff0bn1+9uzPy7Wttj2+v8U6bvdU5Wa1/HAK9g\nA0OnnTP5QtvXK8fpchFwGPDVJOOe+9DsMDi0MQ4HnllVt46z7IHAdVX1qyR70R2dTOZk4IVJntxO\nZL+Njf8QfCBwE3BLkkcBr9rI7UwpyYFJ1n/AX0/3Qf7bNn813Unk9T4FHJbk8e3D8Z10Q0eXAacC\nf5jkgHTfmDoS+C9T7P6BwC3ADUkWAn89Ha9pA2rtpapWA59hgvMpo5Lcl+7fwi+BQ6rqzgm2eRLd\nHwT/mmSPvjVpehgc6q2qflJVqyZY/Grg7Uluphvn/uwU27oQeA3dCdirgJuBa4CJTrxO5o10QXUz\n8M90H1oz5b8CZyW5he4I63VV9dO27K3Ainau5CVVdTrwt8C/0L3GPYCDAKrqWuBA4O/phoQeQ3dy\nebLX/zZgT+BGuuAZ78hvo0xW60Z6O7DdlL26LwS8AHgOXSDe0h5PG6fGFW27Z2zEuRdNg3gjJ80l\n6b63fwPdcNNPp+q/uUlyH7pzHAdX1TeHrkcaj0ccGlySFybZto2nvxf4IXDZsFXNniTPTTK/DQ2t\nPy9z5sBlSRMyODQX7E93UnYtsAQ4qLasQ+EnAT8BrqX7zcUB7avKm5UkTxsZgrrbY+ja1I9DVZKk\nXjzikCT1MvQF02bELrvsUosXLx66DEnapJxzzjnXVtWUv/bfLINj8eLFrFo10bdFJUnjSXL51L0c\nqpIk9WRwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9bJZ/nJcmsri5acOtu/L\njt1vsH1L08EjDklSLzMWHElOTHJNkgtG2t6T5MdJzk/y+STzR5YdnWR1kouTPHek/XmtbXWS5TNV\nryRpw8zkEcdHgeeNaTsN+IOqeizw/4CjAZI8hu6+xr/f1vlgkq2SbAX8E7Av3b2YX9r6SpIGMmPB\nUVXfBq4b0/aNqrqjzZ4JLGrT+wOfrqpft/tMrwb2ao/VVXVpVd0OfLr1lSQNZMhzHC8DvtqmFwJX\njixb09omar+HJEckWZVk1bp162agXEkSDBQcSf4GuAP45PqmcbrVJO33bKw6vqqWVtXSBQumvA+J\nJGkjzfrXcZMsA14A7FN33fB8DbDbSLdFwNo2PVG7JGkAs3rEkeR5wJuAF1XVbSOLVgIHJdkmye7A\nEuA/gbOBJUl2T7I13Qn0lbNZsyTp7mbsiCPJScAzgF2SrAGOofsW1TbAaUkAzqyqP6+qC5N8FvgR\n3RDWkVX127advwC+DmwFnFhVF85UzZKkqc1YcFTVS8dpPmGS/u8A3jFO+1eAr0xjaZKke8FfjkuS\nejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaH\nJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSL/OGLkDa0ixefuog+73s2P0G2a82\nPx5xSJJ6MTgkSb0YHJKkXgwOSVIvMxYcSU5Mck2SC0badkpyWpJL2vOOrT1JjkuyOsn5SfYcWWdZ\n639JkmUzVa8kacPM5BHHR4HnjWlbDpxeVUuA09s8wL7AkvY4AvgQdEEDHAPsDewFHLM+bCRJw5ix\n4KiqbwPXjWneH1jRplcAB4y0f6w6ZwLzkzwEeC5wWlVdV1XXA6dxzzCSJM2i2T7H8eCqugqgPT+o\ntS8Erhzpt6a1TdR+D0mOSLIqyap169ZNe+GSpM5cOTmecdpqkvZ7NlYdX1VLq2rpggULprU4SdJd\nZjs4rm5DULTna1r7GmC3kX6LgLWTtEuSBjLbwbESWP/NqGXAF0faD2nfrnoicGMbyvo68JwkO7aT\n4s9pbZKkgczYtaqSnAQ8A9glyRq6b0cdC3w2yeHAFcCBrftXgOcDq4HbgMMAquq6JH8HnN36vb2q\nxp5wlyTNohkLjqp66QSL9hmnbwFHTrCdE4ETp7E0SdK9MFdOjkuSNhEGhySpF4NDktSLwSFJ6sXg\nkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknqZscuqSxti8fJT\nhy5BUk8ecUiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSehkkOJK8PsmF\nSS5IclKS+yXZPclZSS5J8pkkW7e+27T51W354iFqliR1Zj04kiwEXgssrao/ALYCDgLeDby/qpYA\n1wOHt1UOB66vqocD72/9JEkDGWqoah5w/yTzgG2Bq4BnAie35SuAA9r0/m2etnyfJJnFWiVJI2Y9\nOKrqZ8B7gSvoAuNG4Bzghqq6o3VbAyxs0wuBK9u6d7T+O89mzZKkuwwxVLUj3VHE7sCuwHbAvuN0\nrfWrTLJsdLtHJFmVZNW6deumq1xJ0hhDDFU9C/hpVa2rqt8ApwBPBua3oSuARcDaNr0G2A2gLd8B\nuG7sRqvq+KpaWlVLFyxYMNOvQZK2WEMExxXAE5Ns285V7AP8CPgm8OLWZxnwxTa9ss3Tlp9RVfc4\n4pAkzY4hznGcRXeS+/vAD1sNxwNvAt6QZDXdOYwT2ionADu39jcAy2e7ZknSXQa5A2BVHQMcM6b5\nUmCvcfr+CjhwNuqSJE3NX45LknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBI\nknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoZ5A6Akmbf4uWnDrbvy47d\nb7B9a/pt0BFHkv81Mr3NzJUjSZrrJg2OJEcleRLw4pHm/5jZkiRJc9lUQ1UXAwcCD0vyHeAiYOck\nj6yqi2e8OknSnDPVUNX1wJuB1cAzgONa+/Ik35vBuiRJc9RURxzPA44B9gDeB5wH3FpVh810YZKk\nuWnSI46qenNV7QNcBnyCLmgWJPluki/NQn2SpDlmQ7+O+/WqOhs4O8mrquqpSXaZycIkSXPTBn0d\nt6qOGpk9tLVdu7E7TTI/yclJfpzkoiRPSrJTktOSXNKed2x9k+S4JKuTnJ9kz43dryTp3uv9y/Gq\nOm8a9vu/ga9V1aOAx9F9W2s5cHpVLQFOb/MA+wJL2uMI4EPTsH9J0kaa9UuOJNkeeDpwAkBV3V5V\nNwD7AytatxXAAW16f+Bj1TkTmJ/kIbNctiSpGeJaVQ8D1gEfSXJukg8n2Q54cFVdBdCeH9T6LwSu\nHFl/TWu7myRHJFmVZNW6detm9hVI0hZsiOCYB+wJfKiqngDcyl3DUuPJOG11j4aq46tqaVUtXbBg\nwfRUKkm6hyGCYw2wpqrOavMn0wXJ1euHoNrzNSP9dxtZfxGwdpZqlSSNMevBUVU/B65M8sjWtA/w\nI2AlsKy1LQO+2KZXAoe0b1c9Ebhx/ZCWJGn2DXVZ9dcAn0yyNXApcBhdiH02yeHAFXTXyAL4CvB8\nusue3Nb6SpIGMkhwVNUPgKXjLNpnnL4FHDnjRUmaMUPdC8T7gMwM7wAoSerF4JAk9WJwSJJ6MTgk\nSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReD\nQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSeplsOBIslWSc5N8\nuc3vnuSsJJck+UySrVv7Nm1+dVu+eKiaJUnDHnG8DrhoZP7dwPuraglwPXB4az8cuL6qHg68v/WT\nJA1kkOBIsgjYD/hwmw/wTODk1mUFcECb3r/N05bv0/pLkgYw1BHHPwBHAXe2+Z2BG6rqjja/BljY\nphcCVwK05Te2/neT5Igkq5KsWrdu3UzWLklbtFkPjiQvAK6pqnNGm8fpWhuw7K6GquOramlVLV2w\nYME0VCpJGs+8Afb5FOBFSZ4P3A/Ynu4IZH6See2oYhGwtvVfA+wGrEkyD9gBuG72y5YkwQBHHFV1\ndFUtqqrFwEHAGVV1MPBN4MWt2zLgi216ZZunLT+jqu5xxCFJmh1z6XccbwLekGQ13TmME1r7CcDO\nrf0NwPKB6pMkMcxQ1e9U1beAb7XpS4G9xunzK+DAWS1MkjShuXTEIUnaBBgckqReBh2q0tyxePmp\nQ5cgaRPhEYckqReDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8Eh\nSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJwSJJ6mfXgSLJb\nkm8muSjJhUle19p3SnJakkva846tPUmOS7I6yflJ9pztmiVJdxniiOMO4K+q6tHAE4EjkzwGWA6c\nXlVLgNPbPMC+wJL2OAL40OyXLElab9aDo6quqqrvt+mbgYuAhcD+wIrWbQVwQJveH/hYdc4E5id5\nyCyXLUlqBj3HkWQx8ATgLODBVXUVdOECPKh1WwhcObLamtY2dltHJFmVZNW6detmsmxJ2qINFhxJ\nHgD8C/CXVXXTZF3Haat7NFQdX1VLq2rpggULpqtMSdIYgwRHkvvShcYnq+qU1nz1+iGo9nxNa18D\n7Day+iJg7WzVKkm6uyG+VRXgBOCiqnrfyKKVwLI2vQz44kj7Ie3bVU8Eblw/pCVJmn3zBtjnU4A/\nA36Y5Aet7c3AscBnkxwOXAEc2JZ9BXg+sBq4DThsdsuVJI2a9eCoqu8y/nkLgH3G6V/AkTNalCRp\ng/nLcUlSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSL0NcckSSZsXi5acO\ntu/Ljt1vsH3PNI84JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxa/jziFDfnVQkjaURxySpF48\n4pCkGTDUCMJs/PDQIw5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknrZZIIjyfOSXJxkdZLlQ9cj\nSVuqTeJ3HEm2Av4JeDawBjg7ycqq+tFM7M9fcEvSxDaVI469gNVVdWlV3Q58Gth/4JokaYu0SRxx\nAAuBK0fm1wB7j3ZIcgRwRJu9JcnFwC7AtbNS4abJ92divjcT872Z2ODvTd59r1Z/6IZ02lSCI+O0\n1d1mqo4Hjr/bSsmqqlo6k4Vtynx/JuZ7MzHfm4ltKe/NpjJUtQbYbWR+EbB2oFokaYu2qQTH2cCS\nJLsn2Ro4CFg5cE2StEXaJIaqquqOJH8BfB3YCjixqi7cgFWPn7rLFs33Z2K+NxPzvZnYFvHepKqm\n7iVJUrOpDFVJkuYIg0OS1MtmHRxJtkpybpIvD13LXJLksiQ/TPKDJKuGrmcuSTI/yclJfpzkoiRP\nGrqmuSLJI9u/mfWPm5L85dB1zRVJXp/kwiQXJDkpyf2GrmmmbNbnOJK8AVgKbF9VLxi6nrkiyWXA\n0qryR1xjJFkBfKeqPty+wbdtVd0wdF1zTbsM0M+Avavq8qHrGVqShcB3gcdU1S+TfBb4SlV9dNjK\nZsZme8SRZBGwH/DhoWvRpiHJ9sDTgRMAqup2Q2NC+wA/MTTuZh5w/yTzgG3ZjH9rttkGB/APwFHA\nnUMXMgcV8I0k57RLtajzMGAd8JE2xPnhJNsNXdQcdRBw0tBFzBVV9TPgvcAVwFXAjVX1jWGrmjmb\nZXAkeQFwTVWdM3Qtc9RTqmpPYF/gyCRPH7qgOWIesCfwoap6AnAr4CX8x2hDeC8CPjd0LXNFkh3p\nLry6O7ArsF2SPx22qpmzWQYH8BTgRW0s/9PAM5N8YtiS5o6qWtuerwE+T3f1YXWXtllTVWe1+ZPp\ngkR3ty/w/aq6euhC5pBnAT+tqnVV9RvgFODJA9c0YzbL4Kiqo6tqUVUtpjukPqOqNtv07yPJdkke\nuH4aeA5wwbBVzQ1V9XPgyiSPbE37ADNyz5dN3EtxmGqsK4AnJtk2Sej+7Vw0cE0zZpO45Iim1YOB\nz3f/tpkHfKqqvjZsSXPKa4BPtuGYS4HDBq5nTkmyLd0N1V45dC1zSVWdleRk4PvAHcC5bMaXH9ms\nv44rSZp+m+VQlSRp5hgckqReDA5JUi8GhySpF4NDktSLwSENIMniJH+yEevNT/Lqkfld29dApVlj\ncEjDWAyMGxztInkTmQ/8Ljiqam1VvXh6S5MmZ3Boi5LkkCTnJzkvyceTPDTJ6a3t9CS/1/p9NMlx\nSb6X5NIkLx7ZxlHtfibnJTm2te2R5GvtwpHfSfKoKbZzLPC0dl+L1yc5NMnnknyJ7gKUD2j1fL/t\na/+R9fZo672nHblc0PZ1vyQfaf3PTfLfW/uhSU5p9V2S5O9n5c3W5quqfPjYIh7A7wMXA7u0+Z2A\nLwHL2vzLgC+06Y/SXcTvPsBjgNWtfV/ge3T36QDYqT2fDixp03vTXeZmsu08A/jySG2H0l0ra/32\n5tHdRwZgF2A1ELojlQtG1vvdPPBXwEfa9KPoLoNxv7btS4Ed2vzlwG5D//fwsek+vOSItiTPBE6u\ndgOrqrqu3eHvf7TlHwdG/xr/QlXdCfwoyYNb27PoPpxvG9nGA+guaPe5dikXgG2m2M54Tquq69p0\ngHe2KxffCSyku1zMZJ4KfKDV9eMklwOPaMtOr6obAZL8CHgocOUU25PGZXBoSxK6e5FMZnT5r8es\nO9E27gPcUFWPn2Cb421nPLeOTB8MLAD+qKp+0670PNWtSCfb9mgNv8X/93UveI5DW5LTgZck2Rkg\nyU50w04HteUH093+czLfAF7WLvZHkp2q6ibgp0kObG1J8rgptnMz8MBJlu9Ad0+Z37RzFQ/dgPW+\n3V4DSR4B/B7d0Jw0rQwObTGq6kLgHcC/JTkPeB/wWuCwJOcDfwa8boptfA1YCaxK8gPgjW3RwcDh\nbbsX0t3UZzLnA3e0E+yvH2f5J4GlSVa1bf+47f8XwL8nuSDJe8as80FgqyQ/BD4DHFpVv0aaZl4d\nV5LUi0cckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknr5/1TqS69Z5sPcAAAAAElFTkSu\nQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGjlJREFUeJzt3XvUHXV97/H3RyIiKIRLtBKoUcRb\nT73QHARtLceoLaCGdQotlpZAUaxatVpr0XPqrauWVpe22KM9KCreEKVe8FoRpK1t4RhEEEQWEQJE\nEIJAEFAR/Z4/5hfZPDy3SZ5nz5Pk/Vprrz3zm9/MfPfkyXz2zOw9O1WFJEmzdb+hC5AkbVkMDklS\nLwaHJKkXg0OS1IvBIUnqxeCQJPVicGhBS/K6JO+do2VVkkdtxvxfTLJqLmrZUiV5WpIrktye5LAx\nrveYJF8b1/o0PYNDs5JkbZK7kuwxof2bbYe8bD7WW1VvqaoXzMeyRyU5N8kLJrQdlGTdSC0HV9Wp\ns1jWZgXUAvdm4B+r6kFV9emJE9vfyQ1Jdhppe0GSc6dbaJKHJDktyXVJNiT5jyRPmfvyNRcMDvVx\nFfD8jSNJfhV44KYuLMmiuShqW7IAttnDgUtn6LMIeEXP5T4I+Drwa8BuwKnA55M8qHeFmncGh/r4\nEHD0yPgq4IOjHZIcmuTCJLcluTbJG0emLWvvxo9Lcg1wTms/OsnVSX6Q5C/bu9ZntmlvTPLhCfOv\nSnJNkpuS/K+R5e+f5L+S3Jrk+iT/mGT7uXrxo0clSR6V5F/bu+Obkpze2v+tdb+onc75vdb+wiRr\nktyc5Mwke44s99lJLm/Leldb7sb1HNPefb8jyc3AG5Psk+Sctr1uSvKRJItHlrc2yZ8nuTjJHUlO\nSfLQdqrth0m+kmTXaV7npLUm+S7wSOCz7bU9YIpFvBV49WhNM6mqK6vq7VV1fVX9rKpOBrYHHjNF\njW9N8rUku8x2HZo7Bof6OA/YOcnjkmwH/B7w4Ql97qALl8XAocCLJzkX/pvA44DfSvJ44F3AUcDD\ngF2ApTPU8et0O5QVwOuTPK61/wx4JbAHcGCb/pK+L3KW/gr4MrArsBfwToCqenqb/sR2Ouf0JM8A\n/gb4XbrXeDXwMYB26u8M4LXA7sDlwFMnrOspwJXAQ4C/BtKWtyfddtwbeOOEeX4HeBbwaOC5wBeB\n19Ftm/sBL5/sRU1Xa1XtA1wDPLe9tp9MsW1WA+cCr55i+oySPIkuONZMaL9fkvcATwCeXVUbNnUd\n2nQGh/raeNTxLOA7wPdGJ1bVuVX1rar6eVVdDJxGFxSj3lhVd1TVj4DDgc9W1deq6i7g9cBMN1B7\nU1X9qKouAi4CntjWfUFVnVdVd1fVWuD/TrLu6ZzUjlZuTXIr8Llp+v6U7rTNnlX146qa7sLtUcD7\nquobbWf7WuDAdl3oEODSqvpkVd0NnAR8f8L811XVO9vr+lFVramqs6rqJ1W1Hnj7JK/znVV1Q1V9\nD/h34PyqurCt/1PAkzeh1j5eD7wsyZKe85FkZ7q/szdNCIb70/097UYXXnf2XbbmhsGhvj4E/D5w\nDBNOUwEkeUqSryZZn2QD8Md073JHXTsyvOfoeNsZ/GCGGkZ3rHfSnR8nyaOTfC7J95PcBrxlknVP\n5+VVtXjjA3jONH1fQ/fO//8luTTJH03Td0+6d+4AVNXtdK9xKfd9/QWsmzD/6PbaeCH5Y0m+117n\nh7nv67xhZPhHk4xPde1gulpnraouoQveE/rMl+SBwGeB86rqbyZMfhSwki5Q7uqzXM0tg0O9VNXV\ndBfJDwE+OUmXjwJnAntX1S7AP9HtYO+1mJHh6+lO9QC/2HHsvonlvZvuKGjfqtqZ7tTMxHXPiar6\nflW9sKr2BF4EvCtTf5LqOrqjEwDSfeJod7qjtYmvP6PjG1c3YfxvWtsT2uv8A+budU5Xa19vAF7I\nLEOnXTP5dFvXiybpchlwLPDFJJNe+9B4GBzaFMcBz6iqOyaZ9mDg5qr6cZL96Y5OpnMG8NwkT20X\nst/Epu8EHwzcBtye5LHAizdxOTNKckSSjTv4W+h25D9r4zfQXUTe6KPAsUme1HaOb6E7dbQW+Dzw\nq0kOS/eJqZcCvzTD6h8M3A7cmmQp8Odz8ZpmUWsvVbUGOJ0prqeMSnJ/ur+FHwFHV9XPp1jmaXRv\nCL6SZJ++NWluGBzqraq+W1Wrp5j8EuDNSX5Id5774zMs61LgZXQXYK8HfgjcCEx14XU6r6YLqh8C\n76Hbac2X/w6cn+R2uiOsV1TVVW3aG4FT27WS362qs4G/BP6Z7jXuAxwJUFU3AUcAf0d3SujxdBeX\np3v9bwL2AzbQBc9kR36bZLpaN9GbgZ1m7NV9IOA5wLPpAvH29viNSWo8tS33nE249qI5EH/ISQtJ\nus/t30p3uumqmfpvbZLcj+4ax1FV9dWh65Em4xGHBpfkuUl2bOfT3wZ8C1g7bFXjk+S3kixup4Y2\nXpc5b+CypCkZHFoIVtJdlL0O2Bc4sratQ+EDge8CN9F95+Kw9lHlrUqS3xg5BXWvx9C1qR9PVUmS\nevGIQ5LUy9A3TJsXe+yxRy1btmzoMiRpi3LBBRfcVFUzftt/qwyOZcuWsXr1VJ8WlSRNJsnVM/fy\nVJUkqSeDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqZet8pvj0kyWnfD5wda9\n9sRDB1u3NBc84pAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIv8xYcSd6X5MYkl4y07ZbkrCRXtOddW3uS\nnJRkTZKLk+w3Ms+q1v+KJKvmq15J0uzM5xHHB4DfntB2AnB2Ve0LnN3GAQ4G9m2P44F3Qxc0wBuA\npwD7A2/YGDaSpGHMW3BU1b8BN09oXgmc2oZPBQ4baf9gdc4DFid5GPBbwFlVdXNV3QKcxX3DSJI0\nRuO+xvHQqroeoD0/pLUvBa4d6beutU3Vfh9Jjk+yOsnq9evXz3nhkqTOQrk4nknaapr2+zZWnVxV\ny6tq+ZIlM/7WuiRpE407OG5op6Bozze29nXA3iP99gKum6ZdkjSQcd+r6kxgFXBie/7MSPufJPkY\n3YXwDVV1fZJ/Ad4yckH82cBrx1yz5tGQ94yStGnmLTiSnAYcBOyRZB3dp6NOBD6e5DjgGuCI1v0L\nwCHAGuBO4FiAqro5yV8BX2/93lxVEy+4S5LGaN6Co6qeP8WkFZP0LeClUyznfcD75rA0SdJmWCgX\nxyVJWwiDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqZdx36tK2uYNdX+utSce\nOsh6tfXxiEOS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQ\nJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9DBIcSV6Z\n5NIklyQ5LckOSR6R5PwkVyQ5Pcn2re8D2viaNn3ZEDVLkjpjD44kS4GXA8ur6r8B2wFHAn8LvKOq\n9gVuAY5rsxwH3FJVjwLe0fpJkgYy1KmqRcADkywCdgSuB54BnNGmnwoc1oZXtnHa9BVJMsZaJUkj\nxh4cVfU94G3ANXSBsQG4ALi1qu5u3dYBS9vwUuDaNu/drf/uE5eb5Pgkq5OsXr9+/fy+CEnahg1x\nqmpXuqOIRwB7AjsBB0/StTbOMs20exqqTq6q5VW1fMmSJXNVriRpgiFOVT0TuKqq1lfVT4FPAk8F\nFrdTVwB7Ade14XXA3gBt+i7AzeMtWZK00RDBcQ1wQJId27WKFcC3ga8Ch7c+q4DPtOEz2zht+jlV\ndZ8jDknSeAxxjeN8uovc3wC+1Wo4GfgL4FVJ1tBdwzilzXIKsHtrfxVwwrhrliTdY9HMXeZeVb0B\neMOE5iuB/Sfp+2PgiHHUJUmamd8clyT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ\n6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgc\nkqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpl0GCI8ni\nJGck+U6Sy5IcmGS3JGcluaI979r6JslJSdYkuTjJfkPULEnqDHXE8Q/Al6rqscATgcuAE4Czq2pf\n4Ow2DnAwsG97HA+8e/zlSpI2WjTuFSbZGXg6cAxAVd0F3JVkJXBQ63YqcC7wF8BK4INVVcB57Wjl\nYVV1/ZhL32otO+HzQ5cgaQsyxBHHI4H1wPuTXJjkvUl2Ah66MQza80Na/6XAtSPzr2tt95Lk+CSr\nk6xev379/L4CSdqGjf2Io61zP+BlVXV+kn/gntNSk8kkbXWfhqqTgZMBli9ffp/p0rZuyCPLtSce\nOti6NfeGOOJYB6yrqvPb+Bl0QXJDkocBtOcbR/rvPTL/XsB1Y6pVkjTB2IOjqr4PXJvkMa1pBfBt\n4ExgVWtbBXymDZ8JHN0+XXUAsMHrG5I0nCFOVQG8DPhIku2BK4Fj6ULs40mOA64Bjmh9vwAcAqwB\n7mx9JUkDGSQ4quqbwPJJJq2YpG8BL533oiRJs+I3xyVJvRgckqReZhUcSf73yPAD5q8cSdJCN21w\nJHlNkgOBw0ea/2t+S5IkLWQzXRy/nO7TTY9M8u9095TaPcljquryea9OkrTgzHSq6hbgdXQfhT0I\nOKm1n5DkP+exLknSAjXTEcdvA28A9gHeDlwE3FFVfpdCkrZR0x5xVNXrqmoFsBb4MF3QLEnytSSf\nHUN9kqQFZrZfAPyXqvo68PUkL66qX0+yx3wWJklamGb1cdyqes3I6DGt7ab5KEiStLD1/gJgVV00\nH4VIkrYMfnNcktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXg\nkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoZLDiSbJfkwiSf\na+OPSHJ+kiuSnJ5k+9b+gDa+pk1fNlTNkqRhjzheAVw2Mv63wDuqal/gFuC41n4ccEtVPQp4R+sn\nSRrIIMGRZC/gUOC9bTzAM4AzWpdTgcPa8Mo2Tpu+ovWXJA1gqCOOvwdeA/y8je8O3FpVd7fxdcDS\nNrwUuBagTd/Q+t9LkuOTrE6yev369fNZuyRt08YeHEmeA9xYVReMNk/StWYx7Z6GqpOranlVLV+y\nZMkcVCpJmsyiAdb5NOB5SQ4BdgB2pjsCWZxkUTuq2Au4rvVfB+wNrEuyCNgFuHn8ZUuSYIAjjqp6\nbVXtVVXLgCOBc6rqKOCrwOGt2yrgM234zDZOm35OVd3niEOSNB4L6XscfwG8KskaumsYp7T2U4Dd\nW/urgBMGqk+SxDCnqn6hqs4Fzm3DVwL7T9Lnx8ARYy1MkjSlhXTEIUnaAhgckqReDA5JUi8GhySp\nF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6mXQ26rr3pad\n8PmhS5DmxVB/22tPPHSQ9W7tPOKQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnq\nxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1MvYgyPJ3km+muSyJJcmeUVr3y3JWUmu\naM+7tvYkOSnJmiQXJ9lv3DVLku4xxBHH3cCfVdXjgAOAlyZ5PHACcHZV7Quc3cYBDgb2bY/jgXeP\nv2RJ0kZjD46qur6qvtGGfwhcBiwFVgKntm6nAoe14ZXAB6tzHrA4ycPGXLYkqRn0GkeSZcCTgfOB\nh1bV9dCFC/CQ1m0pcO3IbOta28RlHZ9kdZLV69evn8+yJWmbNlhwJHkQ8M/An1bVbdN1naSt7tNQ\ndXJVLa+q5UuWLJmrMiVJEwwSHEnuTxcaH6mqT7bmGzaegmrPN7b2dcDeI7PvBVw3rlolSfc2xKeq\nApwCXFZVbx+ZdCawqg2vAj4z0n50+3TVAcCGjae0JEnjt2iAdT4N+EPgW0m+2dpeB5wIfDzJccA1\nwBFt2heAQ4A1wJ3AseMtV5I0auzBUVVfY/LrFgArJulfwEvntShJ0qz5zXFJUi8GhySpF4NDktSL\nwSFJ6sXgkCT1YnBIknoxOCRJvRgckqRehvjm+IK37ITPD12CJC1YHnFIknoxOCRJvRgckqReDA5J\nUi8GhySpF4NDktSLwSFJ6sXgkCT14hcAJW21hvoy79oTDx1kvePiEYckqReDQ5LUi8EhSerF4JAk\n9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerFe1VJ0hwb6h5ZMJ77ZHnEIUnqZYsJ\njiS/neTyJGuSnDB0PZK0rdoigiPJdsD/AQ4GHg88P8njh61KkrZNW0RwAPsDa6rqyqq6C/gYsHLg\nmiRpm7SlXBxfClw7Mr4OeMpohyTHA8e30duTXD5hGXsAN81bhZvP+jaP9W0e69s8C6a+/O2kzbOt\n7+GzWceWEhyZpK3uNVJ1MnDylAtIVlfV8rkubK5Y3+axvs1jfZtnW6tvSzlVtQ7Ye2R8L+C6gWqR\npG3alhIcXwf2TfKIJNsDRwJnDlyTJG2TtohTVVV1d5I/Af4F2A54X1Vd2nMxU57GWiCsb/NY3+ax\nvs2zTdWXqpq5lyRJzZZyqkqStEAYHJKkXra64EiyXZILk3xukmnHJFmf5Jvt8YIB6lub5Ftt/asn\nmZ4kJ7Vbq1ycZL8FVt9BSTaMbMPXj7m+xUnOSPKdJJclOXDC9MG23yxqG3rbPWZk3d9McluSP53Q\nZ8jtN5v6BtuGSV6Z5NIklyQ5LckOE6Y/IMnpbdudn2TZuGqbZX1zt/+rqq3qAbwK+CjwuUmmHQP8\n48D1rQX2mGb6IcAX6b67cgBw/gKr76DJtu0Y6zsVeEEb3h5YvFC23yxqG3TbTahlO+D7wMMXyvab\nZX2DbEO6LyFfBTywjX8cOGZCn5cA/9SGjwROX2D1zdn+b6s64kiyF3Ao8N6ha9kMK4EPVuc8YHGS\nhw1d1EKQZGfg6cApAFV1V1XdOqHbINtvlrUtJCuA71bV1RPaF8rf31T1DWkR8MAki4Adue93yVbS\nvXkAOANYkWSyLy8PVd+c2aqCA/h74DXAz6fp8zvtEPyMJHtP02++FPDlJBe026RMNNntVZaOpbLO\nTPUBHJjkoiRfTPIrY6ztkcB64P3tdOR7k+w0oc9Q2282tcFw226iI4HTJmkf+u9vo6nqgwG2YVV9\nD3gbcA1wPbChqr48odsvtl1V3Q1sAHZfQPXBHO3/tprgSPIc4MaqumCabp8FllXVE4CvcM+7g3F6\nWlXtR3en35cmefqE6TPeXmWezVTfN+hOHzwReCfw6THWtgjYD3h3VT0ZuAOYeIv9obbfbGobctv9\nQrov0T4P+MRkkydpG+tn9meob5BtmGRXuiOKRwB7Ajsl+YOJ3SaZdSzbbpb1zdn+b6sJDuBpwPOS\nrKW7e+4zknx4tENV/aCqftJG3wP82nhLhKq6rj3fCHyK7s6/owa9vcpM9VXVbVV1exv+AnD/JHuM\nqbx1wLqqOr+Nn0G3s57YZ4jtN2NtA2+7UQcD36iqGyaZthBu7zNlfQNuw2cCV1XV+qr6KfBJ4KkT\n+vxi27XTRbsAN4+htlnVN5f7v60mOKrqtVW1V1UtozvMPaeq7pW4E87VPg+4bIwlkmSnJA/eOAw8\nG7hkQrczgaPbp1sOoDvkvH6h1Jfklzaet02yP93f0A/GUV9VfR+4NsljWtMK4NsTug2y/WZT25Db\nboLnM/VpoMH+/kZMWd+A2/Aa4IAkO7b1r+C++48zgVVt+HC6fdC4jtZmrG8u939bxC1HNkeSNwOr\nq+pM4OVJngfcTfdO4Jgxl/NQ4FPt734R8NGq+lKSPwaoqn8CvkD3yZY1wJ3AsQusvsOBFye5G/gR\ncOQY/3MAvAz4SDudcSVw7ALafjPVNvS2I8mOwLOAF420LZTtN5v6BtmGVXV+kjPoTpXdDVwInDxh\n/3IK8KEka+j2L0fOd10965uz/Z+3HJEk9bLVnKqSJI2HwSFJ6sXgkCT1YnBIknoxOCRJvRgc0gCS\nLEvy+5sw3+IkLxkZ37N9DFMaG4NDGsYyYNLgaN86nspiuruwAt03/avq8LktTZqewaFtSpKj203e\nLkryoSQPT3J2azs7yS+3fh9I97sU/5nkyiSHjyzjNel+s+SiJCe2tn2SfCndzSH/PcljZ1jOicBv\npPtdhFem+62ETyT5LN1NJh/U6vlGW9fKkfn2afO9tR25XNLWtUOS97f+Fyb5H639mCSfbPVdkeTv\nxrKxtfWai3uz+/CxJTyAXwEup/3eCLAb3Y3fVrXxPwI+3YY/QHeTvfsBjwfWtPaDgf8Edty4jPZ8\nNrBvG34K3e0mplvOQYz8rgTdt3jXjSxvEbBzG96D7pvcoTtSuWRkvl+MA38GvL8NP5buNhQ7tGVf\nSXfvpB2Aq4G9h/738LHlPrb6W45II54BnFFVNwFU1c3pfqXvf7bpHwJG341/uqp+Dnw7yUNb2zPp\nds53jizjQXQ3lPtE7vn5hQfMsJzJnFVVG2+KF+At6e5O/HO6W3ZPNy/Ar9PdMZaq+k6Sq4FHt2ln\nV9UGgCTfBh7OvW+fLs2awaFtSZj5Ntej038yMpyR54nLuB9wa1U9aYplTracydwxMnwUsAT4tar6\nabq7Pu8w6VyzW/ZoDT/D//vaDF7j0LbkbOB3k+wOkGQ3utNOG29GdxTwtRmW8WXgj9rN+EiyW1Xd\nBlyV5IjWliRPnGE5PwQePM30Xeh+X+an7VrFw2cx37+110CSRwO/THdqTppTBoe2GVV1KfDXwL8m\nuQh4O/ByurvYXgz8IfCKGZbxJbrbZ69O8k3g1W3SUcBxbbmX0v2oznQuBu5uF9hfOcn0jwDLk6xu\ny/5OW/8PgP9IckmSt06Y513Adkm+BZxO95vTP0GaY94dV5LUi0cckqReDA5JUi8GhySpF4NDktSL\nwSFJ6sXgkCT1YnBIknr5/4o4RF0/jCAaAAAAAElFTkSuQmCC\n", "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -137,8 +152,10 @@ }, { "cell_type": "code", - "execution_count": 21, - "metadata": {}, + "execution_count": 27, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "do_map3k = pyro.do(mapk_receiver, data={'map3k': 5.678})" @@ -153,7 +170,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 28, "metadata": {}, "outputs": [ { @@ -163,7 +180,7 @@ "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmapk_do_dist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minfer_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdo_map3k\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnoise_marginals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mmapk_marginal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mEmpiricalMarginal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmapk_do_dist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msites\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'mapk'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmapk_do_dist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minfer_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdo_map3k\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnoise_marginals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mmapk_marginal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mEmpiricalMarginal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmapk_do_dist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msites\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'mapk'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m~/Desktop/Git/causal_demon/mapk_demon/inference.py\u001b[0m in \u001b[0;36minfer_dist\u001b[0;34m(prog, n_dist, type)\u001b[0m\n\u001b[1;32m 19\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mtype\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'mcmc'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 20\u001b[0m \u001b[0mhmc_kernel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mHMC\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprog\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstep_size\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.9\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m \u001b[0mposterior\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMCMC\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhmc_kernel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_samples\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwarmup_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m50\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_dist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 22\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mposterior\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/infer/abstract_infer.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 81\u001b[0m \"\"\"\n\u001b[1;32m 82\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_init\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 83\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mtr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlogit\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpoutine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mblock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_traces\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 84\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexec_traces\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 85\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_weights\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlogit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/infer/mcmc/mcmc.py\u001b[0m in \u001b[0;36m_traces\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 35\u001b[0m \u001b[0mlogging_interval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mceil\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_samples\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_samples\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 37\u001b[0;31m \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkernel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 38\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mlogging_interval\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 39\u001b[0m \u001b[0mstage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"WARMUP\"\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m\"SAMPLE\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", From e22069a8b988c5d70636b8c99d4e8a5858270011 Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Sat, 16 Feb 2019 14:21:17 -0500 Subject: [PATCH 10/11] counterfactual on companion model --- mapk_demon/receivers.py | 100 ++++++++++++-- notebooks/Mapk Counterfactual Inference.ipynb | 128 ++++++++++-------- 2 files changed, 161 insertions(+), 67 deletions(-) diff --git a/mapk_demon/receivers.py b/mapk_demon/receivers.py index 05eabfe..083ccf1 100644 --- a/mapk_demon/receivers.py +++ b/mapk_demon/receivers.py @@ -4,7 +4,7 @@ import pyro from pyro import sample -from pyro.distributions import Normal, Uniform +from pyro.distributions import Normal, Uniform, Delta import pyro.optim import torch @@ -29,40 +29,83 @@ def f(mu, N): } -def f_map3k(E1, N_3k, params): +def f_map3k(E1, N_3k, params, mode): total_3k = hyperparameters['t_3k'] alpha_3k = params['alpha_3k'].rsample() nu_3k = params['nu_3k'].rsample() map3k_mu = total_3k * g1(E1 * (alpha_3k/nu_3k)) - map3k = sample("map3k", Normal(f(map3k_mu, N_3k), 1.)) - return map3k + if mode == 'companion': + return sample("map3k", Normal(f(map3k_mu, N_3k), 1.)) + else: + return sample("map3k", Delta(f(map3k_mu, N_3k))) -def f_map2k(map3k, N_2k, params): +def f_map2k(map3k, N_2k, params, mode): total_2k = hyperparameters['t_2k'] alpha_2k = params['alpha_2k'].rsample() nu_2k = params['nu_2k'].rsample() map2k_mu = total_2k * g2(map3k * (alpha_2k / nu_2k)) - map2k = sample("map2k", Normal(f(map2k_mu, N_2k), 1.)) - return map2k + if mode == 'companion': + return sample("map2k", Normal(f(map2k_mu, N_2k), 1.)) + else: + return sample("map2k", Delta(f(map2k_mu, N_2k))) -def f_mapk(map2k, N_k, params): +def f_mapk(map2k, N_k, params, mode): total_k = hyperparameters['t_k'] alpha_k = params['alpha_k'].rsample() nu_k = params['nu_k'].rsample() mapk_mu = total_k * g2(map2k * (alpha_k / nu_k)) - mapk = sample("mapk", Normal(f(mapk_mu, N_k), 1.)) - return mapk + + if mode == 'companion': + return sample("mapk", Normal(f(mapk_mu, N_k), 1.)) + else: + return sample("mapk", Delta(f(mapk_mu, N_k))) + + +def mapk_companion(noise_dists): + mode = 'companion' + # Parameters + al_m = torch.tensor(700.) + al_s = torch.tensor(1.) + + nu_m = torch.tensor(.15) + nu_s = torch.tensor(.05) + + params = { + "alpha_3k": Normal(al_m, al_s), + "alpha_2k": Normal(al_m, al_s), + "alpha_k": Normal(al_m, al_s), + "nu_3k": Normal(nu_m, nu_s), + "nu_2k": Normal(nu_m, nu_s), + "nu_k": Normal(nu_m, nu_s) + } + + with pyro.iarange("model"): + E1 = Uniform(1.5e-5, 10.).rsample() + N_3k = sample('N_3k', noise_dists['N_3k']) + N_2k = sample('N_2k', noise_dists['N_2k']) + N_k = sample('N_k', noise_dists['N_k']) + + map3k = f_map3k(E1, N_3k, params, mode) + map2k = f_map2k(map3k, N_2k, params, mode) + mapk = f_mapk(map2k, N_k, params, mode) + + return { + 'map3k': map3k, + 'map2k': map2k, + 'mapk': mapk + } def mapk_receiver(noise_dists): # Parameters + mode = 'original' al_m = torch.tensor(700.) al_s = torch.tensor(1.) @@ -84,12 +127,43 @@ def mapk_receiver(noise_dists): N_2k = sample('N_2k', noise_dists['N_2k']) N_k = sample('N_k', noise_dists['N_k']) - map3k = f_map3k(E1, N_3k, params) - map2k = f_map2k(map3k,N_2k, params) - mapk = f_mapk(map2k, N_k, params) + map3k = f_map3k(E1, N_3k, params, mode) + map2k = f_map2k(map3k, N_2k, params, mode) + mapk = f_mapk(map2k, N_k, params, mode) return { 'map3k': map3k, 'map2k': map2k, 'mapk': mapk } + + +def mapk_do_receiver(do_value,noise_dists): + mode = 'original' + # Parameters + al_m = torch.tensor(700.) + al_s = torch.tensor(1.) + + nu_m = torch.tensor(.15) + nu_s = torch.tensor(.05) + + params = { + "alpha_3k": Normal(al_m, al_s), + "alpha_2k": Normal(al_m, al_s), + "alpha_k": Normal(al_m, al_s), + "nu_3k": Normal(nu_m, nu_s), + "nu_2k": Normal(nu_m, nu_s), + "nu_k": Normal(nu_m, nu_s) + } + + with pyro.iarange("model"): + E1 = Uniform(1.5e-5, 10.).rsample() + N_3k = sample('N_3k', noise_dists['N_3k']) + N_2k = sample('N_2k', noise_dists['N_2k']) + N_k = sample('N_k', noise_dists['N_k']) + + map3k = do_value + map2k = f_map2k(map3k, N_2k, params, mode) + mapk = f_mapk(map2k, N_k, params, mode) + + return map2k diff --git a/notebooks/Mapk Counterfactual Inference.ipynb b/notebooks/Mapk Counterfactual Inference.ipynb index de5395b..a339259 100644 --- a/notebooks/Mapk Counterfactual Inference.ipynb +++ b/notebooks/Mapk Counterfactual Inference.ipynb @@ -9,10 +9,8 @@ }, { "cell_type": "code", - "execution_count": 2, - "metadata": { - "collapsed": true - }, + "execution_count": 1, + "metadata": {}, "outputs": [], "source": [ "import sys\n", @@ -31,7 +29,7 @@ "from pyro.infer import EmpiricalMarginal\n", "\n", "from mapk_demon.inference import infer_dist\n", - "from mapk_demon.receivers import mapk_receiver\n", + "from mapk_demon.receivers import mapk_receiver, mapk_companion, mapk_do_receiver\n", "\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", @@ -55,25 +53,14 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 3, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "{'map2k': tensor(12.9241), 'map3k': tensor(12.4646), 'mapk': tensor(-0.2236)}" - ] - }, - "execution_count": 16, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "noise_vars = ['N_3k', 'N_2k', 'N_k']\n", "noise_dists = {N: Normal(10., 1.) for N in noise_vars}\n", "\n", - "mapk_receiver(noise_dists)" + "evidence = mapk_receiver(noise_dists)" ] }, { @@ -85,20 +72,18 @@ }, { "cell_type": "code", - "execution_count": 17, - "metadata": { - "collapsed": true - }, + "execution_count": 4, + "metadata": {}, "outputs": [], "source": [ - "conditioned_mapk = pyro.condition(mapk_receiver, data={'map2k': 13.8113, 'map3k': 11.0504, 'mapk': 0.2803})" + "conditioned_mapk = pyro.condition(mapk_companion, data=evidence)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## 2. Infer the noise distributions N_3k, N_2k, N_k\n", + "## 2. Infer the noise distributions N_3k, N_2k, N_k from companion model\n", "\n", "Let there be a prior over Noise distributions. $$N_{3k}, N_{2k}, N_k \\sim Normal\\ (0, 1)$$.\n", "\n", @@ -107,7 +92,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 5, "metadata": { "collapsed": true }, @@ -125,14 +110,14 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 6, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGjlJREFUeJzt3XvUHXV97/H3RyIiKIRLtBKoUcRb\nT73QHARtLceoLaCGdQotlpZAUaxatVpr0XPqrauWVpe22KM9KCreEKVe8FoRpK1t4RhEEEQWEQJE\nEIJAEFAR/Z4/5hfZPDy3SZ5nz5Pk/Vprrz3zm9/MfPfkyXz2zOw9O1WFJEmzdb+hC5AkbVkMDklS\nLwaHJKkXg0OS1IvBIUnqxeCQJPVicGhBS/K6JO+do2VVkkdtxvxfTLJqLmrZUiV5WpIrktye5LAx\nrveYJF8b1/o0PYNDs5JkbZK7kuwxof2bbYe8bD7WW1VvqaoXzMeyRyU5N8kLJrQdlGTdSC0HV9Wp\ns1jWZgXUAvdm4B+r6kFV9emJE9vfyQ1Jdhppe0GSc6dbaJKHJDktyXVJNiT5jyRPmfvyNRcMDvVx\nFfD8jSNJfhV44KYuLMmiuShqW7IAttnDgUtn6LMIeEXP5T4I+Drwa8BuwKnA55M8qHeFmncGh/r4\nEHD0yPgq4IOjHZIcmuTCJLcluTbJG0emLWvvxo9Lcg1wTms/OsnVSX6Q5C/bu9ZntmlvTPLhCfOv\nSnJNkpuS/K+R5e+f5L+S3Jrk+iT/mGT7uXrxo0clSR6V5F/bu+Obkpze2v+tdb+onc75vdb+wiRr\nktyc5Mwke44s99lJLm/Leldb7sb1HNPefb8jyc3AG5Psk+Sctr1uSvKRJItHlrc2yZ8nuTjJHUlO\nSfLQdqrth0m+kmTXaV7npLUm+S7wSOCz7bU9YIpFvBV49WhNM6mqK6vq7VV1fVX9rKpOBrYHHjNF\njW9N8rUku8x2HZo7Bof6OA/YOcnjkmwH/B7w4Ql97qALl8XAocCLJzkX/pvA44DfSvJ44F3AUcDD\ngF2ApTPU8et0O5QVwOuTPK61/wx4JbAHcGCb/pK+L3KW/gr4MrArsBfwToCqenqb/sR2Ouf0JM8A\n/gb4XbrXeDXwMYB26u8M4LXA7sDlwFMnrOspwJXAQ4C/BtKWtyfddtwbeOOEeX4HeBbwaOC5wBeB\n19Ftm/sBL5/sRU1Xa1XtA1wDPLe9tp9MsW1WA+cCr55i+oySPIkuONZMaL9fkvcATwCeXVUbNnUd\n2nQGh/raeNTxLOA7wPdGJ1bVuVX1rar6eVVdDJxGFxSj3lhVd1TVj4DDgc9W1deq6i7g9cBMN1B7\nU1X9qKouAi4CntjWfUFVnVdVd1fVWuD/TrLu6ZzUjlZuTXIr8Llp+v6U7rTNnlX146qa7sLtUcD7\nquobbWf7WuDAdl3oEODSqvpkVd0NnAR8f8L811XVO9vr+lFVramqs6rqJ1W1Hnj7JK/znVV1Q1V9\nD/h34PyqurCt/1PAkzeh1j5eD7wsyZKe85FkZ7q/szdNCIb70/097UYXXnf2XbbmhsGhvj4E/D5w\nDBNOUwEkeUqSryZZn2QD8Md073JHXTsyvOfoeNsZ/GCGGkZ3rHfSnR8nyaOTfC7J95PcBrxlknVP\n5+VVtXjjA3jONH1fQ/fO//8luTTJH03Td0+6d+4AVNXtdK9xKfd9/QWsmzD/6PbaeCH5Y0m+117n\nh7nv67xhZPhHk4xPde1gulpnraouoQveE/rMl+SBwGeB86rqbyZMfhSwki5Q7uqzXM0tg0O9VNXV\ndBfJDwE+OUmXjwJnAntX1S7AP9HtYO+1mJHh6+lO9QC/2HHsvonlvZvuKGjfqtqZ7tTMxHXPiar6\nflW9sKr2BF4EvCtTf5LqOrqjEwDSfeJod7qjtYmvP6PjG1c3YfxvWtsT2uv8A+budU5Xa19vAF7I\nLEOnXTP5dFvXiybpchlwLPDFJJNe+9B4GBzaFMcBz6iqOyaZ9mDg5qr6cZL96Y5OpnMG8NwkT20X\nst/Epu8EHwzcBtye5LHAizdxOTNKckSSjTv4W+h25D9r4zfQXUTe6KPAsUme1HaOb6E7dbQW+Dzw\nq0kOS/eJqZcCvzTD6h8M3A7cmmQp8Odz8ZpmUWsvVbUGOJ0prqeMSnJ/ur+FHwFHV9XPp1jmaXRv\nCL6SZJ++NWluGBzqraq+W1Wrp5j8EuDNSX5Id5774zMs61LgZXQXYK8HfgjcCEx14XU6r6YLqh8C\n76Hbac2X/w6cn+R2uiOsV1TVVW3aG4FT27WS362qs4G/BP6Z7jXuAxwJUFU3AUcAf0d3SujxdBeX\np3v9bwL2AzbQBc9kR36bZLpaN9GbgZ1m7NV9IOA5wLPpAvH29viNSWo8tS33nE249qI5EH/ISQtJ\nus/t30p3uumqmfpvbZLcj+4ax1FV9dWh65Em4xGHBpfkuUl2bOfT3wZ8C1g7bFXjk+S3kixup4Y2\nXpc5b+CypCkZHFoIVtJdlL0O2Bc4sratQ+EDge8CN9F95+Kw9lHlrUqS3xg5BXWvx9C1qR9PVUmS\nevGIQ5LUy9A3TJsXe+yxRy1btmzoMiRpi3LBBRfcVFUzftt/qwyOZcuWsXr1VJ8WlSRNJsnVM/fy\nVJUkqSeDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqZet8pvj0kyWnfD5wda9\n9sRDB1u3NBc84pAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIv8xYcSd6X5MYkl4y07ZbkrCRXtOddW3uS\nnJRkTZKLk+w3Ms+q1v+KJKvmq15J0uzM5xHHB4DfntB2AnB2Ve0LnN3GAQ4G9m2P44F3Qxc0wBuA\npwD7A2/YGDaSpGHMW3BU1b8BN09oXgmc2oZPBQ4baf9gdc4DFid5GPBbwFlVdXNV3QKcxX3DSJI0\nRuO+xvHQqroeoD0/pLUvBa4d6beutU3Vfh9Jjk+yOsnq9evXz3nhkqTOQrk4nknaapr2+zZWnVxV\ny6tq+ZIlM/7WuiRpE407OG5op6Bozze29nXA3iP99gKum6ZdkjSQcd+r6kxgFXBie/7MSPufJPkY\n3YXwDVV1fZJ/Ad4yckH82cBrx1yz5tGQ94yStGnmLTiSnAYcBOyRZB3dp6NOBD6e5DjgGuCI1v0L\nwCHAGuBO4FiAqro5yV8BX2/93lxVEy+4S5LGaN6Co6qeP8WkFZP0LeClUyznfcD75rA0SdJmWCgX\nxyVJWwiDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqZdx36tK2uYNdX+utSce\nOsh6tfXxiEOS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQ\nJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9DBIcSV6Z\n5NIklyQ5LckOSR6R5PwkVyQ5Pcn2re8D2viaNn3ZEDVLkjpjD44kS4GXA8ur6r8B2wFHAn8LvKOq\n9gVuAY5rsxwH3FJVjwLe0fpJkgYy1KmqRcADkywCdgSuB54BnNGmnwoc1oZXtnHa9BVJMsZaJUkj\nxh4cVfU94G3ANXSBsQG4ALi1qu5u3dYBS9vwUuDaNu/drf/uE5eb5Pgkq5OsXr9+/fy+CEnahg1x\nqmpXuqOIRwB7AjsBB0/StTbOMs20exqqTq6q5VW1fMmSJXNVriRpgiFOVT0TuKqq1lfVT4FPAk8F\nFrdTVwB7Ade14XXA3gBt+i7AzeMtWZK00RDBcQ1wQJId27WKFcC3ga8Ch7c+q4DPtOEz2zht+jlV\ndZ8jDknSeAxxjeN8uovc3wC+1Wo4GfgL4FVJ1tBdwzilzXIKsHtrfxVwwrhrliTdY9HMXeZeVb0B\neMOE5iuB/Sfp+2PgiHHUJUmamd8clyT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ\n6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgc\nkqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpl0GCI8ni\nJGck+U6Sy5IcmGS3JGcluaI979r6JslJSdYkuTjJfkPULEnqDHXE8Q/Al6rqscATgcuAE4Czq2pf\n4Ow2DnAwsG97HA+8e/zlSpI2WjTuFSbZGXg6cAxAVd0F3JVkJXBQ63YqcC7wF8BK4INVVcB57Wjl\nYVV1/ZhL32otO+HzQ5cgaQsyxBHHI4H1wPuTXJjkvUl2Ah66MQza80Na/6XAtSPzr2tt95Lk+CSr\nk6xev379/L4CSdqGjf2Io61zP+BlVXV+kn/gntNSk8kkbXWfhqqTgZMBli9ffp/p0rZuyCPLtSce\nOti6NfeGOOJYB6yrqvPb+Bl0QXJDkocBtOcbR/rvPTL/XsB1Y6pVkjTB2IOjqr4PXJvkMa1pBfBt\n4ExgVWtbBXymDZ8JHN0+XXUAsMHrG5I0nCFOVQG8DPhIku2BK4Fj6ULs40mOA64Bjmh9vwAcAqwB\n7mx9JUkDGSQ4quqbwPJJJq2YpG8BL533oiRJs+I3xyVJvRgckqReZhUcSf73yPAD5q8cSdJCN21w\nJHlNkgOBw0ea/2t+S5IkLWQzXRy/nO7TTY9M8u9095TaPcljquryea9OkrTgzHSq6hbgdXQfhT0I\nOKm1n5DkP+exLknSAjXTEcdvA28A9gHeDlwE3FFVfpdCkrZR0x5xVNXrqmoFsBb4MF3QLEnytSSf\nHUN9kqQFZrZfAPyXqvo68PUkL66qX0+yx3wWJklamGb1cdyqes3I6DGt7ab5KEiStLD1/gJgVV00\nH4VIkrYMfnNcktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXg\nkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoZLDiSbJfkwiSf\na+OPSHJ+kiuSnJ5k+9b+gDa+pk1fNlTNkqRhjzheAVw2Mv63wDuqal/gFuC41n4ccEtVPQp4R+sn\nSRrIIMGRZC/gUOC9bTzAM4AzWpdTgcPa8Mo2Tpu+ovWXJA1gqCOOvwdeA/y8je8O3FpVd7fxdcDS\nNrwUuBagTd/Q+t9LkuOTrE6yev369fNZuyRt08YeHEmeA9xYVReMNk/StWYx7Z6GqpOranlVLV+y\nZMkcVCpJmsyiAdb5NOB5SQ4BdgB2pjsCWZxkUTuq2Au4rvVfB+wNrEuyCNgFuHn8ZUuSYIAjjqp6\nbVXtVVXLgCOBc6rqKOCrwOGt2yrgM234zDZOm35OVd3niEOSNB4L6XscfwG8KskaumsYp7T2U4Dd\nW/urgBMGqk+SxDCnqn6hqs4Fzm3DVwL7T9Lnx8ARYy1MkjSlhXTEIUnaAhgckqReDA5JUi8GhySp\nF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6mXQ26rr3pad\n8PmhS5DmxVB/22tPPHSQ9W7tPOKQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnq\nxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1MvYgyPJ3km+muSyJJcmeUVr3y3JWUmu\naM+7tvYkOSnJmiQXJ9lv3DVLku4xxBHH3cCfVdXjgAOAlyZ5PHACcHZV7Quc3cYBDgb2bY/jgXeP\nv2RJ0kZjD46qur6qvtGGfwhcBiwFVgKntm6nAoe14ZXAB6tzHrA4ycPGXLYkqRn0GkeSZcCTgfOB\nh1bV9dCFC/CQ1m0pcO3IbOta28RlHZ9kdZLV69evn8+yJWmbNlhwJHkQ8M/An1bVbdN1naSt7tNQ\ndXJVLa+q5UuWLJmrMiVJEwwSHEnuTxcaH6mqT7bmGzaegmrPN7b2dcDeI7PvBVw3rlolSfc2xKeq\nApwCXFZVbx+ZdCawqg2vAj4z0n50+3TVAcCGjae0JEnjt2iAdT4N+EPgW0m+2dpeB5wIfDzJccA1\nwBFt2heAQ4A1wJ3AseMtV5I0auzBUVVfY/LrFgArJulfwEvntShJ0qz5zXFJUi8GhySpF4NDktSL\nwSFJ6sXgkCT1YnBIknoxOCRJvRgckqRehvjm+IK37ITPD12CJC1YHnFIknoxOCRJvRgckqReDA5J\nUi8GhySpF4NDktSLwSFJ6sXgkCT14hcAJW21hvoy79oTDx1kvePiEYckqReDQ5LUi8EhSerF4JAk\n9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerFe1VJ0hwb6h5ZMJ77ZHnEIUnqZYsJ\njiS/neTyJGuSnDB0PZK0rdoigiPJdsD/AQ4GHg88P8njh61KkrZNW0RwAPsDa6rqyqq6C/gYsHLg\nmiRpm7SlXBxfClw7Mr4OeMpohyTHA8e30duTXD5hGXsAN81bhZvP+jaP9W0e69s8C6a+/O2kzbOt\n7+GzWceWEhyZpK3uNVJ1MnDylAtIVlfV8rkubK5Y3+axvs1jfZtnW6tvSzlVtQ7Ye2R8L+C6gWqR\npG3alhIcXwf2TfKIJNsDRwJnDlyTJG2TtohTVVV1d5I/Af4F2A54X1Vd2nMxU57GWiCsb/NY3+ax\nvs2zTdWXqpq5lyRJzZZyqkqStEAYHJKkXra64EiyXZILk3xukmnHJFmf5Jvt8YIB6lub5Ftt/asn\nmZ4kJ7Vbq1ycZL8FVt9BSTaMbMPXj7m+xUnOSPKdJJclOXDC9MG23yxqG3rbPWZk3d9McluSP53Q\nZ8jtN5v6BtuGSV6Z5NIklyQ5LckOE6Y/IMnpbdudn2TZuGqbZX1zt/+rqq3qAbwK+CjwuUmmHQP8\n48D1rQX2mGb6IcAX6b67cgBw/gKr76DJtu0Y6zsVeEEb3h5YvFC23yxqG3TbTahlO+D7wMMXyvab\nZX2DbEO6LyFfBTywjX8cOGZCn5cA/9SGjwROX2D1zdn+b6s64kiyF3Ao8N6ha9kMK4EPVuc8YHGS\nhw1d1EKQZGfg6cApAFV1V1XdOqHbINtvlrUtJCuA71bV1RPaF8rf31T1DWkR8MAki4Adue93yVbS\nvXkAOANYkWSyLy8PVd+c2aqCA/h74DXAz6fp8zvtEPyMJHtP02++FPDlJBe026RMNNntVZaOpbLO\nTPUBHJjkoiRfTPIrY6ztkcB64P3tdOR7k+w0oc9Q2282tcFw226iI4HTJmkf+u9vo6nqgwG2YVV9\nD3gbcA1wPbChqr48odsvtl1V3Q1sAHZfQPXBHO3/tprgSPIc4MaqumCabp8FllXVE4CvcM+7g3F6\nWlXtR3en35cmefqE6TPeXmWezVTfN+hOHzwReCfw6THWtgjYD3h3VT0ZuAOYeIv9obbfbGobctv9\nQrov0T4P+MRkkydpG+tn9meob5BtmGRXuiOKRwB7Ajsl+YOJ3SaZdSzbbpb1zdn+b6sJDuBpwPOS\nrKW7e+4zknx4tENV/aCqftJG3wP82nhLhKq6rj3fCHyK7s6/owa9vcpM9VXVbVV1exv+AnD/JHuM\nqbx1wLqqOr+Nn0G3s57YZ4jtN2NtA2+7UQcD36iqGyaZthBu7zNlfQNuw2cCV1XV+qr6KfBJ4KkT\n+vxi27XTRbsAN4+htlnVN5f7v60mOKrqtVW1V1UtozvMPaeq7pW4E87VPg+4bIwlkmSnJA/eOAw8\nG7hkQrczgaPbp1sOoDvkvH6h1Jfklzaet02yP93f0A/GUV9VfR+4NsljWtMK4NsTug2y/WZT25Db\nboLnM/VpoMH+/kZMWd+A2/Aa4IAkO7b1r+C++48zgVVt+HC6fdC4jtZmrG8u939bxC1HNkeSNwOr\nq+pM4OVJngfcTfdO4Jgxl/NQ4FPt734R8NGq+lKSPwaoqn8CvkD3yZY1wJ3AsQusvsOBFye5G/gR\ncOQY/3MAvAz4SDudcSVw7ALafjPVNvS2I8mOwLOAF420LZTtN5v6BtmGVXV+kjPoTpXdDVwInDxh\n/3IK8KEka+j2L0fOd10965uz/Z+3HJEk9bLVnKqSJI2HwSFJ6sXgkCT1YnBIknoxOCRJvRgc0gCS\nLEvy+5sw3+IkLxkZ37N9DFMaG4NDGsYyYNLgaN86nspiuruwAt03/avq8LktTZqewaFtSpKj203e\nLkryoSQPT3J2azs7yS+3fh9I97sU/5nkyiSHjyzjNel+s+SiJCe2tn2SfCndzSH/PcljZ1jOicBv\npPtdhFem+62ETyT5LN1NJh/U6vlGW9fKkfn2afO9tR25XNLWtUOS97f+Fyb5H639mCSfbPVdkeTv\nxrKxtfWai3uz+/CxJTyAXwEup/3eCLAb3Y3fVrXxPwI+3YY/QHeTvfsBjwfWtPaDgf8Edty4jPZ8\nNrBvG34K3e0mplvOQYz8rgTdt3jXjSxvEbBzG96D7pvcoTtSuWRkvl+MA38GvL8NP5buNhQ7tGVf\nSXfvpB2Aq4G9h/738LHlPrb6W45II54BnFFVNwFU1c3pfqXvf7bpHwJG341/uqp+Dnw7yUNb2zPp\nds53jizjQXQ3lPtE7vn5hQfMsJzJnFVVG2+KF+At6e5O/HO6W3ZPNy/Ar9PdMZaq+k6Sq4FHt2ln\nV9UGgCTfBh7OvW+fLs2awaFtSZj5Ntej038yMpyR54nLuB9wa1U9aYplTracydwxMnwUsAT4tar6\nabq7Pu8w6VyzW/ZoDT/D//vaDF7j0LbkbOB3k+wOkGQ3utNOG29GdxTwtRmW8WXgj9rN+EiyW1Xd\nBlyV5IjWliRPnGE5PwQePM30Xeh+X+an7VrFw2cx37+110CSRwO/THdqTppTBoe2GVV1KfDXwL8m\nuQh4O/ByurvYXgz8IfCKGZbxJbrbZ69O8k3g1W3SUcBxbbmX0v2oznQuBu5uF9hfOcn0jwDLk6xu\ny/5OW/8PgP9IckmSt06Y513Adkm+BZxO95vTP0GaY94dV5LUi0cckqReDA5JUi8GhySpF4NDktSL\nwSFJ6sXgkCT1YnBIknr5/4o4RF0/jCAaAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGd9JREFUeJzt3Xm0JWV97vHvI+2s0AKtkYbYiu2U\nxAH7glOMESdEbdYNKIaEBlGMczTGtN4bp6woxizHXMlFUdsZJQ4ozqBRY+DaoCCILFpsoAWhEWgE\nnNDf/aPeA9vNmao5Z+9zON/PWnudqrfeqvrtotnPrrf2rp2qQpKk2brVuAuQJC0uBockqReDQ5LU\ni8EhSerF4JAk9WJwSJJ6MTi0oCV5VZL3zNG2Ksm9b8b6X0iybi5qWaySPDLJeUmuSXLACPd7WJJv\njWp/mp7BoVlJsjnJr5PsOtT+vfaCvGo+9ltVb6iqZ8/Htgcl+XqSZw+1PSbJloFa9quqDbPY1s0K\nqAXu9cC/VdWdqurTwwvbv5NLk9xxoO3ZSb4+3UaT3DXJR5NcnGRbkv9Kss/cl6+5YHCojx8Dz5yY\nSfInwO23d2NJls1FUUvJAjhm9wDOnqHPMuAlPbd7J+A7wEOBnYENwIlJ7tS7Qs07g0N9fBA4dGB+\nHfCBwQ5J9k/y3SRXJ7koyWsHlq1q78aPSHIhcHJrPzTJBUl+luQf27vWx7Vlr03yoaH11yW5MMnl\nSf7XwPb3TvLfSa5KckmSf0tym7l68oNnJUnuneQ/27vjy5Mc19q/0bqf0YZzntHan5NkU5IrkpyQ\nZLeB7T4hybltW+9q253Yz2Ht3fdbk1wBvDbJnklObsfr8iQfTrJ8YHubk/x9kjOTXJvk2CR3a0Nt\nP0/y1SR3meZ5Tlprkh8B9wI+257bbafYxJuBlw/WNJOqOr+q3lJVl1TVb6vqGOA2wH2nqPHNSb6V\nZKfZ7kNzx+BQH6cAOya5f5IdgGcAHxrqcy1duCwH9geeN8lY+J8B9weemOQBwLuAQ4C7AzsBK2eo\n41F0Lyj7Aq9Ocv/W/lvgpcCuwMPb8uf3fZKz9E/Al4G7ALsD7wSoqke35Q9qwznHJXks8Ebg6XTP\n8QLgYwBt6O944JXALsC5wCOG9rUPcD5wV+CfgbTt7UZ3HPcAXju0zl8AjwfuAzwV+ALwKrpjcyvg\nxZM9qelqrao9gQuBp7bn9qspjs1G4OvAy6dYPqMkD6YLjk1D7bdK8m7ggcATqmrb9u5D28/gUF8T\nZx2PB34I/GRwYVV9vaq+X1W/q6ozgY/SBcWg11bVtVX1C+BA4LNV9a2q+jXwamCmG6i9rqp+UVVn\nAGcAD2r7Pq2qTqmq66tqM/B/J9n3dN7RzlauSnIV8Llp+v6Gbthmt6r6ZVVNd+H2EOC9VXV6e7F9\nJfDwdl3oycDZVfXJqroeeAfw06H1L66qd7bn9Yuq2lRVX6mqX1XVVuAtkzzPd1bVpVX1E+CbwKlV\n9d22/08BD9mOWvt4NfCiJCt6rkeSHen+nb1uKBhuTffvaWe68Lqu77Y1NwwO9fVB4C+BwxgapgJI\nsk+SryXZmmQb8Dd073IHXTQwvdvgfHsx+NkMNQy+sF5HNz5Okvsk+VySnya5GnjDJPuezouravnE\nA3jKNH1fQffO//8lOTvJs6bpuxvdO3cAquoauue4kps+/wK2DK0/eLwmLiR/LMlP2vP8EDd9npcO\nTP9ikvmprh1MV+usVdVZdMG7vs96SW4PfBY4pareOLT43sBaukD5dZ/tam4ZHOqlqi6gu0j+ZOCT\nk3T5CHACsEdV7QT8O90L7O9tZmD6ErqhHuCGF45dtrO8o+nOglZX1Y50QzPD+54TVfXTqnpOVe0G\nPBd4V6b+JNXFdGcnAKT7xNEudGdrw88/g/MTuxuaf2Nre2B7nn/F3D3P6Wrt6zXAc5hl6LRrJp9u\n+3ruJF3OAQ4HvpBk0msfGg2DQ9vjCOCxVXXtJMvuDFxRVb9Msjfd2cl0jgeemuQR7UL269j+F8E7\nA1cD1yS5H/C87dzOjJIclGTiBf5Kuhfy37b5S+kuIk/4CHB4kge3F8c30A0dbQZOBP4kyQHpPjH1\nAuAPZtj9nYFrgKuSrAT+fi6e0yxq7aWqNgHHMcX1lEFJbk33b+EXwKFV9bsptvlRujcEX02yZ9+a\nNDcMDvVWVT+qqo1TLH4+8PokP6cb5/74DNs6G3gR3QXYS4CfA5cBU114nc7L6YLq58C76V605sv/\nAE5Ncg3dGdZLqurHbdlrgQ3tWsnTq+ok4B+B/6B7jnsCBwNU1eXAQcC/0A0JPYDu4vJ0z/91wF7A\nNrrgmezMb7tMV+t2ej1wxxl7dR8IeArwBLpAvKY9/nSSGje07Z68HddeNAfiDzlpIUn3uf2r6Iab\nfjxT/1uaJLeiu8ZxSFV9bdz1SJPxjENjl+SpSe7QxtP/Ffg+sHm8VY1OkicmWd6Ghiauy5wy5rKk\nKRkcWgjW0l2UvRhYDRxcS+tU+OHAj4DL6b5zcUD7qPItSpI/HRiC+r3HuGtTPw5VSZJ68YxDktTL\nuG+YNi923XXXWrVq1bjLkKRF5bTTTru8qmb8tv8tMjhWrVrFxo1TfVpUkjSZJBfM3MuhKklSTwaH\nJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSL7fIb45LM1m1/sSx7XvzUfuPbd/S\nXPCMQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIv8xYcSd6b5LIkZw207ZzkK0nOa3/v\n0tqT5B1JNiU5M8leA+usa/3PS7JuvuqVJM3OfJ5xvB940lDbeuCkqloNnNTmAfYDVrfHkcDR0AUN\n8BpgH2Bv4DUTYSNJGo95C46q+gZwxVDzWmBDm94AHDDQ/oHqnAIsT3J34InAV6rqiqq6EvgKNw0j\nSdIIjfoax92q6hKA9veurX0lcNFAvy2tbar2m0hyZJKNSTZu3bp1zguXJHUWysXxTNJW07TftLHq\nmKpaU1VrVqxYMafFSZJuNOrguLQNQdH+XtbatwB7DPTbHbh4mnZJ0piM+u64JwDrgKPa388MtL8w\nycfoLoRvq6pLknwJeMPABfEnAK8ccc3SnBrXnXm9K6/myrwFR5KPAo8Bdk2yhe7TUUcBH09yBHAh\ncFDr/nngycAm4DrgcICquiLJPwHfaf1eX1XDF9wlSSM0b8FRVc+cYtG+k/Qt4AVTbOe9wHvnsDRJ\n0s2wUC6OS5IWCYNDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ\n6mXUd8eVNCbjuisveGfeWxrPOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknox\nOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqRexhIcSV6a5Owk\nZyX5aJLbJblnklOTnJfkuCS3aX1v2+Y3teWrxlGzJKkz8uBIshJ4MbCmqv4Y2AE4GHgT8NaqWg1c\nCRzRVjkCuLKq7g28tfWTJI3JuIaqlgG3T7IMuANwCfBY4Pi2fANwQJte2+Zpy/dNkhHWKkkaMPLg\nqKqfAP8KXEgXGNuA04Crqur61m0LsLJNrwQuaute3/rvMrzdJEcm2Zhk49atW+f3SUjSEjaOoaq7\n0J1F3BPYDbgjsN8kXWtilWmW3dhQdUxVramqNStWrJirciVJQ8YxVPU44MdVtbWqfgN8EngEsLwN\nXQHsDlzcprcAewC05TsBV4y2ZEnShHEEx4XAw5LcoV2r2Bf4AfA14MDWZx3wmTZ9QpunLT+5qm5y\nxiFJGo1xXOM4le4i9+nA91sNxwD/ALwsySa6axjHtlWOBXZp7S8D1o+6ZknSjZbN3GXuVdVrgNcM\nNZ8P7D1J318CB42iLknSzPzmuCSpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8G\nhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1\nYnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvSwbx06TLAfe\nA/wxUMCzgHOB44BVwGbg6VV1ZZIAbweeDFwHHFZVp4+hbM2DVetPHHcJknoa1xnH24EvVtX9gAcB\n5wDrgZOqajVwUpsH2A9Y3R5HAkePvlxJ0oSRB0eSHYFHA8cCVNWvq+oqYC2woXXbABzQptcCH6jO\nKcDyJHcfcdmSpGYcZxz3ArYC70vy3STvSXJH4G5VdQlA+3vX1n8lcNHA+lta2+9JcmSSjUk2bt26\ndX6fgSQtYeMIjmXAXsDRVfUQ4FpuHJaaTCZpq5s0VB1TVWuqas2KFSvmplJJ0k2MIzi2AFuq6tQ2\nfzxdkFw6MQTV/l420H+PgfV3By4eUa2SpCEjD46q+ilwUZL7tqZ9gR8AJwDrWts64DNt+gTg0HQe\nBmybGNKSJI3eWD6OC7wI+HCS2wDnA4fThdjHkxwBXAgc1Pp+nu6juJvoPo57+OjLlSRNGEtwVNX3\ngDWTLNp3kr4FvGDei5IkzYrfHJck9TKr4Ejyvwembzt/5UiSFrppgyPJK5I8HDhwoPm/57ckSdJC\nNtM1jnPpLlLfK8k36W4NskuS+1bVufNenSRpwZlpqOpK4FV0n2h6DPCO1r4+ybfnsS5J0gI10xnH\nk4DXAHsCbwHOAK6tKj8SK0lL1LTBUVWvAkhyBvAh4CHAiiTfAq6sqqfOf4mSFrtx3T5/81H7j2W/\nt3Sz/R7Hl6rqO8B3kjyvqh6VZNf5LEyStDDN6uO4VfWKgdnDWtvl81GQJGlh6/0FwKo6Yz4KkSQt\nDn5zXJLUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJw\nSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJwSJJ6GVtwJNkhyXeTfK7N3zPJ\nqUnOS3Jcktu09tu2+U1t+apx1SxJGu8Zx0uAcwbm3wS8tapWA1cCR7T2I4Arq+rewFtbP0nSmIwl\nOJLsDuwPvKfNB3gscHzrsgE4oE2vbfO05fu2/pKkMRjXGcfbgFcAv2vzuwBXVdX1bX4LsLJNrwQu\nAmjLt7X+kqQxGHlwJHkKcFlVnTbYPEnXmsWywe0emWRjko1bt26dg0olSZMZxxnHI4GnJdkMfIxu\niOptwPIky1qf3YGL2/QWYA+Atnwn4IrhjVbVMVW1pqrWrFixYn6fgSQtYSMPjqp6ZVXtXlWrgIOB\nk6vqEOBrwIGt2zrgM236hDZPW35yVd3kjEOSNBoL6Xsc/wC8LMkmumsYx7b2Y4FdWvvLgPVjqk+S\nBCybucv8qaqvA19v0+cDe0/S55fAQSMtbIlZtf7EcZcgaRFZSGcckqRFwOCQJPVicEiSejE4JEm9\nGBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS\n1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSelk2\n7gIkab6sWn/iWPa7+aj9x7LfUfGMQ5LUi8EhSepl5MGRZI8kX0tyTpKzk7ykte+c5CtJzmt/79La\nk+QdSTYlOTPJXqOuWZJ0o3GccVwP/F1V3R94GPCCJA8A1gMnVdVq4KQ2D7AfsLo9jgSOHn3JkqQJ\nIw+Oqrqkqk5v0z8HzgFWAmuBDa3bBuCANr0W+EB1TgGWJ7n7iMuWJDVjvcaRZBXwEOBU4G5VdQl0\n4QLctXVbCVw0sNqW1iZJGoOxBUeSOwH/AfxtVV09XddJ2mqS7R2ZZGOSjVu3bp2rMiVJQ8YSHElu\nTRcaH66qT7bmSyeGoNrfy1r7FmCPgdV3By4e3mZVHVNVa6pqzYoVK+aveEla4sbxqaoAxwLnVNVb\nBhadAKxr0+uAzwy0H9o+XfUwYNvEkJYkafTG8c3xRwJ/DXw/yfda26uAo4CPJzkCuBA4qC37PPBk\nYBNwHXD4aMuVJA0aeXBU1beY/LoFwL6T9C/gBfNalCRp1rxX1QIyrvvqSFIf3nJEktSLwSFJ6sXg\nkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqRe\nDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknpZNu4CJOmWZtX6E8e2781H7T/v+/CMQ5LUi8Eh\nSerF4JAk9eI1jkmMc3xSkhY6zzgkSb0YHJKkXgwOSVIvBockqReDQ5LUy6IJjiRPSnJukk1J1o+7\nHklaqhZFcCTZAfg/wH7AA4BnJnnAeKuSpKVpUQQHsDewqarOr6pfAx8D1o65JklakhbLFwBXAhcN\nzG8B9hnskORI4Mg2e02Sc2/mPncFLr+Z25hvi6FGWBx1LoYawTrn0mKoEXrWmTfdrH3dYzadFktw\nZJK2+r2ZqmOAY+Zsh8nGqlozV9ubD4uhRlgcdS6GGsE659JiqBEWZp2LZahqC7DHwPzuwMVjqkWS\nlrTFEhzfAVYnuWeS2wAHAyeMuSZJWpIWxVBVVV2f5IXAl4AdgPdW1dnzvNs5G/aaR4uhRlgcdS6G\nGsE659JiqBEWYJ2pqpl7SZLULJahKknSAmFwSJJ6WfLBkWSHJN9N8rlJlt02yXHtNienJlk1+gpv\nqGW6Og9LsjXJ99rj2WOqcXOS77caNk6yPEne0Y7nmUn2WoA1PibJtoFj+epR19jqWJ7k+CQ/THJO\nkocPLV8Ix3KmGsd+LJPcd2D/30tydZK/HeqzEI7lbOoc+/GcsCgujs+zlwDnADtOsuwI4MqquneS\ng4E3Ac8YZXEDpqsT4LiqeuEI65nKn1fVVF9W2g9Y3R77AEcz9EXOEZmuRoBvVtVTRlbN5N4OfLGq\nDmyfJLzD0PKFcCxnqhHGfCyr6lzgwXDDrYt+AnxqqNvYj+Us64SF8W9zaZ9xJNkd2B94zxRd1gIb\n2vTxwL5JJvsy4ryaRZ2LxVrgA9U5BVie5O7jLmqhSbIj8GjgWICq+nVVXTXUbazHcpY1LjT7Aj+q\nqguG2hfav8up6lwwlnRwAG8DXgH8borlN9zqpKquB7YBu4ymtN8zU50Af9FOs49Pssc0/eZTAV9O\nclq7BcywyW4ds3Ikld1ophoBHp7kjCRfSPJHoyyuuRewFXhfG558T5I7DvUZ97GcTY0w/mM56GDg\no5O0j/tYDpuqTlggx3PJBkeSpwCXVdVp03WbpG2kn1+eZZ2fBVZV1QOBr3LjWdKoPbKq9qI79X9B\nkkcPLR/78WTmGk8H7lFVDwLeCXx6xPVBN4S8F3B0VT0EuBYY/imBcR/L2dS4EI4lAG0o7WnAJyZb\nPEnbWL6nMEOdC+Z4LtngAB4JPC3JZrq77T42yYeG+txwq5Mky4CdgCtGWSSzqLOqflZVv2qz7wYe\nOtoSb6jj4vb3Mrrx2b2Huoz91jEz1VhVV1fVNW3688Ctk+w6yhrpjtOWqjq1zR9P9yI93Gecx3LG\nGhfIsZywH3B6VV06ybJxH8tBU9a5kI7nkg2OqnplVe1eVavoTg1Prqq/Gup2ArCuTR/Y+oz0nchs\n6hwaj30a3UX0kUpyxyR3npgGngCcNdTtBODQ9imWhwHbquqShVRjkj+YuI6VZG+6/0d+NqoaAarq\np8BFSe7bmvYFfjDUbazHcjY1LoRjOeCZTD38M9ZjOWTKOhfS8fRTVUOSvB7YWFUn0F34+2CSTXRn\nGgePtbgBQ3W+OMnTgOvp6jxsDCXdDfhU+3e9DPhIVX0xyd8AVNW/A58HngxsAq4DDl+ANR4IPC/J\n9cAvgINH/WaheRHw4TZ0cT5w+AI7lrOpcUEcyyR3AB4PPHegbaEdy9nUuSCOJ3jLEUlST0t2qEqS\ntH0MDklSLwaHJKkXg0OS1IvBIUnqxeCQxiDJqiR/uR3rLU/y/IH53ZIcP7fVSdMzOKTxWAVMGhzt\nLgVTWQ7cEBxVdXFVHTi3pUnTMzi0pCQ5tN0M8owkH0xyjyQntbaTkvxh6/f+dL/R8O0k5yc5cGAb\nr0j3mx5nJDmqte2Z5Ivt5onfTHK/GbZzFPCn6X5X4aXpflPlE0k+S3cTxju1ek5v+1o7sN6ebb03\ntzOXs9q+bpfkfa3/d5P8eWs/LMknW33nJfmXkRxs3XJVlQ8fS+IB/BFwLrBrm9+Z7gaR69r8s4BP\nt+n3091o7lbAA4BNrX0/4NvAHSa20f6eBKxu0/vQ3Rpmuu08BvjcQG2H0d0zaWJ7y4Ad2/SudN9q\nDt2ZylkD690wD/wd8L42fT/gQuB2bdvn091r7XbABcAe4/7v4WPxPrzliJaSxwLHV/sRp6q6It2v\n1v3PtvyDwOC78U9X1e+AHyS5W2t7HN2L83UD27gT8AjgE7nx51puO8N2JvOVqpq4iWaAN6S7e+/v\n6G7zPd26AI+iu2sqVfXDJBcA92nLTqqqbQBJfgDcg9+/lbg0awaHlpIw8+2yB5f/amA6A3+Ht3Er\n4KqqevAU25xsO5O5dmD6EGAF8NCq+k26uyPfbpp1Z9r2YA2/xf/3dTN4jUNLyUnA05PsApBkZ7ph\np4mbVx4CfGuGbXwZeFa7IR1Jdq6qq4EfJzmotSXJg2bYzs+BO0+zfCe632H5TbtWcY9ZrPeN9hxI\nch/gD+mG5qQ5ZXBoyaiqs4F/Bv4zyRnAW4AX093V9Uzgr+l+2326bXyR7jbcG5N8D3h5W3QIcETb\n7tl0P0c6nTOB69sF9pdOsvzDwJokG9u2f9j2/zPgv5KcleTNQ+u8C9ghyfeB44DD6sbfaZHmjHfH\nlST14hmHJKkXg0OS1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF7+P082LwRAop0FAAAAAElFTkSu\nQmCC\n", "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -147,18 +132,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 3. do Map3k = 5.678" + "## 3. do Map3k = 0.678" ] }, { "cell_type": "code", - "execution_count": 27, - "metadata": { - "collapsed": true - }, + "execution_count": 12, + "metadata": {}, "outputs": [], "source": [ - "do_map3k = pyro.do(mapk_receiver, data={'map3k': 5.678})" + "do_map3k = pyro.do(mapk_companion, data={'map3k': torch.tensor(50000.678)})" ] }, { @@ -170,31 +153,68 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "# Deriving mean and std by sampling from marginals\n", + "noise_posterior_mean = {\n", + " n: torch.stack([noise_marginals[n]() for _ in range(5000)]).mean()\n", + " for n in noise_vars\n", + "}\n", + "noise_posterior_std = {\n", + " n: torch.stack([noise_marginals[n]() for _ in range(5000)]).std()\n", + " for n in noise_vars\n", + "}\n", + "\n", + "noise_posteriors = {\n", + " n: Normal(noise_posterior_mean[n], noise_posterior_std[n])\n", + " for n in noise_vars\n", + "}\n", + "\n", + "mapk_do_dist = infer_dist(do_map3k, noise_posteriors)\n", + "\n", + "mapk_marginal = EmpiricalMarginal(mapk_do_dist, sites='map2k')" + ] + }, + { + "cell_type": "code", + "execution_count": 14, "metadata": {}, "outputs": [ { - "ename": "RuntimeError", - "evalue": "One of the differentiated Tensors appears to not have been used in the graph. Set allow_unused=True if this is the desired behavior.", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mmapk_do_dist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minfer_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdo_map3k\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnoise_marginals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mmapk_marginal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mEmpiricalMarginal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmapk_do_dist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msites\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'mapk'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Desktop/Git/causal_demon/mapk_demon/inference.py\u001b[0m in \u001b[0;36minfer_dist\u001b[0;34m(prog, n_dist, type)\u001b[0m\n\u001b[1;32m 19\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mtype\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'mcmc'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 20\u001b[0m \u001b[0mhmc_kernel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mHMC\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mprog\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstep_size\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.9\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m \u001b[0mposterior\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMCMC\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhmc_kernel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum_samples\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwarmup_steps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m50\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_dist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 22\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mposterior\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/infer/abstract_infer.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 81\u001b[0m \"\"\"\n\u001b[1;32m 82\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_init\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 83\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mtr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlogit\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpoutine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mblock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_traces\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 84\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexec_traces\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 85\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_weights\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlogit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/infer/mcmc/mcmc.py\u001b[0m in \u001b[0;36m_traces\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 35\u001b[0m \u001b[0mlogging_interval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mceil\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_samples\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 36\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnum_samples\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 37\u001b[0;31m \u001b[0mtrace\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkernel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 38\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mlogging_interval\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 39\u001b[0m \u001b[0mstage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"WARMUP\"\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m<=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwarmup_steps\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m\"SAMPLE\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/infer/mcmc/hmc.py\u001b[0m in \u001b[0;36msample\u001b[0;34m(self, trace)\u001b[0m\n\u001b[1;32m 232\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_potential_energy\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 233\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep_size\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 234\u001b[0;31m self.num_steps)\n\u001b[0m\u001b[1;32m 235\u001b[0m \u001b[0;31m# apply Metropolis correction.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 236\u001b[0m \u001b[0menergy_proposal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_energy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz_new\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr_new\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/ops/integrator.py\u001b[0m in \u001b[0;36mvelocity_verlet\u001b[0;34m(z, r, potential_fn, step_size, num_steps)\u001b[0m\n\u001b[1;32m 22\u001b[0m \u001b[0mz_next\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mz\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0mr_next\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 24\u001b[0;31m \u001b[0mgrads\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpotential_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz_next\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 25\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 26\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum_steps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/pyro/ops/integrator.py\u001b[0m in \u001b[0;36m_grad\u001b[0;34m(potential_fn, z)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[0mnode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrequires_grad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 65\u001b[0m \u001b[0mpotential_energy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpotential_fn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mz\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 66\u001b[0;31m \u001b[0mgrads\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgrad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpotential_energy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mz_nodes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 67\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mnode\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mz_nodes\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0mnode\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrequires_grad\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/anaconda3/lib/python3.6/site-packages/torch/autograd/__init__.py\u001b[0m in \u001b[0;36mgrad\u001b[0;34m(outputs, inputs, grad_outputs, retain_graph, create_graph, only_inputs, allow_unused)\u001b[0m\n\u001b[1;32m 143\u001b[0m return Variable._execution_engine.run_backward(\n\u001b[1;32m 144\u001b[0m \u001b[0moutputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgrad_outputs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mretain_graph\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcreate_graph\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 145\u001b[0;31m inputs, allow_unused)\n\u001b[0m\u001b[1;32m 146\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 147\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mRuntimeError\u001b[0m: One of the differentiated Tensors appears to not have been used in the graph. Set allow_unused=True if this is the desired behavior." - ] + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAGkJJREFUeJzt3Xu4HXV97/H3RygoXgiX6JGEYxDj\nrVarzcFrPR6xCqKG5xxpsbQEpKUXWz3aVoO2Be2jpbWPtthqDxUkKsULtRrFWxq01LZQgghykRIB\nSQAhlItcVIx+zx/z27Lc2dk7A3uvSdjv1/PsZ8385jcz31mE9Vnzm7VmpaqQJGlbPWjoAiRJOxaD\nQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHNouJHlzkvfP0rYqyePux/qfS7JiNmrZUSV5bpIrk9yZ\n5NCh69maJEcl+crQdcw3Bod+QpJrktyTZO9J7V9rL8hL5mK/VfWOqvq1udj2qCRfTvJrk9pekGTj\nSC0HV9WqbdjW/Qqo7dzbgL+uqodV1SeHKiLJI5OckeT6JLcn+dckzxyqHnUMDk3lauBVEzNJfgZ4\nyH3dWJKdZ6Oo+WQ7eM4eA1w6cA0ADwPOB34O2BNYBZyV5GGDVjXPGRyayoeAI0fmVwAfHO2Q5JAk\nFyb5TpINSU4YWbakvRs/Jsm1wNmt/cgk30ryX0n+qJ3dvKgtOyHJhyetvyLJtUluTvKWke0fkOTf\nk9yW5IYkf51kl9k6+NGzkiSPS/LP7d3uzUk+2trPad0vasM5v9Tafz3J+iS3JFmdZJ+R7b44yRVt\nW+9t253Yz1Ht3fS7k9wCnJBk/yRnt+fr5iSnJ1kwsr1rkvxBkouT3JXklCSPakNtdyT5pyR7THOc\nU9aa5JvAY4FPt2PbdYp1e+07yceTfLsd+zlJfnpk2WlJ/jbJmrbuPyd5DEBVXVVV76qqG6rqh1V1\nMrAL8IStHNM7k3wlye7T/1fW/WFwaCrnAo9I8qQkOwG/BHx4Up+76MJlAXAI8FvZciz8fwJPAl6S\n5MnAe4EjgEcDuwOLZqjjeXQvEAcCf5zkSa39h8Drgb2BZ7flv933ILfRnwBfBPYAFgPvAaiq57fl\nT2vDOR9N8kLgT4FfpDvGbwEfAUg39HcmcBywF3AF8JxJ+3omcBXwSODtQNr29qF7HvcFTpi0zv8B\nfgF4PPBy4HPAm+memwcBr53qoKartar2B64FXt6O7ftbeW767PtzwNJ2bF8FTp+0rSPonuu9ga9N\nsXyi7p+lC471k9oflOTvgKcCL66q27dSs2aBwaGtmTjr+AXgG8B1owur6stV9fWq+lFVXQycQRcU\no06oqruq6rvAK4FPV9VXquoe4I+BmW6U9taq+m5VXQRcBDyt7fuCqjq3qjZX1TXA/5ti39M5qZ2t\n3JbkNuAz0/T9Ad2wzT5V9b2qmu5C7BHAqVX11fZiexzw7HTXhV4KXFpVn6iqzcBJwLcnrX99Vb2n\nHdd3q2p9Va2pqu9X1SbgXVMc53uq6saqug74F+C8qrqw7f8fgaffh1q31Tbvu6pOrao72rITgKdN\nOis4q6rOacvf0mrZd3RnSR5B9+/yrZOC4afo/v3tSRd2d/c4Bt0HBoe25kPALwNHMWmYCiDJM5N8\nKcmmJLcDv0n3bnHUhpHpfUbn2//c/zVDDaMvrHfTjXeT5PFJPtOGPr4DvGOKfU/ntVW1YOIPeNk0\nfd9I987/P5JcmuTV0/Tdh+6dOwBVdSfdMS5iy+MvYOOk9Uefr4kLwx9Jcl07zg+z5XHeODL93Snm\nt3YtYLpat9U27TvJTklOTPLNdhzXtD6jxzL63NwJ3NJqpG3jIcCngXOr6k8n1fE4YDldoNzTo37d\nRwaHplRV36K7SP5S4BNTdPl7YDWwb1XtDvwt3QvsT2xmZPoGuqEe4McvBHvdx/LeR3cWtLSqHkE3\nPDJ537Oiqr5dVb9eVfsAvwG8N1v/JNX1dGcnACR5KN0xXseWx5/R+YndTZr/09b21Hacv8LsHed0\ntc62X6Z7YX8R3RDlkondjvT58dlFugvfe7YaaddYPtlq+40ptn85cDTwuSRTXvvQ7DI4NJ1jgBdW\n1V1TLHs4cEtVfS/JAXQvDtM5E3h5kue0C9lv5b6/CD4c+A5wZ5InAr91H7czoySHJZl4gb+V7oX8\nh23+RrqLyBP+Hjg6yc+2F7t30A3fXAOcBfxMkkPTfWLqNcB/m2H3DwfuBG5Lsgj4g9k4pm2odbY9\nHPg+3RnNbm1fk700yfPav40/abVsSPJTdP92vgscWVU/mmoHVXUG3RuIf0qy/xwcg0YYHNqqqvpm\nVa3byuLfBt6W5A666xUfm2FblwK/S3cB9gbgDuAmuheUvn6fLqjuAP4O+Oh92Ma2+h/AeUnupDvD\nel1VXd2WnQCsatdKfrGq1gJ/BPwD3THuDxwOUFU3A4cBf073AvpkYB3TH/9bgWcAt9MFz1RnfvfJ\ndLXOgQ/SDYtdB1xG9+GLyf4eOJ5uiOrn6K7BQPcBgpcBL6YL0Dvb389P3kD77s3bgLN7XqtRT/GH\nnDSENhxxG91w09Uz9X+gSfIgumscR1TVl4auZ0hJTgM2VtUfDl2Lto1nHBqbJC9PslsbT/8L4Ovc\ne6H0AS/JS5IsaENDE9dlpnr3LW3XDA6N03K6C57X032m//CaX6e8zwa+CdxM972HQ9tHlaUdikNV\nkqRePOOQJPUy9I3U5sTee+9dS5YsGboMSdqhXHDBBTdX1cKZ+j0gg2PJkiWsW7e1T5FKkqaS5Fsz\n93KoSpLUk8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReDQ5LUywPym+PSTJasPGuw\nfV9z4iGD7VuaDZ5xSJJ6MTgkSb0YHJKkXgwOSVIvcxYcSU5NclOSS0ba3pnkG0kuTvKPSRaMLDsu\nyfokVyR5yUj7Qa1tfZKVc1WvJGnbzOUZx2nAQZPa1gBPqaqnAv8JHAeQ5MnA4cBPt3Xem2SnJDsB\nfwMcDDwZeFXrK0kayJwFR1WdA9wyqe2LVbW5zZ4LLG7Ty4GPVNX3q+pqYD1wQPtbX1VXVdU9wEda\nX0nSQIa8xvFq4HNtehGwYWTZxta2tfYtJDk2ybok6zZt2jQH5UqSYKAvACZ5C7AZOH2iaYpuxdTB\nVlNts6pOBk4GWLZs2ZR9pO3BUF8+9IuHmi1jD44kK4CXAQdW1cQL/EZg35Fui4Hr2/TW2iVJAxjr\nUFWSg4A3Aa+oqrtHFq0GDk+ya5L9gKXAfwDnA0uT7JdkF7oL6KvHWbMk6SfN2RlHkjOAFwB7J9kI\nHE/3KapdgTVJAM6tqt+sqkuTfAy4jG4I6zVV9cO2nd8BvgDsBJxaVZfOVc2SpJnNWXBU1aumaD5l\nmv5vB94+Rftngc/OYmmSpPvBb45LknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1\nYnBIknoxOCRJvRgckqReDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5J\nUi8GhySpF4NDktSLwSFJ6sXgkCT1MmfBkeTUJDcluWSkbc8ka5Jc2R73aO1JclKS9UkuTvKMkXVW\ntP5XJlkxV/VKkrbNXJ5xnAYcNKltJbC2qpYCa9s8wMHA0vZ3LPA+6IIGOB54JnAAcPxE2EiShjFn\nwVFV5wC3TGpeDqxq06uAQ0faP1idc4EFSR4NvARYU1W3VNWtwBq2DCNJ0hiN+xrHo6rqBoD2+MjW\nvgjYMNJvY2vbWvsWkhybZF2SdZs2bZr1wiVJne3l4nimaKtp2rdsrDq5qpZV1bKFCxfOanGSpHuN\nOzhubENQtMebWvtGYN+RfouB66dplyQNZNzBsRqY+GTUCuBTI+1Htk9XPQu4vQ1lfQF4cZI92kXx\nF7c2SdJAdp6rDSc5A3gBsHeSjXSfjjoR+FiSY4BrgcNa988CLwXWA3cDRwNU1S1J/gQ4v/V7W1VN\nvuAuSRqjOQuOqnrVVhYdOEXfAl6zle2cCpw6i6VJku6H7eXiuCRpB2FwSJJ6MTgkSb0YHJKkXgwO\nSVIvc/apKmlbLFl51tAlSOrJMw5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqRe\nDA5JUi8GhySpF4NDktSLwSFJ6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpl0GCI8nrk1ya\n5JIkZyR5cJL9kpyX5MokH02yS+u7a5tf35YvGaJmSVJn7MGRZBHwWmBZVT0F2Ak4HPgz4N1VtRS4\nFTimrXIMcGtVPQ54d+snSRrIUENVOwMPSbIzsBtwA/BC4My2fBVwaJte3uZpyw9MkjHWKkkaMfbg\nqKrrgL8ArqULjNuBC4Dbqmpz67YRWNSmFwEb2rqbW/+9Jm83ybFJ1iVZt2nTprk9CEmax4YYqtqD\n7ixiP2Af4KHAwVN0rYlVpll2b0PVyVW1rKqWLVy4cLbKlSRNMsRQ1YuAq6tqU1X9APgE8BxgQRu6\nAlgMXN+mNwL7ArTluwO3jLdkSdKEIYLjWuBZSXZr1yoOBC4DvgS8svVZAXyqTa9u87TlZ1fVFmcc\nkqTxGOIax3l0F7m/Cny91XAy8CbgDUnW013DOKWtcgqwV2t/A7By3DVLku6188xdZl9VHQ8cP6n5\nKuCAKfp+DzhsHHVJkmbmN8clSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJwSJJ6MTgk\nSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJwSJJ62abgSPKHI9O7zl05kqTt3bTBkeSN\nSZ7Nvb8FDvDvc1uSJGl7NtNPx15B97Otj03yL8DldL///YSqumLOq5MkbXdmGqq6FXgzsB54AXBS\na1+Z5N/msC5J0nZqpjOOg4Djgf2BdwEXAXdV1dFzXZik2bVk5VmD7fuaEw8ZbN+afdOecVTVm6vq\nQOAa4MN0QbMwyVeSfHoM9UmStjMznXFM+EJVnQ+cn+S3qup5Sfaey8IkSdunbfo4blW9cWT2qNZ2\n81wUJEnavvX+AmBVXXR/d5pkQZIzk3wjyeVJnp1kzyRrklzZHvdofZPkpCTrk1yc5Bn3d/+SpPtu\nqG+O/xXw+ap6IvA0uo/5rgTWVtVSYG2bBzgYWNr+jgXeN/5yJUkTxh4cSR4BPB84BaCq7qmq24Dl\nwKrWbRVwaJteDnywOucCC5I8esxlS5KaIc44HgtsAj6Q5MIk70/yUOBRVXUDQHt8ZOu/CNgwsv7G\n1vYTkhybZF2SdZs2bZrbI5CkeWyI4NgZeAbwvqp6OnAX9w5LTSVTtNUWDVUnV9Wyqlq2cOHC2alU\nkrSFIYJjI7Cxqs5r82fSBcmNE0NQ7fGmkf77jqy/GLh+TLVKkiYZe3BU1beBDUme0JoOBC4DVgMr\nWtsK4FNtejVwZPt01bOA2yeGtCRJ47etXwCcbb8LnJ5kF+Aq4Gi6EPtYkmOAa+lurgjwWeCldPfL\nurv1lSQNZJDgqKqvAcumWHTgFH0LeM2cFyVJ2ib+AqAkqReDQ5LUi8EhSerF4JAk9WJwSJJ6MTgk\nSb0YHJKkXgwOSVIvBockqReDQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvBockqReD\nQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb0YHJKkXgwOSVIvgwVHkp2SXJjkM21+vyTnJbkyyUeT7NLa\nd23z69vyJUPVLEka9ozjdcDlI/N/Bry7qpYCtwLHtPZjgFur6nHAu1s/SdJABgmOJIuBQ4D3t/kA\nLwTObF1WAYe26eVtnrb8wNZfkjSAnQfa718CbwQe3ub3Am6rqs1tfiOwqE0vAjYAVNXmJLe3/jeP\nr9wHviUrzxq6BEk7iLGfcSR5GXBTVV0w2jxF19qGZaPbPTbJuiTrNm3aNAuVSpKmMsRQ1XOBVyS5\nBvgI3RDVXwILkkycAS0Grm/TG4F9Adry3YFbJm+0qk6uqmVVtWzhwoVzewSSNI+NPTiq6riqWlxV\nS4DDgbOr6gjgS8ArW7cVwKfa9Oo2T1t+dlVtccYhSRqP7el7HG8C3pBkPd01jFNa+ynAXq39DcDK\ngeqTJDHcxXEAqurLwJfb9FXAAVP0+R5w2FgLkyRt1fZ0xiFJ2gEYHJKkXgwOSVIvBockqReDQ5LU\ni8EhSerF4JAk9TLo9zgkzQ9D3UTzmhMPGWS/D3SecUiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS\n1IvBIUnqxeCQJPVicEiSejE4JEm9GBySpF4MDklSLwaHJKkXg0OS1IvBIUnqxeCQJPVicEiSehl7\ncCTZN8mXklye5NIkr2vteyZZk+TK9rhHa0+Sk5KsT3JxkmeMu2ZJ0r2GOOPYDPxeVT0JeBbwmiRP\nBlYCa6tqKbC2zQMcDCxtf8cC7xt/yZKkCWMPjqq6oaq+2qbvAC4HFgHLgVWt2yrg0Da9HPhgdc4F\nFiR59JjLliQ1g17jSLIEeDpwHvCoqroBunABHtm6LQI2jKy2sbVN3taxSdYlWbdp06a5LFuS5rXB\ngiPJw4B/AP5vVX1nuq5TtNUWDVUnV9Wyqlq2cOHC2SpTkjTJIMGR5KfoQuP0qvpEa75xYgiqPd7U\n2jcC+46svhi4fly1SpJ+0hCfqgpwCnB5Vb1rZNFqYEWbXgF8aqT9yPbpqmcBt08MaUmSxm/nAfb5\nXOBXga8n+VprezNwIvCxJMcA1wKHtWWfBV4KrAfuBo4eb7mSpFFjD46q+gpTX7cAOHCK/gW8Zk6L\nkiRtM785LknqxeCQJPVicEiSejE4JEm9GBySpF6G+DiutmLJyrOGLkGSZuQZhySpF4NDktSLwSFJ\n6sXgkCT1YnBIknoxOCRJvRgckqReDA5JUi8GhySpF785LukBa8i7MVxz4iGD7XuuecYhSerF4JAk\n9WJwSJJ6MTgkSb14cXwK3t5ckrbOMw5JUi8GhySpF4NDktTLDhMcSQ5KckWS9UlWDl2PJM1XO8TF\n8SQ7AX8D/AKwETg/yeqqumzYyiRpakN9yGYc31jfUc44DgDWV9VVVXUP8BFg+cA1SdK8tEOccQCL\ngA0j8xuBZ452SHIscGybvTPJFWOq7f7aG7h56CIGNt+fg/l+/OBzALP0HOTP7tfqj9mWTjtKcGSK\ntvqJmaqTgZPHU87sSbKuqpYNXceQ5vtzMN+PH3wOYMd6DnaUoaqNwL4j84uB6weqRZLmtR0lOM4H\nlibZL8kuwOHA6oFrkqR5aYcYqqqqzUl+B/gCsBNwalVdOnBZs2WHG16bA/P9OZjvxw8+B7ADPQep\nqpl7SZLU7ChDVZKk7YTBIUnqxeAYUJKdklyY5DND1zKEJAuSnJnkG0kuT/LsoWsatySvT3JpkkuS\nnJHkwUPXNNeSnJrkpiSXjLTtmWRNkivb4x5D1jiXtnL872z/H1yc5B+TLBiyxpkYHMN6HXD50EUM\n6K+Az1fVE4GnMc+eiySLgNcCy6rqKXQf/Dh82KrG4jTgoEltK4G1VbUUWNvmH6hOY8vjXwM8paqe\nCvwncNy4i+rD4BhIksXAIcD7h65lCEkeATwfOAWgqu6pqtuGrWoQOwMPSbIzsBvz4PtJVXUOcMuk\n5uXAqja9Cjh0rEWN0VTHX1VfrKrNbfZcuu+qbbcMjuH8JfBG4EdDFzKQxwKbgA+04br3J3no0EWN\nU1VdB/wFcC1wA3B7VX1x2KoG86iqugGgPT5y4HqG9Grgc0MXMR2DYwBJXgbcVFUXDF3LgHYGngG8\nr6qeDtzFA3t4YgttHH85sB+wD/DQJL8ybFUaUpK3AJuB04euZToGxzCeC7wiyTV0d/p9YZIPD1vS\n2G0ENlbVeW3+TLogmU9eBFxdVZuq6gfAJ4DnDFzTUG5M8miA9njTwPWMXZIVwMuAI2o7/4KdwTGA\nqjquqhZX1RK6i6FnV9W8eqdZVd8GNiR5Qms6EJhvv69yLfCsJLslCd1zMK8+IDBiNbCiTa8APjVg\nLWOX5CDgTcArquruoeuZyQ5xyxE9YP0ucHq7/9hVwNED1zNWVXVekjOBr9INT1zIDnTbifsqyRnA\nC4C9k2wEjgdOBD6W5Bi6QD1suArn1laO/zhgV2BN9x6Cc6vqNwcrcgbeckSS1ItDVZKkXgwOSVIv\nBockqReDQ5LUi8EhSerF4JAGkGRJkl++D+stSPLbI/P7tI/0SmNjcEjDWAJMGRzthodbswD4cXBU\n1fVV9crZLU2ansGheSXJke03Dy5K8qEkj0mytrWtTfLfW7/TkpyU5N+SXJXklSPbeGOSr7dtnNja\n9k/y+SQXJPmXJE+cYTsnAj+f5GvtNzmOSvLxJJ8GvpjkYa2er7Z9LR9Zb/+23jvbmcslbV8PTvKB\n1v/CJP+rtR+V5BOtviuT/PlYnmw9cFWVf/7Niz/gp4ErgL3b/J7Ap4EVbf7VwCfb9GnAx+neXD0Z\nWN/aDwb+DdhtYhvtcS2wtE0/k+42MtNt5wXAZ0ZqO4ru/l0T29sZeESb3htYD4TuTOWSkfV+PA/8\nHvCBNv1Eum9gP7ht+ypg9zb/LWDfof97+Lfj/nnLEc0nLwTOrKqbAarqlvarg/+7Lf8QMPpu/JNV\n9SPgsiSPam0vontxvntkGw+juznhx9vtIqC7fcR025nKmqqa+J2GAO9I8ny6W+8vAqZbF+B5wHta\nXd9I8i3g8W3Z2qq6HSDJZcBjgA0zbE+aksGh+STATPfYGV3+/Unrbm0bDwJuq6qf3co2p9rOVO4a\nmT4CWAj8XFX9oN1JeaaflZ1u26M1/BD/39f94DUOzSdrgV9Mshd0v3NNN+w08XOtRwBfmWEbXwRe\nnWS3iW1U1XeAq5Mc1tqS5GkzbOcO4OHTLN+d7jdbftCuVTxmG9Y7px0DSR4P/He6oTlpVhkcmjeq\n6lLg7cA/J7kIeBfdb34fneRi4Ffpfgd+um18nu4W4OuSfA34/bboCOCYtt1L6X6gaToXA5vbBfbX\nT7H8dGBZknVt299o+/8v4F+TXJLknZPWeS+wU5KvAx8Fjqqq7yPNMu+OK0nqxTMOSVIvBockqReD\nQ5LUi8EhSerF4JAk9WJwSJJ6MTgkSb38f74wl/QOtY4oAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "hist(mapk_marginal, 'map2k')" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "tensor(7.7734)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "mapk_do_dist = infer_dist(do_map3k, noise_marginals)\n", - "mapk_marginal = EmpiricalMarginal(mapk_do_dist, sites = 'mapk')" + "torch.stack([mapk_do_receiver(0.678, noise_posteriors)]*1000).mean()" ] }, { From 487f75abe8b7bc1aa31926c297c6afa062618a00 Mon Sep 17 00:00:00 2001 From: kaushalpaneri Date: Sat, 16 Feb 2019 14:48:39 -0500 Subject: [PATCH 11/11] receivers updated --- mapk_demon/receivers.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/mapk_demon/receivers.py b/mapk_demon/receivers.py index 083ccf1..d1ab6ae 100644 --- a/mapk_demon/receivers.py +++ b/mapk_demon/receivers.py @@ -68,9 +68,14 @@ def f_mapk(map2k, N_k, params, mode): return sample("mapk", Delta(f(mapk_mu, N_k))) -def mapk_companion(noise_dists): - mode = 'companion' +def mapk_receiver(noise_dists): + ''' + Mapk receiver model + :param noise_dists: N_3k, N_2k, N_k + :return: concentration of map3k, map2k and mapk at steady state + ''' # Parameters + mode = 'original' al_m = torch.tensor(700.) al_s = torch.tensor(1.) @@ -103,9 +108,14 @@ def mapk_companion(noise_dists): } -def mapk_receiver(noise_dists): +def mapk_companion(noise_dists): + ''' + Companion model that enables inference + :param noise_dists: N_3k, N_2k, N_k + :return: + ''' + mode = 'companion' # Parameters - mode = 'original' al_m = torch.tensor(700.) al_s = torch.tensor(1.) @@ -139,7 +149,15 @@ def mapk_receiver(noise_dists): def mapk_do_receiver(do_value,noise_dists): + ''' + This model validates the results from doing counterfactual inference on mapk_receiver + :param do_value: model uses this value to do hard intervention on map3k + :param noise_dists: N_3k, N_2k, N_k + :return: + ''' mode = 'original' + + # Parameters al_m = torch.tensor(700.) al_s = torch.tensor(1.)