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) diff --git a/causal_demon/transmitters.py b/causal_demon/transmitters.py index d84eb2c..94dc271 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., @@ -152,26 +153,29 @@ 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. + """ - 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) + + 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 f4db6c9..fa1dc7c 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,140 +32,289 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "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", + "import warnings\n", + "warnings.filterwarnings('ignore')\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", - "\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(\"#\")" + "%matplotlib inline\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Each noise term is modeled with a weakly informed prior." + "### 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": null, + "execution_count": 2, "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}" + "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_importance(model):\n", + " return pyro.infer.Importance(model, num_samples=1000)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "cancer_signaling(noise_prior)" + "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()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 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." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "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)" + "noise_vars = ['N_egf', 'N_igf', 'N_sos', 'N_ras', 'N_pi3k', 'N_akt', 'N_raf', 'N_mek', 'N_erk']\n", + "# 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, + "metadata": {}, + "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)" ] }, { "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", + "Prior to running the experiment, we can simulate the experiment using this model. \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." + "**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([[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", + "conditions = {\n", + " 'egf': torch.tensor(\n", + " [\n", + " [800., 800.],\n", + " [8000., 8000.],\n", + " [80000., 80000.]\n", + " ]\n", + " ),\n", + " 'igf': torch.zeros([3, 2])\n", + "}\n", + "initial_experiment = do(cancer_signaling, data=conditions)\n", + "initial_results = initial_experiment(noise_priors)\n", + "pprint(initial_results)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "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)" ] }, { "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 800?*\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": 8, "metadata": {}, "outputs": [], "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" ] }, { "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": 9, "metadata": {}, "outputs": [], "source": [ - "cancer_dist = infer_dist(cancer_obs, noise_prior)" + "evidence = {'erk': initial_results['erk']}\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": 10, "metadata": {}, "outputs": [], "source": [ - "noise_marginals = {\n", - " n: EmpiricalMarginal(cancer_dist, sites=n)\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", "}" ] @@ -174,42 +323,70 @@ "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" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ - "cancer_do = do(cancer_signaling, data={'igf': tensor(800.)})" + "conditions = {\n", + " 'egf': torch.tensor(\n", + " [\n", + " [800., 800.],\n", + " [8000., 8000.],\n", + " [80000., 80000.],\n", + " ]\n", + " ),\n", + " 'igf': torch.tensor(\n", + " [\n", + " [800., 800.],\n", + " [800., 800.],\n", + " [800., 800.],\n", + " ]\n", + " ),\n", + "}\n", + "counterfactual_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." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ - "cancer_do_dist = infer_dist(cancer_do, noise_marginals)\n", - "erk_cf_marginal = EmpiricalMarginal(cancer_do_dist, sites = 'erk')" + "new_results = counterfactual_experiment(noise_posteriors)" ] }, { "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": [ - "hist(erk_cf_marginal, 'Erk')" + "plot_erk(new_results)" ] }, { @@ -218,7 +395,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" ]