{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Inheritance and body heights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Halerium models consist of Variables, Entities and Graphs. They are the essential [building blocks](./01_building_blocks.ipynb) to create arbitrarily hierarchical models.\n", "\n", "We want to illustrate how this works by implementing a simple inheritance model, in which we look at the dependencies between the height of parents and the height of children." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### We start by importing required packages:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import warnings; warnings.simplefilter('ignore')\n", "import halerium.core as hal\n", "from plots import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Let us create our first entity, a person:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "person = hal.Entity(\"person\")\n", "with person:\n", " hal.Variable(\"height\", shape=[])\n", " hal.Variable(\"genetic_factor\", shape=[], mean=0, variance=75)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have an Entity instance called person. In our simple world a person is fully defined by its genetic and actual height. \n", "We can create a Template from this instance for future usage." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "Person = person.get_template()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Now let's get more specific and create a man and a woman:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "man = Person(\"man\")\n", "with man:\n", " height.mean = 175 + genetic_factor\n", " height.variance = 25.\n", "\n", "woman = Person(\"woman\")\n", "with woman:\n", " height.mean = 167 + genetic_factor\n", " height.variance = 25.\n", " \n", "Man = man.get_template()\n", "Woman = woman.get_template()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, so what happened? We defined both man and woman as a person, but then we modified each of them in a different way.\n", "Essentially we are saying, the typical man is 175cm tall and the typical woman 167cm tall. However, the actual height is modified by genetic factors of this particular person as well as a random noise around the genetic predisposition.\n", "\n", "We also created templates of man and woman." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Now, let's define a Graph in which a man and a woman get a daughter." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "make_daughter = hal.Graph(\"make_daughter\")\n", "with make_daughter:\n", " with inputs:\n", " Man(\"father\")\n", " Woman(\"mother\")\n", " with outputs:\n", " Woman(\"daughter\")\n", " outputs.daughter.genetic_factor.mean = (inputs.father.genetic_factor + inputs.mother.genetic_factor) / 2.\n", " outputs.daughter.genetic_factor.variance = 37.5\n", "\n", "MakeDaughter = make_daughter.get_template()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So what we specified is the following: the Graph make_daughter has two inputs: mother and father. mother is an instace of Woman, father is an instance of Man. The output is a daughter, an instance of Woman.\n", "\n", "If we left it at that there would have been no dependence of the daughter height on her parents. We specified this dependency in the the two last lines of code. We specified that the genetic_height of the daughter fluctuates around the mean genetic height of her parents (a super simple model of combination and mutation).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's repeat that for a son:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "make_son = hal.Graph(\"make_son\")\n", "with make_son:\n", " with inputs:\n", " Man(\"father\")\n", " Woman(\"mother\")\n", " with outputs:\n", " Man(\"son\")\n", " outputs.son.genetic_factor.mean = (inputs.father.genetic_factor + inputs.mother.genetic_factor) / 2.\n", " outputs.son.genetic_factor.variance = 37.5\n", "\n", "MakeSon = make_son.get_template()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Now we can create the Graph for a full family." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's call the parents alice and bob" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "family = hal.Graph(\"family\")\n", "with family:\n", " Woman(\"alice\")\n", " Man(\"bob\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "alice and bob are happy together so let's make two children, harry and sally" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "with family:\n", " MakeSon(\"make_harry\")\n", " MakeDaughter(\"make_sally\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, we have not specified yet, that alice and bob are actually the parents of both of them.\n", "We can do that by linking:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "with family:\n", " hal.link(alice, make_harry.inputs.mother)\n", " hal.link(alice, make_sally.inputs.mother)\n", " hal.link(bob, make_harry.inputs.father)\n", " hal.link(bob, make_sally.inputs.father)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we specified that the father in make_harry was actually bob and not somebody else." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can inspect the family graph in the interactive GUI" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "hal.show(family)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Now we are done and can start to generate random family examples:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start by creating a generative model." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "gen_model = hal.get_generative_model(family)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The generative model is the simplest kind of model, that allows us to create artificial data and/or do forward simulations.\n", "For example we can get example samples:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'alice': [array([155.97561197])],\n", " 'bob': [array([185.9337903])],\n", " 'harry': [array([156.81862308])],\n", " 'sally': [array([159.64632049])]}" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAFuCAYAAABwX0lLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgz0lEQVR4nO3dfbRcdX3v8fdXsHh9am0JrQIpaAGFUunllD54W+2DhdZbo7bacG8tVGtYXtqKtb1LbBW0zartaku97cIWBbe50RBKY0oaG2M0N5CmgSQSQ8IhcAyRJBCTQB4JCcnJ9/6xd06G48nDJOfMzI+8X2vNynx/e8+e7znzcD757b1nIjORJEkqzQu63YAkSdLxMMRIkqQiGWIkSVKRDDGSJKlIhhhJklQkQ4wkSSpSR0NMRJwdEfMjoj8iVkXEB5rxdzb1gYjoG3ab6yNiICJWR8TlnexXkiT1rlM7fH/7gQ9l5jci4mXAsoj4KrASeAfwT60rR8SFwETgIuBVwLyIOD8zBzvctyRJ6jEdnYnJzCcy8xvN9Z1AP3BmZvZn5uoRbjIBuD0z92bmo8AAcFnnOpYkSb2qa8fERMQ5wI8D9x5htTOBdS31+mZMkiSd5Dq9OwmAiHgp8C/AdZm540irjjD2Xd+TEBGTgEkAL3nJSy597WtfOyp9SpKk7lu2bNmWzBw3fLzjISYiXkgdYL6QmTOOsvp64OyW+izg8eErZeYtwC0AfX19uXTp0lHqVpIkdVtEfHuk8U6fnRTArUB/Zv7tMdzkLmBiRJwWEecC5wH3jWWPkiSpDJ2eiXkD8G7ggYhY3ox9BDgN+HtgHDA7IpZn5uWZuSoi7gAepD6z6VrPTJIkSdDhEJOZCxn5OBeALx3mNpOByWPWlCRJKpKf2CtJkopkiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkQwxkiSpSIYYSZJUJEOMJEkqkiFGkiQVyRAjSZKKZIiRJElFMsRIkqQiGWIkSVKRDDGSJKlIhhhJklQkQ4wkSSqSIUaSJBXJECNJkopkiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSitTREBMRZ0fE/Ijoj4hVEfGBZvz7I+KrEfFI8+8rWm5zfUQMRMTqiLi8k/1KkqTe1emZmP3AhzLzdcBPAddGxIXAh4GvZeZ5wNeammbZROAi4Arg5og4pcM9S5KkHtTREJOZT2TmN5rrO4F+4ExgAvD5ZrXPA29rrk8Abs/MvZn5KDAAXNbJniVJUm/q2jExEXEO8OPAvcAPZuYTUAcd4IxmtTOBdS03W9+MDd/WpIhYGhFLN2/ePKZ9S5Kk3tCVEBMRLwX+BbguM3ccadURxvK7BjJvycy+zOwbN27caLUpSZJ6WMdDTES8kDrAfCEzZzTD34mIVzbLXwlsasbXA2e33Pws4PFO9SpJknpXp89OCuBWoD8z/7Zl0V3AVc31q4B/bRmfGBGnRcS5wHnAfZ3qV5Ik9a5TO3x/bwDeDTwQEcubsY8AnwTuiIj3Ao8B7wTIzFURcQfwIPWZTddm5mCHe5YkST2ooyEmMxcy8nEuAL94mNtMBiaPWVOSJKlIfmKvJEkqkiFGkiQVyRAjSZKKZIiRJElFMsRIkqQiGWIkSVKRDDGSJKlIhhhJklQkQ4wkSSqSIUaSJBXJECNJkopkiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkQwxkiSpSIYYSZJUJEOMJEkqkiFGkiQVyRAjSZKKZIiRJElFMsRIkqQiGWIkSVKRDDGSJKlIhhhJklQkQ4wkSSpSR0NMRNwWEZsiYmXL2Osj4j8j4oGImBURL29Zdn1EDETE6oi4vJO9SpKk3tbpmZgKuGLY2GeBD2fmxcCXgD8GiIgLgYnARc1tbo6IUzrXqiRJ6mUdDTGZeTfw1LDhC4C7m+tfBX69uT4BuD0z92bmo8AAcFlHGpUkST2vF46JWQm8tbn+TuDs5vqZwLqW9dY3Y98lIiZFxNKIWLp58+Yxa1SSJPWOXggx7wGujYhlwMuAZ5vxGGHdHGkDmXlLZvZlZt+4cePGqE1JktRLTu12A5n5EPDLABFxPvCWZtF6Ds3KAJwFPN7Z7iRJUq/q+kxMRJzR/PsC4E+Bf2wW3QVMjIjTIuJc4Dzgvu50KUmSek1HZ2IiYhrwJuD0iFgP3AC8NCKubVaZAXwOIDNXRcQdwIPAfuDazBzsZL+SJKl3ReaIh5kUq6+vL5cuXdrtNiRJ0iiJiGWZ2Td8vOu7kyRJko6HIUaSJBXJECNJkopkiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSitTRL4CUpG66+PMXd7sFtXjgqge63YIK50yMJEkqkiFGkiQVyRAjSZKKZIiRJElFMsRIkqQiGWIkSVKRDDGSJKlIhhhJklQkQ4wkSSqSIUaSJBXJECNJkopkiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkInU0xETEbRGxKSJWtoxdEhGLI2J5RCyNiMtall0fEQMRsToiLu9kr5Ikqbd1eiamAq4YNvZXwMcz8xLgY01NRFwITAQuam5zc0Sc0rFOJUlST+toiMnMu4Gnhg8DL2+ufy/weHN9AnB7Zu7NzEeBAeAyJEmSgFO73QBwHfCViPhr6lD1M834mcDilvXWN2PfJSImAZMAxo8fP2aNSpKk3tELB/a+H/hgZp4NfBC4tRmPEdbNkTaQmbdkZl9m9o0bN26M2pQkSb2kF0LMVcCM5vo/c2iX0Xrg7Jb1zuLQriZJknSS64UQ8zjwxub6LwCPNNfvAiZGxGkRcS5wHnBfF/qTJEk9qKPHxETENOBNwOkRsR64AXgf8KmIOBXYQ3NsS2auiog7gAeB/cC1mTnYyX4lSVLv6miIycwrD7Po0sOsPxmYPHYdSZKkUvXC7iRJkqS2GWIkSVKRDDGSJKlIhhhJklQkQ4wkSSqSIUaSJBXJECNJkopkiJEkSUUyxEiSpCIZYiRJUpHaCjERMRgRlx1m2aUR4XcbSZKkjmh3JiaOsOwUIE+gF0mSpGN2TF8AGREv4FCAeUFTt/ovwK8AW0axN0mSpMM6aoiJiBuAjzVlAv9xhNVvHo2mJEmSjuZYdif9P+ATwJ9Rz8bc1tStlz8B3g58YEy6bMOTTz7J8uXLARgcHKSqKlasWAHAvn37qKqKlStXArBnzx6qqqK/vx+A3bt3U1UVq1evBmDXrl1UVcXAwAAA27dvp6oq1qxZA8DWrVupqoq1a9cCsGXLFqqqYt26dQBs2rSJqqrYsGEDABs3bqSqKjZu3AjAhg0bqKqKTZs2AbBu3TqqqmLLlnpCa+3atVRVxdatWwFYs2YNVVWxfft2AAYGBqiqil27dgGwevVqqqpi9+7dAPT391NVFXv27AFg5cqVVFXFvn37AFixYgVVVTE4WB/KtHz5cqqqGvpdLlu2jClTpgzVS5YsYerUqUP14sWLmTZt2lC9aNEipk+fPlQvXLiQO++8c6hesGABM2bMGKrnz5/PzJkzh+p58+Yxa9asoXru3LnMnj17qJ4zZw5z5swZqmfPns3cuXOH6lmzZjFv3ryheubMmcyfP3+onjFjBgsWLBiq77zzThYuXDhUT58+nUWLFg3V06ZNY/HixUP11KlTWbJkyVA9ZcoUli1bNlRXVeVzr5Dn3q7+XexYtuNQvWoXO+4/VO9cuZOdy3ceqlfsZOeKlnr5TnauPFTvuH8Hu1btOlQv28Gu/kP19qXbefqhpw/VS7bz9MMt9b3b2f3I7qF62+Jt7B5oqRdtY/eaQ/XW/9jKM2ufOVTfs5Vnvl3XeSDZes9W9qyrf/e5v6nX1/WBZw/U9eNNvbeu9z6xF4DBPYN1/Z2m3l3Xz256tq6fburNdb1/5/66frKpd9T1vq31Y71v27663tbUW+t6/479wMn33PN97/jf90Zy1JmYzFwALACIiAQ+k5mPH+12kiRJYykyn1/H4vb19eXSpUu73YakHnTx5y/udgtq8cBVD3S7BRUiIpZlZt/w8WM6sHfYht4IXAmMB140bHFm5i8eX4uSJEnHrq0QExHXAJ8GngQeAfYOX2WU+pIkSTqidmdiPgR8EXhPZj47Bv1IkiQdk3Y/7O5M4HMGGEmS1G3thphlwKvHohFJkqR2tBti/gC4LiJ+biyakSRJOlbH8om963judyJ9LzA/InYDW4etnpn5w6PYnyRJ0oiO5cDer+EXO0qSpB5zLJ/Ye3UH+pAkSWpLu8fESJIk9YR2P+zut4+w+ACwHbg/M9efUFeSJElH0e6H3VUcOj6m9dN5W8cORMR04Hf8PBlJkjRW2t2d9Abg28A/AG8EXtv8ezPwGPAW4Hrg7cCNo9alJEnSMO3OxPwRcHtmfqRl7GHgnojYCUzKzLdHxMuB/wl8ZKSNSJIknah2Z2LeTH3K9Ui+Dhz8Buu7qb+iQJIkaUy0G2KeBS49zLJLm+UHt/v08BUi4raI2BQRK1vGpkfE8uayNiKWtyy7PiIGImJ1RFzeZq+SJOl5rN3dSf8MfDwiBoE7gU3AGcA7qY+Bua1Z7xJg9Qi3r6iPp5lycCAzf/Pg9Yj4G+oznIiIC4GJwEXAq4B5EXF+Zg622bMkSXoeajfE/CHwMuCvmkurLwIfaq6vBP5z+I0z8+6IOGekDUdEAO8CfqEZmkB9/M1e4NGIGAAuG2m7kiTp5NNWiMnMZ4DfiohPAD8JvBJ4Arg3Mx9uWW/2cfTys8B3MvORpj4TWNyyfD0eZyNJkhrtzsQA0ASWh4+6YnuuBKa11DHCOiN+h1NETAImAYwfP36U25IkSb3oWL7FejzwRGbua64fUWY+1m4TEXEq8A6ee9DweuDslvos4PHD3OctwC0AfX19flmlJEkngWOZiXkU+GngPmAtR/9G61OOo49fAh4a9nUFdwFfjIi/pT6w97ymB0mSpGMKMe8BvtVy/bhnOiJiGvAm4PSIWA/ckJm3Up+F1LoricxcFRF3AA8C+4FrPTNJkiQddNQQk5mfb7lencidZeaVhxm/+jDjk4HJJ3KfkiTp+andD7sDICJeEBE/GhFvjIiXjHZTkiRJR9N2iImIa4GNwArqrxq4oBmfGRF/MLrtSZIkjaytEBMR7wM+Bcyk/mC61tOg7wF+fdQ6kyRJOoJ2Z2L+EPibzJwEfGnYsodoZmUkSZLGWrsh5lzgK4dZ9jTwfSfUjSRJ0jFqN8RsAc45zLILgA0n1I0kSdIxajfEzAI+FhGvbhnLiDgd+CD1sTJd9eSTT7J8+XIABgcHqaqKFStWALBv3z6qqmLlypUA7Nmzh6qq6O/vB2D37t1UVcXq1fUXcO/atYuqqhgYGABg+/btVFXFmjVrANi6dStVVbF27VoAtmzZQlVVrFu3DoBNmzZRVRUbNtTZbuPGjVRVxcaNGwHYsGEDVVWxadMmANatW0dVVWzZsgWAtWvXUlUVW7duBWDNmjVUVcX27dsBGBgYoKoqdu3aBcDq1aupqordu3cD0N/fT1VV7NmzB4CVK1dSVRX79u0DYMWKFVRVxeBg/fE7y5cvp6qqod/lsmXLmDJl6AvHWbJkCVOnTh2qFy9ezLRphz7eZ9GiRUyfPn2oXrhwIXfeeedQvWDBAmbMmDFUz58/n5kzZw7V8+bNY9asWUP13LlzmT370NdwzZkzhzlz5gzVs2fPZu7cuUP1rFmzmDdv3lA9c+ZM5s+fP1TPmDGDBQsWDNV33nknCxcuHKqnT5/OokWLhupp06axePGhr++aOnUqS5YsGaqnTJnCsmXLhuqqqnzuFfLc29W/ix3LdhyqV+1ix/2H6p0rd7Jz+c5D9Yqd7FzRUi/fyc6Vh+od9+9g16pdh+plO9jVf6jevnQ7Tz/09KF6yXaefrilvnc7ux/ZPVRvW7yN3QMt9aJt7F5zqN76H1t5Zu0zh+p7tvLMt+s6DyRb79nKnnX17z73N/X6uj7w7IG6fryp99b13if2AjC4Z7Cuv9PUu+v62U3P1vXTTb25rvfv3F/XTzb1jrret7V+rPdt21fX25p6a13v37EfOPmee77vHf/73kjaDTF/Cuyl/pbqedQffPd/gH5gEPhEm9uTJEk6LpHZ3gfwRsTLgOuAy4EzgCeBOcBNmbnjCDftiL6+vly6dGm325DUgy7+/MXdbkEtHrjqgW63oEJExLLM7Bs+3va3WGfmTuDPmoskSVJXtB1iIuIq4EpgPPCiYYszM18zGo1JkiQdSVshJiI+Cnyc+piY5dTHx0iSJHVcuzMx7wU+lZkfHItmJEmSjlW7Zyf9APVp1pIkSV3VbohZALx+LBqRJElqx1F3J0VEa9C5DpgREU8CXwaeGr5+Zh4Yte4kSZIO41iOidlP/aF2BwXwucOsm8e4TUmSpBNyLIHjEzw3xEiSJHXdUUNMZt7YgT4kSZLa0u6BvZIkST3BECNJkopkiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkQwxkiSpSB0NMRFxW0RsioiVw8Z/PyJWR8SqiPirlvHrI2KgWXZ5J3uVJEm97dQO318F/AMw5eBARPw8MAH4sczcGxFnNOMXAhOBi4BXAfMi4vzMHOxwz5IkqQd1dCYmM+8Gnho2/H7gk5m5t1lnUzM+Abg9M/dm5qPAAHBZx5qVJEk9rReOiTkf+NmIuDciFkTETzTjZwLrWtZb34xJkiR1fHfSSE4FXgH8FPATwB0R8WogRlg3R9pAREwCJgGMHz9+jNqUJEm9pBdCzHpgRmYmcF9EHABOb8bPblnvLODxkTaQmbcAtwD09fWNGHSk57jxe7vdgYa7cXu3O5BUmF7YnTQT+AWAiDgf+B5gC3AXMDEiTouIc4HzgPu61aQkSeotHZ2JiYhpwJuA0yNiPXADcBtwW3Pa9bPAVc2szKqIuAN4ENgPXOuZSZIk6aCOhpjMvPIwi37rMOtPBiaPXUeSJKlUvbA7SZIkqW2GGEmSVCRDjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkXrhu5MkSRoz/a99XbdbUIvXPdQ/attyJkaSJBXJECNJkopkiJEkSUXymBjgnA/P7nYLarH2k2/pdguSpAI4EyNJkopkiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkQwxkiSpSIYYSZJUJEOMJEkqkiFGkiQVyRAjSZKKZIiRJElFMsRIkqQiGWIkSVKROhpiIuK2iNgUEStbxm6MiA0Rsby5/GrLsusjYiAiVkfE5Z3sVZIk9bZOz8RUwBUjjN+UmZc0ly8DRMSFwETgouY2N0fEKR3rVJIk9bSOhpjMvBt46hhXnwDcnpl7M/NRYAC4bMyakyRJRemVY2J+LyJWNLubXtGMnQmsa1lnfTMmSZLUEyHm08BrgEuAJ4C/acZjhHVzpA1ExKSIWBoRSzdv3jwmTUqSpN7S9RCTmd/JzMHMPAB8hkO7jNYDZ7esehbw+GG2cUtm9mVm37hx48a2YUmS1BO6HmIi4pUt5duBg2cu3QVMjIjTIuJc4Dzgvk73J0mSetOpnbyziJgGvAk4PSLWAzcAb4qIS6h3Fa0FrgHIzFURcQfwILAfuDYzBzvZryRJ6l0dDTGZeeUIw7ceYf3JwOSx60iSJJWq67uTJEmSjochRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkQwxkiSpSIYYSZJUJEOMJEkqkiFGkiQVyRAjSZKKZIiRJElFMsRIkqQiGWIkSVKRDDGSJKlIhhhJklQkQ4wkSSqSIUaSJBXJECNJkopkiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkQwxkiSpSIYYSZJUJEOMJEkqUkdDTETcFhGbImLlCMv+KCIyIk5vGbs+IgYiYnVEXN7JXiVJUm/r9ExMBVwxfDAizgbeDDzWMnYhMBG4qLnNzRFxSmfalCRJva6jISYz7waeGmHRTcD/BrJlbAJwe2buzcxHgQHgsrHvUpIklaDrx8RExFuBDZn5zWGLzgTWtdTrm7GRtjEpIpZGxNLNmzePUaeSJKmXdDXERMSLgT8BPjbS4hHGcoQxMvOWzOzLzL5x48aNZouSJKlHndrl+38NcC7wzYgAOAv4RkRcRj3zcnbLumcBj3e8Q0mS1JO6OhOTmQ9k5hmZeU5mnkMdXP5rZm4E7gImRsRpEXEucB5wXxfblSRJPaTTp1hPA/4TuCAi1kfEew+3bmauAu4AHgTmANdm5mBnOpUkSb2uo7uTMvPKoyw/Z1g9GZg8lj1JkqQydf3sJEmSpONhiJEkSUUyxEiSpCIZYiRJUpEMMZIkqUiGGEmSVCRDjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkQwxkiSpSIYYSZJUJEOMJEkqkiFGkiQVyRAjSZKKZIiRJElFMsRIkqQiGWIkSVKRDDGSJKlIhhhJklQkQ4wkSSqSIUaSJBXJECNJkopkiJEkSUUyxEiSpCIZYiRJUpEiM7vdw6iKiM3At7vdR5ecDmzpdhPqOB/3k4+P+cnpZH7cfzgzxw0ffN6FmJNZRCzNzL5u96HO8nE/+fiYn5x83L+bu5MkSVKRDDGSJKlIhpjnl1u63YC6wsf95ONjfnLycR/GY2IkSVKRnImRJElFMsQUJCKqiFjbUp8TERkRV3evK3VbRFzdPA/OaRlbGxFV97rSsYqIG5vH79Re2pZOnI/H2PMXW7YngJ8GvtXtRiRJ6jRDTMEycy+wuNt9SJI6LyJOa/4ODB8P4IWZ+WwX2uoodyf1gIj4kYj4vxHxaEQ8ExFrIuLTEfGKo9xuxN1JEfHGiPhqRGyPiKcj4psR8d5h67yvGd8TEVsi4taI+P4x+PF0BBFxfkR8KSI2NY/FYxHxzxFxakS8KCJuioiVEbErIjZGxKyIeG2b93Fp8zyZMMKyKiLWR8Qpo/dT6Ti9LiLmR8TuiHgiIj4REUPv0RFxQfNc2da8TyyOiCuOZ1vquHMjYnbzOv52RHzs4ONxrK/zlt3GP9e8R2wD7m2WrY2IqRHxnoh4CHgWeHtEbI6Im4Y307Kttt5LepFP6t7wKmA9cB1wOfAJ4BeBL7e7oeYP1deA7wGuASYAtwE/3LLOJ4GbgXnAW4E/Bq4A/t0/Zh33b8CZwPupH/sPA3upX5unAS8D/hx4S7POi4DFEfFDx3oHmbkMWEL9fBgSEd8HvAv4bGYOnugPohM2k/o1+Tbgi8BHgY8BRMSrgIXA64Hfo37ctgGzI+JX2tmWuuJLwNepH4+ZwMeBq5pl7b7OvwA8CvwG9fvFQT8P/GGz7SuApcDngKsi4kXDtnENsCAzHzrBn6v7MtNLj12od/P9NyCBH28Zr4C1LfU5zTpXN3UAa6mfvC84zLbPAQaBjw0bf0Ozrbd1++c/WS7U34OSwFuPcf1TgBcDO4EPtoxf3WznnJaxtUA1bJ1B6u8fOTj2B8B+4Kxu/y5O5gtwY/P4fXjY+Geax/r7gL9uHqsfGfZ8WA18o51tdfvnPZkuLY/H7wwbfwCYe5jbHO11ftMIt1kL7AZ+aNj4uc3r/t0tYz/WbGdit38/o3FxJqYHRMT3RMRHIuKhiHgG2Afc0yy+oI1NXUA94/LZzDxwmHXeTP2//C80uyxObY6cvxfYAfzc8f0UOg5PAmuATza7984bvkJEvCsi7m2mjvcDTwMvpb3nBcDt1P9zf1/L2DXA7Mxcfxy9a/TdMay+nfqx/lHq1+XizBw4uDDr2bNpwCUR8fI2tqXOmz2sXgmMP1i0+Tr/0mHuY3FmbmwdyMxHga/w3FnYa4DNwIx2foBeZYjpDX9BndinUk8nXga8o1k2fBrwSH6g+fdIf5TOaP4doA5LrZeXt2xDYyzr/xa9mXrm7C+Ah5vjod4PEBG/BkwH+oH/Afwk8BPUb0DtPC/IzD3UU8vvbYLrzwIXAv84Sj+OTtx3DlOfCXw/9dmIw22knoEdfvzckbalzntqWL2X5jV8HK/zkZ4HRxq/GXhDRPxoRLwE+C3gc/k8OejXs5N6w0RgSmb++cGBiHjpcWzn4Fe0H+mN6snm318Gth5huTogM9cAv92cTXDweIebo/48oInAQGZefXD9iHgh9R+04/Fp6n3mE4C3U09Bf+V4e9eo+0HqmbnWGmAD9R/BkY6P+CHqXQPD/0geaVvqLe2+zg/3MfuHG/8y9Wv9GuCb1MffPG++vsCZmN7wYuqZkFa/cxzbeZj6yfq7zR/FkXwVOACMz8ylI1wePY771QnK2nLqkAH1tP+LqaeWW72bep/58dzHt4C51Ady/wbwmSPsdlTnvWtYPRHYRb3rYQHwU/HcDzQ8BfhN4P7M3NnGttRbRvV1PlzzGv+nZpu/B8xr3gueF5yJ6Q1zqI8gf4B6N887gJ9pdyOZmRFxHfW+zq9HxD9ST0m+DjgjM2/IzG9FxF8C/xARF1C/Oe4BzqbetfHZzJw/Gj+Ujiwifgz4FPVU8gD1m9bV1G9oX6c+RultzSmS/wZcSn0w7rYTuNubgX+lDs23ncB2NPre15x2u4T6TLXfBW7MzG3Nc+Bq4KsRcQP1c+N/AedT74I+5m2N+U+hds1h9F/nw91KfcjC64FfH8Xtdp0hpjf8PvV+7clN/WXgSuC+djeUmf8aEW+mPqXy1mb4W8DftazzkYjoB65tLgmsoz41+5Hj+xF0HDYCj1HPvpxFHSYfAP57Zi6LiPupw+V7qKeClwC/xuEP7DsWs6nPYvjy8IMA1XUTgL+nfu1upz7l9s8AMvPxiPhvwF9S7xY8DVgOvCUz57SzLfWczzD6r/PnyMzNEbEAuBi4a7S22wv8FmvpJNIE3LnAL2Xm17rdj6SxF/UHpz4G/F1mfrTb/YwmQ4x0EoiI1wCvBm4C9mbmpV1uSdIYi4hx1KdpfwD4VerPGTrcWUxF8sBe6eTwUeDfqU/t/O0u9yKpM95C/ZljlwFXPd8CDDgTI0mSCuVMjCRJKpIhRpIkFckQI0mSimSIkSRJRTLESJKkIhliJElSkf4/PPwNiCbyuPsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sample = gen_model.get_samples({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,})\n", "\n", "plot_family_heights(sample)\n", "sample" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or the mean values of specific variables in our graph" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'alice': array([167.0568489]),\n", " 'bob': array([174.79392627]),\n", " 'harry': array([174.60622365]),\n", " 'sally': array([166.89633946])}" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "means = gen_model.get_means({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,}, n_samples=1000)\n", "means" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "or the standard_deviations ..." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'alice': array([10.06116075]),\n", " 'bob': array([9.96899764]),\n", " 'harry': array([9.94216823]),\n", " 'sally': array([9.92119949])}" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAFuCAYAAABwX0lLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqeUlEQVR4nO3de3xV9Z3v/9cHwkUFEQVvXCre76U1tTen1rZWz7Qz1I5t8Te9OO2pfUztzLQz5zL2nGk7neNjeuZxZjz9TY+dY6c205toHbRYLIJKsdRASCBAIIDhHhBBrgkhISTf88fahDQNSJDs7CWv5+OxH+zvd6299mdnfffeb9ZtR0oJSZKkvBk00AVIkiSdCEOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKpaKGmIiYEBFzI6I+IlZExF8U+j9WaHdGRHmPx9wXEQ0RsToibitmvZIkqXSVFfn5DgF/lVJaHBEjgZqImAPUAR8F/m/3mSPiamAqcA1wIfBsRFyeUuooct2SJKnEFHVLTErp5ZTS4sL9JqAeGJdSqk8pre7lIVOAaSmltpTSeqABuLF4FUuSpFI1YMfERMRFwFuAhceYbRywuVu7sdAnSZJOccXenQRARIwA/h34ckpp37Fm7aXvd34nISLuAe4BOOOMM2648sorT0qdkiRp4NXU1LyaUhrbs7/oISYihpAFmJ+klKa/xuyNwIRu7fHA1p4zpZQeAh4CKC8vT9XV1SepWkmSNNAiYmNv/cU+OymA7wP1KaV/Oo6HzACmRsSwiJgEXAZU9WeNkiQpH4q9JebdwKeA5RFRW+j7KjAM+GdgLDAzImpTSrellFZExGPASrIzm+71zCRJkgRFDjEppfn0fpwLwBNHecz9wP39VpQkScolr9grSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyqaghJiImRMTciKiPiBUR8ReF/rMjYk5EvFT4d3S3x9wXEQ0RsToibitmvZIkqXQVe0vMIeCvUkpXAe8A7o2Iq4G/Bp5LKV0GPFdoU5g2FbgGuB14MCIGF7lmSZJUgooaYlJKL6eUFhfuNwH1wDhgCvBvhdn+DfhI4f4UYFpKqS2ltB5oAG4sZs2SJKk0DdgxMRFxEfAWYCFwXkrpZciCDnBuYbZxwOZuD2ss9PVc1j0RUR0R1Tt27OjXuiVJUmkYkBATESOAfwe+nFLad6xZe+lLv9OR0kMppfKUUvnYsWNPVpmSJKmEFT3ERMQQsgDzk5TS9EL3KxFxQWH6BcD2Qn8jMKHbw8cDW4tVqyRJKl3FPjspgO8D9Smlf+o2aQbwmcL9zwA/79Y/NSKGRcQk4DKgqlj1SpKk0lVW5Od7N/ApYHlE1Bb6vgp8C3gsIj4HbAI+BpBSWhERjwEryc5sujel1FHkmiVJUgkqaohJKc2n9+NcAN5/lMfcD9zfb0VJkqRc8oq9kiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQplwwxkiQpl4oaYiLi4YjYHhF13freHBGVEbE8Ip6KiDO7TbsvIhoiYnVE3FbMWiVJUmkr9paYCuD2Hn3/Cvx1Suk64AngPwNExNXAVOCawmMejIjBxStVkiSVsqKGmJTSC8CuHt1XAC8U7s8B/qhwfwowLaXUllJaDzQANxalUEmSVPJK4ZiYOuAPC/c/Bkwo3B8HbO42X2Oh73dExD0RUR0R1Tt27Oi3QiVJUukohRDzWeDeiKgBRgIHC/3Ry7yptwWklB5KKZWnlMrHjh3bT2VKkqRSUjbQBaSUVgEfBIiIy4EPFSY1cmSrDMB4YGtxq5MkSaVqwLfERMS5hX8HAf8d+JfCpBnA1IgYFhGTgMuAqoGpUpIklZqibomJiEeA9wJjIqIR+DowIiLuLcwyHfgBQEppRUQ8BqwEDgH3ppQ6ilmvJEkqXZFSr4eZ5FZ5eXmqrq4e6DIkSdJJEhE1KaXynv0DvjtJkiTpRBhiJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLhliJElSLpUNdAGSVCwP1j440CWomy9O/uJAl6CcM8RIkt7Qdvzzdwa6BHUz9s++dNKW5e4kSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS4YYSZKUS0UNMRHxcERsj4i6bn2TI2JBRNRGRHVE3Nht2n0R0RARqyPitmLWKkmSSluxt8RUALf36PsH4G9TSpOBrxXaRMTVwFTgmsJjHoyIwUWrVJIklbSihpiU0gvArp7dwJmF+6OArYX7U4BpKaW2lNJ6oAG4EUmSJKBsoAsAvgw8ExH/iyxUvavQPw5Y0G2+xkLf74iIe4B7ACZOnNhvhUqSpNJRCgf2/inwlZTSBOArwPcL/dHLvKm3BaSUHkoplaeUyseOHdtPZUqSpFJSCiHmM8D0wv2fcWSXUSMwodt84zmyq0mSJJ3iSiHEbAVuLtx/H/BS4f4MYGpEDIuIScBlQNUA1CdJkkpQUY+JiYhHgPcCYyKiEfg68Hng2xFRBrRSOLYlpbQiIh4DVgKHgHtTSh3FrFeSJJWuooaYlNJdR5l0w1Hmvx+4v/8qkiRJeVUKu5MkSZL6zBAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyyRAjSZJyqRR+O2nAPTBnzUCXoG6+cuvlA12CJCkHDDE6Nc39+4GuQD3dct9AVyApZ9ydJEmScskQI0mScqlPISYiOiLixqNMuyEi/G0jSZJUFH3dEhPHmDYYSK+jFkmSpON2XAf2RsQgjgSYQYV2d6cB/wF49STWJkmSdFSvGWIi4uvA1wrNBPzmGLM/eDKKkiRJei3HszvpV8A3gb8j2xrzcKHd/fbfgDuAv+iXKvtg586d1NbWAtDR0UFFRQXLli0DoL29nYqKCurq6gBobW2loqKCVzatBeBg6wGqZk9n++Z1ALQd2E/V7Ons2LIRgAP7m6iaPZ2dL28GoKVpL1Wzp7NrWyMAzXt3UTV7Oru3vwxA0+6dVM2ezp5XXwFg364dVM2ezr5dOwDY8+orVM2eTtPunQDs3v4yVbOn07x3FwC7tjVSNXs6LU17s9f28maqZk/nwP4mAHZs2UjV7Om0HdgPwPbN66iaPZ2DrQcAeGXTWqpmT6f9YBsAL69fQ9Xs6XQcagdg67pVVM2eTmdHdijTlrX1VM2e3vW33LymjkVznuxqb1q9jOpnf97V3lBfy+K5v+hqr1+xmCW/erqrva6umqUvzOpqNyyrYtn82V3tl2oXsPw3z3a11yx+kRWVz3e1Z8+ezcyZM7vas2bNYtasI8ubOXMms2cfWd5TTz3Fs88eWd6TTz7J3Llzu9rTp09n3rx5Xe3H5y1n/vINXe1H5y7jxbqNXe1HnlvKgpWbuto/nrOERasau9o/nL2YmjVbutoVs2qobdgKQEdHJxWzali2NhsL7Yc6qJhVQ936bQC0HmynYlYN9Ru3A9DSepCKWTWs3pSNjeYDbVTMqqFhSzY29ja3UjGrhnVbs7Gxu+kAFbNq2LBtNwCv7t1PxawaNm/fA8D23c1UzKphy45s7Gzb1UTFrBq27crGzpYde6mYVcP23c0AbN6+h4pZNby6NxtLG7btpmJWDbubsrG0busuKmbVsLe5FYCGLTupmFVD84FsbK3etIOKWTW0tB4EoH7jdipm1dB6MBtrdeu3UTGrhvZD2VhbtvZlKmbV0NHRCUBtw1YqKiq6/pY1NTX88Ic/7GovWrSIH//4x13tBQsW8Mgjj3S1X3zxRR599NGu9vz583n88ce72vPmzWP69CNje82iNSydu7SrvWrBKpbPW97Vrq+sp+7XdV3tFb9ZwYrfrOhq1/26jvrK+q728nnLWbVgVVd76dylrFl05PpTtc/V8lL1S13txXMW07Ck4cjrfaaGdbXrutrVv6xm/bL1Xe2qmVVsXHFkbC58aiGbuo3NyhmVbF6dfS51dnRSOaOSLYWx2dHeQeWMSrYWxmZ7WzuVMyrZti4biwcPHKRyRiWvbMg+p1pbWqmcUcn2TdnYPNB0gMoZlbzamG1ob9nXQuWMSnZuzcZm8+5mKmdUsrswFpt2NVE5o5I9hbG499W9VM6oZO+r2Vjcs30PlTMqaSqMxc2bN1NRUcGrr2bL37BhAxUVFezenS1v3bp1VFRUsHdv9viGhgYqKipobs7G7urVq6moqKClpQWA+vp6KioqaG3NxmpdXR3TFi6gvfA5t2LLFqYtXEBHZzb2ljc2Mm3hgiPrbtMmHq1a2NVesnEjP1u06Mi6Wb+e6TXVR9bNunU8ubimq71g7VqeWrKkq/3iSy/xi6W1Xe35a9bw9LIjY2/e6lU8s/zI2JtbX8+cFUfG3vMrV/L8ypVd7Tkr6phbf2TsPbN8OfNWHxl7Ty9byvw1R8beL5bW8uJLR8beU0uWsGDt2q72k4trqFp3ZOxNr6mmev2RsfezRYtYsvHI2Hu0aiFLNx0Ze9MWLmB5Y/a52NHZybSFC1ixJRt77R0dTFu4gPqt2dhrbW/P2oX6W1paqKioYPXq1QA0NzdTUVFBQ0P23ti7dy8VFRWs61ZfT6+5JSalNA+YBxARCfheSmnraz1OkkrNzxuyEN5S1kLnyE42NGzI2kNb6BzSybqG7MNy/7D9MBTWNmQf9vtPy8Ld4Q/X/afvh4A1DdmXRfOIZgbFIFY3FD6MRzYzaPAg6huyD+vmUVl7ZUP2ZdR0VhNlUcaKhiwYNZ3dRBllLG/IvsyazmmiLJWxrCH7D9i+sftYf2g9tQ21WfvcfWw4tIElDdmX5d7z9rKpbRPDG4aTUmLfefvY3LaZYQ3DjrQPZO3Ozk6azmuisaWRoQ1Dj7Sbe7SbsnZHRwfN5zWzZe8WhrQOOdLes4UhLUPoOFRo79rCkOYhHDp0iP3n7Wfrzq2U7SvjUHuhvWMrZXu6tbdv5b/yX/thLetUEim9sY7FLS8vT9XV1a89Yzdesbe0FOWKvV7srvQU4WJ3tz1+W78/h47fM3c+U5Tneel97y/K8+j4XPb8c31+TETUpJTKe/b3+Yq9EXEzcBcwERjeY3JKKTlaJJWkKZdOGegSJJ1EfQoxEfEF4LvATuAloK3nLCepLkmSToqz7rhjoEtQP+nrlpi/An4KfDaldLAf6pEkSToufb3Y3TjgBwYYSZI00PoaYmqAi/ujEEmSpL7oa4j5c+DLEfGe/ihGkiTpeB3PFXs389u/iTQKmBsRLcDuHrOnlNKbTmJ9kiRJvTqeA3ufwx92lCRJJeZ4rth7dxHqkCRJ6pO+HhMjSZJUEvp6sbtPH2NyJ7AXWJJSajzGfJIkSa9bXy92V8GR42O6X523e19nRDwK/InXk5EkSf2lr7uT3g1sBL4D3AxcWfj3QWAT8CHgPuAO4BsnrUpJkqQe+rol5j8B01JKX+3Wtwb4dUQ0AfeklO6IiDOBPwa+2ttCJEmSXq++bom5leyU6948Dxz+BesXyH6iQJIkqV/0NcQcBG44yrQbCtMPL3d/zxki4uGI2B4Rdd36Ho2I2sJtQ0TUdpt2X0Q0RMTqiLitj7VKkqQ3sL7uTvoZ8LcR0QE8DmwHzgU+RnYMzMOF+SYDq3t5fAXZ8TQ/PNyRUvrE4fsR8Y9kZzgREVcDU4FrgAuBZyPi8pRSRx9rliRJb0B9DTF/CYwE/qFw6+6nwF8V7tcBlT0fnFJ6ISIu6m3BERHAx4H3FbqmkB1/0wasj4gG4MbelitJkk49fQoxKaUDwCcj4pvA24ELgJeBhSmlNd3mm3kCtfwe8EpK6aVCexywoNv0RjzORpIkFfR1SwwAhcCy5jVn7Ju7gEe6taOXeXr9DaeIuAe4B2DixIknuSxJklSKjudXrCcCL6eU2gv3jymltKmvRUREGfBRfvug4UZgQrf2eGDrUZ7zIeAhgPLycn+sUpKkU8DxbIlZD7wTqAI28Nq/aD34BOr4ALCqx88VzAB+GhH/RHZg72WFGiRJko4rxHwWWNvt/glv6YiIR4D3AmMiohH4ekrp+2RnIXXflURKaUVEPAasBA4B93pmkiRJOuw1Q0xK6d+63a94PU+WUrrrKP13H6X/fuD+1/OckiTpjamvF7sDICIGRcS1EXFzRJxxsouSJEl6LX0OMRFxL7ANWEb2UwNXFPqfjIg/P7nlSZIk9a5PISYiPg98G3iS7MJ03U+D/jXwRyetMkmSpGPo65aYvwT+MaV0D/BEj2mrKGyVkSRJ6m99DTGTgGeOMm0/cNbrqkaSJOk49TXEvApcdJRpVwBbXlc1kiRJx6mvPzvwFPC1iPgVsLHQlyJiDPAVsmNlBtTOnTupra1l8uTJdHR08KMf/Yi3vvWtXH/99bS3t/OTn/yE8vJyrr32WlpbW5k2bRo/3TiWocNOI3V20Lx3F8NPH8GQYafR2dnB/r27GH76SIYMG05nxyH279vN8DNGMmTokfZpZ5xJ2dBhdBxqp6VpD6eNOJOyId3boygbMpSOQwdpadrL6SNHMbhsKIfaD3KgeS+njzyLwWVDONTexoHmfUfaB9s4sH8fZ5w5mkGDy2g/2Err/qYj7bZWWluaOGPU2QwaNJj2tgO0tjQzYtTZxKDBHGw7QFtLMyPOOoeIQRxsbaHtwP7faY88awxE0HaghYOt+xk5eiwAbQf2c7D1ACNHj8naLc0cbGvtare2NHPoYBsjzjona+9v4lD7wd9uH2pnxKizATiwfx+dhw5xxuF28z46Ozo4Y9ToQnsvnZ2dfOXWywGYPXs27e3tfOhDHwJg1qxZANx+++0AzJw5kyFDhvDBD34QgKeeeorTTjuND3zgAwA8+eSTjBo1iltuuQWA6dOnc84553DzzTcD8Pi85Zx/9khuuu4iAB6du4wJY0fxrmvfBMAjzy1l0gWjecfV2YWqfzxnCVdMGMvbrhwPwA9nL+aai87jhsuzn/SqmFXD5EsvYPKlF9LR0cmP5izhrZddyPWXXED7oQ5+8mwt5VeM49pJ59N6sJ1pzy/j7VdN4Ko3nUtL60Ee+9Vy3nn1RK6YOJbmA208Pq+Om667iEvHncPe5laemL+C91w/iYsvPJvdTQf4+W9W8t7JF3PR+aN5de9+flG5ive/9RImnHsW23c38/TC1dx6w6WMGzuKbbuamFW1httvvJzzzx7Jlh17mVPTwO+//QrOHT2Czdv38NzitXz4nVcyZtQZbNi2m1/VrmPKu69m9MjTWLd1Fy8sW88dN13DqBHDadiyk/nLN3Dnzdcy4rRhrN60g8qVm/j4e6/j9OFDqd+4nYX1m5n6vusZPnQIdeu3Ub16C3/8gckMKRvMsrUvs/ilrXzq1rcwePAgahu2UruxgrvvvhuAmpoaVqxYwac//WkAFi1axOrVq/nkJz8JwIIFC1i/fj133ZVdteHFF19k8+bNfOITnwBg/vz5bNu2jTvvvBOAefPmsXPnTrg4+5xYs2gNB5oP8OZb3gzAqgWraG9r57qbrwOgvrKejkMdXPt71wKw4jcrALjm3dcAUPfrOgaXDeaqd14FwPJ5yxkybAhXvuNKAJbOXcppI07j8rdlY7n2uVrOGHUGl5VfBsDiOYs5c8yZXPqWS7PX+0wNo88bzcWTswKrf1nNOePOYdL1kwComlnFeRedx5uuycbmwqcWcsElFzCxMDYrZ1Qy/orxTLhiAp0dnSycuZCJV05k3OXj6GjvoOqXVbzp6jdx4aUX0t7WTvUz1Uy6dhLnX3w+Bw8cpGZODRdffzHnXXQerS2tLHl2CZdMvoRzJ57LgaYD1M6t5bK3XsaY8WNo2dfC0l8t5fLyyznnwnNo3t3M8l8v58obr2T0+aNp2tVE3fw6rnrHVZx17lnsfXUvK19cydXvuppRY0axZ/se6hfUc+1N2d928+bNPPfcc3z4wx9mzJgxbNiwgV/96ldMmTKF0aNHs27dOl544QXuuOMORo0aRUNDA/Pnz+fOO+9kxIgRrF69msrKSj7+8Y9z+umnU19fz8KFC5k6dSrDhw+nrq6OuQsX8Eflb2PI4MGs2LKF5Y2b+djbbmTwoEEsb2xkxZZGpr79Hdm627SJVdte5hM3vh2AJRs30rB9Ox9729uydbN+PZt27eSjN5Rn62bdOrbu2c1H3ppdcH7B2rXs2LePP3jLW7Kx+dJL7GrZz4ffPDkbm2vWsK/1AL9/fTb25q1eRevBdm67Lht7c+vrOdTZwa3XZH+f51euBOB9V18NwJwVdZQNGswtV2Vj75nlyxk+dAg3X5GNvaeXLeXM4adx0+XZ2PvF0lrOPv0M3nVZNvaeWrKEsWeeyTsuuQSAJxfXcOFZo7nx4mzsTa+pZuLZ51A+KRt7P1u0iEvPPZe3vCkbe49WLeTK8y/gzYWf+Jm2cAHXjBvPdePH09HZyc8WVXHd+AlcM24c7R0d/Hv1It48YSJXXXghre3tPLm4hg984P1cddVVtLS08Nhjj/HOd76TK664gubmZh5//HFuuukmLr30Uvbu3csTTzzBe97zHo6mryHmv5P9ynQdsJDswnf/P3AlsB34Zh+XJw2M2p9wM/sZum8w1A4H4JbYz7C9ZVA7DID3Dd7PaXuOtG8ta+a0XUO62h8c0swZO4dC7VAAbh/azIgdQ6F5KINSob19KDQNZXChPfKVZbB3KENT4vah+xn58nLYPYRhnVn7zJeXw64hnFZoj2qsgx1lnNHRye1DW7L29jJGFNpnbaqDbWWMOlRob6yDrWWcdaiD24ceYPSGFbBlMGcX2mevWwGbBjOmPWuftXYFlA1mbPshbh/ayqiXVkLZIM4/mLVHrFkJgwdxQaF9xuqsPa4ta59WXw+Dgglt7Ywa2sawlVl7Yms7o4e2MXRFPUQwqfUgY4YeZPDyVRBw8YGDnDv0IIOWZe1LDxyk9sIbellRknR0kVLfLsAbESOBLwO3AecCO4FZwAMppX0nu8C+Ki8vT9XV1X16zANzTvZvWer1OLwlpl89cG3/P4f65it1/f4UD9Y+2O/PoeP3xclfLMrz7Pjn7xTleXR8xv7Zl/r8mIioSSmV9+zv869Yp5SagL8r3KR8mvzHA12BJOl16nOIiYjPAHcBE4HhPSanlNIlJ6MwSZKkY+lTiImIvwH+luyYmFqgrR9qkiRJek193RLzOeDbKaWv9EcxkiRJx6uv14k5h+w0a0mSpAHV1xAzD3hzfxQiSZLUF6+5OykiugedLwPTI2In8DSwq+f8KaXOk1adJEnSURzPMTGHyC5qd1gAPzjKvOk4lylJkvS6HE/g+Ca/HWIkKZe+u/S7A12CuinWxe5e/T//pyjPo+NzIhe7O5rXDDEppW+ctGeTJEk6Sfp6YK8kSVJJMMRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcKmqIiYiHI2J7RNT16P+ziFgdESsi4h+69d8XEQ2FabcVs1ZJklTayor8fBXAd4AfHu6IiFuAKcD1KaW2iDi30H81MBW4BrgQeDYiLk8pdRS5ZkmSVIKKuiUmpfQCsKtH958C30optRXm2V7onwJMSym1pZTWAw3AjUUrVpIklbRSOCbmcuD3ImJhRMyLiLcV+scBm7vN11jokyRJKvrupN6UAaOBdwBvAx6LiIuB6GXe1NsCIuIe4B6AiRMn9lOZkiSplJTClphGYHrKVAGdwJhC/4Ru840Htva2gJTSQyml8pRS+dixY/u9YEmSNPBKIcQ8CbwPICIuB4YCrwIzgKkRMSwiJgGXAVUDVaQkSSotRd2dFBGPAO8FxkREI/B14GHg4cJp1weBz6SUErAiIh4DVgKHgHs9M0mSJB1W1BCTUrrrKJM+eZT57wfu77+KJElSXpXC7iRJkqQ+M8RIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcKoXfThpw337upYEuQd185dbL+/9J5n2r/59DfXPLfQNdgaSccUuMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKJUOMJEnKpaKGmIh4OCK2R0Rdt75vRMSWiKgt3H6/27T7IqIhIlZHxG3FrFWSJJW2Ym+JqQBu76X/gZTS5MLtaYCIuBqYClxTeMyDETG4aJVKkqSSVtQQk1J6Adh1nLNPAaallNpSSuuBBuDGfitOkiTlSqkcE/OliFhW2N00utA3DtjcbZ7GQp8kSVJJhJjvApcAk4GXgX8s9Ecv86beFhAR90REdURU79ixo1+KlCRJpWXAQ0xK6ZWUUkdKqRP4Hkd2GTUCE7rNOh7YepRlPJRSKk8plY8dO7Z/C5YkSSVhwENMRFzQrXkHcPjMpRnA1IgYFhGTgMuAqmLXJ0mSSlNZMZ8sIh4B3guMiYhG4OvAeyNiMtmuog3AFwBSSisi4jFgJXAIuDel1FHMeiVJUukqaohJKd3VS/f3jzH//cD9/VeRJEnKqwHfnSRJknQiDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXDDGSJCmXihpiIuLhiNgeEXW9TPtPEZEiYky3vvsioiEiVkfEbcWsVZIklbZib4mpAG7v2RkRE4BbgU3d+q4GpgLXFB7zYEQMLk6ZkiSp1BU1xKSUXgB29TLpAeC/AKlb3xRgWkqpLaW0HmgAbuz/KiVJUh4M+DExEfGHwJaU0tIek8YBm7u1Gwt9vS3jnoiojojqHTt29FOlkiSplAxoiImI04H/Bnytt8m99KVe+kgpPZRSKk8plY8dO/ZklihJkkpU2QA//yXAJGBpRACMBxZHxI1kW14mdJt3PLC16BVKkqSSNKBbYlJKy1NK56aULkopXUQWXN6aUtoGzACmRsSwiJgEXAZUDWC5kiSphBT7FOtHgErgiohojIjPHW3elNIK4DFgJTALuDel1FGcSiVJUqkr6u6klNJdrzH9oh7t+4H7+7MmSZKUTwN+dpIkSdKJMMRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcMsRIkqRcipTSQNdwUkXEDmDjQNcxQMYArw50ESo61/upx3V+ajqV1/ubUkpje3a+4ULMqSwiqlNK5QNdh4rL9X7qcZ2fmlzvv8vdSZIkKZcMMZIkKZcMMW8sDw10ARoQrvdTj+v81OR678FjYiRJUi65JUaSJOWSISZHIqIiIjZ0a18UESki7h64qjTQIuLuwji4qFvfhoioGLiqdLwi4huF9VdWSsvS6+f66H/+YfPtZeCdwNqBLkSSpGIzxORYSqkNWDDQdUiSii8ihhW+B3r2BzAkpXRwAMoqKncnlYCIuDQifhQR6yPiQESsi4jvRsTo13hcr7uTIuLmiJgTEXsjYn9ELI2Iz/WY5/OF/taIeDUivh8RZ/fDy9MxRMTlEfFERGwvrItNEfGziCiLiOER8UBE1EVEc0Rsi4inIuLKPj7HDYVxMqWXaRUR0RgRg0/eq9IJuioi5kZES0S8HBHfjIiuz+iIuKIwVvYUPicWRMTtJ7IsFd2kiJhZeB9vjIivHV4fx/s+77bb+D2Fz4g9wMLCtA0R8eOI+GxErAIOAndExI6IeKBnMd2W1afPklLkoC4NFwKNwJeB24BvAu8Hnu7rggpfVM8BQ4EvAFOAh4E3dZvnW8CDwLPAHwL/Gbgd+KVfZkX3C2Ac8Kdk6/6vgTay9+YwYCTwP4APFeYZDiyIiPOP9wlSSjXAIrLx0CUizgI+DvxrSqnj9b4QvW5Pkr0nPwL8FPgb4GsAEXEhMB94M/AlsvW2B5gZEf+hL8vSgHgCeJ5sfTwJ/C3wmcK0vr7PfwKsB+4k+7w47BbgLwvLvh2oBn4AfCYihvdYxheAeSmlVa/zdQ28lJK3EruR7ea7CUjAW7r1VwAburUvKsxzd6EdwAaywTvoKMu+COgAvtaj/92FZX1koF//qXIj+x2UBPzhcc4/GDgdaAK+0q3/7sJyLurWtwGo6DFPB9nvjxzu+3PgEDB+oP8Wp/IN+EZh/f11j/7vFdb1WcD/KqyrS3uMh9XA4r4sa6Bf76l067Y+/qRH/3Jg9lEe81rv8wd6ecwGoAU4v0f/pML7/lPd+q4vLGfqQP99TsbNLTElICKGRsRXI2JVRBwA2oFfFyZf0YdFXUG2xeVfU0qdR5nnVrL/5f+ksMuirHDk/EJgH/CeE3sVOgE7gXXAtwq79y7rOUNEfDwiFhY2HR8C9gMj6Nu4AJhG9j/3z3fr+wIwM6XUeAK16+R7rEd7Gtm6vpbsfbkgpdRweGLKtp49AkyOiDP7sCwV38we7Tpg4uFGH9/nTxzlORaklLZ170gprQee4be3wn4B2AFM78sLKFWGmNLw92SJ/cdkmxNvBD5amNZzM+CxnFP491hfSucW/m0gC0vdb2d2W4b6Wcr+W3Qr2ZazvwfWFI6H+lOAiPgD4FGgHvj/gLcDbyP7AOrLuCCl1Eq2aflzheD6e8DVwL+cpJej1++Vo7THAWeTnY3Y0zayLbA9j5871rJUfLt6tNsovIdP4H3e2zg4Vv+DwLsj4tqIOAP4JPCD9AY56Nezk0rDVOCHKaX/cbgjIkacwHIO/0T7sT6odhb+/SCw+xjTVQQppXXApwtnExw+3uHByK4HNBVoSCndfXj+iBhC9oV2Ir5Lts98CnAH2SboZ060dp1055FtmeveBthC9iXY2/ER55PtGuj5JXmsZam09PV9frTL7B+t/2my9/oXgKVkx9+8YX6+wC0xpeF0si0h3f3JCSxnDdlg/Y+FL8XezAE6gYkppepebutP4Hn1OqVMLVnIgGyz/+lkm5a7+xTZPvMTeY61wGyyA7nvBL53jN2OKr6P92hPBZrJdj3MA94Rv31Bw8HAJ4AlKaWmPixLpeWkvs97KrzH/29hmV8Cni18FrwhuCWmNMwiO4J8Odluno8C7+rrQlJKKSK+TLav8/mI+BeyTZJXAeemlL6eUlobEf8T+E5EXEH24dgKTCDbtfGvKaW5J+NF6dgi4nrg22SbkhvIPrTuJvtAe57sGKWPFE6R/AVwA9nBuHtex9M+CPycLDQ//DqWo5Pv84XTbheRnan2H4FvpJT2FMbA3cCciPg62dj4InA52S7o415Wv78K9dUsTv77vKfvkx2y8Gbgj07icgecIaY0/BnZfu37C+2ngbuAqr4uKKX084i4leyUyu8XutcC/7vbPF+NiHrg3sItAZvJTs1+6cRegk7ANmAT2daX8WRhcjnw4ZRSTUQsIQuXnyXbFLwI+AOOfmDf8ZhJdhbD0z0PAtSAmwL8M9l7dy/ZKbd/B5BS2hoRNwH/k2y34DCgFvhQSmlWX5alkvM9Tv77/LeklHZExDzgOmDGyVpuKfBXrKVTSCHgzgY+kFJ6bqDrkdT/Irtw6ibgf6eU/mag6zmZDDHSKSAiLgEuBh4A2lJKNwxwSZL6WUSMJTtN+y+A3ye7ztDRzmLKJQ/slU4NfwP8kuzUzk8PcC2SiuNDZNccuxH4zBstwIBbYiRJUk65JUaSJOWSIUaSJOWSIUaSJOWSIUaSJOWSIUaSJOWSIUaSJOXS/wN9xwpoTOTjzAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "stds = gen_model.get_standard_deviations({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,}, n_samples=1000)\n", "\n", "plot_family_heights(means, stds)\n", "stds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Now let's simulate how the family of Alice and Bob would develop if we actually knew both their genetic factor and their height" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's say Alice has a genetic_factor of +5, but still was of average height: 167cm. \n", "Bob has a genetic_factor of +5 as well, but is of average height, too: 175cm \n", "I guess they did not eat their spinach." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "alice_gen_factor = 5.\n", "alice_height = 167.\n", "bob_gen_factor = 5\n", "bob_height = 175." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "family.alice.genetic_factor" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now lets create a generative model conditioned to these data" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "gen_model_wdata = hal.get_generative_model(family, data={family.alice.genetic_factor: alice_gen_factor,\n", " family.alice.height: alice_height,\n", " family.bob.genetic_factor: bob_gen_factor,\n", " family.bob.height: bob_height})" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'alice': array([167.]),\n", " 'bob': array([175.]),\n", " 'harry': array([180.02715988]),\n", " 'sally': array([172.43563186])}" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cond_means = gen_model_wdata.get_means({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,}, n_samples=1000)\n", "cond_means" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'alice': array([0.]),\n", " 'bob': array([0.]),\n", " 'harry': array([7.5037483]),\n", " 'sally': array([7.95827691])}" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cond_stds = gen_model_wdata.get_standard_deviations({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,}, n_samples=1000)\n", "cond_stds" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAFuCAYAAABwX0lLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAApyUlEQVR4nO3deZTV9Z3n/+e7Noq12ApRFkFRBEFRS4yJ7RI10mpCTBvFX9w6mZjTY9Kd7p6ZEzPTSToznmT6THemf52xe0xrTIwtIcZdm7hEMSSKgqKAiKIgICJ7UQXU/pk/7kWvZYEUUrfuV56Pc+6p+/58v/f7fRd3qRff7UZKCUmSpKwp6+0GJEmSDoQhRpIkZZIhRpIkZZIhRpIkZZIhRpIkZZIhRpIkZVJRQ0xEjImIJyJieUQsi4i/yI9/MV93RERdp8fcEBErI2JFRFxQzH4lSVLpqijy+tqAv04pPR8RA4FFEfEosBT4AvB/C2eOiMnALOB44AjgsYg4NqXUXuS+JUlSiSnqlpiU0tsppefz9xuA5cColNLylNKKLh4yE5idUmpOKa0CVgLTi9exJEkqVb12TExEjANOAhbsY7ZRwNqCel1+TJIkHeKKvTsJgIgYAPwa+GZKace+Zu1i7APfkxAR1wHXAfTv3/+U44477qD0KUmSet+iRYs2p5RqO48XPcRERCW5AHNHSunuD5l9HTCmoB4NrO88U0rpZuBmgLq6urRw4cKD1K0kSeptEfFmV+PFPjspgFuA5Smlf9iPh9wPzIqIPhExHjgGeLYne5QkSdlQ7C0xnwKuApZExOL82LeBPsA/AbXAQxGxOKV0QUppWUTMAV4md2bT9Z6ZJEmSoMghJqU0n66PcwG4Zy+PuRG4sceakiRJmeQVeyVJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYVNcRExJiIeCIilkfEsoj4i/z40Ih4NCJey/8cUvCYGyJiZUSsiIgLitmvJEkqXcXeEtMG/HVKaRLwCeD6iJgMfAt4PKV0DPB4viY/bRZwPDADuCkiyovcsyRJKkFFDTEppbdTSs/n7zcAy4FRwEzgZ/nZfgZ8Pn9/JjA7pdScUloFrASmF7NnSZJUmnrtmJiIGAecBCwADkspvQ25oAOMyM82Clhb8LB1+bHOy7ouIhZGxMJNmzb1aN+SJKk09EqIiYgBwK+Bb6aUduxr1i7G0gcGUro5pVSXUqqrra09WG1KkqQSVvQQExGV5ALMHSmlu/PD70TE4fnphwMb8+PrgDEFDx8NrC9Wr5IkqXQV++ykAG4BlqeU/qFg0v3ANfn71wD3FYzPiog+ETEeOAZ4tlj9SpKk0lVR5PV9CrgKWBIRi/Nj3wZ+CMyJiK8Aa4AvAqSUlkXEHOBlcmc2XZ9Sai9yz5IkqQQVNcSklObT9XEuAOfu5TE3Ajf2WFOSJCmTvGKvJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKpKKGmIi4NSI2RsTSgrETI+LpiFgSEQ9ExKCCaTdExMqIWBERFxSzV0mSVNqKvSXmNmBGp7F/Bb6VUpoK3AP8Z4CImAzMAo7PP+amiCgvXquSJKmUFTXEpJSeArZ2Gp4IPJW//yjwJ/n7M4HZKaXmlNIqYCUwvSiNSpKkklcKx8QsBT6Xv/9FYEz+/ihgbcF86/JjHxAR10XEwohYuGnTph5rVJIklY5SCDFfBq6PiEXAQKAlPx5dzJu6WkBK6eaUUl1Kqa62traH2pQkSaWkorcbSCm9AnwGICKOBS7KT1rHe1tlAEYD64vbnSRJKlW9viUmIkbkf5YB/w34l/yk+4FZEdEnIsYDxwDP9k6XkiSp1BR1S0xE3AmcDQyPiHXAd4EBEXF9fpa7gZ8CpJSWRcQc4GWgDbg+pdRezH4lSVLpipS6PMwks+rq6tLChQt7uw1JknSQRMSilFJd5/Fe350kSZJ0IAwxkiQpkwwxkiQpkwwxkiQpkwwxkiQpk3r9YneSJPWkTf/0495uQQVqv/H1g7Yst8RIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMMsRIkqRMKmqIiYhbI2JjRCwtGJsWEc9ExOKIWBgR0wum3RARKyNiRURcUMxeJUlSaaso8vpuA34M/Lxg7O+Av00p/XtEXJivz46IycAs4HjgCOCxiDg2pdRe5J71cfTED3q7A3V2zg293YGkjClqiEkpPRUR4zoPA4Py92uA9fn7M4HZKaVmYFVErASmA08Xo1dJ0sfD9nvu6e0WVKD2G18/aMsq9paYrnwT+E1E/C9yu7c+mR8fBTxTMN+6/NgHRMR1wHUAY8eO7bFGJUlS6SiFEPNnwF+mlH4dEZcBtwDnAdHFvKmrBaSUbgZuBqirq+tyHknSoWnwJZf0dgvqIaVwdtI1wN35+78it8sIcltexhTMN5r3djVJkqRDXCmEmPXAWfn7nwZey9+/H5gVEX0iYjxwDPBsL/QnSZJKUFF3J0XEncDZwPCIWAd8F/gq8I8RUQE0kT+2JaW0LCLmAC8DbcD1npkkSZL2KPbZSVfsZdIpe5n/RuDGnutIkiRlVSnsTpIkSeo2Q4wkScokQ4wkScokQ4wkScokQ4wkScokQ4wkScokQ4wkScokQ4wkScokQ4wkScokQ4wkScqkboWYiGiPiOl7mXZKRPjdRpIkqSi6uyUm9jGtHEgfoRdJkqT9tl9fABkRZbwXYMrydaG+wB8Dmw9ib5IkSXv1oSEmIr4LfCdfJuD3+5j9poPRlCRJ0ofZny0xT+Z/BrkwcwuwrtM8zcDLwIMHrbMDtGXLFhYvXsy0adNob2/n9ttv5+STT+aEE06gtbWVO+64g7q6OqZMmUJTUxOzZ8/mtNNOY9KkSezatYs5c+Zw+umnM3HiRBobG7nrrrs444wzmDBhAvX19dxzzz2ceeaZHHXUUWzbto377ruPs88+m3HjxrF582YefPBBzj33XMaMGcPGjRt5+OGHOf/88xk1ahQbNmxg7ty5zJgxg5EjR/LWW2/x6KOPcuGFFzJixAjWrl3L448/zsUXX8zw4cNZvXo1Tz75JDNnzmTIkCG88cYbPPXUU1xyySXU1NSwcuVK5s+fz6WXXsqAAQNYsWIFTz/9NJdddhn9+vVj+fLlLFiwgFmzZlFdXc3SpUtZuHAhX/rSl6isrOSll17i+eef56qrrqK8vJzFixezePFirr32WgAWLVrEsmXLuPrqqwF47rnnWLFiBVdeeSUAzzzzDKtWreKKK64A4A9/+ANr167l8ssvB2D+/Pls2LCBSy+9FIB58+axZcsWvvCFLwDwxBNPUF9fz+c//3kAHnvsMXbv3s1nP/tZAB555BFaW1u56KKLAJg7dy4AM2bMAOChhx6isrKSz3zmMwA88MAD9O3bl/POOw+Ae++9l5qaGs455xwA7r77boYNG8ZZZ50FwF3zljBy6EDOmDoOgF8+8RJjamv45JQjAbjz8RcZf/gQPjF5LAC/ePQFJo6p5dTjRgPw80ee5/hxh3HKsaMAuG3uIqZNOJxpE46gvb2D2x99gZOPOYITjj6c1rZ27nhsMXUTRzFl/EiaWlqZ/duXOG3SGCYdOYJdTS3MeXIJp08ey8SxtTTubuaueUs5Y+o4JowaRn1jE/fMX8aZJ4znqCOGsq1hN/f9/mXOnnYU40YOYXP9Th58+hXOPfloxowYzMZtjTy8YAXnnzKBUbU1bNjawNxnX2XG9GMZOXQgb22q59FFK7nwtImMGDKAtRu38/jzr3Px6ccxvKY/qzds48nFbzDzU5MZMrAvb6zfylMvreKSM46nZkA1K9/awvwlq7n0rCkM6NuHFWs28fTLa7js7Kn0q65i+ZsbWbB8LbM+fQLVVZUsXbWBhSve4kvnTaOyopyXXn+b519bz1Xnn0R5eRmLV65n8Zu39fhr74F+DwCwfft22traGD58OADbtm2jo6ODYcOGdVlv3boVgKFDh777OVNWVsaQIUO6rDdv3kxFRQWDBw/ust60aRNVVVXU1NQAsHHjRqqrqxk0aFCX9TvvvEO/fv0YOHBgl/WGDRsYMGAAAwYMIKXEO++8w8CBA+nfv/8H6o6ODjZu3MigQYPo16/fB+r29nY2bdpETU0Nffv2pa2tjc2bNzN48GCqq6s/ULe2trJlyxaGDBlCnz59PlC3tLSwdetWhg4dSlVVFc3NzWzbto1hw4bxi8/+oiife08seIY/qTuVyvJylr31FkvWreWLp06nvKyMJevWseytdcw67RMAvLhmDa9seJvLp58GwAtvvsnKjRv54qmnArBw1SrWbN3CF06pA+DZN95g/fZtfP7kU3KvzddfZ9OOHXz2pJNyr83XXmPrrp1cfOK03Gvz1VfZ0bSbC084MffaXPEKTS2tXDB1KgBPLF9OW0c75x8/BYDfvvwyAJ+ePBmAR5ctpaKsnHMmTQLgN0uWUF1VyVkTjwPg4ZdeZFB1X8449lgAHnxxMUP79eeTxxwDwAMvvEDtoEF84uijAbj3+UUcMXgI0486CoC7Fy1k7NBh1I0fD8CvnnuOCSNGcNKRuc/FXz67gONGHs6JY3Ofi7MXPMPxo0YzdfRo2js6+NVzzzJ19BiOHzWK1vZ2fr3wOU4cM5ZJRxxBU2sr9z6/iPPOO7fbf3P35kNDTEppHjAPICIS8JOU0voPe5wklZqF7ywEoH1nO6kjsbp99fvqVW2rAGjb2QaJ9+rGNgAqWiveqwMqWvJ1QxtRFpS3lL+/bi6oywvqHW1ERVDeVFDvDMp3F9S7Cur6NmJ3UL4rV7fWt1LWVPb+urmM8p3lpJRoq2+jvLmcssay9+qWfN2RaNuRr/vso24tp6yqjNSeaGvYR92WaGtsY03bGsoqP1h3tHXQ3tjOmvY1lFWU0dHaQfvOdta2r+2Jp1iHmEjp43Usbl1dXVq4cGFvt6FS98QPersDdXbODT2+iqk/m9rj69D+W3LNkqKsZ9M//bgo69H+qf3G17v9mIhYlFKq6zy+Xwf2dlrQWcAVwFigutPklFI6t9vdSZIkdVO3QkxEfA34Z2AL8Bq5Y2HeN8tB6kuSJGmfursl5q+BfwO+nFJq6YF+JEmS9kt3L3Y3CvipAUaSJPW27oaYRcBRPdGIJElSd3Q3xPw58M2I2PtJ25IkSUWwP1fsXcv7vxOpBngiInYB2zrNnlJKRx7E/iRJkrq0Pwf2Po5f7ChJkkrM/lyx99oi9CEV1+I7ersDdVaEi93p0LT5//yf3m5BBQ7kYnd7091jYiRJkkpCdy92d/U+JncA9cALKaXOXxAplZZpX+rtDiRJH1F3L3Z3G+8dH1N4dd7CsY6I+CXwp15PRpIk9ZTu7k76FPAm8GPgLOC4/M+bgDXARcANwCXA9w5al5IkSZ10d0vMfwJmp5S+XTD2KvC7iGgArkspXRIRg4AvAd/uaiGSJEkfVXe3xJxP7pTrrvwW2PMN1k+R+4oCSZKkHtHdENMCnLKXaafkp+9Z7s7OM0TErRGxMSKWFoz9MiIW52+rI2JxwbQbImJlRKyIiAu62askSfoY6+7upF8BfxsR7cBdwEZgBPBFcsfA3JqfbxqwoovH30bueJqf7xlIKV2+535E/D25M5yIiMnALOB44AjgsYg4NqXU3s2eJUnSx1B3Q8xfAQOBv8vfCv0b8Nf5+0uBpzs/OKX0VESM62rBERHAZcCn80MzyR1/0wysioiVwPSulitJkg493QoxKaXdwJUR8X3gNOBw4G1gQUrp1YL5HjqAXv4IeCel9Fq+HgU8UzB9HR5nI0mS8rq7JQaAfGB59UNn7J4rgDsL6uhini6/wykirgOuAxg7duxBbkuSJJWi/fkW67HA2yml1vz9fUopreluExFRAXyB9x80vA4YU1CPBtbvZZ03AzcD1NXV+WWVkiQdAvZnS8wq4HTgWWA1H/6N1uUH0Md5wCudvq7gfuDfIuIfyB3Ye0y+B0mSpP0KMV8GXi+4f8BbOiLiTuBsYHhErAO+m1K6hdxZSIW7kkgpLYuIOcDLQBtwvWcmSZKkPT40xKSUflZw/7aPsrKU0hV7Gb92L+M3Ajd+lHVKkqSPp+5e7A6AiCiLiCkRcVZE9D/YTUmSJH2YboeYiLge2AC8RO6rBibmx++NiD8/uO1JkiR1rVshJiK+CvwjcC+5C9MVngb9O+BPDlpnkiRJ+9DdLTF/Bfx9Suk64J5O014hv1VGkiSpp3U3xIwHfrOXaTuBwR+pG0mSpP3U3RCzGRi3l2kTgbc+UjeSJEn7qbtfO/AA8J2IeBJ4Mz+WImI48JfkjpXpVVu2bGHx4sVMmzaN9vZ2br/9dk4++WROOOEEWltbueOOO6irq2PKlCk0NTUxe/Zsfr1tDP369aOjo4ONGzdSUzOIvn370d7ezqZNm6ipqaFv3760t7exadNmBg8eTHV1NW1tbWze/F7d2trKli1bGDJkCH369Hm3Hjp0CFVVfWhpaWHr1q0MHTqUqqoqWlqa2bp1G8OGDaOyspLm5ma2bXuvbmpqYvv27QwfPpyKiop369ra4ZSXV7B7927q6+upra2lvLyc3bt3UV+/gxEjRlBWVsauXbvYseO9eufOnTQ0NHDYYYcREe/WI0ceBgSNjY00NjYycuRIABobG9i5cxeHHXYYAA0NDeza9V69Y8cOmpqaGDFixLt1c3MztbW1ANTX19PS0lJQb6e1tY3hw4cDsH37dtra3qu3bdtGR0cHj337YgAeeeQRWltbueiiiwCYO3cuADNmzADgoYceorKyks985jMAPPDAA/Tt25fzzjsPgHvvvZeamhrOOeccAO6++26GDRvGWWedBcBd85YwcuhAzpg6DoBfPvESY2pr+OSUIwG48/EXGX/4ED4xOXeh6l88+gITx9Ry6nGjAfj5I89z/LjDOOXY3Fd63TZ3EdMmHM60CUfQ3t7B7Y++wMnHHMEJRx9Oa1s7dzy2mLqJo5gyfiRNLa3M/u1LnDZpDJOOHMGuphbmPLmE0yePZeLYWhp3N3PXvKWcMXUcE0YNo76xiXvmL+PME8Zz1BFD2dawm/t+/zJnTzuKcSOHsLl+Jw8+/Qrnnnw0Y0YMZuO2Rh5esILzT5nAqNoaNmxtYO6zrzJj+rGMHDqQtzbV8+iilVx42kRGDBnA2o3befz517n49OMYXtOf1Ru28eTiN5j5qckMGdiXN9Zv5amXVnHJGcdTM6CalW9tYf6S1Vx61hQG9O3DijWbePrlNVx29lT6VVex/M2NLFi+llmfPoHqqkqWrtrAwhVv8aXzplFZUc5Lr7/N86+t56rzT6K8vIzFK9ez+M3buPbaawFYtGgRy5Yt4+qrrwbgueeeY8WKFVx55ZUAPPPMM6xatYorrshdteEPf/gDa9eu5fLLLwdg/vz5bNiwgUsvvRSAefPmsWXLlnc/JxqXN9Kxq4NBpwzK1csa6WjpYNBJubphaQO0wcBpA3P1Sw0ADDwhXy9ugAoYOCVX73hhB2VVZQw4fkCuXrSDsn5lDJiUq+sX1lMxoIL+x+VO5qx/rp6Kmgr6H5uvF9RTObSSfsf0y703ntlO1fAq+k3I13/YTtXIKvodlau3/X4b1aOq6Tuub67+3Taqx1bT98i+pI7E9t9vp++4vlSPqSa1JbY/vZ2+4/tSPbqajpYO6hfU0/fovlQfUU1Hcwf1z9bTb0I/+hzeh/amdnY8t4N+x/ajz2F9aN/Vzo5FO+g/sT9VI6po39nOjud30P+4/lTVVtHW0EbD4gb6T+5P1bAq2na00fBiAwOmDKBySCWt21tpXNLIgKkDqBxcSeu2VhqXNjLwxNy/3dq1a3n88ce5+OKLGT58OKtXr+bJJ59k5syZDBkyhDfeeIOnnnqKSy65hJqaGlauXMn8+fO59NJLGTBgACtWrODpp5/msssuo1+/fixfvpwFCxYwa9YsqqurWbp0KffXb+fiQTVURrCiqYmXm5v43KAayiNY3tTEK81NXFIzGIBlTU2sbG5mZk0NAEt272ZVSwufy9cv7t7NutZWLhqUe628sHsXG1rb+ON8vWjXLja3t3HBwFz93K5dbG9v5/yBud93wa6dNLR3cF6+/sPOnTSnDs4ZkKt/v3MnbSlx1oDca+d3OxsB+KP+uXpeYyMVEXyqf+6180RjA32ijE/m68caGhhYXsZp/XL1ow0NDC4v59R+udfObxp2MLy8glPy9b/v2MHIygpO6purH9qxg9GVlZzYN/faur++nvFVVUzN1/fV1zOhTx+Or64G4J767RzXp5pJ1dW0p8T9O+qZ3KeaidXVtKbEgzvqmVLdl2P69KG5o4OHG3bA8uVMmjSJXbt2MWfOHE4//XQmTpxIY2Mjd911F2eccQYTJkygvr6ee+65hzPPPJO96W6I+W/kvmV6KbCA3IXv/n/gOGAj8P1uLq8kvLhuO2WVTaSODjp27SR2pHzdXlDvfq9ugLKKqnfrsgaIiipSexsduz9Yr24IoqKS1N6aqxuDKK8ktbXS0bST1Y1lRHkFqa2l63pnOVFWTkdbC6mwbm0mNe/qot5KlJXR0dpEat7N6l1biSijo6WJ1LKnDjpadpNamj5Qv7l7K0Cubm1i9a7KXN28i9TW3KluYdXOiny9k9TWyhuN5e/V7W3v1U2NpI52Xm8o67Jub2qE1FG8J16SlGmRUvcuwBsRA4FvAhcAI4AtwFzgRymlHQe7we6qq6tLCxcu7NZjxn3rQL50Wz1l9Q8v6vmVPPGDnl+HuuecG3p8FVN/NrXH16H9t+SaJUVZz/LjJhVlPdo/k15Z3u3HRMSilFJd5/Fuf4t1SqkB+O/5myRJUq/odoiJiGuAK4CxQHWnySmldPTBaEySJGlfuhViIuJvgL8ld0zMYqC5B3qSJEn6UN3dEvMV4B9TSn/ZE81IkiTtr+5eJ2YYudOsJUmSelV3Q8w84MSeaESSJKk7PnR3UkQUBp1vAndHxBbgYWBr5/lT8kIfkiSp5+3PMTFt5C5qt0cAP93LvGk/lylJkvSR7E/g+D7vDzGSJEm97kNDTErpe0XoQyqueT/s7Q7UWRGu2Cvp46W7B/ZKkiSVBEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKJEOMJEnKpKKGmIi4NSI2RsTSTuPfiIgVEbEsIv6uYPyGiFiZn3ZBMXuVJEmlraLI67sN+DHw8z0DEXEOMBM4IaXUHBEj8uOTgVnA8cARwGMRcWxKqb3IPUuSpBJU1C0xKaWngK2dhv8M+GFKqTk/z8b8+ExgdkqpOaW0ClgJTC9as5IkqaSVwjExxwJ/FBELImJeRJyaHx8FrC2Yb11+TJIkqei7k7pSAQwBPgGcCsyJiKOA6GLe1NUCIuI64DqAsWPH9lCbkiSplJTClph1wN0p51mgAxieHx9TMN9oYH1XC0gp3ZxSqksp1dXW1vZ4w5IkqfeVQoi5F/g0QEQcC1QBm4H7gVkR0ScixgPHAM/2VpOSJKm0FHV3UkTcCZwNDI+IdcB3gVuBW/OnXbcA16SUErAsIuYALwNtwPWemSRJkvYoaohJKV2xl0lX7mX+G4Ebe64jSZKUVaWwO0mSJKnbDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTDDGSJCmTihpiIuLWiNgYEUsLxr4XEW9FxOL87cKCaTdExMqIWBERFxSzV0mSVNqKvSXmNmBGF+M/SilNy98eBoiIycAs4Pj8Y26KiPKidSpJkkpaUUNMSukpYOt+zj4TmJ1Sak4prQJWAtN7rDlJkpQppXJMzNcj4qX87qYh+bFRwNqCedblxyRJkkoixPwzcDQwDXgb+Pv8eHQxb+pqARFxXUQsjIiFmzZt6pEmJUlSaen1EJNSeiel1J5S6gB+wnu7jNYBYwpmHQ2s38sybk4p1aWU6mpra3u2YUmSVBJ6PcRExOEF5SXAnjOX7gdmRUSfiBgPHAM8W+z+JElSaaoo5soi4k7gbGB4RKwDvgucHRHTyO0qWg18DSCltCwi5gAvA23A9Sml9mL2K0mSSldRQ0xK6Youhm/Zx/w3Ajf2XEeSJCmren13kiRJ0oEwxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwyxEiSpEwqaoiJiFsjYmNELO1i2n+KiBQRwwvGboiIlRGxIiIuKGavkiSptBV7S8xtwIzOgxExBjgfWFMwNhmYBRyff8xNEVFenDYlSVKpK2qISSk9BWztYtKPgP8CpIKxmcDslFJzSmkVsBKY3vNdSpKkLOj1Y2Ii4nPAWymlFztNGgWsLajX5ce6WsZ1EbEwIhZu2rSphzqVJEmlpFdDTET0A/4r8J2uJncxlroYI6V0c0qpLqVUV1tbezBblCRJJaqil9d/NDAeeDEiAEYDz0fEdHJbXsYUzDsaWF/0DiVJUknq1S0xKaUlKaURKaVxKaVx5ILLySmlDcD9wKyI6BMR44FjgGd7sV1JklRCin2K9Z3A08DEiFgXEV/Z27wppWXAHOBlYC5wfUqpvTidSpKkUlfU3UkppSs+ZPq4TvWNwI092ZMkScqmXj87SZIk6UAYYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiYZYiRJUiZFSqm3ezioImIT8GZv99FLhgObe7sJFZ3P+6HH5/zQdCg/70emlGo7D37sQsyhLCIWppTqersPFZfP+6HH5/zQ5PP+Qe5OkiRJmWSIkSRJmWSI+Xi5ubcbUK/weT/0+JwfmnzeO/GYGEmSlEluiZEkSZlkiMmQiLgtIlYX1OMiIkXEtb3XlXpbRFybfx2MKxhbHRG39V5X2l8R8b3881dRSsvSR+fz0fP8h822t4HTgdd7uxFJkorNEJNhKaVm4Jne7kOSVHwR0Sf/d6DzeACVKaWWXmirqNydVAIiYkJE3B4RqyJid0S8ERH/HBFDPuRxXe5OioizIuLRiKiPiJ0R8WJEfKXTPF/NjzdFxOaIuCUihvbAr6d9iIhjI+KeiNiYfy7WRMSvIqIiIqoj4kcRsTQiGiNiQ0Q8EBHHdXMdp+RfJzO7mHZbRKyLiPKD91vpAE2KiCciYldEvB0R34+Idz+jI2Ji/rWyPf858UxEzDiQZanoxkfEQ/n38ZsR8Z09z8f+vs8Ldhufmf+M2A4syE9bHRG/iIgvR8QrQAtwSURsiogfdW6mYFnd+iwpRb6oS8MRwDrgm8AFwPeBc4GHu7ug/B+qx4Eq4GvATOBW4MiCeX4I3AQ8BnwO+M/ADODf/WNWdA8Co4A/I/fcfwtoJvfe7AMMBP4HcFF+nmrgmYgYub8rSCktAp4j93p4V0QMBi4D/jWl1P5RfxF9ZPeSe09+Hvg34G+A7wBExBHAfOBE4OvknrftwEMR8cfdWZZ6xT3Ab8k9H/cCfwtck5/W3ff5HcAq4FJynxd7nAP8VX7ZM4CFwE+BayKiutMyvgbMSym98hF/r96XUvJWYjdyu/nOABJwUsH4bcDqgnpcfp5r83UAq8m9eMv2suxxQDvwnU7jn8ov6/O9/fsfKjdy34OSgM/t5/zlQD+gAfjLgvFr88sZVzC2Grit0zzt5L5/ZM/YnwNtwOje/rc4lG/A9/LP37c6jf8k/1wPBv5X/rma0On1sAJ4vjvL6u3f91C6FTwff9ppfAnwyF4e82Hv8x918ZjVwC5gZKfx8fn3/VUFYyfklzOrt/99DsbNLTElICKqIuLbEfFKROwGWoHf5SdP7MaiJpLb4vKvKaWOvcxzPrn/5d+R32VRkT9yfgGwAzjzwH4LHYAtwBvAD/O7947pPENEXBYRC/KbjtuAncAAuve6AJhN7n/uXy0Y+xrwUEpp3QH0roNvTqd6Nrnnegq59+UzKaWVeyam3NazO4FpETGoG8tS8T3UqV4KjN1TdPN9fs9e1vFMSmlD4UBKaRXwG96/FfZrwCbg7u78AqXKEFMafkAusf+C3ObE6cAX8tM6bwbcl2H5n/v6ozQi/3MlubBUeBtUsAz1sJT7b9H55Lac/QB4NX881J8BRMRngV8Cy4H/DzgNOJXcB1B3XheklJrIbVr+Sj64/hEwGfiXg/Tr6KN7Zy/1KGAoubMRO9tAbgts5+Pn9rUsFd/WTnUz+ffwAbzPu3od7Gv8JuBTETElIvoDVwI/TR+Tg349O6k0zAJ+nlL6H3sGImLAASxnz1e07+uDakv+52eAbfuYriJIKb0BXJ0/m2DP8Q43Re56QLOAlSmla/fMHxGV5P6gHYh/JrfPfCZwCblN0L850N510B1GbstcYQ3wFrk/gl0dHzGS3K6Bzn8k97UslZbuvs/3dpn9vY0/TO69/jXgRXLH33xsvr7ALTGloR+5LSGF/vQAlvMquRfrf8j/UezKo0AHMDaltLCL26oDWK8+opSzmFzIgNxm/37kNi0XuorcPvMDWcfrwCPkDuS+FPjJPnY7qvgu61TPAhrJ7XqYB3wi3n9Bw3LgcuCFlFJDN5al0nJQ3+ed5d/j/ze/zK8Dj+U/Cz4W3BJTGuaSO4J8CbndPF8APtndhaSUUkR8k9y+zt9GxL+Q2yQ5CRiRUvpuSun1iPifwI8jYiK5D8cmYAy5XRv/mlJ64mD8Utq3iDgB+Edym5JXkvvQupbcB9pvyR2j9Pn8KZIPAqeQOxh3+0dY7U3AfeRC860fYTk6+L6aP+32OXJnqv0H4Hsppe3518C1wKMR8V1yr43/CBxLbhf0fi+rx38LdddcDv77vLNbyB2ycCLwJwdxub3OEFMavkFuv/aN+fph4Arg2e4uKKV0X0ScT+6Uylvyw68D/7tgnm9HxHLg+vwtAWvJnZr92oH9CjoAG4A15La+jCYXJpcAF6eUFkXEC+TC5ZfJbQp+Dvgsez+wb388RO4shoc7HwSoXjcT+Cdy7916cqfc/neAlNL6iDgD+J/kdgv2ARYDF6WU5nZnWSo5P+Hgv8/fJ6W0KSLmAVOB+w/WckuB32ItHULyAfcR4LyU0uO93Y+knhe5C6euAf53Sulverufg8kQIx0CIuJo4CjgR0BzSumUXm5JUg+LiFpyp2n/BXAhuesM7e0spkzywF7p0PA3wL+TO7Xz6l7uRVJxXETummPTgWs+bgEG3BIjSZIyyi0xkiQpkwwxkiQpkwwxkiQpkwwxkiQpkwwxkiQpkwwxkiQpk/4fsxLBvqS9OIsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_family_heights(cond_means, cond_stds)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now what has changed? We see that alice's and bob's height are fixed to the linked data. Their children are taller than both of them, since we gave both alice and bob positive genetic factors. \n", "The standard deviation of of the heights of harry and sally was reduced, too. Since we now added information about their parents." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### OK, now let's make this a little bit more interesting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In reality we would not know the genetic_factor of both alice and bob, but just their heights. This time lets make alice and bob taller than average:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "post_model = hal.get_posterior_model(family, data={family.alice.height: 172,\n", " family.bob.height: 180})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This time we used a posterior model. This model can solve any dependency, be it forward, backward or other." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'alice': array([172.]),\n", " 'bob': array([180.]),\n", " 'harry': array([178.75055174]),\n", " 'sally': array([170.74982655])}" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "post_means = post_model.get_means({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,}, n_samples=100)\n", "post_means" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'alice': array([0.]),\n", " 'bob': array([0.]),\n", " 'harry': array([8.89468819]),\n", " 'sally': array([8.36109147])}" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "post_stds = post_model.get_standard_deviations({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,}, n_samples=100)\n", "post_stds" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAFuCAYAAABwX0lLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlWElEQVR4nO3de5hdd13v8feXVsvhJmqDSi+2YikUkHqI9cJRQMVWUSte06PYAlIeT0HBy3koKjftA/qo6AGLFimLnmIv1LS0Boc0kJO21KTJbtJk0mnaaZI2l6a5dDKZyWSSyeR7/lgrM5thkmanM3vv1bxfz7Of7O9vrb32d2atveeTddk7MhNJkqS6eU6nG5AkSToehhhJklRLhhhJklRLhhhJklRLhhhJklRLhhhJklRLbQ0xEXFGRCyOiL6IWBsRf1iN/0ZVH4qIuVMec2VE9EfEuoi4sJ39SpKk7nVym5/vIPDHmXl/RLwQaETEnUAv8KvAvzTPHBHnAfOAVwEvBRZFxMszc7zNfUuSpC7T1j0xmflEZt5f3R8C+oDTMrMvM9dN85CLgRszc39mbgD6gQva17EkSepWHTsnJiLOAn4YWHaU2U4DNjXVm6sxSZJ0gmv34SQAIuIFwL8D78vMPUebdZqxb/mehIi4HLgc4PnPf/7rXvGKV8xIn5IkqfMajcbOzJwzdbztISYivo0ywHwxM+c/zeybgTOa6tOBrVNnysxrgGsA5s6dmytWrJihbiVJUqdFxGPTjbf76qQAPgf0ZebfH8NDbgfmRcQpEXE2cA5w32z2KEmS6qHde2JeD7wNWBMRq6qxDwKnAJ8C5gALImJVZl6YmWsj4mbgQcorm67wyiRJkgRtDjGZeQ/Tn+cCcOsRHnMVcNWsNSVJkmrJT+yVJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm11NYQExFnRMTiiOiLiLUR8YfV+HdFxJ0R8Uj173c2PebKiOiPiHURcWE7+5UkSd2r3XtiDgJ/nJmvBH4MuCIizgM+AHwtM88BvlbVVNPmAa8CLgKujoiT2tyzJEnqQm0NMZn5RGbeX90fAvqA04CLgS9Us30B+JXq/sXAjZm5PzM3AP3ABe3sWZIkdaeOnRMTEWcBPwwsA74nM5+AMugAL6lmOw3Y1PSwzdXY1GVdHhErImLFjh07ZrVvSZLUHToSYiLiBcC/A+/LzD1Hm3WasfyWgcxrMnNuZs6dM2fOTLUpSZK6WNtDTER8G2WA+WJmzq+Gn4yI76umfx+wvRrfDJzR9PDTga3t6lWSJHWvdl+dFMDngL7M/PumSbcDl1b3LwW+3DQ+LyJOiYizgXOA+9rVryRJ6l4nt/n5Xg+8DVgTEauqsQ8CnwBujoh3Ao8DvwGQmWsj4mbgQcorm67IzPE29yxJkrpQW0NMZt7D9Oe5APzMER5zFXDVrDUlSZJqyU/slSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtWSIkSRJtdTWEBMR10bE9ojobRp7bUT8V0SsiYg7IuJFTdOujIj+iFgXERe2s1dJktTd2r0npgAumjL2r8AHMvM1wK3AnwJExHnAPOBV1WOujoiT2teqJEnqZm0NMZl5F/DUlOFzgbuq+3cCv1bdvxi4MTP3Z+YGoB+4oC2NSpKkrtcN58T0Ar9c3f8N4Izq/mnApqb5Nldj3yIiLo+IFRGxYseOHbPWqCRJ6h7dEGLeAVwREQ3ghcCBajymmTenW0BmXpOZczNz7pw5c2apTUmS1E1O7nQDmfkQ8HMAEfFy4C3VpM1M7pUBOB3Y2t7uJElSt+r4npiIeEn173OAPwf+uZp0OzAvIk6JiLOBc4D7OtOlJEnqNm3dExMRNwBvBE6NiM3Ah4EXRMQV1Szzgc8DZObaiLgZeBA4CFyRmePt7FeSJHWvyJz2NJPamjt3bq5YsaLTbUiSpBkSEY3MnDt1vOOHkyRJko6HIUaSJNWSIUaSJNWSIUaSJNVSxz8nRpKk2bTjU5/udAtqMue975mxZbknRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1ZIhRpIk1VJbQ0xEXBsR2yOit2ns/IhYGhGrImJFRFzQNO3KiOiPiHURcWE7e5UkSd2t3XtiCuCiKWN/A3w0M88HPlTVRMR5wDzgVdVjro6Ik9rWqSRJ6mont/PJMvOuiDhr6jDwour+dwBbq/sXAzdm5n5gQ0T0AxcA/9WOXiU9+7y95+2dbkFNPn/R5zvdgmqurSHmCN4HfDUi/pZyz9BPVOOnAUub5ttcjX2LiLgcuBzgzDPPnLVGJdXbiidXdLoFSTOoG07s/X3g/Zl5BvB+4HPVeEwzb063gMy8JjPnZubcOXPmzFKbkiSpm3RDiLkUmF/d/xLlISMo97yc0TTf6UweapIkSSe4bggxW4E3VPd/Gnikun87MC8iTomIs4FzgPs60J8kSepCbT0nJiJuAN4InBoRm4EPA+8C/jEiTgZGqc5tycy1EXEz8CBwELgiM8fb2a+exRZ/vNMdaKo3XdnpDvQstfvWWzvdgprMee97ZmxZ7b466ZIjTHrdEea/Crhq9jqSJEl11Q1XJ0mSNGte/Na3droFzZJuOCdGkiSpZYYYSZJUS4YYSZJUS4YYSZJUS4YYSZJUS4YYSZJUS4YYSZJUS4YYSZJUS4YYSZJUSy2FmIgYj4gLjjDtdRHhdxtJkqS2aHVPTBxl2klAPoNeJEmSjtkxfXdSRDyHyQDznKpu9t+Anwd2zmBvkiRJR/S0ISYiPgx8qCoT+MZRZr96JpqSJEl6OsdyOOn/AR8D/pJyb8y1Vd18+zPgrcAfzkqXLdi1axerVq0CYHx8nKIoWL16NQBjY2MURUFvby8Ao6OjFEVBX18fACMjIxRFwbp16wAYHh6mKAr6+/sBGBwcpCgK1q9fD8DAwABFUbBx40YAdu7cSVEUbNq0CYDt27dTFAVbtmwBYNu2bRRFwbZt2wDYsmULRVGwfft2ADZt2kRRFOzcWe7Q2rhxI0VRMDAwAMD69espioLBwUEA+vv7KYqC4eFhANatW0dRFIyMjADQ19dHURSMjo4C0NvbS1EUjI2NAbB69WqKomB8vDyVadWqVRRFMfG7bDQaXHfddRP18uXLuf766yfqpUuXcsMNN0zU9957LzfddNNEfc8993DLLbdM1EuWLGH+/PkT9eLFi7ntttsm6kWLFnHHHXdM1AsXLmTBggUTdU9PDz09PRP1ggULWLhw4UR9xx13sGjRoon6tttuY/HixRP1/PnzWbJkyUR9y5I13LNm40R90+LV3Nv72ER9w9ceYOmDj0/U19+5kuUPbZ6or1t4P42Ht0zURU+DVf1bARgfP0TR02D1o08AMHZwnKKnQe+Gct2PHhij6GnQ91i57kdGD1D0NFj3+A4Ahvftp+hp0L9lFwCDw6MUPQ3Wb30KgIGhfRQ9DTZuK7eNnYN7KXoabNq+G4DtA8MUPQ227Ci3lW1PDVH0NNj21BAAW3YMUvQ02D5Qbjubtu+m6Gmwc3AvABu3DVD0NBgY2gfA+q1PUfQ0GBwut6X+LbsoehoM79sPwLrHd1D0NBgZPQBA32PbKXoajB4ot7XeDdsoehqMHSy3tdWPPkHR02B8/BAAq/q3tnXbG+4bZk9jz2S9dpg9Kyfrod4hhlYNTdarhxha3VSvGmKod7Les3IPw2uHJ+vGHob7JuvBFYPsfWjvZL18kL0PN9XLBhl5ZGSi3r10NyP9TfW9uxlZP1kPfGOAfRv3TdZ3D7DvsbLOQ8nA3QOMbirXVR6s6s1lfejAobLeWtX7y3r/E+W6HB8dL+snq3qkrA9sL9ft+N6q3lHWB4cOlvWuqt5T1mMD5bof2z1W1rureqCsD+45CLTnfe/GZUsZq97n1m7Zwo3LljJ+qNz21mzezI3Llk78Lh94/HFuum/ZRL3yscf40vLlE/WKDRuY31gxUd+3fj233d+YqJc++ih3rFw5Ud/7yCP8xwOrJup7Hn6Yr6x+YKJesu4hvrpmzUS9uK+PO9f2TtRff/BBvv7ggxP1nWt7WVz9zQL46po1LFn30ET9ldUPcM/DD0/U//HAKu595JGJ+o6VK1n66KMT9W33N7iv+psGML+xghUbNkzUX1q+nJWPTb4v3nTfMh54fPJ98cZlS1mzuXxfHD90iBuXLWVt9TdvbHycG5ctpW9r+b44OjZW1sf5N3c6T7snJjOXAEsAIiKBz2bm1qd7nCRJ0myKzGfXubhz587NFStWPP2MOrEt/ninO9BUb7py1p/iNV94zaw/h47dmkvXPP1MM2DHpz7dlufRsZnz3ve0/JiIaGTm3Knjx3Ri75QFvQG4BDgTeO6UyZmZP9Nyd5IkSS1qKcRExLuBzwC7gEeA/VNnmaG+JEmSjqrVPTF/DPwb8I7MPDAL/UiSJB2TVj/s7jTg8wYYSZLUaa2GmAbwA7PRiCRJUitaDTF/ALwvIn5qNpqRJEk6Vsfyib2b+ObvRPoOYHFEjAADU2bPzPz+GexPkiRpWsdyYu/X8IsdJUlSlzmWT+y9rA19SJIktaTVc2IkSZK6Qqsfdve7R5l8CBgEVmbm5qPM13V+61/+q9MtqMlN7/7xTrcgSaqBVj/srmDy/JjmT+dtHjsUETcBb6/L58ks2/BUp1uQJEktavVw0uuBx4BPA28AXlH9ezXwOPAW4ErgrcBHZqxLSZKkKVrdE/MnwI2Z+cGmsYeBuyNiCLg8M98aES8Cfhv44HQLkSRJeqZa3RPzZspLrqfzdeDwN1jfRfkVBZIkSbOi1T0xB4DXMX2QeV01HcpwtHfqDBFxLfCLwPbMfHU1dhNwbjXLi4HdmXl+Ne1K4J3AOPAHmfnVFvuVprfqi53uQFO96cpOdyCpZloNMV8CPhoR48AtwHbgJcBvUJ4Dc2013/nAumkeX1CeT3Pd4YHM/K3D9yPi7yivcCIizgPmAa8CXgosioiXZ+Z4iz1LkqRnoVZDzB8BLwT+pro1+zfgj6v7vcC3XLecmXdFxFnTLTgiAvhN4KeroYspz7/ZD2yIiH7ggumWK7Xs/N/udAeSpGeopRCTmfuA34mIjwE/Cnwf8ASwLDMfbppvwXH08pPAk5n5SFWfBixtmr4Zz7ORJEmVVvfEAFAFloefdsbWXALc0FTHNPNM+x1OEXE5cDnAmWeeOcNtSZKkbnQs32J9JvBEZo5V948qMx9vtYmIOBn4VcqTgw/bDJzRVJ8ObD3Cc14DXAMwd+5cv6xSkjRh5z/9U6dbUJM5733PjC3rWPbEbAB+HLgP2MjTf6P1ScfRx88CD035uoLbgX+LiL+nPLH3nKoHSZKkYwox7wAebbp/3Hs6IuIG4I3AqRGxGfhwZn6O8iqk5kNJZObaiLgZeBA4CFzhlUmSJOmwpw0xmfmFpvvFM3myzLzkCOOXHWH8KuCqZ/KckiTp2anVT+wFICKeExGvjog3RMTzZ7opSZKkp9NyiImIK4BtwGrKrxo4txq/LSL+YGbbkyRJml5LISYi3gX8I3Ab5QfTNV8GfTfwazPWmSRJ0lG0uifmj4C/y8zLgVunTHuIye9AkiRJmlWthpizgSN9CeNeyi9wlCRJmnWthpidwFlHmHYusOUZdSNJknSMWg0xdwAfiogfaBrLiDgVeD/luTIdtWvXLlatWgXA+Pg4RVGwevVqAMbGxiiKgt7eXgBGR0cpioKxnZsAODQ2yvCaRYztKj9z79CBfWU9UH5Q8KH9e8t697Zy+aPDDK9ZxMHBJ8t6ZLCs9+wo6727y3poZ1kPDzC8ZhHjwwMAHBzaWdZ7d5f1nh1lPTJY1oNPlvXocNn/7m0Mr1nEof17y3pga1kf2FfWuzaX9dhoWe/cxPCaReTBAwAc2LGxrMcPlvX2DWV9qPz4nQNPrmd4zaKJ3+X+bf0M935tsn7iYYbXLp6stzzE3geXTNSjm/vY23dXU72WvQ/dM1k/voaRdfdO1o+tZuThye/z3LdxFSP9yybqhQsXsmDB5Ndw9fT00NPTM1EvWLCAhQsXTtR33HEHixZN9n/bbbexePFkv/Pnz2fJksl+b1myhnvWbJyob1q8mnt7H5uob/jaAyx9cPIDqK+/cyXLH5r8PMbrFt5P4+HJ3F70NFjVX24r4+OHKHoarH70CQDGDo5T9DTo3VBuO6MHxih6GvQ9th2AkdEDFD0N1j1ebjvD+/ZT9DTo37ILgMHhUYqeBuu3PgXAwNA+ip4GG7eV29LOwb0UPQ02bd8NwPaBYYqeBlt2lNvStqeGKHoabHtqCIAtOwYpehpsHyi3rU3bd1P0NNg5WG5bG7cNUPQ0GBgqt631W5+i6GkwOFxuW/1bdlH0NBjetx+AdY/voOhpMDJabmt9j22n6GkwemAMgN4N2yh6GowdLLe11Y8+QdHTYHz8EACr+rdSFMXE77LRaHDddRNfds/y5cu5/vrrJ+qlS5dyww2THy117733ctNNN03U99xzD7fccstEvWTJEubPnz9RD/cNs6exZ7JeO8yelZP1UO8QQ6uGJuvVQwytbqpXDTHUO1nvWbmH4bXDk3VjD8N9k/XgikH2PrR3sl4+yN6Hm+plg4w8MjJR7166m5H+pvre3Yysn6wHvjHAvo37Juu7B9j3WFnnoWTg7gFGN5XrKg9W9eayPnTgUFlvrer9Zb3/iXJdjo+Ol/WTVT1S1ge2l+t2fG9V7yjrg0MHy3pXVe8p67GBct2P7R4r691VPVDWB/eU70ObNm2iKAp27izfJzdu3EhRFAwMlNv2+vXrKYqCwcFyW+7v76coCoaHy9/vunXrKIqCkZHy99PX10dRFIyOlj9fb28vtw7uZizLjzhbNzrKrYO7Ga/qvqo+bO3oKF+ungtgzb593N5UP7BvHwv2TG4rK/eN8J9NdWNkhK8OTdbLR0a4c2hyW1k2spdFTfW9e/eyeHiy/sbevSwZntx27t47zN17J+slw8N8Y+/ktrN4eIh7m+pFQ0MsG5ms7xwaYvnI5Lbz1aE9NJrq/9yzh5X7JusFe/bwwL7Jbev2wUHWNNVfHhxkbfW7Bbh1cDd9VT2eya2Du1lX1WNV/cj+clvaf+hQOX9fHwAjIyMURcG6desAGB4epigK+vv7ARgcHKQoCtavX8+RtBpi/hzYT/kt1YsoP/ju/wB9wDjwsRaXJ0mSdFwis7UP4I2IFwLvAy4EXgLsAnqAT2bmnqM8tC3mzp2bK1asaOkxZ33geL50W7Nl4yfeMvtPsvjjs/8cas2brpz1p3jNF14z68+hY7fm0jVteZ6+V7yyLc+jY/PKh/pafkxENDJz7tTxlr/FOjOHgL+sbpIkSR3RcoiJiEuBS4AzgedOmZyZ+bKZaEySJOloWgoxEfEXwEcpz4lZRXl+jCRJUtu1uifmncA/Zub7Z6MZSZKkY9Xq1UnfTXmZtSRJUke1GmKWAK+djUYkSZJa8bSHkyKiOei8D5gfEbuArwBPTZ0/Mw/NWHeSJElHcCznxByk/FC7wwL4/BHmzWNcpiRJ0jNyLIHjY3xziJEkSeq4pw0xmfmRNvQhSZLUklZP7JUkSeoKhhhJklRLhhhJklRLhhhJklRLhhhJklRLhhhJklRLhhhJklRLhhhJklRLfkWATkxLPtHpDjTVm67sdAeSasY9MZIkqZYMMZIkqZYMMZIkqZYMMZIkqZYMMZIkqZbaGmIi4tqI2B4RvVPG3xsR6yJibUT8TdP4lRHRX027sJ29SpKk7tbuS6wL4NPAdYcHIuJNwMXAD2Xm/oh4STV+HjAPeBXwUmBRRLw8M8fb3LMkSepCbd0Tk5l3AU9NGf594BOZub+aZ3s1fjFwY2buz8wNQD9wQdualSRJXa0bzol5OfCTEbEsIpZExI9U46cBm5rm21yNSZIkdcUn9p4MfCfwY8CPADdHxA8AMc28Od0CIuJy4HKAM888c5balCRJ3aQb9sRsBuZn6T7gEHBqNX5G03ynA1unW0BmXpOZczNz7pw5c2a9YUmS1HndEGJuA34aICJeDnw7sBO4HZgXEadExNnAOcB9nWpSkiR1l7YeToqIG4A3AqdGxGbgw8C1wLXVZdcHgEszM4G1EXEz8CBwELjCK5MkSdJhbQ0xmXnJESb9zhHmvwq4avY6kiRJddUNh5MkSZJaZoiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm11NYQExHXRsT2iOhtGvtIRGyJiFXV7Reapl0ZEf0RsS4iLmxnr5Ikqbu1e09MAVw0zfgnM/P86vYVgIg4D5gHvKp6zNURcVLbOpUkSV2trSEmM+8CnjrG2S8GbszM/Zm5AegHLpi15iRJUq10yzkx74mI1dXhpu+sxk4DNjXNs7kakyRJ6ooQ8xngZcD5wBPA31XjMc28Od0CIuLyiFgRESt27NgxK01KkqTu0vEQk5lPZuZ4Zh4CPsvkIaPNwBlNs54ObD3CMq7JzLmZOXfOnDmz27AkSeoKHQ8xEfF9TeVbgcNXLt0OzIuIUyLibOAc4L529ydJkrrTye18soi4AXgjcGpEbAY+DLwxIs6nPFS0EXg3QGaujYibgQeBg8AVmTnezn4lSVL3amuIycxLphn+3FHmvwq4avY6kiRJddXxw0mSJEnHwxAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqyRAjSZJqqa0hJiKujYjtEdE7zbQ/iYiMiFObxq6MiP6IWBcRF7azV0mS1N3avSemAC6aOhgRZwBvBh5vGjsPmAe8qnrM1RFxUnvalCRJ3a6tISYz7wKemmbSJ4H/DWTT2MXAjZm5PzM3AP3ABbPfpSRJqoOOnxMTEb8MbMnMB6ZMOg3Y1FRvrsamW8blEbEiIlbs2LFjljqVJEndpKMhJiKeB/wZ8KHpJk8zltOMkZnXZObczJw7Z86cmWxRkiR1qZM7/PwvA84GHogIgNOB+yPiAso9L2c0zXs6sLXtHUqSpK7U0T0xmbkmM1+SmWdl5lmUweW/Z+Y24HZgXkScEhFnA+cA93WwXUmS1EXafYn1DcB/AedGxOaIeOeR5s3MtcDNwINAD3BFZo63p1NJktTt2no4KTMveZrpZ02prwKums2eJElSPXX86iRJkqTjYYiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1ZIiRJEm1FJnZ6R5mVETsAB7rdB8dciqws9NNqO1c7yce1/mJ6URe79+fmXOmDj7rQsyJLCJWZObcTveh9nK9n3hc5ycm1/u38nCSJEmqJUOMJEmqJUPMs8s1nW5AHeF6P/G4zk9MrvcpPCdGkiTVkntiJElSLRliaiQiiojY2FSfFREZEZd1rit1WkRcVm0HZzWNbYyIonNd6VhFxEeq9XdyNy1Lz5zrY/b5i623J4AfBx7tdCOSJLWbIabGMnM/sLTTfUiS2i8iTqn+DkwdD+DbMvNAB9pqKw8ndYGI+MGI+L8RsSEi9kXE+oj4TER859M8btrDSRHxhoi4MyIGI2JvRDwQEe+cMs+7qvHRiNgZEZ+LiO+ahR9PRxERL4+IWyNie7UuHo+IL0XEyRHx3Ij4ZET0RsRwRGyLiDsi4hUtPsfrqu3k4mmmFRGxOSJOmrmfSsfplRGxOCJGIuKJiPhYREy8R0fEudW2srt6n1gaERcdz7LUdmdHxILqdfxYRHzo8Po41td502Hjn6reI3YDy6ppGyPi+oh4R0Q8BBwA3hoROyLik1ObaVpWS+8l3ciNuju8FNgMvA+4EPgY8DPAV1pdUPWH6mvAtwPvBi4GrgW+v2meTwBXA4uAXwb+FLgI+E//mLXdfwCnAb9Pue4/AOynfG2eArwQ+CvgLdU8zwWWRsT3HusTZGYDWE65PUyIiBcDvwn8a2aOP9MfRM/YbZSvyV8B/g34C+BDABHxUuAe4LXAeyjX225gQUT8fCvLUkfcCnydcn3cBnwUuLSa1urr/IvABuDXKd8vDnsT8EfVsi8CVgCfBy6NiOdOWca7gSWZ+dAz/Lk6LzO9ddmN8jDf/wAS+OGm8QLY2FSfVc1zWVUHsJFy433OEZZ9FjAOfGjK+OurZf1Kp3/+E+VG+T0oCfzyMc5/EvA8YAh4f9P4ZdVyzmoa2wgUU+YZp/z+kcNjfwAcBE7v9O/iRL4BH6nW3wemjH+2WtcvBv62Wlc/OGV7WAfc38qyOv3znki3pvXx9inja4CFR3jM073OPznNYzYCI8D3Thk/u3rdv61p7Ieq5czr9O9nJm7uiekCEfHtEfHBiHgoIvYBY8Dd1eRzW1jUuZR7XP41Mw8dYZ43U/4v/4vVIYuTqzPnlwF7gJ86vp9Cx2EXsB74RHV475ypM0TEb0bEsmrX8UFgL/ACWtsuAG6k/J/7u5rG3g0syMzNx9G7Zt7NU+obKdf1qylfl0szs//wxCz3nt0AnB8RL2phWWq/BVPqXuDMw0WLr/Nbj/AcSzNzW/NAZm4Avso374V9N7ADmN/KD9CtDDHd4eOUif16yt2JFwC/Wk2buhvwaL67+vdof5ReUv3bTxmWmm8valqGZlmW/y16M+Wes48DD1fnQ/0+QET8EnAT0Af8T+BHgR+hfANqZbsgM0cpdy2/swquPwmcB/zzDP04euaePEJ9GvBdlFcjTrWNcg/s1PPnjrYstd9TU+r9VK/h43idT7cdHG38auD1EfHqiHg+8DvA5/NZctKvVyd1h3nAdZn5V4cHIuIFx7Gcw1/RfrQ3ql3Vvz8HDBxlutogM9cDv1tdTXD4fIero/w8oHlAf2Zednj+iPg2yj9ox+MzlMfMLwbeSrkL+qvH27tm3PdQ7plrrgG2UP4RnO78iO+lPDQw9Y/k0Zal7tLq6/xIH7N/pPGvUL7W3w08QHn+zbPm6wvcE9Mdnke5J6TZ249jOQ9Tbqy/V/1RnM6dwCHgzMxcMc1tw3E8r56hLK2iDBlQ7vZ/HuWu5WZvozxmfjzP8SiwkPJE7l8HPnuUw45qv9+cUs8DhikPPSwBfiy++QMNTwJ+C1iZmUMtLEvdZUZf51NVr/F/qZb5HmBR9V7wrOCemO7QQ3kG+RrKwzy/CvxEqwvJzIyI91Ee6/x6RPwz5S7JVwIvycwPZ+ajEfHXwKcj4lzKN8dR4AzKQxv/mpmLZ+KH0tFFxA8B/0i5K7mf8k3rMso3tK9TnqP0K9Ulkv8BvI7yZNzdz+Bprwa+TBmar30Gy9HMe1d12e1yyivVfg/4SGburraBy4A7I+LDlNvG/wJeTnkI+piXNes/hVrVw8y/zqf6HOUpC68Ffm0Gl9txhpju8F7K49pXVfVXgEuA+1pdUGZ+OSLeTHlJ5eeq4UeBf2ia54MR0QdcUd0S2ER5afYjx/cj6DhsAx6n3PtyOmWYXAP8YmY2ImIlZbh8B+Wu4OXAL3HkE/uOxQLKqxi+MvUkQHXcxcCnKF+7g5SX3P4lQGZujYj/Afw15WHBU4BVwFsys6eVZanrfJaZf51/k8zcERFLgNcAt8/UcruB32ItnUCqgLsQ+NnM/Fqn+5E0+6L84NTHgX/IzL/odD8zyRAjnQAi4mXADwCfBPZn5us63JKkWRYRcygv0/5D4BcoP2foSFcx1ZIn9konhr8A/pPy0s7f7XAvktrjLZSfOXYBcOmzLcCAe2IkSVJNuSdGkiTVkiFGkiTVkiFGkiTVkiFGkiTVkiFGkiTVkiFGkiTV0v8HS3uLbAAHERUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_family_heights(post_means, post_stds)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Interesting... we see that the expected height of harry and sally is taller that average but not quite as tall as their parents. This effect is known as 'regression to the mean'. We can understand it better by looking at the estimated genetic_factors of alice and bob." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'alice genfactor': array([3.75000001]), 'bob genfactor': array([3.75000001])}" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "post_model.get_means({\"alice genfactor\": family.alice.genetic_factor,\n", " \"bob genfactor\": family.bob.genetic_factor,}, n_samples=100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both alice and bob were 5cm taller than the average man or woman. As we can see, 3/4 of this extra height were estimated into the genetic_factor. The rest was explained as a random fluctuation of the height of alice and bob around their genetic predisposition.\n", "\n", "This is because the variance of the genetic factor was exactly 3 times greater that the variance of the actual height around its genetic predisposition (75 vs 25)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consequently, the genetic factor of their children is" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'harry genfactor': array([3.75030801]),\n", " 'sally genfactor': array([3.74999697])}" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "post_model.get_means({\"harry genfactor\": family.make_harry.outputs.son.genetic_factor,\n", " \"sally genfactor\": family.make_sally.outputs.daughter.genetic_factor,}, n_samples=100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And since for them the changes of being taller or shorter than their genetic predisposition are equal, their average tends exactly to their genetic predisposition." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Now lets turn this around. What does the height of the children tell me about the height of the parents." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "post_model2 = hal.get_posterior_model(family, data={family.make_harry.outputs.son.height: 180,\n", " family.make_sally.outputs.daughter.height: 172})" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "{'alice': array([169.72737024]),\n", " 'bob': array([177.72675745]),\n", " 'harry': array([180.]),\n", " 'sally': array([172.])}" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAFuCAYAAABwX0lLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlfElEQVR4nO3de5SddX3v8fdXsHi8tbaEVrkIWkChVnqY0otttRcKradGbbXhnCqoFZfSVqjtWWKroJZV29VKPe3CFgW3FA1QOkAwGCbRGAiZMJnAJJkwJAxJyD2Ty2Tuk0wmv/PH82RmO04uO5nZez/k/VprVvb39zz72d89+9l7PnluO1JKSJIkFc1Lat2AJEnS8TDESJKkQjLESJKkQjLESJKkQjLESJKkQjLESJKkQqpqiImIsyNiYUR0RMTqiPhkPv6+vD4YEQ0T7nNjRHRGxJqIuKKa/UqSpPp1apUf7wDwqZTSUxHxKmB5RMwH2oH3Av9RPnNEXATMAi4GXgcsiIgLUkqjVe5bkiTVmapuiUkpbUspPZXf7gM6gDNTSh0ppTWT3GUmcE9KaV9KaT3QCVxWvY4lSVK9qtkxMRFxLvALwJNHmO1MYFNZvTkfkyRJJ7lq704CICJeCfw3cH1KqfdIs04y9iPfkxAR1wLXArziFa+49E1vetOU9ClJkmpv+fLlu1JKMyaOVz3ERMRLyQLMt1JKjUeZfTNwdll9FrB14kwppduB2wEaGhpSa2vrFHUrSZJqLSJemGy82mcnBXAH0JFS+vIx3GUOMCsiTouI84DzgZbp7FGSJBVDtbfEvA34ALAqItrysc8ApwH/CswA5kZEW0rpipTS6oi4D3iG7Mym6zwzSZIkQZVDTEppMZMf5wLwwGHucwtwy7Q1JUmSCskr9kqSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEKqaoiJiLMjYmFEdETE6oj4ZD7+kxExPyKey/99Tdl9boyIzohYExFXVLNfSZJUv6q9JeYA8KmU0puBXwaui4iLgE8D30spnQ98L6/Jp80CLgauBG6LiFOq3LMkSapDVQ0xKaVtKaWn8tt9QAdwJjAT+GY+2zeBd+e3ZwL3pJT2pZTWA53AZdXsWZIk1aeaHRMTEecCvwA8Cfx0SmkbZEEHOCOf7UxgU9ndNudjE5d1bUS0RkTrzp07p7VvSZJUH2oSYiLilcB/A9enlHqPNOskY+lHBlK6PaXUkFJqmDFjxlS1KUmS6ljVQ0xEvJQswHwrpdSYD++IiNfm018LdOXjm4Gzy+5+FrC1Wr1KkqT6Ve2zkwK4A+hIKX25bNIc4Or89tXAQ2XjsyLitIg4DzgfaKlWv5IkqX6dWuXHexvwAWBVRLTlY58BvgTcFxEfATYC7wNIKa2OiPuAZ8jObLoupTRa5Z4lSVIdqmqISSktZvLjXAB++zD3uQW4ZdqakiRJheQVeyVJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiEZYiRJUiFVNcRExJ0R0RUR7WVjb42I5ohYFREPR8Sry6bdGBGdEbEmIq6oZq+SJKm+VXtLTAm4csLY14FPp5TeAjwA/DVARFwEzAIuzu9zW0ScUr1WJUlSPatqiEkpPQbsmTB8IfBYfns+8If57ZnAPSmlfSml9UAncFlVGpUkSXWvHo6JaQfeld9+H3B2fvtMYFPZfJvzsR8REddGRGtEtO7cuXPaGpUkSfWjHkLMh4HrImI58Cpgfz4ek8ybJltASun2lFJDSqlhxowZ09SmJEmqJ6fWuoGU0rPA7wJExAXAO/NJmxnfKgNwFrC1ut1JkqR6VfMtMRFxRv7vS4C/Bf49nzQHmBURp0XEecD5QEttupQkSfWmqltiImI28A7g9IjYDNwEvDIirstnaQS+AZBSWh0R9wHPAAeA61JKo9XsV5Ik1a9IadLDTAqroaEhtba21roNSZI0RSJieUqpYeJ4zXcnSZIkHQ9DjCRJKiRDjCRJKiRDjCRJKiRDjCRJKqSaX+xOkqrltrbbat2Cynzikk/UugUVnFtiJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIRliJElSIVU1xETEnRHRFRHtZWOXRMTSiGiLiNaIuKxs2o0R0RkRayLiimr2KkmS6lu1t8SUgCsnjP0j8PmU0iXA5/KaiLgImAVcnN/ntog4pWqdSpKkunZqNR8spfRYRJw7cRh4dX77x4Gt+e2ZwD0ppX3A+ojoBC4DmqvRqyTpxeGFD3yw1i2ozOv/864pW1ZVQ8xhXA88GhH/RLZl6Ffz8TOBpWXzbc7HfkREXAtcC3DOOedU3MCt89dWfB9Nnxsuv6DWLUh6ERlctqzWLWia1MOBvR8HbkgpnQ3cANyRj8ck86bJFpBSuj2l1JBSapgxY8Y0tSlJkupJPYSYq4HG/PZ/ke0ygmzLy9ll853F+K4mSZJ0kquHELMVeHt++7eA5/Lbc4BZEXFaRJwHnA+01KA/SZJUh6p6TExEzAbeAZweEZuBm4CPAl+JiFOBYfJjW1JKqyPiPuAZ4ABwXUpptJr9SpKk+lXts5OuOsykSw8z/y3ALdPXkSRJKqp62J0kSZJUMUOMJEkqJEOMJEkqpHq42J0kVcVDnQ/VugWV+cQln6h1Cyo4t8RIkqRCckuMpJPGzJ+dWesWJE0ht8RIkqRCMsRIkqRCMsRIkqRCMsRIkqRCqijERMRoRFx2mGmXRoTfbSRJkqqi0i0xcYRppwDpBHqRJEk6Zsd0inVEvITxAPOSvC73P4DfA3ZNYW+SJEmHddQQExE3AZ/LywQ8cYTZb5uKpiRJko7mWHYn/QD4AvBFsq0xd+Z1+c/fAO8BPjktXVZg9+7dtLW1ATA6OkqpVGLlypUAjIyMUCqVaG9vB2B4eJhSqcSOjc8DsH94iJamRro2rQNg39AALU2N7NzyAgBDA320NDWye9smAAb7emhpamTP9s0A9PfsoaWpke6ubQD0de+mpamRvbt2ANC7ZyctTY307tkJwN5dO2hpaqSvezcA3V3baGlqpL9nDwB7tm+mpamRwb6e7Llt20RLUyNDA30A7NzyAi1NjewbGgCga9M6Wpoa2T88BMCOjc/T0tTIyP59AGxbv5aWpkZGD4wAsHXds7Q0NXJwNDuUacvzHbQ0NY79LjetbWfZ/AfH6o1rVtK6YPyy7Rs62nhq4XfG6vWrn+LpHzwyVq9rb2XFY/PG6s6VLaxc3DRWP9e2lFVPLBir1z61hNXN3x+rm5qamDt37lg9b9485s0bX97cuXNpahpf3sMPP8yCBePLe/DBB1m4cOFY3djYyKJFi8bq+++/n8WLF4/V9957L0uWLBmrZ8+ezdKlS8fqu+++m2XLlo3Vd911F8uXLx+rS6VSxeteR0cHAIODg5RKJdasWQNAf38/pVKJzs5OAHp6eiiVSqxbl62b3d3dlEolNmzYAMCuXbsolUps2pStm11dXZRKJbZs2QLA9u3bKZVKbN++HYAtW7ZQKpXo6uoCYNOmTZRKJXbtyjambtiwgVKpRHd3NwDr1q2jVCrR05Oti52dnZRKJfr7+wFYs2YNpVKJwcFBADo6OiiVSgwPDwPQ3t5OqVRiZCRb91auXEmpVGI0X/fa2toolUpjv8vly5dz1113jdXLli3j7rvvHquXLl3K7Nmzx+olS5Zw7733jtWLFy/m/vvvH6sXLVpEY+P4ur122VpWLFwxVj+79FlWLVo1Vnc0d9D+ePtYvfqJ1ax+YvVY3f54Ox3NHWP1qkWreHbps2P1ioUrWLts7Vjd9r02nmt9bqx+av5TdD7dOf58H13OurZ1Y3Xrd1tZv3L9WN0yt4UXVr8wVj/58JNsfGbjWN08p5lNa7LX/uDoQZrnNLNlbfbaj46M0jynma2dWwEY2TdC85xmtq/L1oX9Q/tpntPMjg3Z59Tw4DDNc5rp2pitG0N9QzTPaWbX5mzdGOwdpHlOM7u3Zp9b/d39NM9ppnt7tq707emjeU4ze7v2AtCzq4fmOc307MrWnb1de2me00zfnuxzrBrr3gM9exlJ2dEOa4aHeaBnL6N53ZHXh6weHuah/LEAVg0NMaesXjE0xNze3rH66aFBvltWLx8c5NG+8XrZ4CDz+/rGX7vBARaU1UsGBljYP14/MTDAovy5ATw+0M/jA+P1ov5+nhgYGKsX9vexpKxe0NfHk4Pj9fy+PpblvxuAR/t6WV5Wf7e3l6eHxuu5vb2sGBoaq+f09LCqrH6op4fV+e8W4IGevXTk9WhKPNCzlzV5PZLXz+3L/gbtO3gwm/84P/cmc9QtMSmlRcAigIhIwNdSSluPdj9JkqTpFCm9uI7FbWhoSK2trRXd59b5a48+k6rmhssvqHULepG6rc093vWkWl8A2fGmN1flcXRs3vxsx9FnmiAilqeUGiaOV/zdSRHxduAq4BzgZRMmp5TSb1fcnSRJUoUqCjER8THgq8Bu4Dlg38RZpqgvSZKkI6p0S8yngG8DH04p7Z+GfiRJko5JpRe7OxP4hgFGkiTVWqUhZjnwhuloRJIkqRKVhpi/AK6PiN+YjmYkSZKO1bFcsXcTP/ydSD8OLIyIQaB7wuwppfT6KexPkiRpUsdyYO/38IsdJUlSnTmWK/ZeU4U+JEmSKlLpMTGSJEl1odKL3X3wCJMPAj3A0ymlzSfUlTTdvvHOWnegiT409+jzSFKZSi92V2L8+Jjyq/OWjx2MiHuBD3k9GdWtFxYffR5JUl2rdHfS24AXgH8D3g68Kf/3NmAj8E7gRuA9wM1T1qUkSdIElW6J+SvgnpTSZ8rG1gKPR0QfcG1K6T0R8Wrg/wCfmWwhkiRJJ6rSLTGXk51yPZnvA4e+wfoxsq8okCRJmhaVhpj9wKWHmXZpPv3QcgcmzhARd0ZEV0S0l43dGxFt+c+GiGgrm3ZjRHRGxJqIuKLCXiVJ0otYpbuT/gv4fESMAvcDXcAZwPvIjoG5M5/vEmDNJPcvkR1Pc9ehgZTSHx+6HRH/THaGExFxETALuBh4HbAgIi5IKY1W2LMkSXoRqjTE/CXwKuAf859y3wY+ld9uB5on3jml9FhEnDvZgiMigPcDv5UPzSQ7/mYfsD4iOoHLJluuJEk6+VQUYlJKQ8CfRMQXgF8CXgtsA55MKa0tm+94Lvjw68COlNJzeX0msLRs+mY8zkaSJOUq3RIDQB5Y1h51xspcBcwuq2OSeSb9DqeIuBa4FuCcc86p+IHvX+61+erJDZdfUOsWJEkFcCzfYn0OsC2lNJLfPqKU0sZKm4iIU4H38sMHDW8Gzi6rzwK2HuYxbwduB2hoaPDLKiVJOgkcy5aY9cCvAC3ABo7+jdanHEcfvwM8O+HrCuYA346IL5Md2Ht+3sOU+6NLz5qOxUqSpGl0LCHmw8DzZbePe0tHRMwG3gGcHhGbgZtSSneQnYVUviuJlNLqiLgPeAY4AFznmUmSJOmQo4aYlNI3y26XTuTBUkpXHWb8msOM3wLcciKPKUmSXpwqvdgdABHxkoj4uYh4e0S8YqqbkiRJOpqKQ0xEXAdsB1aSfdXAhfn4gxHxF1PbniRJ0uQqCjER8VHgK8CDZBemKz8N+nHgD6esM0mSpCOodEvMXwL/nFK6FnhgwrRnybfKSJIkTbdKL3Z3HvDoYaYNAD9xQt1I0jT66oqv1roFlfnEJZ+odQsquEq3xOwCzj3MtAuBLSfUjSRJ0jGqNMQ8DHwuIt5QNpYi4nTgBrJjZWpq9+7dtLW1ATA6OkqpVGLlypUAjIyMUCqVaG9vB2B4eJhSqcSOjdllcPYPD9HS1EjXpnUA7BsaoKWpkZ1bXgBgaKCPlqZGdm/bBMBgXw8tTY3s2Z5do6+/Zw8tTY10d20DoK97Ny1NjezdtQOA3j07aWlqpHfPTgD27tpBS1Mjfd27Aeju2kZLUyP9PXsA2LN9My1NjQz29WTPbdsmWpoaGRroA2DnlhdoaWpk39AAAF2b1tHS1Mj+4SEAdmx8npamRkb27wNg2/q1tDQ1MnpgBICt656lpamRg6PZ5Xe2PN9BS1Pj2O9y09p2ls1/cKzeuGYlrQseGqs3dLTx1MLvjNXrVz/F0z94ZKxe197KisfmjdWdK1tYubhprH6ubSmrnlgwVq99agmrm78/Vjc1NTF37vjXcM2bN49588aXN3fuXJqaxpf38MMPs2DB+PIefPBBFi5cOFY3NjayaNGisfr+Z0ZYvPHAWH1v+whLNo3Xs1eNsHTzeH33yv0s2zJ+qaK7Vuxn+dbxutS2n7btWT16MFFq28/KHVk9MprV7V1ZPXwgqzt2ZvXgSFav2ZXV/fuzunPPQQB6hrN6XXdWdw9l9Ya9Wb1r8CCltv1s6snqroGs3tKb1dv7s3p7f1Zv6c3qroGs3tST1bsGs3rD3qzuHsouC7WuO6t7hrO6c09W9+/P6jW7Rim17WdwJKs7dmb18IGsbu/K6pHRrF65I6tHD2Z12/bsvXrI8uXLueuusS+7Z9myZdx9991j9dKlS5k9e/zSUkuWLOHee+8dqxcvXsz9998/Vi9atIjGxvF1u7+jn97lveP16n56nx6v+9r76GvrG69X9tG3sqxu66OvfbzufbqX/tX94/XyXvo7xuue1h4Gnh0Yr5f1MLC2rH6yh8HnBsfqvUv3MthZVi/Zy+C68br7iW6GNgyN1493M/RCVqeDie7HuxneNJzVB/J6c1Yf3H8wq7fm9b6s3rct+5wYHR7N6h15PZjV+7v2Z/VAXu/M6gN9B7J6d173ZvVId/Y5M7J3JKv35nV3Vh/ozd5bmzZtolQqsWvXLgA2bNhAqVSiu7sbgHXr1lEqlejpyT4HOzs7KZVK9Pdnv981a9ZQKpUYHMx+Px0dHZRKJYaHs+fX3t7OAz17GUn5ujo8zAM9exnN6468PmT18DAP5Y8FsGpoiDll9YqhIeb2jq8rTw8N8t2yevngII/2jdfLBgeZ3ze+rjw5OMCCsnrJwAAL+8frJwYGWNQ/vu48PtDP4wPj9aL+fp4YGF93Fvb3saSsXtDXx5OD4/X8vj6WDY6vO4/29bK8rP5uby9PD43Xc3t7WTE0vm7N6elhVVn9UE8Pq/PfLcADPXvpyOvRlHigZy9r8nokr5/bl61L+w4ezObv6ABgcHCQUqnEmjVrAOjv76dUKtHZ2QlAT08PpVKJdevWcTiVhpi/BfaRfUv1ArIL3/0/oAMYBb5Q4fIkSZKOS6RU2QV4I+JVwPXAFcAZwG5gHnBrSqn3CHetioaGhtTa2lrRfW6dP9XfZakTUZUvgLz5x6f/MVSZm3uOPs8Jess33zLtj6Fjt+rqVVV5nI43vbkqj6Nj8+ZnOyq+T0QsTyk1TByv+FusU0p9wBfzH0mSpJqoOMRExNXAVcA5wMsmTE4ppTdORWOSJElHUlGIiYjPAp8nOyamjez4GEmSpKqrdEvMR4CvpJRumI5mJEmSjlWlZyf9FNlp1pIkSTVVaYhZBLx1OhqRJEmqxFF3J0VEedC5HmiMiN3AI8CeifOnlA5OWXeSJEmHcSzHxBwgu6jdIQF84zDzpmNcpiRJ0gk5lsDxBX44xEiSJNXcUUNMSunmKvQhSZJUkUoP7JUkSaoLhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIhhhJklRIVQ0xEXFnRHRFRPuE8T+PiDURsToi/rFs/MaI6MynXVHNXiVJUn07tcqPVwL+Dbjr0EBE/CYwE/j5lNK+iDgjH78ImAVcDLwOWBARF6SURqvcsyRJqkNV3RKTUnoM2DNh+OPAl1JK+/J5uvLxmcA9KaV9KaX1QCdwWdWalSRJda0ejom5APj1iHgyIhZFxC/m42cCm8rm25yPSZIkVX130mROBV4D/DLwi8B9EfEGICaZN022gIi4FrgW4JxzzpmmNiVJUj2phy0xm4HGlGkBDgKn5+Nnl813FrB1sgWklG5PKTWklBpmzJgx7Q1LkqTaq4ctMQ8CvwX8ICIuAH4M2AXMAb4dEV8mO7D3fKBlOhr4yveem47F6jjdcPkFtW5BklQAVQ0xETEbeAdwekRsBm4C7gTuzE+73g9cnVJKwOqIuA94BjgAXOeZSZIk6ZCqhpiU0lWHmfQnh5n/FuCW6etIkiQVVT0cEyNJklQxQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSokQ4wkSSqkqoaYiLgzIroior1s7OaI2BIRbfnP75dNuzEiOiNiTURcUc1eJUlSfav2lpgScOUk47emlC7Jfx4BiIiLgFnAxfl9bouIU6rWqSRJqmtVDTEppceAPcc4+0zgnpTSvpTSeqATuGzampMkSYVSL8fE/FlErMx3N70mHzsT2FQ2z+Z8TJIkqS5CzFeBNwKXANuAf87HY5J502QLiIhrI6I1Ilp37tw5LU1KkqT6UvMQk1LakVIaTSkdBL7G+C6jzcDZZbOeBWw9zDJuTyk1pJQaZsyYMb0NS5KkulDzEBMRry0r3wMcOnNpDjArIk6LiPOA84GWavcnSZLq06nVfLCImA28Azg9IjYDNwHviIhLyHYVbQA+BpBSWh0R9wHPAAeA61JKo9XsV5Ik1a+qhpiU0lWTDN9xhPlvAW6Zvo4kSVJR1Xx3kiRJ0vEwxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEIyxEiSpEKqaoiJiDsjoisi2ieZ9lcRkSLi9LKxGyOiMyLWRMQV1exVkiTVt2pviSkBV04cjIizgcuBjWVjFwGzgIvz+9wWEadUp01JklTvqhpiUkqPAXsmmXQr8H+BVDY2E7gnpbQvpbQe6AQum/4uJUlSEdT8mJiIeBewJaW0YsKkM4FNZfXmfGyyZVwbEa0R0bpz585p6lSSJNWTmoaYiHg58DfA5yabPMlYmmSMlNLtKaWGlFLDjBkzprJFSZJUp06t8eO/ETgPWBERAGcBT0XEZWRbXs4um/csYGvVO5QkSXWppltiUkqrUkpnpJTOTSmdSxZc/mdKaTswB5gVEadFxHnA+UBLDduVJEl1pNqnWM8GmoELI2JzRHzkcPOmlFYD9wHPAPOA61JKo9XpVJIk1buq7k5KKV11lOnnTqhvAW6Zzp4kSVIx1fzsJEmSpONhiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYVkiJEkSYUUKaVa9zClImIn8EKt+6iR04FdtW5CVefrfvLxNT85ncyv++tTSjMmDr7oQszJLCJaU0oNte5D1eXrfvLxNT85+br/KHcnSZKkQjLESJKkQjLEvLjcXusGVBO+7icfX/OTk6/7BB4TI0mSCsktMZIkqZAMMQUSEaWI2FBWnxsRKSKuqV1XqrWIuCZfD84tG9sQEaXadaVjFRE356/fqfW0LJ04X4/p5y+22LYBvwI8X+tGJEmqNkNMgaWU9gFLa92HJKn6IuK0/O/AxPEAXppS2l+DtqrK3Ul1ICJ+NiL+MyLWR8RQRKyLiK9GxGuOcr9JdydFxNsjYn5E9ETEQESsiIiPTJjno/n4cETsiog7IuInp+Hp6Qgi4oKIeCAiuvLXYmNE/FdEnBoRL4uIWyOiPSL6I2J7RDwcEW+q8DEuzdeTmZNMK0XE5og4ZeqelY7TmyNiYUQMRsS2iPhCRIx9RkfEhfm6sjf/nFgaEVcez7JUdedFxNz8ffxCRHzu0OtxrO/zst3Gv5F/RuwFnsynbYiIuyPiwxHxLLAfeE9E7IyIWyc2U7asij5L6pErdX14HbAZuB64AvgC8NvAI5UuKP9D9T3gx4CPATOBO4HXl83zJeA2YAHwLuCvgSuB7/rHrOq+A5wJfJzstf80sI/svXka8Crg74B35vO8DFgaET9zrA+QUloOLCNbH8ZExE8A7we+nlIaPdEnohP2INl78t3At4HPAp8DiIjXAYuBtwJ/Rva67QXmRsTvVbIs1cQDwPfJXo8Hgc8DV+fTKn2ffwtYD/wR2efFIb8J/GW+7CuBVuAbwNUR8bIJy/gYsCil9OwJPq/aSyn5U2c/ZLv5fg1IwC+UjZeADWX1ufk81+R1ABvIVt6XHGbZ5wKjwOcmjL8tX9a7a/38T5Yfsu9BScC7jnH+U4CXA33ADWXj1+TLObdsbANQmjDPKNn3jxwa+wvgAHBWrX8XJ/MPcHP++n16wvjX8tf6J4B/yl+rn52wPqwBnqpkWbV+vifTT9nr8aEJ46uApsPc52jv81snuc8GYBD4mQnj5+Xv+w+Ujf18vpxZtf79TMWPW2LqQET8WER8JiKejYghYAR4PJ98YQWLupBsi8vXU0oHDzPP5WT/y/9Wvsvi1PzI+SeBXuA3ju9Z6DjsBtYBX8p3750/cYaIeH9EPJlvOj4ADACvpLL1AuAesv+5f7Rs7GPA3JTS5uPoXVPvvgn1PWSv9c+RvS+XppQ6D01M2daz2cAlEfHqCpal6ps7oW4HzjlUVPg+f+Awj7E0pbS9fCCltB54lB/eCvsxYCfQWMkTqFeGmPrw92SJ/W6yzYmXAe/Np03cDHgkP5X/e6Q/Smfk/3aShaXyn1eXLUPTLGX/LbqcbMvZ3wNr8+OhPg4QEX8A3At0AP8b+CXgF8k+gCpZL0gpDZNtWv5IHlx/HbgI+Pcpejo6cTsOU58J/CTZ2YgTbSfbAjvx+LkjLUvVt2dCvY/8PXwc7/PJ1oMjjd8GvC0ifi4iXgH8CfCN9CI56Nezk+rDLOCulNLfHRqIiFcex3IOfUX7kT6oduf//i7QfYTpqoKU0jrgg/nZBIeOd7gtsusBzQI6U0rXHJo/Il5K9gfteHyVbJ/5TOA9ZJugHz3e3jXlfppsy1x5DbCF7I/gZMdH/AzZroGJfySPtCzVl0rf54e7zP7hxh8he69/DFhBdvzNi+brC9wSUx9eTrYlpNyHjmM5a8lW1j/N/yhOZj5wEDgnpdQ6yc/643hcnaCUaSMLGZBt9n852ablch8g22d+PI/xPNBEdiD3HwFfO8JuR1Xf+yfUs4B+sl0Pi4Bfjh++oOEpwB8DT6eU+ipYlurLlL7PJ8rf4/+RL/PPgAX5Z8GLglti6sM8siPIV5Ht5nkv8KuVLiSllCLierJ9nd+PiH8n2yT5ZuCMlNJNKaXnI+IfgH+LiAvJPhyHgbPJdm18PaW0cCqelI4sIn4e+ArZpuROsg+ta8g+0L5PdozSu/NTJL8DXEp2MO7eE3jY24CHyELznSewHE29j+an3S4jO1PtT4GbU0p783XgGmB+RNxEtm58AriAbBf0MS9r2p+FKjWPqX+fT3QH2SELbwX+cAqXW3OGmPrw52T7tW/J60eAq4CWSheUUnooIi4nO6Xyjnz4eeBfyub5TER0ANflPwnYRHZq9nPH9xR0HLYDG8m2vpxFFiZXAf8rpbQ8Ip4mC5cfJtsUvAz4Aw5/YN+xmEt2FsMjEw8CVM3NBP6V7L3bQ3bK7RcBUkpbI+LXgH8g2y14GtAGvDOlNK+SZanufI2pf5//kJTSzohYBLwFmDNVy60Hfou1dBLJA24T8Dsppe/Vuh9J0y+yC6duBP4lpfTZWvczlQwx0kkgIt4IvAG4FdiXUrq0xi1JmmYRMYPsNO1PAr9Pdp2hw53FVEge2CudHD4LfJfs1M4P1rgXSdXxTrJrjl0GXP1iCzDglhhJklRQbomRJEmFZIiRJEmFZIiRJEmFZIiRJEmFZIiRJEmFZIiRJEmF9P8BEsOGixZPQ+YAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "post_means2 = post_model2.get_means({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,}, n_samples=100)\n", "post_stds2 = post_model2.get_standard_deviations({\"alice\": family.alice.height,\n", " \"bob\": family.bob.height,\n", " \"harry\": family.make_harry.outputs.son.height,\n", " \"sally\": family.make_sally.outputs.daughter.height,}, n_samples=100)\n", "plot_family_heights(post_means2, post_stds2)\n", "post_means2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see the regression to the mean also works backwards." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }