patched heritability and varE calculations

cjGO · Jun 14, 2024 · 103b2df · 103b2df
1 parent dfff446
commit 103b2df
Showing 1 changed file with 137 additions and 35 deletions.
diff --git a/nbs/x01_chewc.ipynb b/nbs/x01_chewc.ipynb
@@ -3,7 +3,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "4ab6a3d9",
+   "id": "3444796d",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -43,7 +43,7 @@
     "        random_effects *= scaling_factors\n",
     "        self.effects = random_effects\n",
     "        \n",
-    "        self.intercept = target_mean - founder_mean\n",
+    "        self.intercept = founder_mean - target_mean\n",
     "        \n",
     "    def __matmul__(self,other):\n",
     "        if isinstance(other,Individual):\n",
@@ -52,24 +52,7 @@
     "            return np.sum(np.array([self @ ind + self.intercept for ind in other.individuals]), axis=(1,2))\n",
     "\n",
     "        \n",
-    "class Trait:\n",
-    "    def __init__(self, genome, founder_population, target_mean, target_variance):\n",
-    "        #store attributes\n",
-    "        self.target_mean = target_mean\n",
-    "        self.target_variance = target_variance\n",
-    "        \n",
-    "        #sample initial random effects\n",
-    "        initial_effects = np.random.randn(g.n_chr * g.n_loci)\n",
-    "        initial_effects -= initial_effects.mean()\n",
-    "        #calculate the founder_population mean and var given these effects\n",
-    "        pop_scores = np.sum(founder_population.get_haplo(),axis=1).reshape(founder_population.size, g.n_chr*g.n_loci)\n",
-    "        founder_scores = np.sum(initial_effects*pop_scores,axis=1)\n",
-    "        founder_mean, founder_var = founder_scores.mean(), founder_scores.var()\n",
-    "        scaling_factors = np.sqrt(self.target_variance / founder_var)\n",
-    "        initial_effects *= scaling_factors\n",
-    "        scaled_effects = initial_effects # ADDS A FIXED EFFECT SOMEHOW????\n",
-    "        self.effects = scaled_effects\n",
-    "        self.intercept = target_mean - founder_mean\n",
+    "\n",
     "\n",
     "\n",
     "        \n",
@@ -93,9 +76,6 @@
     "    def get_pheno(self):\n",
     "        return np.array([x.fitness for x in self.individuals])\n",
     "    \n",
-    "    def trial(self, h2):\n",
-    "        [x.phenotype(h2) for x in self.individuals]\n",
-    "\n",
     "    def __getitem__(self, index):\n",
     "        return self.individuals[index]\n",
     "\n",
@@ -191,6 +171,27 @@
     "        shuffled_haplotypes[:, i, :] = np.array(shuffled_chromosome).reshape(1, n_loci)  # Store the shuffled chromosome\n",
     "    return shuffled_haplotypes[0,:,:]\n",
     "\n",
+    "\n",
+    "\n",
+    "class Trait:\n",
+    "    def __init__(self, genome, founder_population, target_mean, target_variance):\n",
+    "        #store attributes\n",
+    "        self.target_mean = target_mean\n",
+    "        self.target_variance = target_variance\n",
+    "        \n",
+    "        #sample initial random effects\n",
+    "        initial_effects = np.random.randn(g.n_chr * g.n_loci)\n",
+    "        initial_effects -= initial_effects.mean()\n",
+    "        #calculate the founder_population mean and var given these effects\n",
+    "        pop_scores = np.sum(founder_population.get_haplo(),axis=1).reshape(founder_population.size, g.n_chr*g.n_loci)\n",
+    "        founder_scores = np.sum(initial_effects*pop_scores,axis=1)\n",
+    "        founder_mean, founder_var = founder_scores.mean(), founder_scores.var()\n",
+    "        scaling_factors = np.sqrt(self.target_variance / founder_var)\n",
+    "        initial_effects *= scaling_factors\n",
+    "        scaled_effects = initial_effects # ADDS A FIXED EFFECT SOMEHOW????\n",
+    "        self.effects = scaled_effects\n",
+    "        self.intercept = target_mean - founder_mean\n",
+    "\n",
     "@patch\n",
     "def __mul__(self:Individual, partner):\n",
     "    if isinstance(partner,Individual):\n",
@@ -200,37 +201,138 @@
     "        progeny = Individual(self.genome, progeny_haplo, self.id, partner.id,source=source, chewc = chewc)\n",
     "        return progeny\n",
     "\n",
-    "@patch\n",
-    "def phenotype(self:Individual, h2):\n",
-    "    breeding_value = chewc.trait @ self\n",
-    "    genetic_variance = np.var(breeding_value, ddof=1)\n",
-    "    environmental_variance = (1 - h2) / h2 * genetic_variance\n",
-    "    phenotype_value = breeding_value + np.random.normal(0, np.sqrt(environmental_variance))\n",
-    "    self.fitness = np.sum(phenotype_value)\n",
     "\n",
     "    \n",
     "@patch\n",
     "def __matmul__(self:Trait,other):\n",
     "    if isinstance(other,Individual):\n",
     "#             print(f' intercept {self.intercept}')\n",
-    "        return self.effects * np.sum(other.haplotype,axis=0).flatten()\n",
+    "        breeding_value = self.effects * np.sum(other.haplotype,axis=0).flatten()\n",
+    "        return breeding_value\n",
     "    else:\n",
     "        print('ffff')\n",
-    "        \n",
+    "@patch\n",
+    "def phenotype(self:Individual, h2, environmental_variance):\n",
+    "    \"\"\"\n",
+    "    Calculate the phenotype for the individual.\n",
+    "\n",
+    "    Args:\n",
+    "        h2 (float): Heritability of the trait.\n",
+    "        environmental_variance (float): The calculated environmental variance for the population.\n",
+    "    \"\"\"\n",
+    "    breeding_value = chewc.trait @ self\n",
+    "    phenotype_value = breeding_value + np.random.normal(0, np.sqrt(environmental_variance))\n",
+    "    self.fitness = np.sum(phenotype_value)\n",
+    "    \n",
+    "@patch\n",
+    "def trial(self:Population, h2):\n",
+    "    \"\"\"\n",
+    "    Simulate phenotypes for all individuals in the population.\n",
+    "\n",
+    "    Args:\n",
+    "        h2 (float): Heritability of the trait.\n",
+    "    \"\"\"\n",
+    "    # 1. Calculate environmental variance at the population level\n",
+    "    population_genetic_variance = np.var(chewc.population.get_pheno(), ddof=1)\n",
+    "    environmental_variance = (1 - h2) / h2 * population_genetic_variance\n",
+    "\n",
+    "    # 2. Phenotype each individual, passing the environmental variance\n",
+    "    for individual in self.individuals:\n",
+    "        individual.phenotype(h2, environmental_variance) \n",
+    "\n",
+    "    def __getitem__(self, index):\n",
+    "        return self.individuals[index]\n",
+    "\n",
+    "    def __repr__(self):\n",
+    "        return f'Population of size: {self.size}'\n",
+    "    \n",
     "    \n",
     "#Define the Simulation Parameters\n",
     "g = Genome(3, 1000)\n",
     "population = Population(g, size=999)\n",
-    "trait = Trait(g, population,0,1)\n",
+    "trait = Trait(g, population,0,1000)\n",
     "\n",
     "#Plug them into ChewC\n",
     "chewc = ChewC()\n",
     "chewc.trait = trait\n",
     "chewc.population = population\n",
-    "chewc.genome = g\n",
+    "chewc.genome = g"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ff8aa8cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chewc.population.trial(1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1b32e029",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "-23.519044978017305"
+      ]
+     },
+     "execution_count": null,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "chewc.population[0].fitness"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fc80ef41",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(array([  8.,  35.,  89., 206., 234., 232., 123.,  53.,  13.,   6.]),\n",
+       " array([-94.98156687, -74.78008894, -54.57861101, -34.37713308,\n",
+       "        -14.17565514,   6.02582279,  26.22730072,  46.42877865,\n",
+       "         66.63025659,  86.83173452, 107.03321245]),\n",
+       " <BarContainer object of 10 artists>)"
+      ]
+     },
+     "execution_count": null,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAh0UlEQVR4nO3df1DUdeLH8ReorJDuEvJLChS1/JE/Mk2kzOrkQvQ0L+9GPaa08/TysKvwzLhKT+c6POu0ybG8ZgrupsyrudIy0yFMrQtNOcnTkFEPQ8MF04MVS0B5f/9o3G+bmILgvpeej5kd5fN5s74/81b36Wf38zHIGGMEAABgkWB/TwAAAOC7CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1mnv7wk0R0NDg8rLy9W5c2cFBQX5ezoAAOASGGN08uRJxcXFKTj4+8+RBGSglJeXKz4+3t/TAAAAzXD48GFde+213zsmIAOlc+fOkr45QKfT6efZAACAS+HxeBQfH+99Hf8+ARko597WcTqdBAoAAAHmUj6ewYdkAQCAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgnfb+ngCA1tX9sXf9PYUmO7R4rL+nAMDPOIMCAACsQ6AAAADr8BYPAOvwthQAzqAAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA67f09ASCQdH/sXX9PAQB+EDiDAgAArEOgAAAA6xAoAADAOgQKAACwTpMCJTs7WzfffLM6d+6s6OhoTZgwQSUlJT5jTp8+rYyMDHXp0kWdOnXSxIkTVVFR4TOmrKxMY8eOVVhYmKKjozV37lydOXPm8o8GAAC0CU0KlC1btigjI0Pbtm1TXl6e6uvrddddd+nUqVPeMY888ojeeecdvfHGG9qyZYvKy8t1zz33ePefPXtWY8eOVV1dnT7++GP97W9/U25urubPn99yRwUAAAJakDHGNPebjx07pujoaG3ZskUjR45UdXW1oqKitGrVKv3sZz+TJO3bt099+/ZVQUGBhg8frvfee08/+clPVF5erpiYGEnSypUrNW/ePB07dkwhISEX/XU9Ho9cLpeqq6vldDqbO32gybjMGBdyaPFYf08BsF5TXr8v6zMo1dXVkqSIiAhJUmFhoerr65WSkuId06dPHyUkJKigoECSVFBQoAEDBnjjRJJSU1Pl8Xi0d+/eRn+d2tpaeTwenwcAAGi7mh0oDQ0Nevjhh3Xrrbeqf//+kiS3262QkBCFh4f7jI2JiZHb7faO+XacnNt/bl9jsrOz5XK5vI/4+PjmThsAAASAZgdKRkaG9uzZo9WrV7fkfBqVlZWl6upq7+Pw4cOt/msCAAD/adat7mfPnq1169Zp69atuvbaa73bY2NjVVdXp6qqKp+zKBUVFYqNjfWO+eSTT3ye79xVPufGfJfD4ZDD4WjOVAEAQABq0hkUY4xmz56tt956S5s2bVJiYqLP/iFDhqhDhw7Kz8/3bispKVFZWZmSk5MlScnJyfrPf/6jyspK75i8vDw5nU7169fvco4FAAC0EU06g5KRkaFVq1Zp7dq16ty5s/czIy6XS6GhoXK5XJo+fboyMzMVEREhp9OpBx98UMnJyRo+fLgk6a677lK/fv107733asmSJXK73XriiSeUkZHBWRIAACCpiYHywgsvSJLuuOMOn+05OTmaNm2aJGnZsmUKDg7WxIkTVVtbq9TUVD3//PPese3atdO6des0a9YsJScn66qrrtLUqVO1aNGiyzsSAADQZlzWfVD8hfugwF+4DwouhPugABd3xe6DAgAA0BoIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1mlyoGzdulXjxo1TXFycgoKCtGbNGp/906ZNU1BQkM9j9OjRPmNOnDih9PR0OZ1OhYeHa/r06aqpqbmsAwEAAG1HkwPl1KlTGjRokFasWHHBMaNHj9bRo0e9j9dee81nf3p6uvbu3au8vDytW7dOW7du1cyZM5s+ewAA0Ca1b+o3pKWlKS0t7XvHOBwOxcbGNrqvuLhYGzZs0I4dOzR06FBJ0vLlyzVmzBg988wziouLa+qUAABAG9Mqn0HZvHmzoqOj1bt3b82aNUvHjx/37isoKFB4eLg3TiQpJSVFwcHB2r59e2tMBwAABJgmn0G5mNGjR+uee+5RYmKiDh48qN///vdKS0tTQUGB2rVrJ7fbrejoaN9JtG+viIgIud3uRp+ztrZWtbW13q89Hk9LTxsAAFikxQNl8uTJ3p8PGDBAAwcOVM+ePbV582aNGjWqWc+ZnZ2thQsXttQUAQCA5Vr9MuMePXooMjJSBw4ckCTFxsaqsrLSZ8yZM2d04sSJC35uJSsrS9XV1d7H4cOHW3vaAADAj1o9UI4cOaLjx4+ra9eukqTk5GRVVVWpsLDQO2bTpk1qaGhQUlJSo8/hcDjkdDp9HgAAoO1q8ls8NTU13rMhklRaWqqioiJFREQoIiJCCxcu1MSJExUbG6uDBw/q0UcfVa9evZSamipJ6tu3r0aPHq0ZM2Zo5cqVqq+v1+zZszV58mSu4AEAAJKacQZl586dGjx4sAYPHixJyszM1ODBgzV//ny1a9dOu3fv1vjx43X99ddr+vTpGjJkiD788EM5HA7vc7z66qvq06ePRo0apTFjxmjEiBF68cUXW+6oAABAQGvyGZQ77rhDxpgL7t+4ceNFnyMiIkKrVq1q6i8NAAB+IPi/eAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGCd9v6eAAC0Bd0fe9ffU2iyQ4vH+nsKwAVxBgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1mlyoGzdulXjxo1TXFycgoKCtGbNGp/9xhjNnz9fXbt2VWhoqFJSUrR//36fMSdOnFB6erqcTqfCw8M1ffp01dTUXNaBAACAtqPJgXLq1CkNGjRIK1asaHT/kiVL9Nxzz2nlypXavn27rrrqKqWmpur06dPeMenp6dq7d6/y8vK0bt06bd26VTNnzmz+UQAAgDalfVO/IS0tTWlpaY3uM8bo2Wef1RNPPKG7775bkvT3v/9dMTExWrNmjSZPnqzi4mJt2LBBO3bs0NChQyVJy5cv15gxY/TMM88oLi7uMg4HAAC0BS36GZTS0lK53W6lpKR4t7lcLiUlJamgoECSVFBQoPDwcG+cSFJKSoqCg4O1ffv2Rp+3trZWHo/H5wEAANquFg0Ut9stSYqJifHZHhMT493ndrsVHR3ts799+/aKiIjwjvmu7OxsuVwu7yM+Pr4lpw0AACwTEFfxZGVlqbq62vs4fPiwv6cEAABaUYsGSmxsrCSpoqLCZ3tFRYV3X2xsrCorK332nzlzRidOnPCO+S6HwyGn0+nzAAAAbVeLBkpiYqJiY2OVn5/v3ebxeLR9+3YlJydLkpKTk1VVVaXCwkLvmE2bNqmhoUFJSUktOR0AABCgmnwVT01NjQ4cOOD9urS0VEVFRYqIiFBCQoIefvhh/fGPf9R1112nxMREPfnkk4qLi9OECRMkSX379tXo0aM1Y8YMrVy5UvX19Zo9e7YmT57MFTwAAEBSMwJl586duvPOO71fZ2ZmSpKmTp2q3NxcPfroozp16pRmzpypqqoqjRgxQhs2bFDHjh293/Pqq69q9uzZGjVqlIKDgzVx4kQ999xzLXA4AACgLQgyxhh/T6KpPB6PXC6Xqqur+TwKrqjuj73r7ykALebQ4rH+ngJ+YJry+h0QV/EAAIAfFgIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYJ0m3+oeaCnclRUAcCGcQQEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAddq39BP+4Q9/0MKFC3229e7dW/v27ZMknT59WnPmzNHq1atVW1ur1NRUPf/884qJiWnpqQAAvkf3x9719xSa7NDisf6eAq6QVjmDcsMNN+jo0aPex0cffeTd98gjj+idd97RG2+8oS1btqi8vFz33HNPa0wDAAAEqBY/gyJJ7du3V2xs7Hnbq6ur9dJLL2nVqlX60Y9+JEnKyclR3759tW3bNg0fPrw1pgMAAAJMq5xB2b9/v+Li4tSjRw+lp6errKxMklRYWKj6+nqlpKR4x/bp00cJCQkqKCi44PPV1tbK4/H4PAAAQNvV4oGSlJSk3NxcbdiwQS+88IJKS0t122236eTJk3K73QoJCVF4eLjP98TExMjtdl/wObOzs+VyubyP+Pj4lp42AACwSIu/xZOWlub9+cCBA5WUlKRu3brp9ddfV2hoaLOeMysrS5mZmd6vPR4PkQIAQBvW6pcZh4eH6/rrr9eBAwcUGxururo6VVVV+YypqKho9DMr5zgcDjmdTp8HAABou1o9UGpqanTw4EF17dpVQ4YMUYcOHZSfn+/dX1JSorKyMiUnJ7f2VAAAQIBo8bd4fve732ncuHHq1q2bysvLtWDBArVr105TpkyRy+XS9OnTlZmZqYiICDmdTj344INKTk7mCh4AAODV4oFy5MgRTZkyRcePH1dUVJRGjBihbdu2KSoqSpK0bNkyBQcHa+LEiT43agMAADgnyBhj/D2JpvJ4PHK5XKqurubzKAEsEO9iCcC/uJNsYGvK6zf/Fw8AALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDrt/T0BtIzuj73r7ykAANBiOIMCAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsE57f08AAIBL1f2xd/09hSY7tHisv6cQkDiDAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrcJlxIwLxMjYAANoSzqAAAADrECgAAMA6vMUDAEArCtSPDfj7DricQQEAANYhUAAAgHUIFAAAYB0CBQAAWMevgbJixQp1795dHTt2VFJSkj755BN/TgcAAFjCb4Hyj3/8Q5mZmVqwYIH+/e9/a9CgQUpNTVVlZaW/pgQAACzht0BZunSpZsyYofvvv1/9+vXTypUrFRYWppdfftlfUwIAAJbwy31Q6urqVFhYqKysLO+24OBgpaSkqKCg4LzxtbW1qq2t9X5dXV0tSfJ4PK0yv4bar1rleQEACBSt8Rp77jmNMRcd65dA+fLLL3X27FnFxMT4bI+JidG+ffvOG5+dna2FCxeetz0+Pr7V5ggAwA+Z69nWe+6TJ0/K5XJ975iAuJNsVlaWMjMzvV83NDToxIkT6tKli4KCgr73ez0ej+Lj43X48GE5nc7WnuoVx/EFNo4vcLXlY5M4vkBn6/EZY3Ty5EnFxcVddKxfAiUyMlLt2rVTRUWFz/aKigrFxsaeN97hcMjhcPhsCw8Pb9Kv6XQ6rVqklsbxBTaOL3C15WOTOL5AZ+PxXezMyTl++ZBsSEiIhgwZovz8fO+2hoYG5efnKzk52R9TAgAAFvHbWzyZmZmaOnWqhg4dqmHDhunZZ5/VqVOndP/99/trSgAAwBJ+C5RJkybp2LFjmj9/vtxut2688UZt2LDhvA/OXi6Hw6EFCxac9xZRW8HxBTaOL3C15WOTOL5A1xaOL8hcyrU+AAAAVxD/Fw8AALAOgQIAAKxDoAAAAOsQKAAAwDptKlCeeuop3XLLLQoLC7vgjdzKyso0duxYhYWFKTo6WnPnztWZM2d8xmzevFk33XSTHA6HevXqpdzc3NaffBNt3rxZQUFBjT527NghSTp06FCj+7dt2+bn2V+a7t27nzf3xYsX+4zZvXu3brvtNnXs2FHx8fFasmSJn2bbNIcOHdL06dOVmJio0NBQ9ezZUwsWLFBdXZ3PmEBevxUrVqh79+7q2LGjkpKS9Mknn/h7Ss2SnZ2tm2++WZ07d1Z0dLQmTJigkpISnzF33HHHeev0wAMP+GnGTfOHP/zhvLn36dPHu//06dPKyMhQly5d1KlTJ02cOPG8m2zaqrG/Q4KCgpSRkSEp8NZt69atGjdunOLi4hQUFKQ1a9b47DfGaP78+eratatCQ0OVkpKi/fv3+4w5ceKE0tPT5XQ6FR4erunTp6umpuYKHkUTmDZk/vz5ZunSpSYzM9O4XK7z9p85c8b079/fpKSkmF27dpn169ebyMhIk5WV5R3z3//+14SFhZnMzEzz2WefmeXLl5t27dqZDRs2XMEjubja2lpz9OhRn8evfvUrk5iYaBoaGowxxpSWlhpJ5v333/cZV1dX5+fZX5pu3bqZRYsW+cy9pqbGu7+6utrExMSY9PR0s2fPHvPaa6+Z0NBQ89e//tWPs7407733npk2bZrZuHGjOXjwoFm7dq2Jjo42c+bM8Y4J5PVbvXq1CQkJMS+//LLZu3evmTFjhgkPDzcVFRX+nlqTpaammpycHLNnzx5TVFRkxowZYxISEnx+L95+++1mxowZPutUXV3tx1lfugULFpgbbrjBZ+7Hjh3z7n/ggQdMfHy8yc/PNzt37jTDhw83t9xyix9nfOkqKyt9jisvL89IMh988IExJvDWbf369ebxxx83b775ppFk3nrrLZ/9ixcvNi6Xy6xZs8Z8+umnZvz48SYxMdF8/fXX3jGjR482gwYNMtu2bTMffvih6dWrl5kyZcoVPpJL06YC5ZycnJxGA2X9+vUmODjYuN1u77YXXnjBOJ1OU1tba4wx5tFHHzU33HCDz/dNmjTJpKamtuqcL1ddXZ2JiooyixYt8m479wK3a9cu/03sMnTr1s0sW7bsgvuff/55c/XVV3vXzhhj5s2bZ3r37n0FZtfylixZYhITE71fB/L6DRs2zGRkZHi/Pnv2rImLizPZ2dl+nFXLqKysNJLMli1bvNtuv/1289BDD/lvUpdhwYIFZtCgQY3uq6qqMh06dDBvvPGGd1txcbGRZAoKCq7QDFvOQw89ZHr27On9R1wgr9t3A6WhocHExsaap59+2rutqqrKOBwO89prrxljjPnss8+MJLNjxw7vmPfee88EBQWZL7744orN/VK1qbd4LqagoEADBgzwuRlcamqqPB6P9u7d6x2TkpLi832pqakqKCi4onNtqrffflvHjx9v9E6848ePV3R0tEaMGKG3337bD7NrvsWLF6tLly4aPHiwnn76aZ+34woKCjRy5EiFhIR4t6WmpqqkpET/+9///DHdy1JdXa2IiIjztgfa+tXV1amwsNDnz1FwcLBSUlKs/3N0KaqrqyXpvLV69dVXFRkZqf79+ysrK0tfffWVP6bXLPv371dcXJx69Oih9PR0lZWVSZIKCwtVX1/vs5Z9+vRRQkJCwK1lXV2dXnnlFf3yl7/0+U9mA3ndvq20tFRut9tnrVwul5KSkrxrVVBQoPDwcA0dOtQ7JiUlRcHBwdq+ffsVn/PFBMT/ZtxS3G73eXeqPfe12+3+3jEej0dff/21QkNDr8xkm+ill15Samqqrr32Wu+2Tp066S9/+YtuvfVWBQcH65///KcmTJigNWvWaPz48X6c7aX57W9/q5tuukkRERH6+OOPlZWVpaNHj2rp0qWSvlmrxMREn+/59npeffXVV3zOzXXgwAEtX75czzzzjHdboK7fl19+qbNnzzb652jfvn1+mlXLaGho0MMPP6xbb71V/fv3927/xS9+oW7duikuLk67d+/WvHnzVFJSojfffNOPs700SUlJys3NVe/evXX06FEtXLhQt912m/bs2SO3262QkJDzPtMXExPj/TszUKxZs0ZVVVWaNm2ad1sgr9t3nVuPxv7cffv1LTo62md/+/btFRERYeV6Wh8ojz32mP785z9/75ji4mKfD3UFsuYc75EjR7Rx40a9/vrrPuMiIyOVmZnp/frmm29WeXm5nn76ab+9wDXl+L4994EDByokJES//vWvlZ2dbe3tm5uzfl988YVGjx6tn//855oxY4Z3u43r90OXkZGhPXv26KOPPvLZPnPmTO/PBwwYoK5du2rUqFE6ePCgevbseaWn2SRpaWnenw8cOFBJSUnq1q2bXn/9dWv/QdYcL730ktLS0hQXF+fdFsjr9kNgfaDMmTPHp3gb06NHj0t6rtjY2POuJDj3afTY2Fjvj9/9hHpFRYWcTucV+cPanOPNyclRly5dLulFKykpSXl5eZczxctyOeuZlJSkM2fO6NChQ+rdu/cF10r6//W80pp6fOXl5brzzjt1yy236MUXX7zo8/t7/S5FZGSk2rVr1+ja+GtdWsLs2bO1bt06bd261edMZWOSkpIkfXNmLNBe6MLDw3X99dfrwIED+vGPf6y6ujpVVVX5nEUJtLX8/PPP9f7771/0zEggr9u59aioqFDXrl292ysqKnTjjTd6x1RWVvp835kzZ3TixAkr19P6QImKilJUVFSLPFdycrKeeuopVVZWek9z5eXlyel0ql+/ft4x69ev9/m+vLw8JScnt8gcLqapx2uMUU5Oju677z516NDhouOLiop8fvNeaZeznkVFRQoODvauXXJysh5//HHV19d7jz0vL0+9e/f229s7TTm+L774QnfeeaeGDBminJwcBQdf/CNh/l6/SxESEqIhQ4YoPz9fEyZMkPTNWyP5+fmaPXu2fyfXDMYYPfjgg3rrrbe0efPm895WbExRUZEkWb9WjampqdHBgwd17733asiQIerQoYPy8/M1ceJESVJJSYnKysqu2N+JLSEnJ0fR0dEaO3bs944L5HVLTExUbGys8vPzvUHi8Xi0fft2zZo1S9I3f2dWVVWpsLBQQ4YMkSRt2rRJDQ0N3jizir8/pduSPv/8c7Nr1y6zcOFC06lTJ7Nr1y6za9cuc/LkSWPM/19mfNddd5mioiKzYcMGExUV1ehlxnPnzjXFxcVmxYoVVl5mfM77779vJJni4uLz9uXm5ppVq1aZ4uJiU1xcbJ566ikTHBxsXn75ZT/MtGk+/vhjs2zZMlNUVGQOHjxoXnnlFRMVFWXuu+8+75iqqioTExNj7r33XrNnzx6zevVqExYWFhCXGR85csT06tXLjBo1yhw5csTnMsdzAnn9Vq9ebRwOh8nNzTWfffaZmTlzpgkPD/e5gi5QzJo1y7hcLrN582afdfrqq6+MMcYcOHDALFq0yOzcudOUlpaatWvXmh49epiRI0f6eeaXZs6cOWbz5s2mtLTU/Otf/zIpKSkmMjLSVFZWGmO+ucw4ISHBbNq0yezcudMkJyeb5ORkP8/60p09e9YkJCSYefPm+WwPxHU7efKk93VNklm6dKnZtWuX+fzzz40x31xmHB4ebtauXWt2795t7r777kYvMx48eLDZvn27+eijj8x1113HZcZXwtSpU42k8x7nrnk3xphDhw6ZtLQ0ExoaaiIjI82cOXNMfX29z/N88MEH5sYbbzQhISGmR48eJicn58oeSBNMmTLlgvckyM3NNX379jVhYWHG6XSaYcOG+VwuaLPCwkKTlJRkXC6X6dixo+nbt6/505/+ZE6fPu0z7tNPPzUjRowwDofDXHPNNWbx4sV+mnHT5OTkNPp79dv/Zgjk9TPGmOXLl5uEhAQTEhJihg0bZrZt2+bvKTXLhdbp3N8LZWVlZuTIkSYiIsI4HA7Tq1cvM3fuXKvvp/FtkyZNMl27djUhISHmmmuuMZMmTTIHDhzw7v/666/Nb37zG3P11VebsLAw89Of/tQnpG23ceNGI8mUlJT4bA/Edfvggw8a/b04depUY8w3lxo/+eSTJiYmxjgcDjNq1Kjzjvv48eNmypQpplOnTsbpdJr777/f+4942wQZY8wVO10DAABwCX5Q90EBAACBgUABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgnf8Duc0E6sIXUbMAAAAASUVORK5CYII=",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
     "\n",
-    "chewc.population.trial(h2=1)"
+    "plt.hist(chewc.population.get_pheno())"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2127a115",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {