From a30c7025e9e735303316fed50d90cc3828f3101b Mon Sep 17 00:00:00 2001 From: Robert Ness Date: Sat, 27 Oct 2018 11:33:17 -0400 Subject: [PATCH 1/4] remove inference file --- causal_demon/inference.py | 21 --------------------- 1 file changed, 21 deletions(-) delete mode 100644 causal_demon/inference.py diff --git a/causal_demon/inference.py b/causal_demon/inference.py deleted file mode 100644 index 37abf35..0000000 --- a/causal_demon/inference.py +++ /dev/null @@ -1,21 +0,0 @@ -import pyro - - -def infer_dist(prog, n_dist): - """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. - - TODO(Kuashal): Do an implementation with stochastic variational inference. - http://pyro.ai/examples/svi_part_i.html - - - `prog`: the subroutine encoding the SCM. - `n_dist`: a dictionary containing distributions for each - noise object. - """ - return pyro.infer.Importance(prog, num_samples=5000).run(n_dist) From 0530b0e918765c62da37efd5f26f0887342d99fd Mon Sep 17 00:00:00 2001 From: Robert Ness Date: Sat, 27 Oct 2018 17:05:39 -0400 Subject: [PATCH 2/4] seems like it should work --- causal_demon/transmitters.py | 56 +++++-- notebooks/biochemical_model.ipynb | 256 ++++++++++++++++++++++++------ 2 files changed, 255 insertions(+), 57 deletions(-) diff --git a/causal_demon/transmitters.py b/causal_demon/transmitters.py index d84eb2c..c66bd63 100644 --- a/causal_demon/transmitters.py +++ b/causal_demon/transmitters.py @@ -4,6 +4,7 @@ import pyro from pyro import sample +import torch from torch import tensor """ @@ -17,7 +18,7 @@ marginals from conditioned programs will sometimes fail when this distribution is used. """ -Delta = partial(pyro.distributions.Normal, scale=1e-10) +Delta = partial(pyro.distributions.Normal, scale=1e-6) constants = { "sos_tot": 120000., @@ -140,7 +141,7 @@ def f_erk(mek, N_erk): return sample('erk', Delta(erk_mu + N_erk)) -def cancer_signaling(noise_dists): +def cancer_signaling(noise_dists, sample_shape=torch.Size([1])): """Model of Biochemical Signal Transduction in Lung Cancer This is the model of the EGFR and IGF1R pathway in lung cancer by Bianconi @@ -152,16 +153,49 @@ def cancer_signaling(noise_dists): equation modeling approach. This model assumes the ODEs have been solved for steady state, and models the steady state with structural causal models. For more details on this math, see the bianconi_math document. + + The following illustrates how you would simulate an experiment with this + model. + + Suppose there were 3 conditions: EGF is varied from low (800), + medium (2000), to high (8000) concentrations. Each condition had 4 + replicates. + + from pyro.distributions import LogNormal + + noise_vars = [ + 'N_egf', 'N_igf', 'N_sos', + 'N_ras', 'N_pi3k', 'N_akt', + 'N_raf', 'N_mek', 'N_erk' + ] + + noise_priors = {N: LogNormal(0, .001) for N in noise_vars} + + # replicates on the rows, conditions on the columns + # + conditions = { + 'egf': torch.tensor( + [[800., 8000., 80000.], + [800., 8000., 80000.], + [800., 8000., 80000.], + [800., 8000., 80000.]] + ), + 'igf': torch.zeros([4, 3]) + } + experiment = pyro.do(cancer_signaling, data=conditions) + # Get some Erk values + experiment(noise_priors, [4, 3])['erk'] """ - N_egf = sample('N_egf', noise_dists['N_egf']) - N_igf = sample('N_igf', noise_dists['N_igf']) - N_sos = sample('N_sos', noise_dists['N_sos']) - N_ras = sample('N_ras', noise_dists['N_ras']) - N_pi3k = sample('N_pi3k', noise_dists['N_pi3k']) - N_akt = sample('N_akt', noise_dists['N_akt']) - N_raf = sample('N_raf', noise_dists['N_raf']) - N_mek = sample('N_mek', noise_dists['N_mek']) - N_erk = sample('N_erk', noise_dists['N_erk']) + + N_egf = sample('N_egf', noise_dists['N_egf'], sample_shape=sample_shape) + N_igf = sample('N_igf', noise_dists['N_igf'], sample_shape=sample_shape) + N_sos = sample('N_sos', noise_dists['N_sos'], sample_shape=sample_shape) + N_ras = sample('N_ras', noise_dists['N_ras'], sample_shape=sample_shape) + N_pi3k = sample('N_pi3k', noise_dists['N_pi3k'], sample_shape=sample_shape) + N_akt = sample('N_akt', noise_dists['N_akt'], sample_shape=sample_shape) + N_raf = sample('N_raf', noise_dists['N_raf'], sample_shape=sample_shape) + N_mek = sample('N_mek', noise_dists['N_mek'], sample_shape=sample_shape) + N_erk = sample('N_erk', noise_dists['N_erk'], sample_shape=sample_shape) egf = f_egf(N_egf) igf = f_igf(N_igf) diff --git a/notebooks/biochemical_model.ipynb b/notebooks/biochemical_model.ipynb index f4db6c9..7dd5fbd 100644 --- a/notebooks/biochemical_model.ipynb +++ b/notebooks/biochemical_model.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Counterfactual Modeling of aa Biochemical Network\n", + "# Counterfactual Modeling with a Biochemical Network\n", "\n", "This workflow demonstrates counterfactual modeling using a [dynamic model of biochemical signaling in lung cancer](https://wwwdev.ebi.ac.uk/biomodels/BIOMD0000000427).\n", "\n", @@ -32,23 +32,56 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "from math import exp\n", "\n", + "from pprint import pprint\n", + "import pyro\n", "from pyro.distributions import LogNormal\n", "from pyro import condition, do, infer, sample\n", "from pyro.infer import EmpiricalMarginal\n", - "from torch import tensor\n", + "from pyro.infer.mcmc import MCMC\n", + "from pyro.infer.mcmc.nuts import NUTS\n", + "import torch\n", "\n", - "from causal_demon.inference import infer_dist\n", "from causal_demon.transmitters import cancer_signaling\n", "\n", "from matplotlib import pyplot as plt\n", - "%matplotlib inline\n", + "%matplotlib inline\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Helper functions for inference\n", + "\n", + "Since the inference algorithm is not the focus of this tutorial, I hide the choice of algorithm in an abstraction." + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "def infer_noise(model):\n", + " nuts_kernel = NUTS(model, adapt_step_size=True)\n", + " return MCMC(nuts_kernel, num_samples=500, warmup_steps=300)\n", "\n", + "def infer_noise2(model):\n", + " return pyro.infer.Importance(model, num_samples=1000)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ "def hist(marginal, name):\n", " plt.hist([marginal() for _ in range(5000)])\n", " plt.title(\"Marginal Histogram of {}\".format(name))\n", @@ -60,121 +93,237 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Each noise term is modeled with a weakly informed prior." + "# Simulating \"wild type\" data from the model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\"Wild type\" refers to data forward-generated from the model. \n", + "\n", + "Pass in noise variables that capture the variation in the model. Each noise term is modeled with a weakly informed prior. The `independent` method reshapes the distribution such that data replicates will be independent." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "noise_vars = ['N_egf', 'N_igf', 'N_sos', 'N_ras', 'N_pi3k', 'N_akt', 'N_raf', 'N_mek', 'N_erk']\n", - "noise_prior = {N: LogNormal(0, 10) for N in noise_vars}" + "noise_priors = {N: LogNormal(0, 10).independent() for N in noise_vars}" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "cancer_signaling(noise_prior)" + "`sample_shape` captures the dimensions of the data. Here we simulate 4 observations." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "{'egf': tensor([ 1.3451e-04, 8.5265e-01, -5.9698e-07, 1.2893e+05]),\n", + " 'igf': tensor([1.1811e+09, 8.2959e-05, 2.6991e-05, 4.0145e+07]),\n", + " 'sos': tensor([1.1957e+05, 2.6486e-01, 3.8057e-03, 1.1350e+05]),\n", + " 'ras': tensor([55184.7930, 1.4579, 77630.7500, 53698.1367]),\n", + " 'pi3k': tensor([ 120010.9844, 1179.5409, 4540044.5000, 119999.7891]),\n", + " 'akt': tensor([405095.7500, 29134.3945, 598950.0625, 437960.0312]),\n", + " 'raf': tensor([2171.0269, 0.1378, 1444.7947, 1316.2133]),\n", + " 'mek': tensor([77864.0000, 4.9157, 47464.5547, 43546.7773]),\n", + " 'erk': tensor([428376.6875, 94.5480, 362289.4375, 349577.0000])}" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "# Experimental use only, please ignore\n", - "\n", - "#noise_priors = {N: LogNormal(0, .001) for N in noise_vars}\n", - "#evidence = {'egf': tensor(800.), 'igf': tensor(2.)}\n", - "#scm_obs = condition(scm, data=evidence)\n", - "#scm_obs(noise_priors)" + "cancer_signaling(noise_priors, sample_shape=torch.Size([4]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Counterfactual Inference\n", + "# Simulating an intervention experiment with the do-operator\n", + "\n", + "Suppose we have an *intervention query*:\n", + "\n", + "### Intervention Query: *What would happen if we blocked IGF and varied the amount of EGF?*\n", + "\n", + "We could answer this query with an experiment:\n", + "* 3 conditions: EGF is varied from low (set to 800), medium (8000), to high (80000) concentrations\n", + "* 2 replicates of each condition.\n", + "* IGF is blocked (set to 0).\n", "\n", - "The goal is to observe the system under natural (or perhaps experimental) conditions and then use this to make counterfactual predictions - i.e. use observations to inform inferences on what the system's behavior would have been if the observations had been different." + "Prior to running the experiment, we can simulate the experiment using this model. \n", + "\n", + "**Why simulate experiments?** Interventions with the `do` operator allow you to simulate the results of an experiment prior to spending the resources to conduct the experiment.\n", + "\n", + "The interventions on IFG and EGF in this example are one set from a very large set of possible interventions. Which sent of interventions is the best? You can answer this by defining a cost function on the experimental results, applying that function to simulated data, and iteratively searching for a set of interventions that optimizes that cost function. " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'akt': tensor([[391338.7188, 404718.7812, 407006.2500],\n", + " [392281.4375, 403664.8438, 573941.7500]]),\n", + " 'egf': tensor([[ 800., 8000., 80000.],\n", + " [ 800., 8000., 80000.]]),\n", + " 'erk': tensor([[5.5875e+03, 5.6264e+05, 2.0091e+04],\n", + " [9.7106e+05, 9.5360e+08, 5.2224e+05]]),\n", + " 'igf': tensor([[0., 0., 0.],\n", + " [0., 0., 0.]]),\n", + " 'mek': tensor([[ 293.2357, 469678.3750, 956.1735],\n", + " [ 52.5269, 1588.2506, 209497.5625]]),\n", + " 'pi3k': tensor([[ 108291.1875, 119667.8672, 119870.2812],\n", + " [ 109045.3203, 118715.5391, 1271766.0000]]),\n", + " 'raf': tensor([[ 8.2239, 60607.8164, 11.9922],\n", + " [ 1.4725, 44.3995, 9023.0674]]),\n", + " 'ras': tensor([[2.2908e+02, 4.7263e+06, 1.1677e+02],\n", + " [5.5069e+01, 1.1748e+03, 4.3400e+01]]),\n", + " 'sos': tensor([[ 0.5087, 5.0809, 136.9033],\n", + " [ 51.3112, 6.9473, 50.8112]])}\n" + ] + } + ], + "source": [ + "conditions = {\n", + " 'egf': torch.tensor(\n", + " [\n", + " [800., 8000., 80000.],\n", + " [800., 8000., 80000.],\n", + " ]\n", + " ),\n", + " 'igf': torch.zeros([2, 3])\n", + "}\n", + "initial_experiment = do(cancer_signaling, data=conditions)\n", + "initial_results = initial_experiment(noise_priors, torch.Size([2, 3]))\n", + "pprint(initial_results)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Conditioning on data from previous observations\n", + "# Simulating a follow-up experiment using counterfactual inference\n", + "\n", + "Suppose that after that initial experiment, we want to try a new experiment under a different set of conditions. We could repeat the simulation we just ran under the new conditions, but we can improve the quality of the simulation by using the old data.\n", "\n", - "A scientist might observe values for each of the variables (or a subset thereof) in previous experiments.\n", + "However, the old data was collected under different experimental conditions than in the proposed experiment. This is addressed by posing a counterfactual query to the model: \n", "\n", - "Suppose that in this experiment Igf was blocked. Egf and Erk were observed to have concentration values of 800. in these settings.\n", + "### Counterfactual Query: *Given the observations of Erk after fixing IGF to 0, what would Erk have been if IGF were set to 8000?*\n", "\n", - "**Counterfactual query**: What would Erk levels be if there had Igf concentration also been 800?" + "The following illustrates how to use the counterfactual simulation algorithm to answer this query." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "1. Condition program on EGF being 800, IGF being 0, and Erk being 800, and infer the conditional distribution." + "### Conditioning on previous observations\n", + "\n", + "Firstly, we pretend the simulated data from the previous section is real, and that we only collected data for Erk." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "2" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "evidence = {'egf': tensor(800.), 'igf': tensor(0.), 'erk': tensor(800.)}\n", - "cancer_obs = condition(cancer_signaling, data=evidence)" + "erk_data = initial_results['erk']\n", + "erk_data.shape[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "2. Infer an observational distribution." + "## Step 1: Condition the model of the initial experiment on the observed data" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ - "cancer_dist = infer_dist(cancer_obs, noise_prior)" + "with pyro.iarange(\"erk\", erk_data.shape[0]) as ind:\n", + " evidence = {'erk': initial_results['erk'][ind, :]}\n", + " initial_experiment_conditioned = condition(initial_experiment, data=evidence)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "3. Do posterior inference on the noise variables, and obtain a marginal distribution for each variable." + "## Step 2: Update the noise priors given the observed data." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 30, "metadata": {}, - "outputs": [], + "outputs": [ + { + "ename": "TypeError", + "evalue": "log_prob() got an unexpected keyword argument 'sample_shape'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m~/Dropbox/code/ppl_causal/venv/lib/python3.6/site-packages/pyro/poutine/trace_struct.py\u001b[0m in \u001b[0;36mlog_prob_sum\u001b[0;34m(self, site_filter)\u001b[0m\n\u001b[1;32m 228\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 229\u001b[0;31m \u001b[0msite_log_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"log_prob_sum\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 230\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyError\u001b[0m: 'log_prob_sum'", + "\nDuring handling of the above exception, another exception occurred:\n", + "\u001b[0;31mTypeError\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[0mmodel_posterior\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minfer_noise2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minitial_experiment_conditioned\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnoise_priors\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;31m#noise_posteriors = {\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;31m# n: EmpiricalMarginal(model_posterior, sites=n)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;31m# for n in noise_vars\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;31m#}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Dropbox/code/ppl_causal/venv/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~/Dropbox/code/ppl_causal/venv/lib/python3.6/site-packages/pyro/infer/importance.py\u001b[0m in \u001b[0;36m_traces\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 43\u001b[0m model_trace = poutine.trace(\n\u001b[1;32m 44\u001b[0m poutine.replay(self.model, trace=guide_trace)).get_trace(*args, **kwargs)\n\u001b[0;32m---> 45\u001b[0;31m \u001b[0mlog_weight\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel_trace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_prob_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mguide_trace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_prob_sum\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 46\u001b[0m \u001b[0;32myield\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmodel_trace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlog_weight\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Dropbox/code/ppl_causal/venv/lib/python3.6/site-packages/pyro/poutine/trace_struct.py\u001b[0m in \u001b[0;36mlog_prob_sum\u001b[0;34m(self, site_filter)\u001b[0m\n\u001b[1;32m 230\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 231\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"args\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"kwargs\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 232\u001b[0;31m \u001b[0msite_log_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"fn\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_prob\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"value\"\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[0m\n\u001b[0m\u001b[1;32m 233\u001b[0m \u001b[0msite_log_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscale_tensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msite_log_p\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"scale\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 234\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"log_prob_sum\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msite_log_p\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mTypeError\u001b[0m: log_prob() got an unexpected keyword argument 'sample_shape'" + ] + } + ], "source": [ - "noise_marginals = {\n", - " n: EmpiricalMarginal(cancer_dist, sites=n)\n", - " for n in noise_vars\n", - "}" + "model_posterior = infer_noise2(initial_experiment_conditioned).run(noise_priors, sample_shape=torch.Size([2, 3]))\n", + "#noise_posteriors = {\n", + "# n: EmpiricalMarginal(model_posterior, sites=n)\n", + "# for n in noise_vars\n", + "#}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "4. Apply do-operator to original program to obtain intervention program." + "## Step 3. Apply do-operator to original model to obtain a model of the new experiment" ] }, { @@ -183,14 +332,28 @@ "metadata": {}, "outputs": [], "source": [ - "cancer_do = do(cancer_signaling, data={'igf': tensor(800.)})" + "conditions = {\n", + " 'egf': torch.tensor(\n", + " [\n", + " [800., 8000., 80000.],\n", + " [800., 8000., 80000.],\n", + " ]\n", + " ),\n", + " 'igf': torch.tensor(\n", + " [\n", + " [8000., 8000., 8000.],\n", + " [8000., 8000., 8000.],\n", + " ]\n", + " ),\n", + "}\n", + "new_experiment = do(cancer_signaling, data=conditions)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "5. Pass updated noise marginals to intervention program, and obtain counterfactual distribution on Erk." + "### 4. Pass the updated noise marginals to the model of the new experiment." ] }, { @@ -199,8 +362,7 @@ "metadata": {}, "outputs": [], "source": [ - "cancer_do_dist = infer_dist(cancer_do, noise_marginals)\n", - "erk_cf_marginal = EmpiricalMarginal(cancer_do_dist, sites = 'erk')" + "noise_marginals['N_egf'].sample()" ] }, { @@ -209,6 +371,8 @@ "metadata": {}, "outputs": [], "source": [ + "cancer_do_dist = infer_dist(cancer_do, noise_marginals)\n", + "erk_cf_marginal = EmpiricalMarginal(cancer_do_dist, sites = 'erk')\n", "hist(erk_cf_marginal, 'Erk')" ] }, From c2d1e8413e93cdd4d88f4bb62de680ddcc74eb1f Mon Sep 17 00:00:00 2001 From: Robert Ness Date: Sat, 27 Oct 2018 19:14:17 -0400 Subject: [PATCH 3/4] working with multiple dimensions --- causal_demon/transmitters.py | 72 ++++------- notebooks/biochemical_model.ipynb | 192 ++++++++++-------------------- 2 files changed, 84 insertions(+), 180 deletions(-) diff --git a/causal_demon/transmitters.py b/causal_demon/transmitters.py index c66bd63..94dc271 100644 --- a/causal_demon/transmitters.py +++ b/causal_demon/transmitters.py @@ -141,7 +141,7 @@ def f_erk(mek, N_erk): return sample('erk', Delta(erk_mu + N_erk)) -def cancer_signaling(noise_dists, sample_shape=torch.Size([1])): +def cancer_signaling(noise_dists): """Model of Biochemical Signal Transduction in Lung Cancer This is the model of the EGFR and IGF1R pathway in lung cancer by Bianconi @@ -154,58 +154,28 @@ def cancer_signaling(noise_dists, sample_shape=torch.Size([1])): for steady state, and models the steady state with structural causal models. For more details on this math, see the bianconi_math document. - The following illustrates how you would simulate an experiment with this - model. - - Suppose there were 3 conditions: EGF is varied from low (800), - medium (2000), to high (8000) concentrations. Each condition had 4 - replicates. - - from pyro.distributions import LogNormal - - noise_vars = [ - 'N_egf', 'N_igf', 'N_sos', - 'N_ras', 'N_pi3k', 'N_akt', - 'N_raf', 'N_mek', 'N_erk' - ] - - noise_priors = {N: LogNormal(0, .001) for N in noise_vars} - - # replicates on the rows, conditions on the columns - # - conditions = { - 'egf': torch.tensor( - [[800., 8000., 80000.], - [800., 8000., 80000.], - [800., 8000., 80000.], - [800., 8000., 80000.]] - ), - 'igf': torch.zeros([4, 3]) - } - experiment = pyro.do(cancer_signaling, data=conditions) - # Get some Erk values - experiment(noise_priors, [4, 3])['erk'] """ - N_egf = sample('N_egf', noise_dists['N_egf'], sample_shape=sample_shape) - N_igf = sample('N_igf', noise_dists['N_igf'], sample_shape=sample_shape) - N_sos = sample('N_sos', noise_dists['N_sos'], sample_shape=sample_shape) - N_ras = sample('N_ras', noise_dists['N_ras'], sample_shape=sample_shape) - N_pi3k = sample('N_pi3k', noise_dists['N_pi3k'], sample_shape=sample_shape) - N_akt = sample('N_akt', noise_dists['N_akt'], sample_shape=sample_shape) - N_raf = sample('N_raf', noise_dists['N_raf'], sample_shape=sample_shape) - N_mek = sample('N_mek', noise_dists['N_mek'], sample_shape=sample_shape) - N_erk = sample('N_erk', noise_dists['N_erk'], sample_shape=sample_shape) - - egf = f_egf(N_egf) - igf = f_igf(N_igf) - sos = f_sos(egf, igf, N_sos) - ras = f_ras(sos, N_ras) - pi3k = f_pi3k(ras, egf, igf, N_pi3k) - akt = f_akt(pi3k, N_akt) - raf = f_raf(ras, akt, N_raf) - mek = f_mek(raf, N_mek) - erk = f_erk(mek, N_erk) + with pyro.iarange("model"): + N_egf = sample('N_egf', noise_dists['N_egf']) + N_igf = sample('N_igf', noise_dists['N_igf']) + N_sos = sample('N_sos', noise_dists['N_sos']) + N_ras = sample('N_ras', noise_dists['N_ras']) + N_pi3k = sample('N_pi3k', noise_dists['N_pi3k']) + N_akt = sample('N_akt', noise_dists['N_akt']) + N_raf = sample('N_raf', noise_dists['N_raf']) + N_mek = sample('N_mek', noise_dists['N_mek']) + N_erk = sample('N_erk', noise_dists['N_erk']) + + egf = f_egf(N_egf) + igf = f_igf(N_igf) + sos = f_sos(egf, igf, N_sos) + ras = f_ras(sos, N_ras) + pi3k = f_pi3k(ras, egf, igf, N_pi3k) + akt = f_akt(pi3k, N_akt) + raf = f_raf(ras, akt, N_raf) + mek = f_mek(raf, N_mek) + erk = f_erk(mek, N_erk) return { 'egf': egf, diff --git a/notebooks/biochemical_model.ipynb b/notebooks/biochemical_model.ipynb index 7dd5fbd..f302465 100644 --- a/notebooks/biochemical_model.ipynb +++ b/notebooks/biochemical_model.ipynb @@ -32,7 +32,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -46,6 +46,8 @@ "from pyro.infer.mcmc import MCMC\n", "from pyro.infer.mcmc.nuts import NUTS\n", "import torch\n", + "import warnings\n", + "warnings.filterwarnings('ignore')\n", "\n", "from causal_demon.transmitters import cancer_signaling\n", "\n", @@ -64,29 +66,30 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "def infer_noise(model):\n", + "def infer_mcmc(model):\n", " nuts_kernel = NUTS(model, adapt_step_size=True)\n", " return MCMC(nuts_kernel, num_samples=500, warmup_steps=300)\n", "\n", - "def infer_noise2(model):\n", + "def infer_importance(model):\n", " return pyro.infer.Importance(model, num_samples=1000)" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "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(\"#\")" + "def plot_erk(data):\n", + " plt.boxplot(data['erk'])\n", + " plt.title(\"Empirical Distribution of Erk\")\n", + " plt.xlabel(\"Condition\")\n", + " plt.ylabel(\"Concentration\")\n", + " plt.show()" ] }, { @@ -102,52 +105,28 @@ "source": [ "\"Wild type\" refers to data forward-generated from the model. \n", "\n", - "Pass in noise variables that capture the variation in the model. Each noise term is modeled with a weakly informed prior. The `independent` method reshapes the distribution such that data replicates will be independent." + "Pass in noise variables that capture the variation in the model. Each noise term is modeled with a weakly informed prior." ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "noise_vars = ['N_egf', 'N_igf', 'N_sos', 'N_ras', 'N_pi3k', 'N_akt', 'N_raf', 'N_mek', 'N_erk']\n", - "noise_priors = {N: LogNormal(0, 10).independent() for N in noise_vars}" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "`sample_shape` captures the dimensions of the data. Here we simulate 4 observations." + "# 4 observations\n", + "mean = torch.zeros(4)\n", + "noise_priors = {N: LogNormal(mean, 10) for N in noise_vars}" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "{'egf': tensor([ 1.3451e-04, 8.5265e-01, -5.9698e-07, 1.2893e+05]),\n", - " 'igf': tensor([1.1811e+09, 8.2959e-05, 2.6991e-05, 4.0145e+07]),\n", - " 'sos': tensor([1.1957e+05, 2.6486e-01, 3.8057e-03, 1.1350e+05]),\n", - " 'ras': tensor([55184.7930, 1.4579, 77630.7500, 53698.1367]),\n", - " 'pi3k': tensor([ 120010.9844, 1179.5409, 4540044.5000, 119999.7891]),\n", - " 'akt': tensor([405095.7500, 29134.3945, 598950.0625, 437960.0312]),\n", - " 'raf': tensor([2171.0269, 0.1378, 1444.7947, 1316.2133]),\n", - " 'mek': tensor([77864.0000, 4.9157, 47464.5547, 43546.7773]),\n", - " 'erk': tensor([428376.6875, 94.5480, 362289.4375, 349577.0000])}" - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "cancer_signaling(noise_priors, sample_shape=torch.Size([4]))" + "cancer_signaling(noise_priors)" ] }, { @@ -174,49 +153,36 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "{'akt': tensor([[391338.7188, 404718.7812, 407006.2500],\n", - " [392281.4375, 403664.8438, 573941.7500]]),\n", - " 'egf': tensor([[ 800., 8000., 80000.],\n", - " [ 800., 8000., 80000.]]),\n", - " 'erk': tensor([[5.5875e+03, 5.6264e+05, 2.0091e+04],\n", - " [9.7106e+05, 9.5360e+08, 5.2224e+05]]),\n", - " 'igf': tensor([[0., 0., 0.],\n", - " [0., 0., 0.]]),\n", - " 'mek': tensor([[ 293.2357, 469678.3750, 956.1735],\n", - " [ 52.5269, 1588.2506, 209497.5625]]),\n", - " 'pi3k': tensor([[ 108291.1875, 119667.8672, 119870.2812],\n", - " [ 109045.3203, 118715.5391, 1271766.0000]]),\n", - " 'raf': tensor([[ 8.2239, 60607.8164, 11.9922],\n", - " [ 1.4725, 44.3995, 9023.0674]]),\n", - " 'ras': tensor([[2.2908e+02, 4.7263e+06, 1.1677e+02],\n", - " [5.5069e+01, 1.1748e+03, 4.3400e+01]]),\n", - " 'sos': tensor([[ 0.5087, 5.0809, 136.9033],\n", - " [ 51.3112, 6.9473, 50.8112]])}\n" - ] - } - ], + "outputs": [], "source": [ + "mean = torch.zeros([3, 2])\n", + "noise_priors = {N: LogNormal(mean, 1) for N in noise_vars}\n", "conditions = {\n", " 'egf': torch.tensor(\n", " [\n", - " [800., 8000., 80000.],\n", - " [800., 8000., 80000.],\n", + " [800., 800.],\n", + " [8000., 8000.],\n", + " [80000., 80000.]\n", " ]\n", " ),\n", - " 'igf': torch.zeros([2, 3])\n", + " 'igf': torch.zeros([3, 2])\n", "}\n", "initial_experiment = do(cancer_signaling, data=conditions)\n", - "initial_results = initial_experiment(noise_priors, torch.Size([2, 3]))\n", + "initial_results = initial_experiment(noise_priors)\n", "pprint(initial_results)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plot_erk(initial_results)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -227,7 +193,7 @@ "\n", "However, the old data was collected under different experimental conditions than in the proposed experiment. This is addressed by posing a counterfactual query to the model: \n", "\n", - "### Counterfactual Query: *Given the observations of Erk after fixing IGF to 0, what would Erk have been if IGF were set to 8000?*\n", + "### Counterfactual Query: *Given the observations of Erk after fixing IGF to 0, what would Erk have been if IGF were set to 800?*\n", "\n", "The following illustrates how to use the counterfactual simulation algorithm to answer this query." ] @@ -243,23 +209,11 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "2" - ] - }, - "execution_count": 14, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ - "erk_data = initial_results['erk']\n", - "erk_data.shape[0]" + "erk_data = initial_results['erk']\n" ] }, { @@ -271,13 +225,12 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "with pyro.iarange(\"erk\", erk_data.shape[0]) as ind:\n", - " evidence = {'erk': initial_results['erk'][ind, :]}\n", - " initial_experiment_conditioned = condition(initial_experiment, data=evidence)" + "evidence = {'erk': initial_results['erk']}\n", + "initial_experiment_conditioned = condition(initial_experiment, data=evidence)" ] }, { @@ -289,34 +242,15 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "TypeError", - "evalue": "log_prob() got an unexpected keyword argument 'sample_shape'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m~/Dropbox/code/ppl_causal/venv/lib/python3.6/site-packages/pyro/poutine/trace_struct.py\u001b[0m in \u001b[0;36mlog_prob_sum\u001b[0;34m(self, site_filter)\u001b[0m\n\u001b[1;32m 228\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 229\u001b[0;31m \u001b[0msite_log_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"log_prob_sum\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 230\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mKeyError\u001b[0m: 'log_prob_sum'", - "\nDuring handling of the above exception, another exception occurred:\n", - "\u001b[0;31mTypeError\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[0mmodel_posterior\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minfer_noise2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minitial_experiment_conditioned\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnoise_priors\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;31m#noise_posteriors = {\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;31m# n: EmpiricalMarginal(model_posterior, sites=n)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;31m# for n in noise_vars\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;31m#}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Dropbox/code/ppl_causal/venv/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~/Dropbox/code/ppl_causal/venv/lib/python3.6/site-packages/pyro/infer/importance.py\u001b[0m in \u001b[0;36m_traces\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 43\u001b[0m model_trace = poutine.trace(\n\u001b[1;32m 44\u001b[0m poutine.replay(self.model, trace=guide_trace)).get_trace(*args, **kwargs)\n\u001b[0;32m---> 45\u001b[0;31m \u001b[0mlog_weight\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodel_trace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_prob_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mguide_trace\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_prob_sum\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 46\u001b[0m \u001b[0;32myield\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmodel_trace\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlog_weight\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/Dropbox/code/ppl_causal/venv/lib/python3.6/site-packages/pyro/poutine/trace_struct.py\u001b[0m in \u001b[0;36mlog_prob_sum\u001b[0;34m(self, site_filter)\u001b[0m\n\u001b[1;32m 230\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 231\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"args\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"kwargs\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 232\u001b[0;31m \u001b[0msite_log_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"fn\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog_prob\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"value\"\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[0m\n\u001b[0m\u001b[1;32m 233\u001b[0m \u001b[0msite_log_p\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscale_tensor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msite_log_p\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"scale\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 234\u001b[0m \u001b[0msite\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"log_prob_sum\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msite_log_p\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mTypeError\u001b[0m: log_prob() got an unexpected keyword argument 'sample_shape'" - ] - } - ], + "outputs": [], "source": [ - "model_posterior = infer_noise2(initial_experiment_conditioned).run(noise_priors, sample_shape=torch.Size([2, 3]))\n", - "#noise_posteriors = {\n", - "# n: EmpiricalMarginal(model_posterior, sites=n)\n", - "# for n in noise_vars\n", - "#}" + "model_posterior = infer_mcmc(initial_experiment_conditioned).run(noise_priors)\n", + "noise_posteriors = {\n", + " n: EmpiricalMarginal(model_posterior, sites=n)\n", + " for n in noise_vars\n", + "}" ] }, { @@ -335,18 +269,20 @@ "conditions = {\n", " 'egf': torch.tensor(\n", " [\n", - " [800., 8000., 80000.],\n", - " [800., 8000., 80000.],\n", + " [800., 800.],\n", + " [8000., 8000.],\n", + " [80000., 80000.],\n", " ]\n", " ),\n", " 'igf': torch.tensor(\n", " [\n", - " [8000., 8000., 8000.],\n", - " [8000., 8000., 8000.],\n", + " [800., 800.],\n", + " [800., 800.],\n", + " [800., 800.],\n", " ]\n", " ),\n", "}\n", - "new_experiment = do(cancer_signaling, data=conditions)" + "counterfactual_experiment = do(cancer_signaling, data=conditions)" ] }, { @@ -362,7 +298,7 @@ "metadata": {}, "outputs": [], "source": [ - "noise_marginals['N_egf'].sample()" + "new_results = counterfactual_experiment(noise_posteriors)" ] }, { @@ -371,9 +307,7 @@ "metadata": {}, "outputs": [], "source": [ - "cancer_do_dist = infer_dist(cancer_do, noise_marginals)\n", - "erk_cf_marginal = EmpiricalMarginal(cancer_do_dist, sites = 'erk')\n", - "hist(erk_cf_marginal, 'Erk')" + "plot_erk(new_results)" ] }, { @@ -382,7 +316,7 @@ "source": [ "## Extentions\n", "\n", - "* Condition on IID data\n", + "* ~~Condition on IID data~~\n", "* Use SVI\n", "* Extend the model to accomadate uncertainty in parameters" ] From 1720e3249c5eb19664b12a4a2de7df13d60d8d96 Mon Sep 17 00:00:00 2001 From: Robert Ness Date: Mon, 21 Jan 2019 14:42:39 -0500 Subject: [PATCH 4/4] mod --- notebooks/biochemical_model.ipynb | 117 +++++++++++++++++++++++++----- 1 file changed, 98 insertions(+), 19 deletions(-) diff --git a/notebooks/biochemical_model.ipynb b/notebooks/biochemical_model.ipynb index f302465..fa1dc7c 100644 --- a/notebooks/biochemical_model.ipynb +++ b/notebooks/biochemical_model.ipynb @@ -32,7 +32,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -66,7 +66,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -80,7 +80,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -110,7 +110,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -122,9 +122,28 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "{'egf': tensor([2.7954e+01, 4.8299e-05, 1.5163e+03, 9.9427e-03]),\n", + " 'igf': tensor([9.4691e-07, 3.3414e-01, 3.6192e-07, 4.3954e-01]),\n", + " 'sos': tensor([6.6813e+05, 9.2959e-03, 9.6350e-01, 1.4639e-02]),\n", + " 'ras': tensor([3.1628e+05, 3.5068e+08, 8.2278e-01, 1.3921e-02]),\n", + " 'pi3k': tensor([113766.5859, 120093.1250, 113535.1094, 619.6009]),\n", + " 'akt': tensor([397999.4688, 405185.7500, 397726.4688, 6370.0747]),\n", + " 'raf': tensor([7.7653e+03, 1.1844e+05, 2.1608e-02, 8.0440e+02]),\n", + " 'mek': tensor([1.8952e+05, 5.3718e+05, 1.7275e+00, 2.7386e+04]),\n", + " 'erk': tensor([5.1520e+05, 5.7953e+05, 1.4578e+08, 2.8049e+05])}" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "cancer_signaling(noise_priors)" ] @@ -153,12 +172,46 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'akt': tensor([[391321.7812, 391321.9688],\n", + " [403665.8125, 403666.0000],\n", + " [404943.1250, 404943.2812]]),\n", + " 'egf': tensor([[ 800., 800.],\n", + " [ 8000., 8000.],\n", + " [80000., 80000.]]),\n", + " 'erk': tensor([[ 696.6604, 773.3032],\n", + " [ 715.8817, 864.9516],\n", + " [1529.1610, 1453.6689]]),\n", + " 'igf': tensor([[0., 0.],\n", + " [0., 0.],\n", + " [0., 0.]]),\n", + " 'mek': tensor([[36.2059, 40.2062],\n", + " [37.2172, 44.9751],\n", + " [79.6481, 75.7078]]),\n", + " 'pi3k': tensor([[108277.4531, 108277.5469],\n", + " [118715.4531, 118715.6094],\n", + " [119870.9766, 119871.1562]]),\n", + " 'raf': tensor([[0.9847, 1.0942],\n", + " [1.0185, 1.2332],\n", + " [2.2002, 2.0968]]),\n", + " 'ras': tensor([[ 2.3053, 2.2083],\n", + " [ 6.3968, 6.3944],\n", + " [45.1486, 45.3369]]),\n", + " 'sos': tensor([[ 1.5161, 1.4281],\n", + " [ 6.1106, 6.0218],\n", + " [51.8502, 51.8491]])}\n" + ] + } + ], "source": [ "mean = torch.zeros([3, 2])\n", - "noise_priors = {N: LogNormal(mean, 1) for N in noise_vars}\n", + "noise_priors = {N: LogNormal(mean, .1) for N in noise_vars}\n", "conditions = {\n", " 'egf': torch.tensor(\n", " [\n", @@ -176,9 +229,22 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAGMBJREFUeJzt3Xm4ZHV95/H3R1CWiNItHZZmaQMIg8QFG1wwjIoi4AJJVESigEhrxhjNZNwdWdTRjBnRqNHgiKBRQMUFkYmiIsQIYoPsoraI0ijS2M3mQli+88c5F8rrXer0vXWrbvf79Tz1dNXvnPqdb1X1U597fr8656SqkCSpXw8YdgGSpPnF4JAkdWJwSJI6MTgkSZ0YHJKkTgwOSVInBodGSpI7kvzJFMs/nOR/znAbT0mycoZ9/L8kh8+kj56+/izJD3oeX5fk6bPRd9vfVUmeMlv99bnNJPlYkjVJLpqlPmf1fdHa23DYBWj0JbkO2BK4p6f55Kr6m9neVlU9eJrlr5jtbY6XpIDfAAXcCVwKnFhVp/fUcUCHvnauqhWTrVNV/w7sMqOi79/eycDKqnpLT/+PnI2+O3oy8Axg26r69fiFSY4APgr8dtyiR1TVzwdfnmbC4FC/nlNVXxtmAUk2qKp7pl9zVjy6qlYk2QI4APhAkl2r6rjZ3EiSDavq7tnsc0TsAFw3UWj0uKCqnjxdR+vwezRvOVSlGUlyRJL/SHJCkluSXJvkSW379Ulu6h3SSXJyO9x0TpLbk5yXZIee5ZVkp551P5Tk7CS/Bp7atr29Z/2Dklya5LYkP06yf9t+ZJLvt9u4NsnL1+b1VdXNVfUJ4K+BNyZ5WNv/N5O8rL2/U/s6bk1yc5LT2/bz224ua4fgDhkbJkvy+iQ3Ah+bZOhszyRXt0M9H0uycc/7/a1xn0G1NSwDDgNe127vS+3y+4Z4kmyU5L1Jft7e3ptko3bZWG1/335uv0hy5GTvTZJtkpyZZHWSFUmObtuPAv4v8MS2js5h29b8+iSXA79OsuG45f8lyU+SHNq1b82cwaHZ8HjgcuBhwKeA04A9gZ2Av6L5a713COow4G3AFjTDQJ+cou8XAe8ANgPGf2HuBXwceC2wObAPcF27+Cbg2cBDgCOBE5LssbYvEPgizR76XhMsexvwVWABsC3wfoCq2qdd/uiqenDPUNdWwEKav8qXTbK9w4BnAjsCjwDeMsl696mqE2ney//dbu85E6z2ZuAJwGOAR7evp7fvrYCHAouBo4APJlkwySZPA1YC2wDPA/5XkqdV1UeBV9DsUTy4qo6ZrvZJHAo8C9i8d4+j/Ry/Aryqqk5dy741AwaH+vWFdo9i7HZ0z7KfVNXH2mGk04HtgOOr6s6q+irwnzQhMubLVXV+Vd1J80X2xCTbTbLdL1bVf1TVvVX1u3HLjgJOqqpz2uU3VNU1AFX15ar6cTXOo/li/7O1ffFVdRdwM80X/nh30YTANlX1u6r61gTr9LoXOKZ9f8aP8Y/5QFVdX1WraYJztv6yPozms7mpqlYBxwEv7ll+V7v8rqo6G7iDCeZf2s9rb+D17Wu+lGYv4yUdannCuP9TPx63/J/a96D3Pfoz4EzgJVV1VodtaRYZHOrXwVW1ec/tIz3Lftlz/7cAVTW+rXeP4/qxO1V1B7Ca5q/WiVw/STs0ATX+ywaAJAckubAdRrkFOJBmD2etJHkgsKitdbzXAQEuSvMLppdO092qCUJwvN7X/VMmf3+62qbtb7K+fzVuPuE3/P5n19vP6qq6fVxfizvUcuG4/1M7jls+0Wf/CuDbVfXNDtvRLDM4NAz37V20Q1gLgcl+STPV6ZuvpxnK+T3tmP0ZwD8CW1bV5sDZNF/ua+sg4G7gD35aWlU3VtXRVbUN8HLgn8fmaSbRzympe/fAtuf+9+fXwKZjC5Js1bHvn9PsHU3Udxc/BxYm2WxcXzesRV+Tmei1vALYPskJs7gddWRwaBgOTPLkJA+imR+4sKqm2rOYzEeBI5Psm+QBSRYn2RV4ELARsAq4O8kBwH5rU2iShUkOAz4I/ENV/WqCdZ6fZNv24RqaL7x728e/BCY9LmUKr0yybZKFNMN5Y/MjlwGPTPKYdsL82HHPm257pwJvSbIozS/G3gr8a9fi2s/r28A7k2yc5FE0Q4ed++rodmB/YJ8k7xrwtjQJg0P9+lL7C5mx2+dn0NengGNohn0eRzOB3llVXUQ78Q3cCpwH7NAOn/wt8GmaL/IX0YyLd3FZkjuAFcDLgL+rqrdOsu6ewHfa9c8EXl1V17bLjgVOacfwX9Bh+5+imZe5lmY47u0AVfVD4Hjga8CPGPeDAZow3a3d3hcm6PftwHKaHzNcAVwy1vdaOBRYQrP38XmaeZsuP9l+4rj/U3ck2XO6J1XVLTTHiByQ5G1rU7hmJl7ISXMpExygJml+cY9DktSJwSFJ6sShKklSJ+5xSJI6WSdPcrjFFlvUkiVLhl2GJM0rF1988c1VtWi69dbJ4FiyZAnLly8fdhmSNK8k+en0azlUJUnqyOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTg0OS1Mk6eQCgJK2NZCYXibzfun4OQINDklrTfeEnWedDoR8OVUmSOjE4JEmdGBySpE4MDknrjYULF5JkrW/AjJ6fhIULFw75XZg5J8clrTfWrFkz9Mnt2frl1jC5xyFJ6sTgkCR14lCVpPVGHfMQOPahw69hnjM4JK03ctxtIzHHUccOtYQZc6hKktSJwSFJ6sTgkCR14hyHpPXKsI+jWLBgwVC3PxsMDknrjZlOjHt23IZDVZKkTgwOSVInDlVJUquf+Y9+1lnXh7MMDklqretf+LPFoSpJUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqROBhYcSU5KclOSKydY9vdJKskW7eMk+ackK5JcnmSPnnUPT/Kj9nb4oOqVJPVnkHscJwP7j29Msh2wH/CznuYDgJ3b2zLgQ+26C4FjgMcDewHHJJn/5ySWpHlsYMFRVecDqydYdALwOqD32P6DgI9X40Jg8yRbA88Ezqmq1VW1BjiHCcJIkjR35nSOI8lBwA1Vddm4RYuB63ser2zbJmufqO9lSZYnWb5q1apZrFqS1GvOgiPJpsCbgLcOov+qOrGqllbV0kWLFg1iE5Ik5naPY0fg4cBlSa4DtgUuSbIVcAOwXc+627Ztk7VLkoZkzoKjqq6oqj+uqiVVtYRm2GmPqroROBN4SfvrqicAt1bVL4CvAPslWdBOiu/XtkmShmSQP8c9FbgA2CXJyiRHTbH62cC1wArgI8B/A6iq1cDbgO+2t+PbNknSkGRdvHDJ0qVLa/ny5cMuQ5LmlSQXV9XS6dbzyHFJUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjrZsJ+VkiwCjgaW9D6nql46mLIkSaOqr+AAvgj8O/A14J7BlSNJGnX9BsemVfX6gVYiSZoX+p3jOCvJgQOtRJI0L/QbHK+mCY/fJbm9vd02yMIkSaOpr6Gqqtps0IVIkuaHfuc4SPJcYJ/24Ter6qzBlCRJGmV9DVUleRfNcNXV7e3VSd45yMIkSaOp3zmOA4FnVNVJVXUSsD/wrKmekOSkJDclubKn7d1JrklyeZLPJ9m8Z9kbk6xI8oMkz+xp379tW5HkDd1eniRptnU5cnzznvsP7WP9k2kCptc5wO5V9Sjgh8AbAZLsBrwQeGT7nH9OskGSDYAPAgcAuwGHtutKkoak3zmOdwLfS3IuEJq5jin/+q+q85MsGdf21Z6HFwLPa+8fBJxWVXcCP0myAtirXbaiqq4FSHJau+7VfdYtSZpl/f6q6tQk3wT2bJteX1U3znDbLwVOb+8vpgmSMSvbNoDrx7U/fqLOkiwDlgFsv/32MyxNkjSZKYeqkuza/rsHsDXNF/dKYJu2ba0keTNwN/DJte1jvKo6saqWVtXSRYsWzVa3kqRxptvj+O80f8X/nwmWFfC0rhtMcgTwbGDfqqq2+QZgu57Vtm3bmKJdkjQEUwZHVS1r7x5QVb/rXZZk464bS7I/8Drgv1bVb3oWnQl8Ksl7gG2AnYGLaOZTdk7ycJrAeCHwoq7blSTNnn5/VfXtPtvuk+RU4AJglyQrkxwFfADYDDgnyaVJPgxQVVcBn6aZ9P434JVVdU9V3Q38DfAV4PvAp9t1JUlDMuUeR5KtaCapN0nyWJo9AICHAJtO9dyqOnSC5o9Osf47gHdM0H42cPZU25IkzZ3p5jieCRxBM7fwnp7224E3DagmSdIIm26O4xTglCR/WVVnzFFNkqQR1u9xHGckeRbNkd0b97QfP6jCJEmjqd+THH4YOAR4Fc08x/OBHQZYlyRpRPX7q6onVdVLgDVVdRzwROARgytLkjSq+g2OsWM4fpNkG+AumiPJJUnrmX5Pcvil9hTo7wYuoTlq/CMDq0qSNLKmDY4kDwC+XlW3AGckOQvYuKpuHXh1kqSRM+1QVVXdS3NNjLHHdxoakrT+6neO4+tJ/jJJpl9VkrQu6zc4Xg58BrgzyW1Jbk9y2wDrkiSNqH4PANxs0IVIkuaHfg8A/Ho/bZKkdd90Z8fdmOYsuFskWcDvnx138aRPlCSts6Ybqno58BqaiytdzP3BcRvNtTUkSeuZ6c6O+z7gfUleVVXvn6OaJEkjrN/J8fcneRKwpPc5VfXxAdUlSRpRfQVHkk8AOwKXAve0zQUYHJK0nun3XFVLgd2qqgZZjCRp9PV7AOCVwFaDLESSND/0u8exBXB1kouAO8caq+q5A6lKkjSy+g2OYwdZhCRp/uj3V1XnJdkB2LmqvpZkU2CDwZYmSRpF/Z5y5Gjgs8C/tE2LgS8MqihJ0ujqd3L8lcDeNEeMU1U/Av54UEVJkkZXv8FxZ1X959iDJBvSHMchSVrP9Bsc5yV5E7BJkmfQXJvjS4MrS5I0qvoNjjcAq4AraE58eDbwlkEVJUkaXf3+HHcT4KSq+ghAkg3att8MqjBJ0mjq+5rjNEExZhPga7NfjiRp1PUbHBtX1R1jD9r7mw6mJEnSKOs3OH6dZI+xB0keB/x2MCVJkkZZv3McrwE+k+TnNFcB3Ao4ZGBVSZJGVr+nHPlukl2BXdqmH1TVXYMrS5I0qvrd4wDYk/uvALhHEq8AKEnrIa8AKEnqxCsASpI68QqAkqROvAKgJKmTgV0BMMlJwLOBm6pq97ZtIXA6zST7dcALqmpNkgDvAw6kOY3JEVV1Sfucw7n/vFhvr6pTutYiSZo9fQ1VVdV5wDXAZu3t+23bVE4G9h/X9gbg61W1M81pTN7Qth8A7NzelgEfgvuC5hjg8cBewDFJFvRTsyRpMPq9AuALgIuA5wMvAL6T5HlTPaeqzgdWj2s+CBjbYzgFOLin/ePVuBDYPMnWwDOBc6pqdVWtAc7hD8NIkjSH+h2qejOwZ1XdBJBkEc1JDj/bcXtbVtUv2vs3Alu29xcD1/est7Jtm6z9DyRZRrO3wvbbb9+xLElSv/r9VdUDxkKj9asOz51Q+9PeWft5b1WdWFVLq2rpokWLZqtbSdI4/X75/1uSryQ5IskRwJdpLubU1S/bISjaf8fC6AZgu571tm3bJmuXJA3JlMGRZKcke1fVa4F/AR7V3i4ATlyL7Z0JHN7ePxz4Yk/7S9J4AnBrO6T1FWC/JAvaSfH92jZJ0pBMN8fxXuCNAFX1OeBzAEn+tF32nMmemORU4CnAFklW0vw66l3Ap5McBfyUZqIdmr2XA4EVND/HPbLd5uokbwO+2653fFWNn3CXJM2h6YJjy6q6YnxjVV2RZMlUT6yqQydZtO8E6xbwykn6OQk4aZo6JUlzZLo5js2nWLbJFMskSeuo6YJjeZKjxzcmeRlw8WBKkiSNsumGql4DfD7JYdwfFEuBBwF/PsjCJEmjacrgqKpfAk9K8lRg97b5y1X1jYFXJkkaSf1eOvZc4NwB1yJJmgdmdPS3JGn9Y3BIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIkjoxOCRJnRgckqRODA5JUicbDrsAaV2SZFb6qapZ6UcaBINDmkX9fOEnMRg0rzlUJUnqxOCQJHVicEiSOjE4JEmdGBySpE4MDklSJwaHJKkTg0PqYOHChSSZ0Q2YcR8LFy4c8juh9ZkHAEodrFmzZiQO3putI9SlteEehySpk6EER5K/S3JVkiuTnJpk4yQPT/KdJCuSnJ7kQe26G7WPV7TLlwyjZklSY86DI8li4G+BpVW1O7AB8ELgH4ATqmonYA1wVPuUo4A1bfsJ7XqSpCEZ1lDVhsAmSTYENgV+ATwN+Gy7/BTg4Pb+Qe1j2uX7xgFeSRqaOQ+OqroB+EfgZzSBcStwMXBLVd3drrYSWNzeXwxc3z737nb9h43vN8myJMuTLF+1atVgX4QkrceGMVS1gGYv4uHANsAfAfvPtN+qOrGqllbV0kWLFs20O0nSJIYxVPV04CdVtaqq7gI+B+wNbN4OXQFsC9zQ3r8B2A6gXf5Q4FdzW7IkacwwjuP4GfCEJJsCvwX2BZYD5wLPA04DDge+2K5/Zvv4gnb5N2oUfkiv9VId8xA49qHDLqOpQxqSDOM7OMlxwCHA3cD3gJfRzGWcBixs2/6qqu5MsjHwCeCxwGrghVV17VT9L126tJYvXz7AV6D11ahcvW9U6tC6JcnFVbV02vXWxf98BocGZVS+sEelDq1b+g0OjxyXJHXiuapG0GwdpuJfpJIGweAYQdN94TtMIWmYHKqSJHVicEiSOjE4hmCmFwMCLwQkaXic4xiCUbgYkOeJXHuj8N4tWLBg2CVoPWZwSB3MRuD74wbNdw5VSZI6MTgkSZ04VDUEo3CiPE+SJ2ltGRxDkONuG/oYdxLq2KGWIGmecqhKktSJwSFJ6sShqiEZ9rEAHgcwGP1+rtOtN+yhTGkqBscQzPRLweMARpefi9YHDlVJkjoxOCRJnRgckqROnOMYQf1MsPazjuPtkgbB4BhBfuFLGmUOVUmSOjE4JEmdGBySpE4MDklSJwaHJKkTg0OS1InBIUnqxOCQJHWSdfFgsySrgJ8Ou44B2gK4edhFaK35+c1f6/pnt0NVLZpupXUyONZ1SZZX1dJh16G14+c3f/nZNRyqkiR1YnBIkjoxOOanE4ddgGbEz2/+8rPDOQ5JUkfucUiSOjE4JEmdGBzzSJKTktyU5Mph16JukmyX5NwkVye5Ksmrh12T+pdk4yQXJbms/fyOG3ZNw+QcxzySZB/gDuDjVbX7sOtR/5JsDWxdVZck2Qy4GDi4qq4ecmnqQ5prNf9RVd2R5IHAt4BXV9WFQy5tKNzjmEeq6nxg9bDrUHdV9YuquqS9fzvwfWDxcKtSv6pxR/vwge1tvf2r2+CQ5liSJcBjge8MtxJ1kWSDJJcCNwHnVNV6+/kZHNIcSvJg4AzgNVV127DrUf+q6p6qegywLbBXkvV2uNjgkOZIOzZ+BvDJqvrcsOvR2qmqW4Bzgf2HXcuwGBzSHGgnVz8KfL+q3jPsetRNkkVJNm/vbwI8A7hmuFUNj8ExjyQ5FbgA2CXJyiRHDbsm9W1v4MXA05Jc2t4OHHZR6tvWwLlJLge+SzPHcdaQaxoaf44rSerEPQ5JUicGhySpE4NDktSJwSFJ6sTgkCR1YnBIk0iyVZLTkvw4ycVJzk7yiBn2+ZQkZ7X3n5vkDe39g5Ps1rPe8UmePrNXIA3GhsMuQBpF7QF7nwdOqaoXtm2PBrYEfjgb26iqM4Ez24cHA2cBV7fL3job25AGwT0OaWJPBe6qqg+PNVTVZcC3krw7yZVJrkhyCNy3J/HNJJ9Nck2ST7bhQ5L927ZLgL8Y6y/JEUk+kORJwHOBd7cHBu6Y5OQkz2vX2zfJ99rtnZRko7b9uiTHJbmkXbbrnL07Wq8ZHNLEdqe5ZsZ4fwE8Bng08HSaL/ut22WPBV4D7Ab8CbB3ko2BjwDPAR4HbDW+w6r6Ns2ex2ur6jFV9eOxZe3zTwYOqao/pRkl+Ouep99cVXsAHwL+x1q/WqkDg0Pq5snAqe2ZUn8JnAfs2S67qKpWVtW9wKXAEmBX4CdV9aNqTtPwrx23t0v7/LHhsVOAfXqWj50s8eJ2e9LAGRzSxK6i2UPo4s6e+/cwN3OIY9ucq+1JBoc0iW8AGyVZNtaQ5FHALcAh7UV9FtH89X/RFP1cAyxJsmP7+NBJ1rsd2GyC9h+0z9+pffximr0caWgMDmkC7bDSnwNPb3+OexXwTuBTwOXAZTTh8rqqunGKfn4HLAO+3E6O3zTJqqcBr20nwXcc9/wjgc8kuQK4F/jwJH1Ic8Kz40qSOnGPQ5LUicEhSerE4JAkdWJwSJI6MTgkSZ0YHJKkTgwOSVIn/x+gmEHMHeaD7wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "plot_erk(initial_results)" ] @@ -209,7 +275,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -225,7 +291,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -242,11 +308,11 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ - "model_posterior = infer_mcmc(initial_experiment_conditioned).run(noise_priors)\n", + "model_posterior = infer_importance(initial_experiment_conditioned).run(noise_priors)\n", "noise_posteriors = {\n", " n: EmpiricalMarginal(model_posterior, sites=n)\n", " for n in noise_vars\n", @@ -262,7 +328,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -294,7 +360,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ @@ -303,9 +369,22 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAHa1JREFUeJzt3XucXGWd5/HPVxDQ8UIyaQWSYBgMuOAoYoGKyuIFCIwaxxswriAqERdccGdVUNdwcVbXy3gDceMSAzOYiNeJiIPgKoyjGDoISBCkBdw0oAmGq2jk8p0/ztNQln2p093VVd35vl+velH1O895zq+qQv36PM+5yDYRERHteky3E4iIiOklhSMiImpJ4YiIiFpSOCIiopYUjoiIqCWFIyIiaknhiJ4i6T5JfzXK8s9L+p8T3MYBkgYn2Md3JB01kT6a+nqxpBuaXt8i6eWT0Xfpb52kAyarvza3KUlflHSnpDWT1Oekfi4xflt3O4HofZJuAZ4KPNQUXmH7+Mnelu0njLH82MneZitJBu4HDGwGrgKW2f5yUx6H1Ohroe2BkdrY/jdg9wkl/ej2VgCDtj/Q1P+ek9F3TS8CDgTm2f5d60JJbwbOBn7fsmg327d1Pr2YiBSOaNcrbV/SzQQkbWX7obFbTopn2x6QNAc4BDhD0jNsnzqZG5G0te0HJ7PPHvE04JbhikaTH9t+0VgdzeDPaNrKUFVMiKQ3S/p3SZ+UdJekmyTtV+LrJW1oHtKRtKIMN10s6V5Jl0p6WtNyS3p6U9uzJF0o6XfAS0rsQ03tF0u6StI9kn4paVGJHy3p52UbN0l6+3jen+07bP8T8A7gZEl/Wfr/gaS3ledPL+/jbkl3SPpyiV9Wurm6DMEdNjRMJum9kn4NfHGEobN9JF1Xhnq+KGm7ps/7hy3fgUsOS4A3Au8p2/tWWf7IEI+kbSV9StJt5fEpSduWZUO5/X353m6XdPRIn42knSStlrRJ0oCkY0r8rcD/BV5Q8qhdbEvO75V0DfA7SVu3LP9Pkm6WdETdvmPiUjhiMjwPuAb4S+BLwCpgH+DpwH+h+mu9eQjqjcDpwByqYaDzRun774B/AJ4ItP5g7gucC7wb2B7YH7ilLN4AvAJ4EnA08ElJe4/3DQL/QrWHvu8wy04HvgvMAuYBnwWwvX9Z/mzbT2ga6toBmE31V/mSEbb3RuBgYFdgN+ADI7R7hO1lVJ/lR8v2XjlMs/cDzwf2Ap5d3k9z3zsATwbmAm8FzpQ0a4RNrgIGgZ2A1wH/S9JLbZ8NHEu1R/EE20vHyn0ERwB/A2zfvMdRvseLgHfaXjnOvmMCUjiiXd8sexRDj2Oalt1s+4tlGOnLwHzgNNubbX8X+CNVERnybduX2d5M9UP2AknzR9juv9j+d9sP2/5Dy7K3AsttX1yW32r7egDb37b9S1cupfphf/F437ztB4A7qH7wWz1AVQR2sv0H2z8cpk2zh4Gl5fNpHeMfcobt9bY3URXOyfrL+o1U380G2xuBU4E3NS1/oCx/wPaFwH0MM/9Svq8XAu8t7/kqqr2MI2vk8vyWf1O/bFn+mfIZNH9GLwZWA0favqDGtmISpXBEu15te/umxxealv2m6fnvAWy3xpr3ONYPPbF9H7CJ6q/W4awfIQ5VgWr9sQFA0iGSLi/DKHcBh1Lt4YyLpMcCfSXXVu8BBKxRdQTTW8bobuMwRbBV8/v+FSN/PnXtVPobqe/ftswn3M+ffnfN/WyyfW9LX3Nr5HJ5y7+pXVuWD/fdHwv8yPYPamwnJlkKR3TDI3sXZQhrNjDSkTSjXb55PdVQzp8oY/ZfAz4OPNX29sCFVD/u47UYeBD4s0NLbf/a9jG2dwLeDnxuaJ5mBO1ckrp5D2xnHv18fgc8fmiBpB1q9n0b1d7RcH3XcRswW9ITW/q6dRx9jWS493IssLOkT07idqKmFI7ohkMlvUjSNlTzA5fbHm3PYiRnA0dLepmkx0iaK+kZwDbAtsBG4EFJhwAHjSdRSbMlvRE4E/jftn87TJvXS5pXXt5J9YP3cHn9G2DE81JGcZykeZJmUw3nDc2PXA3sKWmvMmF+Sst6Y21vJfABSX2qjhj7IPDPdZMr39ePgA9L2k7Ss6iGDmv3VdO9wCJgf0kf6fC2YgQpHNGub5UjZIYe35hAX18CllIN+zyXagK9NttrKBPfwN3ApcDTyvDJfwPOp/oh/zuqcfE6rpZ0HzAAvA14l+0PjtB2H+Anpf1q4ATbN5VlpwDnlDH8N9TY/peo5mVuohqO+xCA7V8ApwGXADfScsAAVTHdo2zvm8P0+yGgn+pghp8BVw71PQ5HAAuo9j6+QTVvU+eQ7Re0/Ju6T9I+Y61k+y6qc0QOkXT6eBKPiVFu5BRTScOcoBYR00v2OCIiopYUjoiIqCVDVRERUUv2OCIiopYZeZHDOXPmeMGCBd1OIyJiWlm7du0dtvvGajcjC8eCBQvo7+/vdhoREdOKpF+N3SpDVRERUVMKR0RE1JLCERERtaRwRERELSkcERFRSwpHRETUksIRERG1pHBEREQtM/IEwIiI8ZAmcpPIR830awCmcEREFGP94Eua8UWhHRmqioiIWlI4IiKilhSOiIioJYUjIrYYs2fPRtK4H8CE1pfE7Nmzu/wpTFwmxyNii3HnnXd2fXJ7so7c6qbscURERC0pHBERUUsKR0RE1NKxOQ5Jy4FXABtsP7PE9gI+D2wHPAj8V9trVA36fRo4FLgfeLPtK8s6RwEfKN1+yPY5nco5ImY2L30SnPLk7ucwzXVycnwFcAZwblPso8Cptr8j6dDy+gDgEGBheTwPOAt4nqTZwFKgARhYK2m17Ts7mHdEzFA69Z6emBz3KV1NYcI6NlRl+zJgU2sYGCq3TwZuK88XA+e6cjmwvaQdgYOBi21vKsXiYmBRp3KOiIixTfXhuCcCF0n6OFXR2q/E5wLrm9oNlthI8YiI6JKpnhx/B/Au2/OBdwFnT1bHkpZI6pfUv3HjxsnqNiJmmImewDfRx6xZs7r9EUzYVBeOo4Cvl+dfAfYtz28F5je1m1diI8X/jO1lthu2G319fZOadETMDLYn9JiMPjZtah3Bn36munDcBvzn8vylwI3l+WrgSFWeD9xt+3bgIuAgSbMkzQIOKrGIiOiSTh6Ou5LqiKk5kgapjo46Bvi0pK2BPwBLSvMLqQ7FHaA6HPdoANubJJ0OXFHanWZ7+pfriIhpTN0+NK0TGo2G+/v7u51GRMwwM/1GTpLW2m6M1S4XOYyIKNq5AGE7bWZycYEUjoiIR8z0H/zJkmtVRURELSkcERFRSwpHRETUksIRERG1pHBEREQtKRwREVFLCkdERNSSwhEREbWkcERERC0pHBERUUsKR0RE1JLCERERtaRwRERELSkcERFRS8cKh6TlkjZIurYp9mVJV5XHLZKualp2sqQBSTdIOrgpvqjEBiSd1Kl8IyKiPZ28H8cK4Azg3KGA7cOGnkv6BHB3eb4HcDiwJ7ATcImk3UrTM4EDgUHgCkmrbV/XwbwjImIUHSscti+TtGC4ZapuofUG4KUltBhYZXszcLOkAWDfsmzA9k1lvVWlbQpHRESXdGuO48XAb2zfWF7PBdY3LR8ssZHif0bSEkn9kvo3btzYgZQjIgK6VziOAFZOZoe2l9lu2G709fVNZtcREdFkyu85Lmlr4DXAc5vCtwLzm17PKzFGiUdERBd0Y4/j5cD1tgebYquBwyVtK2kXYCGwBrgCWChpF0nbUE2gr57yjCMi4hGdPBx3JfBjYHdJg5LeWhYdTsswle11wPlUk97/Chxn+yHbDwLHAxcBPwfOL20jIqJLZLvbOUy6RqPh/v7+bqcRETGtSFpruzFWu5w5HhERtaRwRERELSkcERFRSwpHRETUksIRERG1pHBEREQtKRwREVFLCkdERNSSwhEREbWkcERERC0pHBERUUsKR0RE1JLCERERtaRwRERELSkcERFRSwpHRETU0sk7AC6XtEHStS3xd0q6XtI6SR9tip8saUDSDZIOboovKrEBSSd1Kt+IiGjP1h3sewVwBnDuUEDSS4DFwLNtb5b0lBLfg+qWsnsCOwGXSNqtrHYmcCAwCFwhabXt6zqYd0REjKJjhcP2ZZIWtITfAXzE9ubSZkOJLwZWlfjNkgaAfcuyAds3AUhaVdqmcEREdMlUz3HsBrxY0k8kXSppnxKfC6xvajdYYiPF/4ykJZL6JfVv3LixA6lHRARMfeHYGpgNPB94N3C+JE1Gx7aX2W7YbvT19U1GlxERMYxOznEMZxD4um0DayQ9DMwBbgXmN7WbV2KMEo+IiC6Y6j2ObwIvASiT39sAdwCrgcMlbStpF2AhsAa4AlgoaRdJ21BNoK+e4pwjIqJJW3sckvqAY4AFzevYfsso66wEDgDmSBoElgLLgeXlEN0/AkeVvY91ks6nmvR+EDjO9kOln+OBi4CtgOW219V8jxERMYlU/W6P0Uj6EfBvwFrgoaG47a91LrXxazQa7u/v73YaERHTiqS1thtjtWt3juPxtt87wZwiImIGaHeO4wJJh3Y0k4iImBbaLRwnUBWPP0i6tzzu6WRiERHRm9oaqrL9xE4nEhER00Pb53FIehWwf3n5A9sXdCaliIjoZW0NVUn6CNVw1XXlcYKkD3cysYiI6E3t7nEcCuxl+2EASecAPwVO7lRiERHRm+qcOb590/MnT3YiERExPbS7x/Fh4KeSvg+Iaq4jN1WKiNgCtXtU1UpJPwCGLoP+Xtu/7lhWERHRs0YdqpL0jPLfvYEdqa5uOwjsVGIREbGFGWuP478DS4BPDLPMwEsnPaOIiOhpoxYO20vK00Ns/6F5maTtOpZVRET0rHaPqvpRm7GIiJjhRt3jkLQD1T2+HyfpOVRHVAE8CXh8h3OLiIgeNNYex8HAx6lu2fqPVHMdn6Ca+3jfaCtKWi5pQ7lp01DsFEm3SrqqPA5tWnaypAFJN0g6uCm+qMQGJOUQ4IiILhtrjuMc4BxJrx3HTZtWAGcA57bEP2n7480BSXtQ3RZ2T2An4JJya1mAM4EDqY7mukLSatvX1cwlIiImSbvncXxN0t9Q/bBv1xQ/bZR1LpO0oM08FgOrbG8GbpY0AOxblg3YvglA0qrSNoUjIqJL2r3I4eeBw4B3Us1zvB542ji3ebyka8pQ1qwSmwusb2ozWGIjxSMiokvaPapqP9tHAnfaPhV4AbDbGOsM5yxgV2Av4HaGPz9kXCQtkdQvqX/jxo2T1W1ERLRot3AMncNxv6SdgAeoziSvxfZvbD9UrrL7BR4djroVmN/UdF6JjRQfru9lthu2G319fXVTi4iINrVbOL4laXvgY8CVwC3Al+puTFJzsflbYOiIq9XA4ZK2lbQLsBBYA1wBLJS0i6RtqCbQV9fdbkRETJ4xJ8clPQb4nu27gK9JugDYzvbdY6y3EjgAmCNpEFgKHCBpL6rLldwCvB3A9jpJ51NNej8IHGf7odLP8cBFwFbActvrxvNGIyJicsj22I2kn9p+zhTkMykajYb7+/u7nUZExLQiaa3txljt2h2q+p6k10rS2E0jImIma7dwvB34CrBZ0j2S7pV0TwfzioiIHtXuCYBP7HQiERExPbR7AuD32olFRMTMN9bVcbejugrunHKWd/PVcXMGd0TEFmisoaq3AydSXXhwLY8WjnuoLmAYERFbmLGujvtp4NOS3mn7s1OUU0RE9LB2J8c/K2k/YEHzOrZbL5keEREzXFuFQ9I/UV2c8CrgoRI2f36vjYiImOHaKhxAA9jD7ZxmHhERM1q7JwBeC+zQyUQiImJ6aHePYw5wnaQ1wOahoO1XdSSriIjoWe0WjlM6mUREREwf7R5VdamkpwELbV8i6fFUlzmPiIgtTLuXHDkG+Crwf0poLvDNTiUVERG9q93J8eOAF1KdMY7tG4GndCqpiIjoXe0Wjs22/zj0QtLWVOdxjEjSckkbJF07zLK/l2RJc8prSfqMpAFJ10jau6ntUZJuLI+j2sw3IiI6pN3Ccamk9wGPk3Qg1b05vjXGOiuARa1BSfOBg4D/3xQ+hOo+4wuBJcBZpe1sqlvOPg/YF1haLrYYERFd0m7hOAnYCPyM6sKHFwIfGG0F25cBm4ZZ9EngPfzpHsti4FxXLge2l7QjcDBwse1Ntu8ELmaYYhQREVOn3cNxHwcst/0FAElbldj9dTYmaTFwq+2rW+5COxdY3/R6sMRGig/X9xKqvRV23nnnOmlFREQNbd9znKpQDHkccEmdDZVDeN8HfLDOeu2yvcx2w3ajr6+vE5uIiAjaLxzb2b5v6EV5/via29oV2AW4WtItwDzgSkk7ALcC85vaziuxkeIREdEl7RaO37Uc6fRc4Pd1NmT7Z7afYnuB7QVUw0572/41sBo4shxd9Xzgbtu3AxcBB0maVSbFDyqxiIjoknbnOE4EviLpNqq7AO4AHDbaCpJWAgdQ3XZ2EFhq++wRml8IHAoMUM2bHA1ge5Ok04ErSrvTbA834R4REVNE7V4pXdJjgd3LyxtsP9CxrCao0Wi4v7+/22lEREwrktbabozVrt09DoB9ePQOgHtLyh0AIyK2QLkDYERE1JI7AEZERC25A2BERNSSOwBGREQtuQNgRETUUucOgE+lOrIKYI3tDZ1LKyIielW7dwB8A7AGeD3wBuAnkl7XycQiIqI3tTtU9X5gn6G9DEl9VBc5/GqnEouIiN7U7lFVj2kZmvptjXUjImIGaXeP418lXQSsLK8Po7q+VEREbGFGLRySng481fa7Jb0GeFFZ9GPgvE4nFxERvWesPY5PAScD2P468HUASX9dlr2yo9lFRETPGWue4qm2f9YaLLEFHckoIiJ62liFY/tRlj1ulGURETFDjVU4+iUd0xqU9DZgbWdSioiIXjZW4TgROFrSDyR9ojwuBd4KnDDaipKWS9og6dqm2OmSrpF0laTvStqpxCXpM5IGyvLm29QeJenG8jhq/G81IiImw6iFw/ZvbO8HnArcUh6n2n5BuVf4aFYAi1piH7P9LNt7ARcAHyzxQ4CF5bEEOAtA0mxgKfA8YF9gabn3eEREdEm716r6PvD9Oh3bvkzSgpbYPU0v/4LqZlAAi4Fzy/0+Lpe0vaQdqe5ZfvHQfcYlXUxVjFYSERFdUefWsZNC0j8ARwJ3Ay8p4bnA+qZmgyU2Uny4fpdQ7a2w8847T27SERHxiCm/bIjt99ueT3UC4fGT2O8y2w3bjb6+vsnqNiIiWnTzelPnAa8tz28F5jctm1diI8UjIqJLprRwSFrY9HIxcH15vho4shxd9Xzgbtu3AxcBB0maVSbFDyqxiIjoko7NcUhaSTW5PUfSINXRUYdK2h14GPgVcGxpfiFwKDAA3A8cDWB7k6TTgStKu9OGJsojIqI7VB3INLM0Gg339/d3O42IiGlF0lrbjbHa5Z4aERFRSwpHRETUksIRERG1pHBEREQtKRwREVFLCkdERNSSwhEREbWkcERERC0pHBERUUsKR0RE1JLCERERtaRwRERELSkcERFRSwpHRETUksIRERG1dKxwSFouaYOka5tiH5N0vaRrJH1D0vZNy06WNCDpBkkHN8UXldiApJM6lW9ERLSnk3scK4BFLbGLgWfafhbwC+BkAEl7AIcDe5Z1PidpK0lbAWcChwB7AEeUthER0SUdKxy2LwM2tcS+a/vB8vJyYF55vhhYZXuz7ZupbiG7b3kM2L7J9h+BVaVtRER0STfnON4CfKc8nwusb1o2WGIjxSMioku6UjgkvR94EDhvEvtcIqlfUv/GjRsnq9uIiGix9VRvUNKbgVcAL7PtEr4VmN/UbF6JMUr8T9heBiwDaDQaHq5NRKdJmpR+Hv1fI6L3TOkeh6RFwHuAV9m+v2nRauBwSdtK2gVYCKwBrgAWStpF0jZUE+irpzLniDpsj/lop11EL+vYHoeklcABwBxJg8BSqqOotgUuLn+ZXW77WNvrJJ0PXEc1hHWc7YdKP8cDFwFbActtr+tUzhERMTbNxL9uGo2G+/v7u51GxLAkZa8iepKktbYbY7XLmeMREVFLCkdERNSSwhFRw+zZs5E0oQcw4T5mz57d5U8itmRTfjhuxHR255139sT8xGQd9hsxHtnjiIiIWlI4IiKilhSOiIioJYUjIiJqyeR4D8r1jiKil6Vw9KCxfvBz5nH3eOmT4JQndzuNKo+ILknhiKhBp97TE0VbEj6l21nElipzHBERUUsKR0RE1JLCERERtaRwdMFEr3cEudZRRHRPJse7oBeud5RrHY1fL3x2s2bN6nYKsQXr2B6HpOWSNki6tin2eknrJD0sqdHS/mRJA5JukHRwU3xRiQ1IOqlT+Ua0o51bw07GrWPHemzatKnLn0RsyTo5VLUCWNQSuxZ4DXBZc1DSHlT3E9+zrPM5SVtJ2go4EzgE2AM4orSNiIgu6dhQle3LJC1oif0cht3VXwyssr0ZuFnSALBvWTZg+6ay3qrS9rpO5R0REaPrlcnxucD6pteDJTZS/M9IWiKpX1L/xo0bO5ZoRMSWbsZMjtteBiwDaDQa3T+1dxS9cNmKXLIiIsarVwrHrcD8ptfzSoxR4tNWL1y2Ipes6Ix2j7gaq123/31EjKZXhqpWA4dL2lbSLsBCYA1wBbBQ0i6StqGaQF/dxTwjRjUZR12laESv69geh6SVwAHAHEmDwFJgE/BZoA/4tqSrbB9se52k86kmvR8EjrP9UOnneOAiYCtgue11nco5IiLGppn4102j0XB/f3+30xhRL1wWvRdyiIjeImmt7cZY7XpljmOL0+2zj3PmcUSMVwpHF0z0L/3sLUREN/XK5HhEREwTKRwREVFLCkdERNSSwhEREbVkcrwHtXPEVTttMoEeEZ2QwtGD8oMfEb0sQ1UREVFLCkdERNSSwhEREbWkcERERC0pHBERUUsKR0RE1JLCERERtaRwRERELTPyRk6SNgK/6nYeHTQHuKPbScS45fubvmb6d/c0231jNZqRhWOmk9Tfzl26ojfl+5u+8t1VMlQVERG1pHBEREQtKRzT07JuJxATku9v+sp3R+Y4IiKipuxxRERELSkcERFRSwrHNCJpuaQNkq7tdi5Rj6T5kr4v6TpJ6ySd0O2con2StpO0RtLV5fs7tds5dVPmOKYRSfsD9wHn2n5mt/OJ9knaEdjR9pWSngisBV5t+7oupxZtUHWv5r+wfZ+kxwI/BE6wfXmXU+uK7HFMI7YvAzZ1O4+oz/bttq8sz+8Ffg7M7W5W0S5X7isvH1seW+xf3SkcEVNM0gLgOcBPuptJ1CFpK0lXARuAi21vsd9fCkfEFJL0BOBrwIm27+l2PtE+2w/Z3guYB+wraYsdLk7hiJgiZWz8a8B5tr/e7XxifGzfBXwfWNTtXLolhSNiCpTJ1bOBn9v+x27nE/VI6pO0fXn+OOBA4PruZtU9KRzTiKSVwI+B3SUNSnprt3OKtr0QeBPwUklXlceh3U4q2rYj8H1J1wBXUM1xXNDlnLomh+NGREQt2eOIiIhaUjgiIqKWFI6IiKglhSMiImpJ4YiIiFpSOCJGIGkHSask/VLSWkkXStptgn0eIOmC8vxVkk4qz18taY+mdqdJevnE3kFEZ2zd7QQielE5Ye8bwDm2Dy+xZwNPBX4xGduwvRpYXV6+GrgAuK4s++BkbCOiE7LHETG8lwAP2P78UMD21cAPJX1M0rWSfibpMHhkT+IHkr4q6XpJ55Xig6RFJXYl8Jqh/iS9WdIZkvYDXgV8rJwYuKukFZJeV9q9TNJPy/aWS9q2xG+RdKqkK8uyZ0zZpxNbtBSOiOE9k+qeGa1eA+wFPBt4OdWP/Y5l2XOAE4E9gL8CXihpO+ALwCuB5wI7tHZo+0dUex7vtr2X7V8OLSvrrwAOs/3XVKME72ha/Q7bewNnAf9j3O82ooYUjoh6XgSsLFdK/Q1wKbBPWbbG9qDth4GrgAXAM4Cbbd/o6jIN/1xze7uX9YeGx84B9m9aPnSxxLVlexEdl8IRMbx1VHsIdWxuev4QUzOHOLTNqdpeRApHxAj+H7CtpCVDAUnPAu4CDis39emj+ut/zSj9XA8skLRreX3ECO3uBZ44TPyGsv7Ty+s3Ue3lRHRNCkfEMMqw0t8CLy+H464DPgx8CbgGuJqquLzH9q9H6ecPwBLg22VyfMMITVcB7y6T4Lu2rH808BVJPwMeBj4/Qh8RUyJXx42IiFqyxxEREbWkcERERC0pHBERUUsKR0RE1JLCERERtaRwRERELSkcERFRy38AcmUpW5VDDxsAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "plot_erk(new_results)" ]