diff --git a/.gitignore b/.gitignore index 4096078..49c06e6 100644 --- a/.gitignore +++ b/.gitignore @@ -97,4 +97,72 @@ venv.bak/ *.pdf # avoid uploading data -notebooks/data/* \ No newline at end of file +notebooks/data/* + +# avoid PyCharm files +*.iml +*.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/__init__.py b/mapk_demon/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/mapk_demon/inference.py b/mapk_demon/inference.py new file mode 100644 index 0000000..aa3575a --- /dev/null +++ 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 new file mode 100644 index 0000000..d1ab6ae --- /dev/null +++ b/mapk_demon/receivers.py @@ -0,0 +1,187 @@ +from __future__ import division + +from functools import partial + +import pyro +from pyro import sample +from pyro.distributions import Normal, Uniform, Delta + +import pyro.optim +import torch + + +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 + + +hyperparameters = { + "t_3k": 1.2, + "t_2k": 1.2, + "t_k": 0.003 +} + + +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)) + + 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, 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)) + + 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, 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)) + + 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_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.) + + 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_companion(noise_dists): + ''' + Companion model that enables inference + :param noise_dists: N_3k, N_2k, N_k + :return: + ''' + 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_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.) + + 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/mapk_demon/svi_test.py b/mapk_demon/svi_test.py new file mode 100644 index 0000000..1fe3d84 --- /dev/null +++ b/mapk_demon/svi_test.py @@ -0,0 +1,142 @@ +''' +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 + +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) diff --git a/mapk_demon/transmitters.py b/mapk_demon/transmitters.py new file mode 100644 index 0000000..24b2149 --- /dev/null +++ b/mapk_demon/transmitters.py @@ -0,0 +1,93 @@ +from __future__ import division + +from functools import partial + +import pyro +from pyro import sample +from pyro.distributions import Normal, Uniform +from torch import tensor + + +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 + + +hyperparameters = { + "t_3k": 1.2, + "t_2k": 1.2, + "t_k": 0.003 +} + +params = { + "alpha_3k": 500, + "alpha_2k": 500, + "alpha_k": 500, + "nu_3k": .15, + "nu_2k": .15, + "nu_k": .15 +} + + +def f_map3k(E1, N_3k): + total_3k = hyperparameters['t_3k'] + alpha_3k = params['alpha_3k'] + nu_3k = params['nu_3k'] + + map3k_mu = total_3k * g1(E1 * (alpha_3k/nu_3k)) + + return f(map3k_mu, N_3k) + + +def f_map2k(map3k, N_2k): + total_2k = hyperparameters['t_2k'] + alpha_2k = params['alpha_2k'] + nu_2k = params['nu_2k'] + + map2k_mu = total_2k * g2(map3k * (alpha_2k / nu_2k)) + + return f(map2k_mu, N_2k) + + +def f_mapk(map2k, N_k): + total_k = hyperparameters['t_k'] + alpha_k = params['alpha_k'] + nu_k = params['nu_k'] + + mapk_mu = total_k * g2(map2k * (alpha_k / nu_k)) + + return f(mapk_mu, N_k) + + +def mapk_signaling(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) + map2k = f_map2k(map3k,N_2k) + mapk = f_mapk(map2k, N_k) + + return { + 'map3k': map3k, + 'map2k': map2k, + 'mapk': mapk + } + +if __name__ == "__main__": + + # Noise Distributions + noise_vars = ['N_3k', 'N_2k', 'N_k'] + noise_prior = {N: Normal(0, 1) for N in noise_vars} + + print(mapk_signaling(noise_dists=noise_prior)) diff --git a/notebooks/Mapk Counterfactual Inference.ipynb b/notebooks/Mapk Counterfactual Inference.ipynb new file mode 100644 index 0000000..a339259 --- /dev/null +++ b/notebooks/Mapk Counterfactual Inference.ipynb @@ -0,0 +1,251 @@ +{ + "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, mapk_companion, mapk_do_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": "markdown", + "metadata": {}, + "source": [ + "Let $$N_{3k}, N_{2k}, N_k \\sim Normal\\ (10, 1)$$." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "noise_vars = ['N_3k', 'N_2k', 'N_k']\n", + "noise_dists = {N: Normal(10., 1.) for N in noise_vars}\n", + "\n", + "evidence = mapk_receiver(noise_dists)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Condition the model on observed map3k, map2k, mapk" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "conditioned_mapk = pyro.condition(mapk_companion, data=evidence)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 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", + "I have used HMC and MCMC for inference." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": true + }, + "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": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "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": {}, + "output_type": "display_data" + } + ], + "source": [ + "hist(noise_marginals['N_2k'], 'N_2k')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. do Map3k = 0.678" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "do_map3k = pyro.do(mapk_companion, data={'map3k': torch.tensor(50000.678)})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 4. Infer Mapk" + ] + }, + { + "cell_type": "code", + "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": [ + { + "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": [ + "torch.stack([mapk_do_receiver(0.678, noise_posteriors)]*1000).mean()" + ] + }, + { + "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