{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The Emitter-Detector Problem"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove-cell"
]
},
"source": [
"Think Bayes, Second Edition\n",
"\n",
"Copyright 2021 Allen B. Downey\n",
"\n",
"License: [Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"# If we're running on Colab, install empiricaldist\n",
"# https://pypi.org/project/empiricaldist/\n",
"\n",
"import sys\n",
"IN_COLAB = 'google.colab' in sys.modules\n",
"\n",
"if IN_COLAB:\n",
" !pip install empiricaldist"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"tags": [
"remove-cell"
]
},
"outputs": [],
"source": [
"# Get utils.py\n",
"\n",
"from os.path import basename, exists\n",
"\n",
"def download(url):\n",
" filename = basename(url)\n",
" if not exists(filename):\n",
" from urllib.request import urlretrieve\n",
" local, _ = urlretrieve(url, filename)\n",
" print('Downloaded ' + local)\n",
" \n",
"download('https://github.com/AllenDowney/ThinkBayes2/raw/master/soln/utils.py')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Click here to run this notebook on Colab](https://colab.research.google.com/github/AllenDowney/ThinkBayes2/blob/master/examples/radiation.ipynb)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modeling a radiation sensor\n",
"\n",
"Here's an example from Jaynes, *Probability Theory*, page 168:\n",
"\n",
"> We have a radioactive source ... which is emitting particles of some sort ... There is a rate $p$, in particles per second, at which a radioactive nucleus sends particles through our counter; and each particle passing through produces counts at the rate $\\theta$. From measuring the number {c1 , c2 , . . .} of counts in different seconds, what can we say about the numbers {n1 , n2 , . . .} actually passing through the counter in each second, and\n",
"what can we say about the strength of the source?\n",
"\n",
"I presented a [version of this problem](https://www.greenteapress.com/thinkbayes/html/thinkbayes015.html#sec130) in the first edition of *Think Bayes*, but I don't think I explained it well, and my solution was a bit of a mess.\n",
"In the second edition, I use more NumPy and SciPy, which makes it possible to express the solution more clearly and concisely, so let me give it another try.\n",
"\n",
"As a model of the radioactive source, Jaynes suggests we imagine \"$N$ nuclei, each of which has independently the probability $r$ of sending a particle through our counter in any one second\".\n",
"If $N$ is large and $r$ is small, the number or particles emitted in a given second is well modeled by a Poisson distribution with parameter $s = N r$, where $s$ is the strength of the source. \n",
"\n",
"As a model of the sensor, we'll assume that \"each particle passing through the counter\n",
"has independently the probability $\\phi$ of making a count\".\n",
"So if we know the actual number of particles, $n$, and the efficiency of the sensor, $\\phi$, the distribution of the count is $\\mathrm{Binomial}(n, \\phi)$.\n",
"\n",
"With that, we are ready to solve the problem, but first, an aside: I am not sure why Jaynes states the problem in terms of $p$ and $\\theta$, and then solves it in terms of $s$ and $\\phi$.\n",
"It might have been an oversight, or there might be subtle distinction he intended to draw the reader's attention to.\n",
"The book is full of dire warnings about distinctions like this, but in this case I don't see an explanation.\n",
"\n",
"Anyway, following Jaynes, I'll start with a uniform prior for $s$, over a range of values wide enough to cover the region where the likelihood of the data is non-negligible."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0.0
\n",
"
1
\n",
"
\n",
"
\n",
"
3.5
\n",
"
1
\n",
"
\n",
"
\n",
"
7.0
\n",
"
1
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"0.0 1\n",
"3.5 1\n",
"7.0 1\n",
"Name: , dtype: int64"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"from empiricaldist import Pmf\n",
"\n",
"ss = np.linspace(0, 350, 101)\n",
"prior_s = Pmf(1, ss)\n",
"prior_s.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For each value of $s$, the distribution of $n$ is Poisson, so we can form the joint prior of $s$ and $n$ using the `poisson` function from SciPy.\n",
"I'll use a range of values for $n$ that, again, covers the region where the likelihood of the data is non-negligible."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(350, 101)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.stats import poisson\n",
"\n",
"ns = np.arange(0, 350)\n",
"S, N = np.meshgrid(ss, ns)\n",
"ps = poisson(S).pmf(N)\n",
"ps.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is an array with one row for each value of $n$ and one column for each value of $s$.\n",
"To get the prior probability for each pair, we multiply each row by the prior probabilities of $s$.\n",
"The following function encapsulates this computation and puts the result in a Pandas `DataFrame` that represents the joint prior."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"\n",
"def make_joint(prior_s, ns):\n",
" ss = prior_s.qs\n",
" S, N = np.meshgrid(ss, ns)\n",
" ps = poisson(S).pmf(N) * prior_s.ps\n",
" joint = pd.DataFrame(ps, index=ns, columns=ss)\n",
" joint.index.name = 'n'\n",
" joint.columns.name = 's'\n",
" return joint"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the joint prior: "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
s
\n",
"
0.0
\n",
"
3.5
\n",
"
7.0
\n",
"
10.5
\n",
"
14.0
\n",
"
17.5
\n",
"
21.0
\n",
"
24.5
\n",
"
28.0
\n",
"
31.5
\n",
"
...
\n",
"
318.5
\n",
"
322.0
\n",
"
325.5
\n",
"
329.0
\n",
"
332.5
\n",
"
336.0
\n",
"
339.5
\n",
"
343.0
\n",
"
346.5
\n",
"
350.0
\n",
"
\n",
"
\n",
"
n
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
1.0
\n",
"
0.030197
\n",
"
0.000912
\n",
"
0.000028
\n",
"
8.315287e-07
\n",
"
2.510999e-08
\n",
"
7.582560e-10
\n",
"
2.289735e-11
\n",
"
6.914400e-13
\n",
"
2.087968e-14
\n",
"
...
\n",
"
4.755624e-139
\n",
"
1.436074e-140
\n",
"
4.336568e-142
\n",
"
1.309530e-143
\n",
"
3.954438e-145
\n",
"
1.194137e-146
\n",
"
3.605981e-148
\n",
"
1.088912e-149
\n",
"
3.288229e-151
\n",
"
9.929590e-153
\n",
"
\n",
"
\n",
"
1
\n",
"
0.0
\n",
"
0.105691
\n",
"
0.006383
\n",
"
0.000289
\n",
"
1.164140e-05
\n",
"
4.394249e-07
\n",
"
1.592338e-08
\n",
"
5.609850e-10
\n",
"
1.936032e-11
\n",
"
6.577099e-13
\n",
"
...
\n",
"
1.514666e-136
\n",
"
4.624158e-138
\n",
"
1.411553e-139
\n",
"
4.308354e-141
\n",
"
1.314851e-142
\n",
"
4.012300e-144
\n",
"
1.224230e-145
\n",
"
3.734968e-147
\n",
"
1.139371e-148
\n",
"
3.475357e-150
\n",
"
\n",
"
\n",
"
2
\n",
"
0.0
\n",
"
0.184959
\n",
"
0.022341
\n",
"
0.001518
\n",
"
8.148981e-05
\n",
"
3.844967e-06
\n",
"
1.671955e-07
\n",
"
6.872067e-09
\n",
"
2.710445e-10
\n",
"
1.035893e-11
\n",
"
...
\n",
"
2.412106e-134
\n",
"
7.444895e-136
\n",
"
2.297302e-137
\n",
"
7.087242e-139
\n",
"
2.185939e-140
\n",
"
6.740663e-142
\n",
"
2.078131e-143
\n",
"
6.405469e-145
\n",
"
1.973961e-146
\n",
"
6.081874e-148
\n",
"
\n",
"
\n",
"
3
\n",
"
0.0
\n",
"
0.215785
\n",
"
0.052129
\n",
"
0.005313
\n",
"
3.802858e-04
\n",
"
2.242898e-05
\n",
"
1.170368e-06
\n",
"
5.612188e-08
\n",
"
2.529749e-09
\n",
"
1.087688e-10
\n",
"
...
\n",
"
2.560853e-132
\n",
"
7.990854e-134
\n",
"
2.492573e-135
\n",
"
7.772342e-137
\n",
"
2.422749e-138
\n",
"
7.549543e-140
\n",
"
2.351752e-141
\n",
"
7.323587e-143
\n",
"
2.279925e-144
\n",
"
7.095520e-146
\n",
"
\n",
"
\n",
"
4
\n",
"
0.0
\n",
"
0.188812
\n",
"
0.091226
\n",
"
0.013946
\n",
"
1.331000e-03
\n",
"
9.812677e-05
\n",
"
6.144433e-06
\n",
"
3.437465e-07
\n",
"
1.770824e-08
\n",
"
8.565541e-10
\n",
"
...
\n",
"
2.039079e-130
\n",
"
6.432637e-132
\n",
"
2.028331e-133
\n",
"
6.392751e-135
\n",
"
2.013910e-136
\n",
"
6.341616e-138
\n",
"
1.996049e-139
\n",
"
6.279975e-141
\n",
"
1.974985e-142
\n",
"
6.208580e-144
\n",
"
\n",
" \n",
"
\n",
"
5 rows × 101 columns
\n",
"
"
],
"text/plain": [
"s 0.0 3.5 7.0 10.5 14.0 17.5 \\\n",
"n \n",
"0 1.0 0.030197 0.000912 0.000028 8.315287e-07 2.510999e-08 \n",
"1 0.0 0.105691 0.006383 0.000289 1.164140e-05 4.394249e-07 \n",
"2 0.0 0.184959 0.022341 0.001518 8.148981e-05 3.844967e-06 \n",
"3 0.0 0.215785 0.052129 0.005313 3.802858e-04 2.242898e-05 \n",
"4 0.0 0.188812 0.091226 0.013946 1.331000e-03 9.812677e-05 \n",
"\n",
"s 21.0 24.5 28.0 31.5 ... 318.5 \\\n",
"n ... \n",
"0 7.582560e-10 2.289735e-11 6.914400e-13 2.087968e-14 ... 4.755624e-139 \n",
"1 1.592338e-08 5.609850e-10 1.936032e-11 6.577099e-13 ... 1.514666e-136 \n",
"2 1.671955e-07 6.872067e-09 2.710445e-10 1.035893e-11 ... 2.412106e-134 \n",
"3 1.170368e-06 5.612188e-08 2.529749e-09 1.087688e-10 ... 2.560853e-132 \n",
"4 6.144433e-06 3.437465e-07 1.770824e-08 8.565541e-10 ... 2.039079e-130 \n",
"\n",
"s 322.0 325.5 329.0 332.5 336.0 \\\n",
"n \n",
"0 1.436074e-140 4.336568e-142 1.309530e-143 3.954438e-145 1.194137e-146 \n",
"1 4.624158e-138 1.411553e-139 4.308354e-141 1.314851e-142 4.012300e-144 \n",
"2 7.444895e-136 2.297302e-137 7.087242e-139 2.185939e-140 6.740663e-142 \n",
"3 7.990854e-134 2.492573e-135 7.772342e-137 2.422749e-138 7.549543e-140 \n",
"4 6.432637e-132 2.028331e-133 6.392751e-135 2.013910e-136 6.341616e-138 \n",
"\n",
"s 339.5 343.0 346.5 350.0 \n",
"n \n",
"0 3.605981e-148 1.088912e-149 3.288229e-151 9.929590e-153 \n",
"1 1.224230e-145 3.734968e-147 1.139371e-148 3.475357e-150 \n",
"2 2.078131e-143 6.405469e-145 1.973961e-146 6.081874e-148 \n",
"3 2.351752e-141 7.323587e-143 2.279925e-144 7.095520e-146 \n",
"4 1.996049e-139 6.279975e-141 1.974985e-142 6.208580e-144 \n",
"\n",
"[5 rows x 101 columns]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"joint = make_joint(prior_s, ns)\n",
"joint.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we're ready to compute the likelihood of the data.\n",
"In this problem, it depends only on $n$, regardless of $s$, so we only have to compute it once for each value of $n$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(350,)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from scipy.stats import binom\n",
"\n",
"phi = 0.1\n",
"c = 10\n",
"likelihood = binom(ns, phi).pmf(c)\n",
"likelihood.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is an array of likelihoods, one for each value of $n$.\n",
"To do the Bayesian update, we need to multiply each column in the prior by this array of likelihoods.\n",
"We can do that using the `multiply` method with the `axis` argument."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"posterior_n = marginal(posterior, 1)\n",
"posterior_n.plot()\n",
"posterior_n.mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The posterior mean of $n$ is close to 109, which is consistent with Equation 6.116. "
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(99, 100.0)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_n.idxmax(), c/phi"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The MAP is 99, which is one less than the analytic result in Equation 6.113, which is 100.\n",
"It looks like the posterior probabilities for 99 and 100 are the same, but the floating-point results differ slightly."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5.065392549852277e-16"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_n[99] - posterior_n[100]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Jeffreys prior\n",
"\n",
"Instead of a uniform prior for $s$, we can use a Jeffreys prior, in which the prior probability for each value of $s$ is proportional to $1/s$.\n",
"This has the advantage of \"invariance under certain changes of parameters\", which is \"the only correct way to express complete ignorance of a scale parameter.\" \n",
"However, Jaynes suggests that it is not clear \"whether $s$ can properly be regarded as a scale parameter in this problem.\"\n",
"\n",
"Nevertheless, he suggests we try it and see what happens.\n",
"Here's the Jeffreys prior for $s$."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
probs
\n",
"
\n",
" \n",
" \n",
"
\n",
"
3.5
\n",
"
0.285714
\n",
"
\n",
"
\n",
"
7.0
\n",
"
0.142857
\n",
"
\n",
"
\n",
"
10.5
\n",
"
0.095238
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
"3.5 0.285714\n",
"7.0 0.142857\n",
"10.5 0.095238\n",
"Name: , dtype: float64"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"prior_jeff = Pmf(1/ss[1:], ss[1:])\n",
"prior_jeff.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use it to compute the joint prior of $s$ and $n$, and update it with `c1`."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"joint_jeff = make_joint(prior_jeff, ns)\n",
"posterior_jeff = update(joint_jeff, phi, c1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the marginal posterior distribution of $n$:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"99.99995605790188"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqo0lEQVR4nO3deXRU93338fdXowVJSAhtIJBAQsiAAIcQGQPOAnZsg+Nj4uQ0NmlrN8kJpbZzujfOyWma9unTpumSEzdeaidu4pzUS54mMXXIY/tJjHdsMGbfJMQigdAG2kHb/J4/ZmRPZC0jtNxZPq9zdDRz7+/OfOdi6zP3d3/3d805h4iIxJ8ErwsQERFvKABEROKUAkBEJE4pAERE4pQCQEQkTiV6XcBY5ObmuuLiYq/LEBGJKu+8806Tcy5v8PKoCoDi4mJ2797tdRkiIlHFzE4PtVxdQCIicUoBICISpxQAIiJxKqwAMLMNZnbMzKrM7P4h1puZPRBcv9/MVoase9zMGszs4DCv/Rdm5sws98o/hoiIjNWoAWBmPuBBYCNQDmw2s/JBzTYCZcGfLcDDIet+CGwY5rWLgBuBM2MtXERExiecI4BVQJVzrto51wM8BWwa1GYT8IQL2AlkmVkBgHPuFeDCMK/9HeCvAM1IJyIyxcIJgLlATcjz2uCysbb5LWZ2G3DWObcvjBpERGSChXMdgA2xbPA39nDavN/YLA34OnDTqG9utoVAtxLz5s0brblMoCN1bew41khvv5+K4pmsWZCD2VD/1CISjcIJgFqgKOR5IXDuCtqEKgVKgH3BPyiFwB4zW+WcOx/a0Dn3KPAoQEVFhbqKpkB3Xz/f+MUhnt5d81vLVxVn893NKyiYkepRZSIykcLpAtoFlJlZiZklA3cC2wa12QbcFRwNtBpodc7VDfeCzrkDzrl851yxc66YQICsHPzHX6Zeb7+fLz/xDk/vruEPP7GAPX99I4f/7mb+4fblHK5r49MPvk7NhS6vyxSRCTBqADjn+oD7gOeBI8AzzrlDZrbVzLYGm20HqoEq4DHgnoHtzexJ4E1gkZnVmtmXJvgzyAT61q+O8srxRv7xM8v52sYlZKcnk5acyOevncdPt66hu8/P3Y+/TfvlXq9LFZFxsmi6JWRFRYXTXECTZ2d1M3c+upO718znbzctG7bN5x/bye0fLuRfP/ehKa5QRK6Emb3jnKsYvFxXAgsA/X7HN7cdoig7lfs3Lhm23eoFOWz9RCn/vaeWt6qbp7BCEZloCgAB4Nm9Zzl6vp2vbVxCarJvxLZfub6MuVmp/PWzB+nt909RhSIy0RQAgnOOR1+pZtGsDDYumz1q+9RkH9+8bSnH6zt48m1dxC0SrRQAwmtVTRw9386XPlYS9jj/Ty7Jp2L+TB7ZcYKePh0FiEQjBYDw6CvV5GWksGnFnLC3MTPuvX4h51ov84t3z05idSIyWRQAce5UUyevVjZx95r5pCSO3Pc/2Lqr8lg2N5OHdlTR74+e0WQiEqAAiHO/2HsWM/jsRwrHvK2Zcc+6hZxq7uI3RxsmoToRmUwKgDjmnOPZvee4tiT7iqd3uLF8FrMyU/jJW0PeclREIpgCII7tr23lZFMnn14x4sStI0ryJXDHNfN4+XijpogQiTIKgDj27N5zJPsS2Li8YFyvs3lVEQb8l4aEikQVBUCccs6x/UAd6xblMSM1aVyvVTAjlRuWzOKnu2vp04VhIlFDARCnDte1cb7tMp8snzUhr/fZlYU0dXTz+glNDyESLRQAceql4Kid9YvyJ+T11i/OI3Naoq4JEIkiCoA49ZujDXyocAZ5GSkT8nopiT4+dXUBzx86T1dP34S8pohMLgVAHGru6ObdmhbWL56Yb/8DPr1iLl09/bx4uH5CX1dEJocCIA69fLwR5+D6CQ6Aa4qzmTNjGtv2jnQ3UBGJFAqAOPTK8UZypyezbM6MCX3dhATj5mWzebWqiY5udQOJRDoFQJxxzvFmdTOrF+SQkBDezJ9jcfPS2fT0+Xn5WOOEv7aITCwFQJw51dxFfVs3a0pzJuX1rynOJic9mecPnZ+U1xeRiaMAiDNvBsfpr1kwOQHgSzA+uWQWLx1toLuvf1LeQ0QmRlgBYGYbzOyYmVWZ2f1DrDczeyC4fr+ZrQxZ97iZNZjZwUHb/LOZHQ22/7mZZY3708io3qxuJj8jhZLc9El7j5uXzaK9u483dFGYSEQbNQDMzAc8CGwEyoHNZlY+qNlGoCz4swV4OGTdD4ENQ7z0i8Ay59zVwHHga2MtXsbGOcfO6mbWlOaEfeevK7G2NJf0ZB8vqBtIJKKFcwSwCqhyzlU753qAp4BNg9psAp5wATuBLDMrAHDOvQJcGPyizrkXnHMDQ0V2AmOfkF7G5ERjJ43t3ZPW/TNgWpKP9YvzefFwvW4UIxLBwgmAuUBNyPPa4LKxthnJF4FfDbXCzLaY2W4z293YqJEl47HrVCCHV5VkT/p73bx0Nk0dPbx75uKkv5eIXJlwAmCovoLBX+vCaTP0i5t9HegDfjLUeufco865CudcRV5eXjgvKcPYc/oi2enJk9r/P+DjV+XhSzBeOqY7hYlEqnACoBYoCnleCAy+1DOcNh9gZncDtwK/65xTX8Eke+fMRVbOy5rU/v8BM1KT+Mi8mezQ9QAiESucANgFlJlZiZklA3cC2wa12QbcFRwNtBpodc7VjfSiZrYB+Cpwm3NOt5KaZC1dPVQ3dvLheTOn7D3XLc7j0Lk26tsuT9l7ikj4Rg2A4Ina+4DngSPAM865Q2a21cy2BpttB6qBKuAx4J6B7c3sSeBNYJGZ1ZrZl4KrvgdkAC+a2V4ze2SiPpR80LtnWgBYOZUBcFVgriFdFSwSmRLDaeSc207gj3zoskdCHjvg3mG23TzM8oXhlynjtefMRXwJxoeKJnb+n5EsKchgVmYKO4438LlrikbfQESmlK4EjhPvnL7I4tkZpCWHlfkTwsxYd1U+rx5vole3ihSJOAqAONDvd+yraZnS7p8B6xbl0d7dx57TGg4qEmkUAHHgRGMHnT39rCjKmvL3vq4sl8QE4yWdBxCJOAqAOLC/thVgSvv/B2ROS+Ij82eyQ9cDiEQcBUAcOFDbQnqyj5Lc6Z68/8evyuPo+XYa27s9eX8RGZoCIA7sP9vK0jkz8E3CDWDC8dGFuQC8caLJk/cXkaEpAGJcX7+fw+faWF449d0/A5bNnUHmtEReq1QAiEQSBUCMq2zooLvPz9UeBoAvwVhbmsvrVU1oxg+RyKEAiHEHgieAl8/1LgAAPlqWy7nWy5xs6vS0DhF5nwIgxu0/28L0lESKcyZ/BtCRDJwHeL1K3UAikUIBEOMOnG1j6ZxMEjw6ATxgfk4ac7NSeU0BIBIxFAAxrN/vOHa+jaVzvO3+gcC0EB9dmMsbJ5p1lzCRCKEAiGGnmju53OtnSUGG16UAgfMA7Zf7OHC21etSRAQFQEw7UtcGwJKCTI8rCVhbGrgXsc4DiEQGBUAMO1LXRmKCUTbLmyuAB8uZnkJ5QSavVmpeIJFIoACIYUfq2inNm05Kos/rUt7z0bJc9pxuoaunz+tSROKeAiCGHalri5j+/wHXLcylp9/PrlOaHlrEawqAGNXS1UNd6+WI6f8fsKo4mySfaV4gkQigAIhRhyPsBPCA1GQfH543kzeqmr0uRSTuhRUAZrbBzI6ZWZWZ3T/EejOzB4Lr95vZypB1j5tZg5kdHLRNtpm9aGaVwd9Tf7uqGHakrh2IvACAwGigg+daae3q9boUkbg2agCYmQ94ENgIlAObzax8ULONQFnwZwvwcMi6HwIbhnjp+4FfO+fKgF8Hn8sEOVLXRu70FPIyUrwu5QOuW5iLc7DzpI4CRLwUzhHAKqDKOVftnOsBngI2DWqzCXjCBewEssysAMA59wpwYYjX3QT8KPj4R8Cnr6B+GUYkngAe8KHCLFKTfLx5QgEg4qVwAmAuUBPyvDa4bKxtBpvlnKsDCP7OD6MWCUNvv5/K+g7KI7D7ByA5MYFrSrJ1QZiIx8IJgKFmERs8mUs4ba6ImW0xs91mtruxURcQhaO6sZOefn9E9v8PWFuaQ2VDBw3tl70uRSRuhRMAtUBRyPNC4NwVtBmsfqCbKPh7yLuGO+cedc5VOOcq8vLywihXIm0KiKFcVxqYHlrdQCLeCScAdgFlZlZiZsnAncC2QW22AXcFRwOtBloHundGsA24O/j4buDZMdQtIzhe305igrEgz9t7AIykfE4mmdMSNRxUxEOjBoBzrg+4D3geOAI845w7ZGZbzWxrsNl2oBqoAh4D7hnY3syeBN4EFplZrZl9KbjqW8CNZlYJ3Bh8LhOgsqGDktx0knyRe5mHL8FYvSCHN6p1HkDEK4nhNHLObSfwRz502SMhjx1w7zDbbh5meTNwQ9iVStiqGjoidgRQqLWlObxwuJ6aC10UZad5XY5I3Incr4hyRS739nO6uZOF+ZEfANct1HkAES8pAGLMyaZO/A7K8iNjCuiRLMyfTu70FF7XvEAinlAAxJjKhg6AiLkHwEjMjLWlObxxoplAL6KITCUFQIypauggwaAkN3JHAIVaW5pDY3s3Jxo7vC5FJO4oAGJMVUM7xTnpEXUTmJEMnAd4XcNBRaacAiDGVNZ3sDAK+v8HFGWnUTgzVfcHEPGAAiCG9Pb7OdnUGVUBAIFuoJ3VF+j36zyAyFRSAMSQ082d9PldVJwADrW2NJfWS73vTWEhIlNDARBDKuuDI4Ci4BqAUGtLcwA0O6jIFFMAxJDKhg7MoDQvuo4A8jOnsTB/Om/ogjCRKaUAiCGVDR0UzkwlNTk6RgCFWluaw65TF+jp83tdikjcUADEkMr69qjr/hmwtjSXrp5+9te2eF2KSNxQAMSIvn4/1U2dUTEFxFBWL8jGTNcDiEwlBUCMqLl4iZ4+f9QNAR2QlZbM0jmZuh5AZAopAGJEZX07AGWzorMLCALdQO+eaeFST7/XpYjEBQVAjBiYBC5ajwAgcCK4p9/P7tMXvC5FJC4oAGLEiYYO5syYxvSUsO7xE5GuKc4mMcE0HFRkiigAYkRlQwelUfztHyA9JZEVRVkKAJEpogCIAX6/o6qhI2qHgIZaW5rDgdoW2i73el2KSMxTAMSAsy2XuNTbH3VzAA1l7cJc/A7eqtZ5AJHJFlYAmNkGMztmZlVmdv8Q683MHgiu329mK0fb1sxWmNlOM9trZrvNbNXEfKT4UzVwF7Ao7wIC+PC8LFISEzQcVGQKjBoAZuYDHgQ2AuXAZjMrH9RsI1AW/NkCPBzGtt8G/tY5twL4RvC5XIHKhsAQ0GgeATQgJdHHNcXZulG8yBQI5whgFVDlnKt2zvUATwGbBrXZBDzhAnYCWWZWMMq2DsgMPp4BnBvnZ4lblfUd5GWkkJWW7HUpE2JNaQ5Hz7fT1NHtdSkiMS2cAJgL1IQ8rw0uC6fNSNv+CfDPZlYD/AvwtaHe3My2BLuIdjc2NoZRbvypbOiIie6fAQO3idRRgMjkCicAbIhlg2/dNFybkbb9I+BPnXNFwJ8CPxjqzZ1zjzrnKpxzFXl5eWGUG1+cGxgBFDsBsGxOJhkpiRoOKjLJwgmAWqAo5HkhH+yuGa7NSNveDfws+PinBLqLZIzOt12mo7uPhVE8BcRgib4Erl2QzZs6ESwyqcIJgF1AmZmVmFkycCewbVCbbcBdwdFAq4FW51zdKNueAz4RfHw9UDnOzxKX3r8LWOwcAQCsKc3lVHMXZ1sueV2KSMwadd4A51yfmd0HPA/4gMedc4fMbGtw/SPAduAWoAroAr4w0rbBl/4y8F0zSwQuExg9JGNUGUNDQENdtzBwm8g3qpr4nYqiUVqLyJUIa+IY59x2An/kQ5c9EvLYAfeGu21w+WvAR8ZSrHxQVUM7M9OSyJme4nUpE+qq/Axy0pN580SzAkBkkuhK4CgXK1NADJaQYKwuzeGNE80Evl+IyERTAEQx5xzH6ztYGANTQAzlutJczrddprqp0+tSRGKSAiCKNXX00HqpN+b6/wesLQ2eB9BwUJFJoQCIYgNTQMRiFxDA/Jw05syYpuGgIpNEARDF3psELka7gMyMNaW5vHmiGb9f5wFEJpoCIIpV1neQMS2R/IzYGgEU6rqFOVzs6uXI+TavSxGJOQqAKFbZ0E5Z/nTMhppxIzasCZ4H0LxAIhNPARDFYnUIaKiCGaksyE3XiWCRSaAAiFIXOnto6uiJ2f7/UGtKc3irupnefr/XpYjEFAVAlBo4ARwLN4EZzXULc+ns6Wd/bavXpYjEFAVAlIqlu4CNZvWCgfMAGg4qMpEUAFGqsr6DtGQfc2akel3KpMtOT2ZJQSavV+k8gMhEUgBEqaqGDhbmTychIXZHAIVaW5rDO2cucrm33+tSRGKGAiBKDQRAvLhuYQ49fX52n7rodSkiMUMBEIXaLvdyvu1yzA8BDXVtSQ7JvgRePt7gdSkiMUMBEIWqYvQmMCNJT0lkVUk2O441el2KSMxQAEShqvrYngNoOOsW5VHZ0EHNhS6vSxGJCQqAKFTZ0E5KYgKFM9O8LmVKrVuUD8CO4zoKEJkICoAoVNnQQWnedHxxMgJoQGleOkXZqbx8TOcBRCZCWAFgZhvM7JiZVZnZ/UOsNzN7ILh+v5mtDGdbM/tKcN0hM/v2+D9OfKis74i77h8ITA+9flE+r1c1azioyAQYNQDMzAc8CGwEyoHNZlY+qNlGoCz4swV4eLRtzWw9sAm42jm3FPiXifhAsa6zu4+zLZfi6gRwqHWL8rjU28/bJy94XYpI1AvnCGAVUOWcq3bO9QBPEfjDHWoT8IQL2AlkmVnBKNv+EfAt51w3gHNOx/VhiKc5gIayZkEuyYkJGg0kMgHCCYC5QE3I89rgsnDajLTtVcDHzOwtM3vZzK4ZS+Hx6nh98DaQs+LnGoBQqck+1izIYYfOA4iMWzgBMNSZxsH35xuuzUjbJgIzgdXAXwLP2BB3NjGzLWa228x2NzbqW19lQwfJiQnMz46vEUCh1i/Ko7qpk9PNnV6XIhLVwgmAWqAo5HkhcC7MNiNtWwv8LNht9DbgB3IHv7lz7lHnXIVzriIvLy+McmPb8fp2SvOmk+iL3wFc7w0HVTeQyLiE81dkF1BmZiVmlgzcCWwb1GYbcFdwNNBqoNU5VzfKtr8Argcws6uAZEDz/Y6isr6Dq+JwBFCo4tx0SnLT+c1RdQOJjMeoAeCc6wPuA54HjgDPOOcOmdlWM9sabLYdqAaqgMeAe0baNrjN48ACMztI4OTw3c65wV1LEqIjOALoqjjt/w91/eJ83jzRTEd3n9eliEStxHAaOee2E/gjH7rskZDHDrg33G2Dy3uA3xtLsfGucuAEcJyOAAp1U/ksfvDaSV4+1sinri7wuhyRqBS/HclRqDI4B5COAOAj82eSnZ7Mi4fPe12KSNRSAESRY/WBOYCK4ngE0IBEXwI3LM7n10cbdLN4kSukAIgix+vbWZgff3MADeempbNpv9zHW9W6KljkSigAokhgBJC6fwZ8dGEu05ISeEHdQCJXRAEQJVovBe4CpgB4X2qyj4+X5fHi4Xo0gExk7BQAUaKqITACKN6vARjspqWzqWu9zMGzbV6XIhJ1FABR4rhGAA3phsX5JBjqBhK5AgqAKHG8vp3UJB9zs1K9LiWizExPZlVJNi8cqve6FJGoowCIEgM3gUnQCKAPuLF8Nsfq2znVpMnhRMZCARAljte3U5av7p+h3Lx0FgDbD9Z5XIlIdFEARIHWrl4a2rt1AngYhTPTWFGUxXP7FAAiY6EAiALH3xsBpCOA4dx6dQGH69qobuzwuhSRqKEAiALv3wVMRwDDGZgQ7pf7dRQgEi4FQBQ4fr6d9GSNABpJwYxUKubP5DkFgEjYFABR4Mj5dhbNzmCIO2ZKiFuvLuBYfft702aLyMgUABHOOceRujaWFGR6XUrEu2V5AWbwPzoKEAmLAiDCnW25RPvlPgVAGPIzp3FtSTbP7TunuYFEwqAAiHBH6gLdGQqA8Hx6xVyqmzrZX9vqdSkiEU8BEOGO1LVhBotnawhoOG65uoDkxAR+tqfW61JEIp4CIMIdqWtjfnYa6Slh3b457mVOS+Km8lls23eOnj7dKUxkJGEFgJltMLNjZlZlZvcPsd7M7IHg+v1mtnIM2/6FmTkzyx3fR4lNOgE8dp9dWcjFrl52HGvwuhSRiDZqAJiZD3gQ2AiUA5vNrHxQs41AWfBnC/BwONuaWRFwI3Bm3J8kBnV293H6QpcCYIw+VpZL7vRkfrbnrNeliES0cI4AVgFVzrlq51wP8BSwaVCbTcATLmAnkGVmBWFs+x3grwAN2RjC0fPtOKcTwGOV6Etg04q5/PpoPS1dPV6XIxKxwgmAuUBNyPPa4LJw2gy7rZndBpx1zu0b6c3NbIuZ7Taz3Y2NjWGUGzsO1wXucrWkQCeAx+ozK+fS2+/4n33nvC5FJGKFEwBDXX46+Bv7cG2GXG5macDXgW+M9ubOuUedcxXOuYq8vLxRi40lR+rayJyWqCkgrkB5QSblBZk8+XaNrgkQGUY4AVALFIU8LwQGf60ars1wy0uBEmCfmZ0KLt9jZrPHUnysO1LXxuKCTE0BcQXMjM9fO4/DdW3s0zUBIkMKJwB2AWVmVmJmycCdwLZBbbYBdwVHA60GWp1zdcNt65w74JzLd84VO+eKCQTFSuecbuwa5Pc7jp1vp1z9/1ds04o5pCX7ePItjTEQGcqoAeCc6wPuA54HjgDPOOcOmdlWM9sabLYdqAaqgMeAe0badsI/RQw6faGLrp5+BcA4ZExL4rYPzWHbvnO0Xe71uhyRiBPW1UXOue0E/siHLnsk5LED7g132yHaFIdTRzw5dC7QbVE+RwEwHptXzeOpXTU8++5Zfn9NsdfliEQUXQkcoQ7UtpLsS9BdwMbp6sIZLJ2TyU/eOqOTwSKDKAAi1P7aVpYUZJCcqH+i8Rg4GXz0fDvvnL7odTkiEUV/XSKQ3+84eLaV5YUzvC4lJtz+4bnMSE3iB6+d9LoUkYiiAIhAp5o7ae/uY/lcBcBESEtOZPOqeTx/6Dw1F7q8LkckYigAItCBs4ETwMvnZnlbSAy5e+18Esz40RunvC5FJGIoACLQgdpWUhITKJs13etSYkbBjFRuWV7A07tq6Oju87ockYigAIhA+8+2Uj4nkySf/nkm0hc/WkJ7dx8/3V0zemOROKC/MBGm3+84dLZV/f+TYEVRFhXzZ/L9V0/S26+bxYgoACLM8fp2Onv6WVGU5XUpMene9Qs523KJn7+rewWIKAAizJ4zgbHqK+fN9LiS2LRuUR7L5mby0EtV9Pt1YZjENwVAhNlzuoXs9GTm56R5XUpMMjPuW1/GqeYuntuvewVIfFMARJh3z1xk5byZmgJ6Et1UPotFszL43m+q8OsoQOKYAiCCXOzsobqpk5Xzs7wuJaYlJBj3Xr+QyoYOth+s87ocEc8oACLIuzXq/58qn1pewKJZGfzL88c0IkjilgIgguw53YIvwbhacwBNOl+C8dWNizjV3MVTb+uGMRKfFAAR5J3TF1lSkEFacli3aZBxWr8on2tLsvnuryvp1NXBEocUABGiu6+fPWcucm1JjtelxA0z4/6Ni2nq6OGxV6u9LkdkyikAIsTeMy109/lZvUABMJU+PG8mtyyfzaOvVFPXesnrckSmlAIgQuysvoAZrCrO9rqUuPO1jUvo9zv+/rkjXpciMqXCCgAz22Bmx8ysyszuH2K9mdkDwfX7zWzlaNua2T+b2dFg+5+bWdaEfKIotbO6mfKCTGakJXldStwpyk7j3vUL+eWBOl6tbPS6HJEpM2oAmJkPeBDYCJQDm82sfFCzjUBZ8GcL8HAY274ILHPOXQ0cB7427k8TpS73Bvr/1f3jnS0fX0BxThp/8+whuvv6vS5HZEqEcwSwCqhyzlU753qAp4BNg9psAp5wATuBLDMrGGlb59wLzrmBoRc7gcIJ+DxRaV+N+v+9Ni3JxzdvW0p1Uyf/8bJOCEt8CCcA5gKhE6jXBpeF0yacbQG+CPxqqDc3sy1mttvMdjc2xubh+RsnmtX/HwHWLcrn1qsL+PffVHL4XJvX5YhMunACYKhJaQZPoDJcm1G3NbOvA33AT4Z6c+fco865CudcRV5eXhjlRp8dxxtZUZSl/v8I8HebljEjNZk/e2YvPX26QlhiWzgBUAsUhTwvBAZPozhcmxG3NbO7gVuB33XOxeWsXM0d3eyvbWHdVflelyJAdnoy//iZ5Rw9386//6bS63JEJlU4AbALKDOzEjNLBu4Etg1qsw24KzgaaDXQ6pyrG2lbM9sAfBW4zTnXNUGfJ+q8UtmIc7B+cWwe3USjG8tn8dmVhTy04wS7Tl3wuhyRSTNqAARP1N4HPA8cAZ5xzh0ys61mtjXYbDtQDVQBjwH3jLRtcJvvARnAi2a218wembiPFT12HGskd3oyy+Zo/p9I8je3lVM4M5X7/msPTR3dXpcjMiksmnpeKioq3O7du70uY8L0+x0Vf/8i6xfl8293rPC6HBnk0LlWPvPQG1QUz+SJL16LL0H3aJDoZGbvOOcqBi/XlcAe2ltzkYtdvXxikbp/ItHSOTP4X5uW8XpVM9958bjX5YhMOE076aFf7j9PcmIC1y/WCeBI9blrithz5iLfe6mKktx0PvuRuL1cRWKQAsAjfr/jVwfr+HhZHhnTNPwzkv3dpmWcudDF/T/bT0HWNNaW5npdksiEUBeQR96tuUhd62VuvbrA61JkFMmJCTz8ex+hJDedP/zxOxw73+51SSITQgHgkef215GcmMANS9T9Ew1mpCbx+B9cQ1qyj9/9/k4q6xUCEv0UAB7o9zu2H6hj3VXq/okmhTPT+K8vr8bM2PzYW1Q1dHhdksi4KAA8sONYA/Vt3Xxm5VDTIkkkK82bzpNfXg3A5sd2cuhcq8cViVw5BYAHfvLWGfIyUrhhySyvS5ErsDB/Ok9tuZakBOOO/9ipewhI1FIATLGzLZfYcayBOyqKSPJp90erhfkZ/Oye6yicmcoX/nMXT+8643VJImOmv0BT7Om3z+CAO1cVjdpWItvsGdP46dY1rCnN4av/fYCv/p/9XO7VzWQkeigAplBHdx9P7DzN9YvyKZyZ5nU5MgEypiXxwy+s4ivXL+Tp3TV85qE3qGrQCCGJDgqAKfTEm6do6erlKzeUeV2KTCBfgvHnNy3iP//gGupaL3HLA6/xyMsn6PdHzzxbEp8UAFOks7uPx16p5hNX5bGiKMvrcmQSrF+czwt/+gnWL8rjW786ymceep13z1z0uiyRYSkApsj3Xz3Jxa5e/viT+vYfy/IyUnjk9z7Cd+9cwbnWy9z+0Bv82dN7OddyyevSRD5AcwFNgerGDh7cUcWtVxewct5Mr8uRSWZmbFoxlxuWzOKhl6r4/qsneW5/HZ+7ppB71i1kTlaq1yWKALofwKRzzvH5x97i0LlW/t+ff4L8jGlelyRT7GzLJR56qYpndtfgHNy8bDZ3rZ7PqpJszHSPAZl8w90PQEcAk+yhHSd4s7qZf7h9uf74x6m5Wan879uXc8/6hfznayd5ZncNv9xfx6JZGXzumiI+tbyA2TP034ZMPR0BTKIXDp1ny4/f4dMr5vCdO1bo254AcKmnn//Zd44f7zzNgbOBqSSuKZ7JLcsLWLcon+KcNP23IhNquCMABcAkeeHQeb7y5Lssnp3B03+4hmlJPq9Lkgh0orGDX+6v47n95zheH5hcrig7lY+V5XFtSTYrirKYl61AkPFRAEyRvn4/P3jtJP/0f4+yvDCLx++uIGd6itdlSRQ42dTJa5WNvFLZxJsnmuno7gMgOz2ZDxXOYHlhFmX50ymbNZ2S3HRSEvWlQsIzrnMAZrYB+C7gA77vnPvWoPUWXH8L0AX8gXNuz0jbmlk28DRQDJwCPueci9pB036/4+XKRv71hWMcPNvGzUtn8Z07VpCWrNMsEp6S3HRKctP5/TXF9PX7OV7fwd6aFvbWXGRvTQsvH29k4NqyBIN52WkUZacxZ0Yqc7JSmZM1jblZqeRlpJCdnkxWWrJuZC8jGvUIwMx8wHHgRqAW2AVsds4dDmlzC/AVAgFwLfBd59y1I21rZt8GLjjnvmVm9wMznXNfHamWSDkCcM7R3t3HuZZLVNZ3sOvUBV461kDNhUvMykzhr28t51PLC3TYLhPqcm8/J5s6qWzooKqhgxONHZy9eIlzLZdoaO/+QHszyEpNIjs9mZz0FDJTE0lPSSQtOZHpKT7SUxJJTw4sS0/xkZLoIyUxgSRfAsmJgZ8kn5GSmECyz/fe86TEBHxmJJiRkAAJZoHnCpuINZ4jgFVAlXOuOvhCTwGbgMMhbTYBT7hAmuw0sywzKyDw7X64bTcB64Lb/wjYAYwYAFfqgV9X8uzeszgAB47AH/HAb3A4/P5A28HL3XvtAxv3+R3tl/t+6zL/1CQf1y7I5i9vXsyGpbNJTtT1dTLxpiX5WFKQyZKCzA+s6+7r53zrZc62XKKpo4cLHd1c6OrlQmc3Fzp7aO7o4VzLZTp7+ujs7qOzu59LkzBxXYIFpsawYCgEHgeWJQyERkhODHxHMmzQ84H1vx0q760PfY1RtrX3G4742pHuH25fzqqS7Al9zXACYC5QE/K8lsC3/NHazB1l21nOuToA51ydmQ15b0Qz2wJsAZg3b14Y5X7QrMwUFs/OBAv845tZ8HfIcwv8h/T+spDnBgQf+8yYkZpEVloS+ZnTWJCbzlWzMvRHXzyVkuhjfk4683PSw96m3+/o7Omjq7ufju4+uvv66enz09vv6Onz09PfT0+fo6ffH1z+/m+/c/T7we8cfr+j3zn8LtAV6nfB5/7Asv7gsoFtAl+pBr5UhfwevJyh1xPSafF+GzfMNkOvD32NaJGeMvHnfMIJgKFicvDuG65NONuOyDn3KPAoBLqAxrLtgDuumccd11xZeIjEKl+CkTktiUzdljRuhfO1tRYInby+EDgXZpuRtq0PdhMR/N0QftkiIjJe4QTALqDMzErMLBm4E9g2qM024C4LWA20Brt3Rtp2G3B38PHdwLPj/CwiIjIGo3YBOef6zOw+4HkCQzkfd84dMrOtwfWPANsJjACqIjAM9AsjbRt86W8Bz5jZl4AzwO9M6CcTEZER6UIwEZEYN9wwUA1dERGJUwoAEZE4pQAQEYlTCgARkTgVVSeBzawROH2Fm+cCTRNYzmSLpnqjqVZQvZMpmmqF6Kp3PLXOd87lDV4YVQEwHma2e6iz4JEqmuqNplpB9U6maKoVoqveyahVXUAiInFKASAiEqfiKQAe9bqAMYqmeqOpVlC9kymaaoXoqnfCa42bcwAiIvLb4ukIQEREQigARETiVFwEgJltMLNjZlYVvP9wRDGzU2Z2wMz2mtnu4LJsM3vRzCqDv2d6WN/jZtZgZgdDlg1bn5l9Lbivj5nZzRFS7zfN7GxwH+8N3sfa83rNrMjMXjKzI2Z2yMz+OLg84vbvCLVG6r6dZmZvm9m+YL1/G1wecft2lHonb/8652L6h8A01CeABUAysA8o97quQTWeAnIHLfs2cH/w8f3AP3lY38eBlcDB0eoDyoP7OAUoCe57XwTU+03gL4Zo62m9QAGwMvg4AzgerCni9u8ItUbqvjVgevBxEvAWsDoS9+0o9U7a/o2HI4D3bmrvnOsBBm5MH+k2AT8KPv4R8GmvCnHOvQJcGLR4uPo2AU8557qdcycJ3CNi1VTUOWCYeofjab3OuTrn3J7g43bgCIF7aUfc/h2h1uF4vW+dc64j+DQp+OOIwH07Sr3DGXe98RAAw92wPpI44AUze8fMtgSXzXKBu6oR/J3vWXVDG66+SN7f95nZ/mAX0cBhf8TUa2bFwIcJfPOL6P07qFaI0H1rZj4z20vglrMvOuciet8OUy9M0v6NhwAY943pp8B1zrmVwEbgXjP7uNcFjUOk7u+HgVJgBVAH/GtweUTUa2bTgf8G/sQ51zZS0yGWTWm9Q9QasfvWOdfvnFtB4H7kq8xs2QjNI7XeSdu/8RAA4dzU3lPOuXPB3w3AzwkcxtWbWQFA8HeDdxUOabj6InJ/O+fqg/9z+YHHeP9Q2fN6zSyJwB/UnzjnfhZcHJH7d6haI3nfDnDOtQA7gA1E6L4NFVrvZO7feAiAcG5q7xkzSzezjIHHwE3AQQI13h1sdjfwrDcVDmu4+rYBd5pZipmVAGXA2x7U91sG/ocPup3APgaP6zUzA34AHHHO/VvIqojbv8PVGsH7Ns/MsoKPU4FPAkeJwH07Ur2Tun+n6gy3lz8Eblh/nMBZ8q97Xc+g2hYQOJO/Dzg0UB+QA/waqAz+zvawxicJHHr2EvjW8aWR6gO+HtzXx4CNEVLvj4EDwP7g/zgFkVAv8FECh+37gb3Bn1sicf+OUGuk7turgXeDdR0EvhFcHnH7dpR6J23/aioIEZE4FQ9dQCIiMgQFgIhInFIAiIjEKQWAiEicUgCIiMQpBYCISJxSAIiIxCkFgMg4mFlxcH78x4JzuL8QvIpTJOIpAETGrwx40Dm3FGgBPuttOSLhUQCIjN9J59ze4ON3gGLvShEJnwJAZPy6Qx73A4leFSIyFgoAEZE4pQAQEYlTmg1URCRO6QhARCROKQBEROKUAkBEJE4pAERE4pQCQEQkTikARETilAJARCRO/X/aplZbIwTZNAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"posterior_n = marginal(posterior_jeff, 1)\n",
"posterior_n.plot()\n",
"posterior_n.mean()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"91"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_n.idxmax()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The posterior mean is close to 100 and the MAP is 91; both are consistent with the results in Equation 6.122."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Robot A\n",
"\n",
"Now we get to what I think is the most interesting part of this example, which is to take into account a second observation under two models of the scenario:\n",
"\n",
"> Two robots, [A and B], have different prior information about the source of the particles.\n",
"The source is hidden in another room which A and B are not allowed to enter.\n",
"A has no knowledge at all about the source of particles; for all [it] knows, ... the other room might be full of little [people] who run back and forth, holding first one radioactive source, then another, up to the exit window. \n",
">\n",
"> B has one additional qualitative fact: [it] knows that the source is a radioactive sample of long lifetime, in a fixed position.\n",
"\n",
"In other words, B has reason to believe that the source strength $s$ is constant from one interval to the next, while A admits the possibility that $s$ is different for each interval.\n",
"\n",
"The following figure, from Jaynes, represents these models graphically (Jaynes calls them \"logical situations\" because he seems to be allergic to the word \"model\").\n",
"\n",
"\n",
"\n",
"For A, the \"different intervals are logically independent\", so the update with `c2 = 16` starts with the same prior."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"c2 = 16\n",
"posterior2 = update(joint, phi, c2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the posterior marginal distribution of `n2`."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"168.947980523708"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAp5ElEQVR4nO3deXxV9Z3/8dcnN3tCEiAJCSGQAAFkUYSAKFqV1grYilo71S7a1bFqZ9pOO7XT+XXamV87Haedzji1+tNOFzu2ajelFbWudUUIyh6WGBICSUhCNiBk//7+yI3GmOWGLOcu7+fjkUfuck7u+x5I3vd877nfY845REQk8kR5HUBERLyhAhARiVAqABGRCKUCEBGJUCoAEZEIFe11gJFIT093eXl5XscQEQkp27Ztq3POZfS/PaQKIC8vj6KiIq9jiIiEFDMrH+h2DQGJiEQoFYCISIRSAYiIRCgVgIhIhFIBiIhEKBWAiEiEUgGIiESokPocgMh46ezqZsuhenYebaKlvYuctHhWz01nxuREr6OJjBsVgES8Tbuq+M5jxRxtPP2u+96/cBrfunIR09MSPEgmMr5UABKxWju6+KdH9/BQUQWLpqfwjSvO4sKCdJJjoymtO8XGHZX85MVSLv/hC3znmiVcec50ryOLjCkVgESkjq5ubvvV6zxdXMMtl8zhS5fNI8b39lticzOT+fJl8/jQshz+7uEd/M2v36CxpZ0bzs/zLrTIGNObwBJxnHP83cM7eLq4hn+5ajF/v3bBO/749zVrahIPfO483nfWNL756B5+uXnAKVVEQpIKQCLOz14uY+OOSr56+Xw+sWrWsMvHRfu4++PLeO+CTL61cQ+vvFk3ASlFxp8KQCLK7qNNfO/xfbzvrExuuWROwOvF+KL4z+uWkp+exK0PvD7gG8YioSagAjCztWa238xKzOz2Ae43M7vTf/9OM1vW576fmlmNme3ut84UM3vKzA76v08e/dMRGVxnVzdf+c0O0hJjuOPaczCzEa0/KT6G+24opL2zm7//7Q6cc+OUVGRiDFsAZuYD7gLWAQuB681sYb/F1gEF/q+bgLv73PdzYO0AP/p24BnnXAHwjP+6yLh54LXD7Ks+wT9vWMSUpNgz+hn56Un8wxVn8XLJcR547fAYJxSZWIHsAawESpxzpc65duBBYEO/ZTYA97sem4E0M8sGcM69ANQP8HM3AL/wX/4FcNUZ5BcJSP2pdn7w5/1cODedyxdljepnfXTlTC6cm853NxVT3dQ6RglFJl4gBZADVPS5fsR/20iX6W+ac64KwP89c6CFzOwmMysys6La2toA4oq8253PHORUexf/9MGFIx766c/M+O7VS+jsdvzbE/vGKKHIxAukAAb6bek/+BnIMmfEOXevc67QOVeYkfGuU1qKDKumuZVfbTnMtctmUDBt0pj8zJlTE/ncRfn84Y2jbCtvGJOfKTLRAimAI0Bun+szgMozWKa/Y73DRP7vNQFkERmxe/5SSle349ZL547pz73lkrlkTorjX/60V28IS0gKpAC2AgVmlm9mscB1wMZ+y2wEbvAfDbQKaOod3hnCRuBG/+UbgUdHkFskILUn2njgtXKuPjeHmVPHdmK3pLhovnzZPLZXNPLsPr1+kdAzbAE45zqB24AngWLgYefcHjO72cxu9i+2CSgFSoD7gFt61zezXwOvAvPN7IiZfcZ/1/eAy8zsIHCZ/7rImHrgtXLau7pHdMz/SHxo+QxmTknkB38+QHe39gIktAQ0F5BzbhM9f+T73nZPn8sOuHWQda8f5PbjwHsDTioyQu2d3fzv5sNcMi+D2RnJ4/IYMb4ovvi+Ar788A6e2FPN+iXZ4/I4IuNBnwSWsLVpVxV1J9v45Or8cX2cDUtzmJ2exI+fL9F7ARJSVAAStn72ShmzM5K4aG76uD6OL8r47EWz2X20mVdLj4/rY4mMJRWAhKW9lc3sqGjkE6tmERU1uuP+A3HNshzSk2O594XScX8skbGiApCw9JttFcT6orhq6XCfRxwb8TE+bjw/j+f317K/+sSEPKbIaKkAJOy0d3bzyBtHuWzhNCaf4Zw/Z+Ljq2aREOPTXoCEDBWAhJ1nio/R0NLBtYUzJvRxJyfF8pEVuWzccZRjzZojSIKfCkDCzu9eP8K0lDjeUzDxU4d8anUeHV2OB7dUDL+wiMdUABJWmk538JcDtXzw7On4JuDN3/5mTU3iooJ0Htx6mM6u7gl/fJGRUAFIWHlq7zE6uhwfOGe6Zxk+vmoWVU2tmh5Cgp4KQMLKYzsryUlL4JwZqZ5leO+CTLJS4vlfnTBGgpwKQMJGU0sHLx6s4wNnZ496zv/RiPZFcd3KXF44UEv58VOe5RAZjgpAwsaTe6vp7HZccbb38/Fct2ImvijjV9oLkCCmApCw8djOKnKnJLAkx7vhn15ZqfG8d0Emv3v9CB16M1iClApAwkLDqXZeLqnjiiXTPR3+6euvCnOpO9nO8/t1KlMJTioACQtPFR/rGf4JoumYL56fQXpyLL/dps8ESHBSAUhYeLa4huzUeBbnpHgd5S0x/rmInimu4fjJNq/jiLyLCkBCXltnFy8erGXNgsygGf7pdW3hDDq7HRt3DHeKbJGJpwKQkPdaaT2n2rt471mZXkd5lwVZKSzJSeU3RUe8jiLyLioACXnP7qshPiaKC+aM74lfztS1y2ewt6qZPZVNXkcReQcVgIQ05xzP7DvG6jnpxMf4vI4zoCvPmU6Mz3h0u4aBJLioACSkHaw5SUX9adYE4fBPr8lJsVw8L4ON2yvp7tY5gyV4qAAkpD1T3DPh2poFwVsAAFcuzaG6uZUtZfVeRxF5iwpAQtqz+46xaHoK2akJXkcZ0vvOyiQx1qdhIAkqKgAJWU2nO9hW3hD0r/4BEmOjuWzhNB7fXUV7p6aGkOCgApCQ9eqbdXQ7uMiDM3+diQ1Lp9PY0sGLBzU1hAQHFYCErBcO1pEU6+PcmWleRwnIRQUZTE6M0TCQBA0VgISslw7Wcf6cqcT4QuO/cYwvivVLsnlq7zFOtXV6HUdEBSChqfz4KQ7Xt4TM8E+vDUtzON3RxdPFx7yOIqICkND04sE6AC4sCM5P/w6mcNZkpqfGaxhIgkJABWBma81sv5mVmNntA9xvZnan//6dZrZsuHXNbKmZbTaz7WZWZGYrx+YpSSR46WAd01PjmZ2e5HWUEYmKMj5wznRePFhL0+kOr+NIhBu2AMzMB9wFrAMWAteb2cJ+i60DCvxfNwF3B7DuHcC3nXNLgW/6r4sMq6vb8cqbdVxUkBF0s38GYv2SbDq6HE/t1TCQeCuQPYCVQIlzrtQ51w48CGzot8wG4H7XYzOQZmbZw6zrgN7J21MB7RNLQHYeaaS5tTPkhn96nTMjlZy0BDbtqvI6ikS4QAogB+h7SqMj/tsCWWaodb8I/LuZVQDfB74+0IOb2U3+IaKi2lodPy094/9msHpuaBaAmbF+SZaGgcRzgRTAQPvY/We0GmyZodb9PPAl51wu8CXgfwZ6cOfcvc65QudcYUZGaB3xIePj1TePszA7hSlJsV5HOWO9w0BPaxhIPBRIARwBcvtcn8G7h2sGW2aodW8Efu+//Bt6hotEhtTa0cXrhxs4f/ZUr6OMytLcNA0DiecCKYCtQIGZ5ZtZLHAdsLHfMhuBG/xHA60CmpxzVcOsWwlc7L+8Bjg4yuciEWBHRSNtnd2sCvECMDPWLc7ixYN1NLdqGEi8MWwBOOc6gduAJ4Fi4GHn3B4zu9nMbvYvtgkoBUqA+4BbhlrXv87ngB+Y2Q7gu/QcPSQypM2l9ZjBivwpXkcZtfVnZ9Pe1a1hIPFMdCALOec20fNHvu9t9/S57IBbA13Xf/tLwPKRhBXZXHqcRdNTSE2I8TrKqJ2bm8b01Hg27arimmUzvI4jEUifBJaQ0drRxbbDDazKD+3hn15mxrol2bxwQMNA4g0VgISM7RWNtIfB+H9f65f0DAM9o7mBxAMqAAkZm0uPh834f69zc9PITo3nsZ3VXkeRCKQCkJARTuP/vaKijHWLs3nhYC0nNAwkE0wFICGh5/j/xrAZ/+/rirOzaO/s5tl9NV5HkQijApCQEI7j/73OzZ3MtJQ4Ht+lYSCZWCoACQmbS48TFWbj/716h4Ge21+jM4XJhFIBSEjoGf9PDavx/77WLs6irbOb5/drwkOZOCoACXpvjf/PDr9X/71W5E0hPTmWTbs1N5BMHBWABL1dR5to7+xmRV74FoAvyrh8URbP7auhtaPL6zgSIVQAEvS2ltUDUBjGBQCwbnE2Le1d/OWAhoFkYqgAJOhtK2tgTkZSSM//H4jzZk9hcmIMj2uKaJkgKgAJat3djqLyBgpnhferf4AYXxTvX5jFM8U1tHVqGEjGnwpAgtqbtSdpOt1BYd5kr6NMiHVLsjjR1snLJXVeR5EIoAKQoLa1rAEI//H/XhfMSWdSfDSb9KEwmQAqAAlqReX1pCfHkjc10esoEyI2OorLFk7jz3uqae/s9jqOhDkVgAS1orIGls+ajJl5HWXCrF+cTXNrJ6+WHvc6ioQ5FYAErZrmVg7Xt4T18f8DubAgneS4aB0NJONOBSBBq6i8Z/x/+azIeAO4V3yMjzULMvnz3mN0dmkYSMaPCkCCVlFZA/ExUSyanup1lAm3fkkW9afa2XKo3usoEsZUABK0isrrOWdGGrHRkfff9OJ5mSTE+DQ3kIyryPvNkpDQ0t7JnsrmiBv/75UQ2zMM9MTuY3R1O6/jSJhSAUhQ2n64ka5ux/II+QDYQNYuzqLuZBvb/O+FiIw1FYAEpaLyBsxg2czILYBLF2QSFx3FJh0NJONEBSBBaWtZPfOnTQrbE8AEIjkumovnZfDE7mq6NQwk40AFIEGnq9vxxuHGiJn/Zyjrl2RT3dzKGxWNXkeRMKQCkKCzr7qZk22dETED6HDWnJVJjM94QkcDyThQAUjQ6X3TU3sAkBIfw0UFGWzaVY1zGgaSsRVQAZjZWjPbb2YlZnb7APebmd3pv3+nmS0LZF0z+4L/vj1mdsfon46Eg61lDWSlxJOTluB1lKCwbnEWRxtPs+tok9dRJMwMWwBm5gPuAtYBC4HrzWxhv8XWAQX+r5uAu4db18wuBTYAZzvnFgHfH4snJKFvW1k9hXmRNQHcUC5bOI3oKNMU0TLmAtkDWAmUOOdKnXPtwIP0/OHuawNwv+uxGUgzs+xh1v088D3nXBuAc65mDJ6PhLijjaepbGqlMMLm/xlKWmIs58+ZyuO7qzQMJGMqkALIASr6XD/ivy2QZYZadx5wkZm9ZmZ/MbMVIwku4akoQk4AP1Lrl2RTfryF4qoTXkeRMBJIAQy0H97/Zchgywy1bjQwGVgFfBV42AbY5zezm8ysyMyKamtrA4groayorIGkWB8LsiZ5HSWovH/hNKIMHtfRQDKGAimAI0Bun+szgMoAlxlq3SPA7/3DRluAbiC9/4M75+51zhU65wozMjICiCuhbGtZPctmTSbapwPU+pqaHMd5+VN5bJeGgWTsBPJbthUoMLN8M4sFrgM29ltmI3CD/2igVUCTc65qmHUfAdYAmNk8IBbQmbAjWHNrB/uPnYi4+f8DtX5JFqW1pzhYc9LrKBImhi0A51wncBvwJFAMPOyc22NmN5vZzf7FNgGlQAlwH3DLUOv61/kpMNvMdtPz5vCNTi9tItrr5Q04R8TOADqcyxdlYQaP62ggGSPRgSzknNtEzx/5vrfd0+eyA24NdF3/7e3Ax0cSVsLbtvIGfFHG0tw0r6MEpcyUeFbMmsLju6v42/cVeB1HwoAGWiVobC2rZ2F2CklxAb0uiUhrF2exr/oEb9ZqGEhGTwUgQaGjq5vtFY0a/x/G2sVZADyxW8NAMnoqAAkKeyqbae3o1vj/MKanJXDuzDQdDipjQgUgQeHtD4BpD2A46xdns/toM4ePt3gdRUKcCkCCQlFZA7lTEpiWEu91lKDXOwykvQAZLRWAeM45R1F5g+b/D1DulESW5KSySe8DyCipAMRz5cdbqDvZpuGfEVi3JIsdFY0cbTztdRQJYSoA8VxR7wlgtAcQsHWLswEdDSSjowIQzxWV1ZMSH01BZrLXUUJGfnoSC7Im8fguvQ8gZ04FIJ4rKm9g+azJREXpBDAjsX5JNkXlDVQ3tXodRUKUCkA8VX+qnZKak6zI1/DPSF1xds8w0J929p+cVyQwKgDxVO/x//oA2MjNyUhmSU4qj25XAciZUQGIp4rKG4j1RbEkJ9XrKCFpw9Lp7DrapLmB5IyoAMRTW8vqOXtGKvExPq+jhKQPnjMdM7QXIGdEBSCeOd3exe6jTRr/H4VpKfGcP3sqj24/qjOFyYipAMQz2ysa6ehyrNAHwEblqqU5lB9vYceRJq+jSIhRAYhnet8AXj5TewCjcfniLGJ9UTzyxlGvo0iIUQGIZ7aWNzB/2iRSE2O8jhLSUhNiWLMgkz/trKKzq9vrOBJCVADiia5ux+vlDazI1/DPWLjq3OnUnWzjlTePex1FQogKQDyxr7qZk22dOv5/jFwyP5NJ8dE8sl3DQBI4FYB4Yuuh3hPAqADGQnyMj3WLs3hydzWtHV1ex5EQoQIQT2wtb2B6ajw5aQleRwkbG5bmcKq9i6eLj3kdRUKECkAmnHOOorJ6vfofY6tmT2VaShx/eF3DQBIYFYBMuCMNpznW3KYPgI0xX5Rx9bkzeP5ALTUnNEOoDE8FIBNu61sTwOkIoLF27fIZdHU7fSZAAqICkAm3tayeSfHRzMuc5HWUsDM3M5lzZ6bx221HNDWEDEsFIBNua1kDhToBzLj58PJcDhw7yU5NDSHDUAHIhNIJYMbfB87JJi46it9uO+J1FAlyKgCZUFsO6QQw4y0lPoa1i7N4dPtRfSZAhhRQAZjZWjPbb2YlZnb7APebmd3pv3+nmS0bwbpfMTNnZumjeyoSCjaXHic+JopzZqR5HSWsXbt8Bs2tnfpMgAxp2AIwMx9wF7AOWAhcb2YL+y22Dijwf90E3B3IumaWC1wGHB71M5GQ8NqhegpnTSE2Wjuf4+mCOelMT43XMJAMKZDfwpVAiXOu1DnXDjwIbOi3zAbgftdjM5BmZtkBrPtD4O8BHa4QARpb2tlX3cyq2Rr+GW++KOOaZTN44UAtlY2nvY4jQSqQAsgBKvpcP+K/LZBlBl3XzK4Ejjrndgz14GZ2k5kVmVlRbW1tAHElWL12qB7nej6xKuPvIytyccBDWyuGXVYiUyAFMNCxev1fsQ+2zIC3m1ki8A3gm8M9uHPuXudcoXOuMCMjY9iwErx6x//P1vj/hMidksh7CjJ4aGuFzhMgAwqkAI4AuX2uzwD6n4F6sGUGu30OkA/sMLMy/+2vm1nWSMJLaNlcqvH/ifbR82ZS3dzKs/tqvI4iQSiQ38StQIGZ5ZtZLHAdsLHfMhuBG/xHA60CmpxzVYOt65zb5ZzLdM7lOefy6CmKZc656rF6YhJcesf/z9Px/xPqvQsymZYSx6+26DgLebdhC8A51wncBjwJFAMPO+f2mNnNZnazf7FNQClQAtwH3DLUumP+LCTovTX+P0fj/xMp2hfFRwpz+cuBWirqW7yOI0EmOpCFnHOb6Pkj3/e2e/pcdsCtga47wDJ5geSQ0PVaab1//D/V6ygR5yMrZ/Kj50p4aGsFX7l8vtdxJIhoMFYmxObS4yyfNZm4aJ/XUSJOTloCl87P5KGiCjr0ZrD0oQKQcdfY0k5xdTOr8jX845WPnjeT2hNtPL1XnwyWt6kAZNxt0fi/5y6Zn0lOWgK/eLXM6ygSRFQAMu42a/zfc74o48YLZrG5tJ49lZomWnqoAGTcvfJmHYWzpmj832MfKZxJQoyPn79c5nUUCRIqABlXNSda2Vd9ggsLNNmr11ITY7h2+Qwe3V5J3ck2r+NIEFAByLh6uaQOgAvnqgCCwSdX59He1c2vXtMHw0QFIOPspYPHmZwYw8LsFK+jCDAnI5lL5mfwy83ltHXqZDGRTgUg48Y5x0sltVwwN13n/w0in16dT+2JNh7bWeV1FPGYCkDGzZu1JznW3MZFGv4JKhcVpDM3M5mfvnyIng/xS6RSAci4efFgz/j/ahVAUDEzPr06n91Hm3n1zeNexxEPqQBk3Lx0sI68qYnkTkn0Oor0c82yHDInxXHX8yVeRxEPqQBkXHR0dbO59Lhe/Qep+Bgfn7toNi+XHGd7RaPXccQjKgAZF28cbuRUe5cO/wxi1583k9SEGH78nPYCIpUKQMbF8/tr8EUZq/UBsKCVHBfNJy/I4897j3Hg2Amv44gHVAAyLp7bX0vhrMmkxMd4HUWG8MkL8kiM9XHP8296HUU8oAKQMVfVdJriqmYuXZDpdRQZxuSkWD66ciaP7qjUGcMikApAxtzz+2sBWKMCCAmfvWg2PjN+rCOCIo4KQMbcc/tqyElLoCAz2esoEoCs1HiuX5nLb4qOUH78lNdxZAKpAGRMtXV28VJJHZfMz8BM0z+EilsvnUu0z/ivZw56HUUmkApAxtTWQw20tHdp+CfEZKbEc+P5eTzyxlFKanREUKRQAciYenZfDbHRUZyv0z+GnL++eA4JMT5++JT2AiKFCkDGjHOO5/bXcF7+FBJjo72OIyM0JSmWz1yYz2O7qth1RKeNjAQqABkzB2tOcqjuFO9flOV1FDlDn3vPbKYmxfJ/H9urmUIjgApAxsyTu6sxg8sXTvM6ipyhSfExfPGyebx2qJ6ni2u8jiPjTAUgY+aJPdUsmzmZzJR4r6PIKFy/Ipc5GUn866ZiOrq6vY4j40gFIGOior6FPZXNrNXwT8iL9kXxD+vPorTulM4dHOZUADImntxTDcDlKoCwsGZBJhfMmcp/PHWA+lPtXseRcaICkDHxxO5qzspOYeZUnfwlHJgZ37pyEafaOrnjiX1ex5FxElABmNlaM9tvZiVmdvsA95uZ3em/f6eZLRtuXTP7dzPb51/+D2aWNibPSCZcTXMr2w43aPgnzMybNolPX5jPg1sreP1wg9dxZBwMWwBm5gPuAtYBC4HrzWxhv8XWAQX+r5uAuwNY9ylgsXPubOAA8PVRPxvxxJN7j+EcrF2sAgg3f/PeArJS4vnmo7vp6tZhoeEmkD2AlUCJc67UOdcOPAhs6LfMBuB+12MzkGZm2UOt65z7s3Ou07/+ZmDGGDwf8cAft1cyJyOJedM0+Vu4SY6L5h8/cBa7jzbzi1fKvI4jYyyQAsgBKvpcP+K/LZBlAlkX4NPA4wM9uJndZGZFZlZUW1sbQFyZSBX1LWwpq+fqc3M0+VuYumJJNpfOz+Dfn9zP4eM6Z0A4CaQABvqt7r8vONgyw65rZt8AOoEHBnpw59y9zrlC51xhRkZGAHFlIm3cUQnAhqUD9bqEAzPju9csITrK+NrvduoTwmEkkAI4AuT2uT4DqAxwmSHXNbMbgQ8AH3P6XxVynHP84Y2jrMibTO4UHf0TzrJTE/iHK87i1dLj/HpLxfArSEgIpAC2AgVmlm9mscB1wMZ+y2wEbvAfDbQKaHLOVQ21rpmtBb4GXOmc035lCNpT2UxJzUmuOlev/iPBdStyWT13Kt95bK+GgsLEsAXgf6P2NuBJoBh42Dm3x8xuNrOb/YttAkqBEuA+4Jah1vWv8yNgEvCUmW03s3vG7mnJRHjkjaPE+IwrlmR7HUUmgJlxx7XnEBVl/O1Db2iaiDBgoTTyUlhY6IqKiryOIUBXt+P8f32Gc3LTuO+GQq/jyAT6445KvvDrN/ibNXP58vvnex1HAmBm25xz7/pF1SeB5Yz85UANNSfauEbDPxHng+dM50PLZvCj50rYcqje6zgyCioAOSO/eq2C9OQ43qepnyPStzcsYuaURL7w69epOdHqdRw5QyoAGbGqptM8u+8YHy6cQYxP/4UiUXJcNHd/fDlNpzu47Vd6PyBU6bdXRuzhrUfodnD9ipleRxEPnZWdwveuOZsth+r5t8c1YVwo0olbZUQ6urp5cOthLipI18yfwlXn5rC9opGfvHSIBdkpXLtcM7qEEu0ByIg8vruaqqZWPnlBntdRJEh844qzWD13Kl///U5eebPO6zgyAioACZhzjv956RCz05O4dH6m13EkSMT4ovjxx5aTNzWJm3+5jZKaE15HkgCpACRgrx9uYEdFI59anUdUlCZ+k7elJsTws0+tIDbax40/3Upl42mvI0kAVAASsHv+UkpqQgwf0jivDGDG5ER+/qkVNJ/u4GM/eU2Hh4YAFYAEpLiqmaf2HuNTq/NIjNWxAzKwxTmp/PzTK6huauUTP9lCg84nHNRUABKQHz1bQnJcNJ+6IN/rKBLkls+awk9uLOTQ8VNcd+9m7QkEMRWADOvAsRNs2l3FJy/IIzUxxus4EgJWz03nZ59cweH6Fj7y/zZzVO8JBCUVgAzrjif2kRwbzWcu1Kt/Cdzquen872dXUneyjQ/f/Qr7q3V0ULBRAciQNpce5+niGj5/6RwmJ8V6HUdCzPJZU3jwplV0djuuvfsVXjig07oGExWADKq72/HdTcVkp8bz6dV69S9nZtH0VB65dTU5kxP41M+3cv+rZTqtZJBQAcigHiqqYOeRJr56+XziY3xex5EQNj0tgd9+/gIumZfBNx/dw5ce2s6ptk6vY0U8FYAMqPZEG/+6qZjz8qdwteb8lzGQHBfNfTcU8pX3z2Pjjkquuutl9lY2ex0roqkAZED//Ke9tHZ0852rl2CmT/3K2IiKMm5bU8AvP3MeDS0dbLjrJf77mYOaTtojKgB5l0e3H+WPOyq5bc1c5mYmex1HwtDquek89aX3sHZxNj946gDX/FhHCXlBBSDvUFHfwj/+YTfLZ03mlkvmeB1HwtjkpFj++/pz+fHHlnG08TTr73yRb23cQ1NLh9fRIoYKQN5yur2LWx54HYD//MhSonW2L5kA65dk8/SXL+a6Fbnc/2oZl3z/OX65uVzDQhNAv+EC9Ez1/NXf7mB3ZRM//MhScqfoZC8ycaYkxfKdq5fwxy9cyLxpk/g/j+xmzQ+e58Eth2nvVBGMFxWA4FzP8f5/2lnF19Yu0InexTOLpqfy4E2r+MkNhUxOjOX23+/i0u8/z/2vlnFSh42OOQulD2QUFha6oqIir2OEFecc//HUAf772RJuPH8W37pykY76kaDgnOP5A7Xc+cxB3jjcSHJcNB9alsMnzs/TwQkjZGbbnHOF77pdBRC5Oru6+fYf9/LLzeVctyKX7169RCd6kaDjnGN7RSO/fLWcP+2sor2rm6W5aVy1dDofOGc66clxXkcMeioAeYf6U+18+eHtPL+/lr++eDZfu3yB/vhL0Ks72cbvth3hke2VFFc144syVs2ewpoF01izIJP89CSvIwYlFYC85fn9NXztdztpONXBP125kI+dN8vrSCIjduDYCR554yhP7T3GwZqTAOSnJ3HJ/AzOy5/C8llTyJikvQNQAQhQWnuSO57YzxN7qpmTkcSd15/LoumpXscSGbWK+hae21/DM8U1bC49Tpv/yKH89CQKZ03m7BmpnJWdwoLsFJLjIu+MdiqACNXd7dhcepxfb63gsZ2VxEZH8YU1BXz2onziojXBm4Sf9s5udh1toqisnq1lDWwrr6ehz4fLZk1N5KysFOZNS2bm1CTypiYya2oS6cmxYXsAxKgKwMzWAv8F+ICfOOe+1+9+89+/HmgBPumce32odc1sCvAQkAeUAX/lnGsYKocKIDCn2jrZcqiel0vqeHx3NUcbTzMpPprrV87kcxfN1m6xRBTnHFVNreytbKa4qpni6maKq05QfvwU3X3+/CXF+sidkkhmSjxZKXFMS4l/6ytzUhyTE2NJTYxhUlx0yL1fdsYFYGY+4ABwGXAE2Apc75zb22eZ9cAX6CmA84D/cs6dN9S6ZnYHUO+c+56Z3Q5Mds59bagskV4AzjlaO7o51d5JS1sXza0dHGtu5VhzG9XNrbxZc5L9x05wqO4UXd2O2Ogozp89lQ8tn8H7F07TlM4ifbR3dnO08TRlx09RXneK8voWKupPU3OileqmVupOtr2jIHpFGaQlxpKWEENqYgxpCTGkJMSQGBtNYqyPpFgfCbHRJMX5SIjxkRQXTUKsj7joKOKio4jxvf319nUj1n851hc15gUzWAEEMhi2EihxzpX6f9CDwAZgb59lNgD3u5422WxmaWaWTc+r+8HW3QBc4l//F8DzwJAFcKbufOYgG3dUvnUSirf+Td07vr3rfvfW/e6d1/v9pwh4vX7390kywDLv/JmdXY5T7Z3veuxeZpA7OZH5WZNYtziL8/KnUpg3WX/0RQYRGx1FfnpSz5FD8999f2dXN3Un2znW3ErNiTaaTnfQ2NJOY0sHjad7vjed7qDuZDuldadoae+ipa2Tlo6uQX9PAxVlEB0VRVQU+MzwRRk//thyLixIH90P7ieQAsgBKvpcP0LPq/zhlskZZt1pzrkqAOdclZllDvTgZnYTcBPAzJkzA4j7bpmT4pg/bZL/B77j21tjfm9fH/r+t9e3QZYf5P5+P6D/ekOuC/iiokiO85EYF/3WK4zkuGimpcSRlRpPRnKc5u4RGUPRviiyUuPJSo0f0Xq9e+ot7Z20tHdxqr2TU21dtHd209HV89Xe2U17VzcdXe6t6x1dPbe1d3bT2eXoco6u7re/slLHfug2kAIYaF+kf78Ntkwg6w7JOXcvcC/0DAGNZN1e162cyXUrz6w8RERGwsxIiPWREOtjqtdhhhHIS8YjQG6f6zOAygCXGWrdY/5hIvzfawKPLSIioxVIAWwFCsws38xigeuAjf2W2QjcYD1WAU3+4Z2h1t0I3Oi/fCPw6Cifi4iIjMCwQ0DOuU4zuw14kp5DOX/qnNtjZjf7778H2ETPEUAl9BwG+qmh1vX/6O8BD5vZZ4DDwIfH9JmJiMiQ9EEwEZEwN9hhoDpsREQkQqkAREQilApARCRCqQBERCJUSL0JbGa1QPkZrp4O1I1hnPEWSnlDKSso73gKpawQWnlHk3WWcy6j/40hVQCjYWZFA70LHqxCKW8oZQXlHU+hlBVCK+94ZNUQkIhIhFIBiIhEqEgqgHu9DjBCoZQ3lLKC8o6nUMoKoZV3zLNGzHsAIiLyTpG0ByAiIn2oAEREIlREFICZrTWz/WZW4j//cFAxszIz22Vm282syH/bFDN7yswO+r9P9jDfT82sxsx297lt0Hxm9nX/tt5vZpcHSd5vmdlR/zbe7j+Pted5zSzXzJ4zs2Iz22Nmf+u/Pei27xBZg3XbxpvZFjPb4c/7bf/tQbdth8k7ftvXORfWX/RMQ/0mMBuIBXYAC73O1S9jGZDe77Y7gNv9l28H/s3DfO8BlgG7h8sHLPRv4zgg37/tfUGQ91vAVwZY1tO8QDawzH95EnDAnynotu8QWYN12xqQ7L8cA7wGrArGbTtM3nHbvpGwB/DWSe2dc+1A74npg90G4Bf+y78ArvIqiHPuBaC+382D5dsAPOica3POHaLnHBErJyJnr0HyDsbTvM65Kufc6/7LJ4Bies6lHXTbd4isg/F62zrn3En/1Rj/lyMIt+0weQcz6ryRUACDnbA+mDjgz2a2zcxu8t82zfWcVQ3/90zP0g1ssHzBvL1vM7Od/iGi3t3+oMlrZnnAufS88gvq7dsvKwTptjUzn5ltp+eUs08554J62w6SF8Zp+0ZCAYz6xPQTYLVzbhmwDrjVzN7jdaBRCNbtfTcwB1gKVAE/8N8eFHnNLBn4HfBF51zzUIsOcNuE5h0ga9BuW+dcl3NuKT3nI19pZouHWDxY847b9o2EAgjkpPaecs5V+r/XAH+gZzfumJllA/i/13iXcECD5QvK7e2cO+b/5eoG7uPtXWXP85pZDD1/UB9wzv3ef3NQbt+Bsgbztu3lnGsEngfWEqTbtq++ecdz+0ZCAQRyUnvPmFmSmU3qvQy8H9hNT8Yb/YvdCDzqTcJBDZZvI3CdmcWZWT5QAGzxIN879P7C+11NzzYGj/OamQH/AxQ75/6jz11Bt30HyxrE2zbDzNL8lxOA9wH7CMJtO1Tecd2+E/UOt5df9Jyw/gA975J/w+s8/bLNpued/B3Ant58wFTgGeCg//sUDzP+mp5dzw56XnV8Zqh8wDf823o/sC5I8v4S2AXs9P/iZAdDXuBCenbbdwLb/V/rg3H7DpE1WLft2cAb/ly7gW/6bw+6bTtM3nHbvpoKQkQkQkXCEJCIiAxABSAiEqFUACIiEUoFICISoVQAIiIRSgUgIhKhVAAiIhFKBSAyCmaW558f/z7/HO5/9n+KUyToqQBERq8AuMs5twhoBD7kbRyRwKgAREbvkHNuu//yNiDPuygigVMBiIxeW5/LXUC0V0FERkIFICISoVQAIiIRSrOBiohEKO0BiIhEKBWAiEiEUgGIiEQoFYCISIRSAYiIRCgVgIhIhFIBiIhEqP8PCB5ZjWUGb/UAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"posterior_n2 = marginal(posterior2, 1)\n",
"posterior_n2.plot()\n",
"posterior_n2.mean()"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"160"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_n2.idxmax()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The posterior mean is close to 169, which is consistent with the result in Equation 6.124.\n",
"The MAP is 160, which is consistent with 6.123."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Robot B\n",
"\n",
"For B, the \"logical situation\" is different. If we consider $s$ to be constant, we can -- and should! -- take the information from the first update into account when we perform the second update.\n",
"We can do that by using the posterior distribution of $s$ from the first update to form the joint prior for the second update, like this:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"137.499999315101"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArG0lEQVR4nO3de3SV9Z3v8fc3OzcSApuQADEEghBAVAwYkWpbb2MLtlM6vep0RnuZMrTaOZ2z2hm7eqZn5pyz5vQ6F9e4tNo61d6s0+mFntqx1latWpQgBIkQCOEWAuQCuUKu+3f+2E90G3eSnWQnz758Xmvtlb2f5/c8+7sf2Pnkuf1+5pxDRETST4bfBYiIiD8UACIiaUoBICKSphQAIiJpSgEgIpKmMv0uYCKKiopceXm532WIiCSVXbt2tTrnikdOT6oAKC8vp7q62u8yRESSipkdizZdh4BERNJUTAFgZpvMrM7M6s3s7ijzzczu8ebvNbP1EfMeMrNmM9sXZbnPeOutNbOvTu2jiIjIRIwbAGYWAO4FNgNrgNvMbM2IZpuBCu+xFbgvYt53gE1R1nsDsAVY65y7FPj6JOoXEZFJimUPYANQ75xrcM71A48S/sUdaQvwiAvbAQTNrATAOfcscDbKej8FfNk51+e1a57shxARkYmLJQBKgRMRrxu9aRNtM9JK4G1m9qKZPWNmV0VrZGZbzazazKpbWlpiKFdERGIRSwBYlGkje5CLpc1ImcA8YCPweeAxM3vTepxzDzjnqpxzVcXFb7qKSUREJimWAGgEyiJeLwaaJtEm2np/4h02egkIAUUx1CMiInEQSwDsBCrMbJmZZQO3AttHtNkO3O5dDbQR6HDOnRpnvT8DbgQws5VANtA6keIlvYRCjmcPtvBPTx7kkT8cpbW7z++SRJLauDeCOecGzewu4AkgADzknKs1s23e/PuBx4FbgHrgPPCx4eXN7IfA9UCRmTUC/9M5923gIeAh7/LQfuAOp8EJZBSdvQN85ge7eebg6+eBvv5EHd/4UCU3r1noY2UiycuS6XduVVWV053A6ad/MMSff/tFdh07x/941yXcumEJR9t6+Jsf7+XVpk4evL2KG1Yv8LtMkYRlZrucc1Ujp+tOYEl49zx1iBePnOXrH7yCj167jNysAKsXzeH7f3E1qxYV8Nkf7eFk+wW/yxRJOgoASWj1zV3c98xhPnDlYt677o1XFhfkZnHvn65nYCjE32+v9alCkeSlAJCE9pX/qiMvK8AXNq+OOr+8KJ/P3FjBk6+e4feHdJ+IyEQoACRhHTjdyZOvnuEv3nYx82fnjNru428tpzQ4i3968iDJdE5LxG8KAElYDz57hLzsAHdcs3TMdjmZAe68YQW7j7ezoyFaryMiEo0CQBJSV+8Av3yliS2VpQTzssdt/771pQTzsvjujqPTX5xIilAASEL61Sun6R0I8YErF8fUPjcrwIeqynii9gynO3qnuTqR1KAAkIT045cbubgon/VLgjEv82dXLyXkHD946fj0FSaSQhQAknCOt53npSNnef+Vi4nSP+ColszP4/qVxfxo53GGQjoZLDIeBYAknO01JwH4k3Xj9Sj+Zu+/cjFnOvt46YhOBouMRwEgCec3+5u5oizIRcFZE172xtULyMsO8Iu943VGKyIKAEkoLV191DS2c9Mk+/bJy87kpksW8qtXTjEwFIpzdSKpRQEgCeXpumacC/8lP1l/vLaEc+cHeOFwWxwrE0k9CgBJKE/tb2bRnFwuvWjOpNdx3apiCnIz+X81OgwkMhYFgCSMvsEhfn+ohRsvWTChq39GyskMcMOqBfyurpmQrgYSGZUCQBLGrmPn6Okf4oZVU+/b/8bVC2jt7mfvyY44VCaSmmIKADPbZGZ1ZlZvZndHmW9mdo83f6+ZrY+Y95CZNXsjf0Vb9+fMzJmZxgNOczsazpJhcPXFhVNe13Uri8kw+O3+M3GoTCQ1jRsAZhYA7gU2A2uA28xszYhmm4EK77EVuC9i3neATaOsuwy4GdCtm8KOhjYuK53LnNysKa9rXn42Vy6dx1MHmuNQmUhqimUPYANQ75xrcM71A48CW0a02QI84sJ2AEEzKwFwzj0LjHZXzj8DfwPoQG2a6x0YYs/xdjZePD9u67xx9UJqmzrVN5DIKGIJgFLgRMTrRm/aRNu8gZm9BzjpnKsZp91WM6s2s+qWFg34kapePnaO/qEQG+Nw+GfYTZeEzyU8Xae9AJFoYgmAaJdjjPyLPZY2rzc2ywO+CHxpvDd3zj3gnKtyzlUVFxeP11yS1I6GNjIMriqPXwBULJjNwjk5PFffGrd1iqSSWAKgESiLeL0YGHmBdSxtIi0HlgE1ZnbUa/+ymS2KoR5JQTsaznJ56VwK4nD8f5iZce2KIl443KbLQUWiiCUAdgIVZrbMzLKBW4HtI9psB273rgbaCHQ4506NtkLn3CvOuQXOuXLnXDnhAFnvnDs9uY8hyaxvcIg9J9rZsCx+f/0Pe+uKIs729LP/dGfc1y2S7MYNAOfcIHAX8ASwH3jMOVdrZtvMbJvX7HGgAagHHgQ+Pby8mf0Q+AOwyswazewTcf4MkuRqmzrpHwpx5dJ5cV/3tSvCVxc/r8NAIm+SGUsj59zjhH/JR067P+K5A+4cZdnbYlh/eSx1SGrafbwdgHVL4h8AC+fkUrFgNs/Vt7H17cvjvn6RZKY7gcV3u4+f46K5uSyckzst6792RREvHWmjb3BoWtYvkqwUAOK73cfbp+Wv/2HXLJ9P70CIvY3qFkIkkgJAfNXc2cvJ9gusm8DYvxM1fGmpRgkTeSMFgPhq94l2gGkNgHn52axcOFsBIDKCAkB89fLxc2QFjEsvmjut73NVeSG7jp3TYPEiERQA4qu9Jzq4pGQOuVmBaX2fDcsK6e4bZP8p3Q8gMkwBIL5xzlHb1DHtf/0Dr91kpsNAIq9TAIhvGs9doLN3kMtKJz/8Y6xK5s6irHCWAkAkggJAfLPPG61rJvYAIHweYOfRs4TvWxQRBYD4prapk0CGsXpRwYy839XLCmnr6edwS8+MvJ9IolMAiG9qmzpYUTx72k8ADxu+H2DnUR0GEgEFgPhoX1Mnl87A8f9hy4ryKZqdo/MAIh4FgPiiubOXlq6+GTv+D+HxAaqWzmPXsXMz9p4iiUwBIL6obQpfj3/ZRTO3BwCwfmmQ42fP09rdN6PvK5KIFADii9qm8BVAa2Y4AIY7ndvjdUEtks4UAOKLfSc7WTo/L65DQMbisovmkplh7D6hw0AiMQWAmW0yszozqzezu6PMNzO7x5u/18zWR8x7yMyazWzfiGW+ZmYHvPY/NbPglD+NJI3aUx1cNoPH/4fNyg5wSckcXj7WPuPvLZJoxg0AMwsA9wKbgTXAbWa2ZkSzzUCF99gK3Bcx7zvApiirfhK4zDm3FjgIfGGixUty6jg/wImzF2b88M+w9UuC1DS2q2M4SXux7AFsAOqdcw3OuX7gUWDLiDZbgEdc2A4gaGYlAM65Z4E3XXfnnPu1N94wwA5g8WQ/hCSX2lPDdwD7EwDrlszjfP8QB890+fL+IokilgAoBU5EvG70pk20zVg+DvxqAu0lidWdDv/i9WsPYHjsgd06ESxpLpYAsCjTRu47x9Im+srNvggMAt8fZf5WM6s2s+qWlpZYVikJru50F/PysiienePL+y8pzKMwP5vdx3UiWNJbLAHQCJRFvF4MNE2izZuY2R3Au4GPuFF66HLOPeCcq3LOVRUXF8dQriS6ujNdrFpUgFm0vxumn5mxrizIywoASXOxBMBOoMLMlplZNnArsH1Em+3A7d7VQBuBDufcqbFWamabgL8F3uOcOz+J2iUJhUKOg6e7WLVwZjqAG836pfM43NJDx/kBX+sQ8dO4AeCdqL0LeALYDzzmnKs1s21mts1r9jjQANQDDwKfHl7ezH4I/AFYZWaNZvYJb9a/AQXAk2a2x8zuj9eHksR1sv0CPf1DrFrkz/H/YevKggDsaWz3tQ4RP2XG0sg59zjhX/KR0+6PeO6AO0dZ9rZRpq+IvUxJFcMngFfNUBfQo1lbFsQMdh8/x3UrdWhR0pPuBJYZVedderly4Wxf65idk8mqhQW8rCuBJI0pAGRG1Z3uojQ4a8a7gIhm3ZJ57Dl+jpBuCJM0pQCQGVV3umvGRgAbz7olQTp7B2lo7fa7FBFfKABkxvQPhjjc0s3KRAmA4RPBJzr8LUTEJwoAmTFHWnsYDLmE2QNYXjyb2TmZ7FHPoJKmFAAyY14/AZwYAZCRYaxdPJc9J9r9LkXEFwoAmTF1pzvJzDCWF/t7BVCkyrIgB0510Tsw5HcpIjNOASAzpu50FxcX55OdmTj/7SrLggyG3GsjlImkk8T5JkrKqzvTlTCHf4ZVeieC1TOopCMFgMyI8/2DNJ67kHABsGBOLhfNzaWmUXsAkn4UADIjGlp6cA4qFiTO8f9hlUuCuhJI0pICQGZEfXP4ZqsKn7uAiKayLMiJsxdo6+7zuxSRGaUAkBlxqLmLzAxj6fx8v0t5kysWBwF0OaikHQWAzIhDZ7opL8onK5B4/+UuXzyXQIZRowCQNJN430ZJSfUt3axIoOv/I+VlZ7JyYQG7FQCSZhQAMu36B0McazufkMf/h1WWBak50a6eQSWtKABk2h1t62Eo5FiRgFcADassm0tn7yBH23r8LkVkxsQUAGa2yczqzKzezO6OMt/M7B5v/l4zWx8x7yEzazazfSOWKTSzJ83skPdz3tQ/jiSiQ2fCVwAldgCE//vpRLCkk3EDwMwCwL3AZmANcJuZrRnRbDNQ4T22AvdFzPsOsCnKqu8GnnLOVQBPea8lBR1q7sKMhOoDaKQVC2aTnx1QAEhaiWUPYANQ75xrcM71A48CW0a02QI84sJ2AEEzKwFwzj0LnI2y3i3Aw97zh4H3TqJ+SQL1zd2UzcsjNyvgdymjCmQYl6tnUEkzsQRAKXAi4nWjN22ibUZa6Jw7BeD9XBCtkZltNbNqM6tuaWmJoVxJNPXN3Ql5B/BIlWXz2H+qUz2DStqIJQAsyrSRl0rE0mZSnHMPOOeqnHNVxcXF8VilzKDBoRANrT0Jffx/WGVZkIEhx6unOv0uRWRGxBIAjUBZxOvFQNMk2ox0ZvgwkfezOYZaJMmcOHeB/sFQUgTAuiVBAPaoZ1BJE7EEwE6gwsyWmVk2cCuwfUSb7cDt3tVAG4GO4cM7Y9gO3OE9vwP4+QTqliRxyBsFrCLBegGNZuGcXBbNyaWmsd3vUkRmxLgB4JwbBO4CngD2A48552rNbJuZbfOaPQ40APXAg8Cnh5c3sx8CfwBWmVmjmX3Cm/Vl4GYzOwTc7L2WFFPfEr4EdHlx4vUBFE1lWVAngiVtZMbSyDn3OOFf8pHT7o947oA7R1n2tlGmtwE3xVypJKX6M92UzM2lIDfL71JiUrkkyH/VnuZsTz+F+dl+lyMyrXQnsEyr+pbupDj+P2y4Z1AdBpJ0oACQaRMKOeqbkysA1i6eS4bpRLCkBwWATJumjguc7x9KqgDIzwn3DKrzAJIOFAAybV4bBWxB4l8BFKmyLEhNYzvhU1siqUsBINPm9QBInj0AgCvKgrSfH+BY23m/SxGZVgoAmTb1zd3Mz89mXpJdTVNZFgTUM6ikPgWATJtDSXYCeNjKhQXkqWdQSQMKAJkWzjkOnelKygAIZBiXlapnUEl9CgCZFi3dfXT2Dibd8f9h68qCvNrUSd+gegaV1KUAkGlR/9ooYMl1BdCwyrIg/UMh9p/q8rsUkWmjAJBpcWj4CqAEHgh+LJWv9Qx6zt9CRKaRAkCmRX1zNwW5mSwoyPG7lElZNCeXBQU51DR2+F2KyLRRAMi0ONTcRcWC2ZhFGyso8ZmZegaVlKcAkGlR35wco4CNpXJJkCOtPbSf7/e7FJFpoQCQuDvX009rd1/SdQExUuVrPYPqMJCkJgWAxN3wIDDJvgdw+eK5mHoGlRQWUwCY2SYzqzOzejO7O8p8M7N7vPl7zWz9eMuaWaWZ7TCzPWZWbWYb4vORxG/DfQAlewAU5GZRsWA2e07oSiBJTeMGgJkFgHuBzcAa4DYzWzOi2WagwntsBe6LYdmvAv/gnKsEvuS9lhRw6Ew3s7IClAZn+V3KlIV7Bu1Qz6CSkmLZA9gA1DvnGpxz/cCjwJYRbbYAj7iwHUDQzErGWdYBc7znc4GmKX4WSRCHmrtYviCfjIzkvAIo0hVlQc729HPi7AW/SxGJu1gCoBQ4EfG60ZsWS5uxlv0s8DUzOwF8HfhCtDc3s63eIaLqlpaWGMoVvx1u7k76E8DDhnsG3a3DQJKCYgmAaH/GjdwfHq3NWMt+Cvhr51wZ8NfAt6O9uXPuAedclXOuqri4OIZyxU9dvQM0dfQm/fH/YasWFpCblaH7ASQlxRIAjUBZxOvFvPlwzWhtxlr2DuAn3vP/IHy4SJLc4ZYeIPlPAA/LDGRweelcahQAkoJiCYCdQIWZLTOzbOBWYPuINtuB272rgTYCHc65U+Ms2wRc5z2/ETg0xc8iCSBZRwEbS2VZkH1NnfQPhvwuRSSuxg0A59wgcBfwBLAfeMw5V2tm28xsm9fscaABqAceBD491rLeMp8EvmFmNcA/Er56SJLcoeYusgMZLCnM87uUuKksm0f/YIgDpzv9LkUkrjJjaeSce5zwL/nIafdHPHfAnbEu601/DrhyIsVK4qs/082yonwyA6lzj+Fwz6A1J9pZ690dLJIKUudbKgmhvqWbFUnaBfRoLpqbS9HsHHbrPICkGAWAxE3vwBDHz55nRXFqBYB6BpVUpQCQuDnc0o1zyTsIzFjWLQnS0NJDx/kBv0sRiRsFgMTN61cApcZNYJGu8I797z3Z7msdIvGkAJC4qW/uJsOgvCh1rgAatrYs3DPobvUMKilEASBxc+hMN+Xz88nJDPhdStzN8XoGrT6mLiEkdSgAJG7qW7pT5g7gaDYsK+TlY+cYCqlnUEkNCgCJi/7BEEdbk38YyLFcVV5Id98g+0/phjBJDQoAiYtjbT0MhlxKXgE0bMOyQgBeOnLW50pE4kMBIHGRylcADSuZO4vF82YpACRlKAAkLg55AXBxcb7PlUyvDcsK2Xn0rEYIk5SgAJC4ONTcTWlwFnnZMXUvlbQ2lBfS1tNPQ2uP36WITJkCQOKi7nQnqxel7uGfYVd55wF26jCQpAAFgExZ/2CIhpYeVqVBAFxclE/R7GxeOqoAkOSnAJApa2jtZjDk0iIAzIyrygt1IlhSggJApqzudBdAWgQAhO8HaDx3gVMdF/wuRWRKYgoAM9tkZnVmVm9md0eZb2Z2jzd/r5mtj2VZM/uMN6/WzL469Y8jfqg73UVmhnFxUereAxBJ9wNIqhg3AMwsANwLbAbWALeZ2ZoRzTYDFd5jK3DfeMua2Q3AFmCtc+5S4Ovx+EAy8w6e6eLi4nyyM9Njh/KSkjnMzslkp84DSJKL5Ru7Aah3zjU45/qBRwn/4o60BXjEhe0AgmZWMs6ynwK+7JzrA3DONcfh84gPDpzuYuXC9Dj8AxDIMK5cOk97AJL0YgmAUuBExOtGb1osbcZadiXwNjN70cyeMbOror25mW01s2ozq25paYmhXJlJ3X2DNJ67kBaXgEbasKyQg2e6OdvT73cpIpMWSwBYlGkjb4Mcrc1Yy2YC84CNwOeBx8zsTe2dcw8456qcc1XFxcUxlCsz6eCZ8AngdNoDAHjL8vkA/OFwm8+ViExeLAHQCJRFvF4MNMXYZqxlG4GfeIeNXgJCQFHspUsiOOhdAbR60RyfK5lZa0vnUpCTyfOHW/0uRWTSYgmAnUCFmS0zs2zgVmD7iDbbgdu9q4E2Ah3OuVPjLPsz4EYAM1sJZAP6NiWZA6e7yMsOsHjeLL9LmVGZgQyuvriQF+r1X1aS17gdtzjnBs3sLuAJIAA85JyrNbNt3vz7gceBW4B64DzwsbGW9Vb9EPCQme0D+oE7nHrYSjoHz3RRsbCAjIxoR/tS2zXLi/jN/mZOtl+gNJheASipIaaeu5xzjxP+JR857f6I5w64M9Zlven9wJ9NpFhJPHWnu7jpkgV+l+GLa1eEj1g+X9/Kh6rKxmktknjS48JtmRat3X209fSzKs2O/w9buXA2RbNzdBhIkpYCQCbttS4g0uwKoGFmxjXL5/P84TaNDyBJSQEgk1bb1AHAmovScw8A4NoV82np6nttRDSRZKIAkEmrbeqkZG4uhfnZfpfim2uWv34eQCTZKABk0mqbOrk0jf/6BygrzGNJYR7PKQAkCSkAZFIu9A/R0NLNmovm+l2K765bWcwLh9voGxzyuxSRCVEAyKQcON1JyMGakvTeAwC4flUx5/uH2HnknN+liEyIAkAmpbapEyDtDwFBuF+g7MwMnq5Th7aSXBQAMimvnupk7qystOsCIpq87EyuXlbI7xQAkmQUADIptU2drCmZQ5QOXNPS9asWcLilhxNnz/tdikjMFAAyYYNDIQ6c6kzr6/9HumFVuKtyHQaSZKIAkAk70tpD32BIx/8jLCvKZ0lhHk/XadAiSR4KAJmw108A6xLQYWbG9auKef5wK70DuhxUkoMCQCastqmD7MwMlhfn+11KQrlh1QJ6B0LsaNAoYZIcFAAyYbVNnaxeVEBmQP99Ir1l+XzysgP8+tUzfpciEhN9g2VCnHNeFxA6/DNSblaAG1Yt4Ne1ZxgKqXdQSXwxBYCZbTKzOjOrN7O7o8w3M7vHm7/XzNZPYNnPmZkzM40HnASOtp2n48IAlWUKgGjecelCWrv72H1cdwVL4hs3AMwsANwLbAbWALeZ2ZoRzTYDFd5jK3BfLMuaWRlwM3B8yp9EZkTNiXYArigL+lpHorpx9QKyAxk8UXva71JExhXLHsAGoN451+AN4/gosGVEmy3AIy5sBxA0s5IYlv1n4G8A7S8niT0n2pmVFWBF8Wy/S0lIBblZXLNiPk/UntEgMZLwYgmAUuBExOtGb1osbUZd1szeA5x0ztWM9eZmttXMqs2suqVF11j7bW9jO5eXztUJ4DFsunQRx8+eZ/+pLr9LERlTLN/iaPf6j/zTZrQ2UaebWR7wReBL4725c+4B51yVc66quLh43GJl+gwMhdjX1MkVOv4/pj9as5AMQ4eBJOHFEgCNQFnE68VAU4xtRpu+HFgG1JjZUW/6y2a2aCLFy8yqO91F/2BIx//HUTQ7h6ryQgWAJLxYAmAnUGFmy8wsG7gV2D6izXbgdu9qoI1Ah3Pu1GjLOudecc4tcM6VO+fKCQfFeuecvjEJbM/wCeDFQV/rSAbvuryEA6e7qDutw0CSuMYNAOfcIHAX8ASwH3jMOVdrZtvMbJvX7HGgAagHHgQ+Pdaycf8UMiP2nGinMD9bXUDH4F1rSwhkGD/bc9LvUkRGlRlLI+fc44R/yUdOuz/iuQPujHXZKG3KY6lD/LXr2DnWL5mnLqBjUDQ7h7dXFPHz3Sf5/DtWkZGhbSaJR5dySExau/s40trDVeXz/C4labx3XSlNHb3sPHrW71JEolIASEyqj4bvbK1SAMTs5jULycsO6DCQJCwFgMRk17GzZGdmcFmpLgGNVV52Ju+8dBG/3HuKvkF1ES2JRwEgMak+do4rFs8lJzPgdylJZUvlRXT2DvK7A7qJURKPAkDG1TswxL6THVy5tNDvUpLOW1cUUVyQw39Unxi/scgMUwDIuGpOtDMw5KhaquP/E5UZyODDVWX8rq6Zk+0X/C5H5A0UADKuF4+cxUwngCfr1g1lOOBHL6nTW0ksCgAZ1/P1rVx60RyCedl+l5KUFs/L4/qVxTy68wQDQyG/yxF5jQJAxnShf4jdx9u5ZrnG65mKj1y9lOauPp7a3+x3KSKvUQDImHYdO0f/UIi3LJ/vdylJ7YbVC7hobi7ff/GY36WIvEYBIGN64XArmRnGhnJdATQVgQzjw1ct4feHWmlo6fa7HBFAASDjeP5wG5VlQfJzYuo2Ssbwp1cvITszgwd/3+B3KSKAAkDG0HFhgFca27lGh3/iorgghw9euZj/3HWS5s5ev8sRUQDI6J471ErIwdtWaiS2eNn69osZDIV46PmjfpciogCQ0f32QDPBvCzWaQSwuFk6P5/Nl5fw/R3H6Owd8LscSXMKAIkqFHI8XdfMdSuLNQB8nH3quuV09Q3ygxd1Y5j4S99siaqmsZ22nn5uXL3A71JSzmWlc3lbRREPPttAT9+g3+VIGospAMxsk5nVmVm9md0dZb6Z2T3e/L1mtn68Zc3sa2Z2wGv/UzMLxuUTSVz87kAzGQbX6fj/tPjrm1fS1tPPvz9/xO9SJI2NGwBmFgDuBTYDa4DbzGzNiGabgQrvsRW4L4ZlnwQuc86tBQ4CX5jyp5G4eepAM+uXzFP3D9Nk/ZJ5/NElC/nmsw20n+/3uxxJU7HsAWwA6p1zDc65fuBRYMuINluAR1zYDiBoZiVjLeuc+7U3aDzADmBxHD6PxMHR1h5qmzrZdNkiv0tJaZ9750q6+wa593f1fpciaSqWACgFIjszb/SmxdImlmUBPg78Ktqbm9lWM6s2s+qWFg2qMRN++copAG65vMTnSlLb6kVz+OCVi/nOC0d1d7D4IpYAsCjTXIxtxl3WzL4IDALfj/bmzrkHnHNVzrmq4mIdj54Jv6hp4sql87goOMvvUlLe5965ipzMAP/nl/v9LkXSUCwB0AiURbxeDDTF2GbMZc3sDuDdwEeccyNDRXxQ39zNgdNdvEt//c+IBQW5/NVNK/jtgWaefPWM3+VImoklAHYCFWa2zMyygVuB7SPabAdu964G2gh0OOdOjbWsmW0C/hZ4j3PufJw+j0zRL2qaMIN3rVUAzJSPXrOM1YsK+Luf7dPNYTKjxg0A70TtXcATwH7gMedcrZltM7NtXrPHgQagHngQ+PRYy3rL/BtQADxpZnvM7P74fSyZjKGQ48e7Grlm+XwWzsn1u5y0kZ2ZwZffv5bmrl6+8qsDfpcjaSSmLh6dc48T/iUfOe3+iOcOuDPWZb3pKyZUqUy73x9q4WT7Bb5wy2q/S0k7lWVBPn7tMr713BFuXrOQ61fpBjyZfroTWF7z6EsnmJ+fzTvW6PJPP3zunatYtbCAz/1HDS1dfX6XI2lAASAANHf18pv9Z/jAlYvJztR/Cz/kZgW457Z1dPUO8t8f28Ogxg+WaaZvugDw3T8cY8g5bt2wxO9S0tqqRQX8ry2X8vtDrXztiTq/y5EUp2GehK7eAR5+4SjvXLOIZUX5fpeT9j581RJeOdnBN59tYNWiAt63XjfJy/TQHoDwgxeP09k7yKdvWO53KeL50rsvZePFhfzNj/fyzEHdAS/TQwGQ5i70D/Gt547w1hVFrF0c9Lsc8WRnZvDA7VVULCzgU9/bxYsNbX6XJClIAZDmvvnsYVq6+virmyr8LkVGmJObxcMfu4qSubnc8e8vaU9A4k4BkMaa2i9w/zOHedflJWxYVuh3ORLFgjm5/Ogv38LFRbP55MPV/Ne+U36XJClEAZDG/u+vDhBycPdm3fiVyIpm5/DDT27k0tI5bPvey/zLbw4SCqnrLJk6BUCa+kVNE7+oaeLO61dQVpjndzkyjrl5Wfzwkxt53/pS/uU3h/jkI9UaSEamTAGQhk62X+CLP32FdUuC3Kkrf5JGblaAb3zwCv73lkt55mALN//zszokJFOiAEgz3X2D/OV3qxkKOf7lw5VkBvRfIJmYGX/+lnJ+due1FM/OYdv3Xmbbd3dxrK3H79IkCenbn0b6B0N86nu72H+qi3/7yHqWztdNX8nqstK5/Pyua/n8O1fx9MFmbvrGM/zdz/ZxprPX79IkiVgyjcNSVVXlqqur/S4jKXVcGOBT39vFC4fb+OoH1vKhqrLxF5Kk0NzZyz2/PcSjL4VHX33X2hI+ek05lWVBzKINyifpxsx2Oeeq3jRdAZD6aps6+OyjezjS2sNX3r+W91+prgVS0fG283znhaM8Vn2C7r5Blhfns6WylFsuL2F5cb7CII0pANLQuZ5+vvlsA9/6fQPBvGzuubWSa1YU+V2WTLPuvkF+vuckP9/TxEtHzgJQGpzFtSvm89aKYtYvCVIanKVASCNTCgBv+MZ/BQLAt5xzXx4x37z5twDngY86514ea1kzKwR+BJQDR4EPOefOjVWHAmB8QyFHTWM7P9t9kv/c1UhP/xDvW1/Kl969hmBett/lyQw71XGBp/Y389yhVl443Epn7yAA8/KyuKx0LpeUzKF8fj7l8/NYWpRPyZxcMjIUDKlm0gFgZgHgIHAz4UHedwK3OedejWhzC/AZwgFwNfCvzrmrx1rWzL4KnHXOfdnM7gbmOef+dqxaFABhfYNDdPUO0tU7yNmePo61nedY23lqGtvZdfQcXX2D5GRmsOmyRXz6+hWsWlTgd8mSAIZCjtqmDmoaO9jX2MErJzuob+6mP2LcgayAUTQ7x3tkh38W5DAnN4vZOQHyczKZ7T3yczLJzQqQnZlBVsDIDmR4z4cfpr2MBDFaAMTSHfQGoN451+Ct6FFgC/BqRJstwCPe0JA7zCxoZiWE/7ofbdktwPXe8g8DTxMeJD7u7nnqENtrmhgOu9ciz73hB5Fh+Pq04dfuja9H5ObIdY9c7o3TRq7DjTL/zeu8MDBE/+CbBwoxgxXFs/njyovYUF7IjZcsYE5u1pvaSfoKZBhrFwff0OnfUMhxurOXY609HG07z/Gz52np6qO1u4+W7j72n+qiraePgaHJHSrOChgZNvyADDPMICMjPM0IX9o6PC/DvNcZYITbjmasaBkreMaMpJl+vwn4x/ddzlXl8e2yJZYAKAVORLxuJPxX/nhtSsdZdqFz7hSAc+6UmUUdBNXMtgJbAZYsmdxgJQsKcli10Psr2N7w47V/uNdfR7z3KG1eX4e9YZmR6xg5/43PJ7bs8PzcrAAFuZnMmZVFQW4mwbxslhTmURqcRW5WIPoGEBlFIMMoDc6iNDiLa0YZpds5R+9AiO6+Qbr7Bunxfnb3DtI3GGJgKET/YIj+ofDz1187BoZChJzDOQiFHCGH9zriufceoVD4dch5r8c4OjFWHI11UGPs5Sb3fmPNdGMvOSGzpuH7HUsARAuwkZ9qtDaxLDsm59wDwAMQPgQ0kWWH3bphiUa6EpkkM2NWdoBZ2QGKC3L8LkfiKJYbwRqByIvGFwNNMbYZa9kz3mEivJ/NsZctIiJTFUsA7AQqzGyZmWUDtwLbR7TZDtxuYRuBDu/wzljLbgfu8J7fAfx8ip9FREQmYNxDQM65QTO7C3iC8KWcDznnas1smzf/fuBxwlcA1RO+DPRjYy3rrfrLwGNm9gngOPDBuH4yEREZk24EExFJcaNdBqrO4ERE0pQCQEQkTSkARETSlAJARCRNJdVJYDNrAY5NcvEioDWO5Uy3ZKo3mWoF1TudkqlWSK56p1LrUudc8ciJSRUAU2Fm1dHOgieqZKo3mWoF1TudkqlWSK56p6NWHQISEUlTCgARkTSVTgHwgN8FTFAy1ZtMtYLqnU7JVCskV71xrzVtzgGIiMgbpdMegIiIRFAAiIikqbQIADPbZGZ1ZlbvjT+cUMzsqJm9YmZ7zKzam1ZoZk+a2SHv5zwf63vIzJrNbF/EtFHrM7MveNu6zszemSD1/r2ZnfS28R5vHGvf6zWzMjP7nZntN7NaM/tv3vSE275j1Jqo2zbXzF4ysxqv3n/wpifcth2n3unbvs4bni1VH4S7oT4MXAxkAzXAGr/rGlHjUaBoxLSvAnd7z+8GvuJjfW8H1gP7xqsPWONt4xxgmbftAwlQ798Dn4vS1td6gRJgvfe8ADjo1ZRw23eMWhN12xow23ueBbwIbEzEbTtOvdO2fdNhD+C1Qe2dc/3A8MD0iW4L8LD3/GHgvX4V4px7Fjg7YvJo9W0BHnXO9TnnjhAeI2LDTNQ5bJR6R+Nrvc65U865l73nXcB+wmNpJ9z2HaPW0fi9bZ1zrtt7meU9HAm4bcepdzRTrjcdAmC0AesTiQN+bWa7zGyrN22hC4+qhvdzgW/VRTdafYm8ve8ys73eIaLh3f6EqdfMyoF1hP/yS+jtO6JWSNBta2YBM9tDeMjZJ51zCb1tR6kXpmn7pkMATHlg+hlwrXNuPbAZuNPM3u53QVOQqNv7PmA5UAmcAr7hTU+Ies1sNvCfwGedc51jNY0ybUbrjVJrwm5b59yQc66S8HjkG8zssjGaJ2q907Z90yEAYhnU3lfOuSbvZzPwU8K7cWfMrATA+9nsX4VRjVZfQm5v59wZ78sVAh7k9V1l3+s1syzCv1C/75z7iTc5IbdvtFoTedsOc861A08Dm0jQbRspst7p3L7pEACxDGrvGzPLN7OC4efAO4B9hGu8w2t2B/Bzfyoc1Wj1bQduNbMcM1sGVAAv+VDfGwx/4T1/Qngbg8/1mpkB3wb2O+f+KWJWwm3f0WpN4G1bbGZB7/ks4I+AAyTgth2r3mndvjN1htvPB+EB6w8SPkv+Rb/rGVHbxYTP5NcAtcP1AfOBp4BD3s9CH2v8IeFdzwHCf3V8Yqz6gC9627oO2Jwg9X4XeAXY631xShKhXuCthHfb9wJ7vMctibh9x6g1UbftWmC3V9c+4Eve9ITbtuPUO23bV11BiIikqXQ4BCQiIlEoAERE0pQCQEQkTSkARETSlAJARCRNKQBERNKUAkBEJE0pAESmwMzKvf7xH/T6cP+1dxenSMJTAIhMXQVwr3PuUqAdeL+/5YjERgEgMnVHnHN7vOe7gHL/ShGJnQJAZOr6Ip4PAZl+FSIyEQoAEZE0pQAQEUlT6g1URCRNaQ9ARCRNKQBERNKUAkBEJE0pAERE0pQCQEQkTSkARETSlAJARCRN/X8kh6AKjI0Z+wAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"joint = make_joint(posterior_s, ns)\n",
"posterior = update(joint, phi, c2)\n",
"posterior_n = marginal(posterior, 1)\n",
"posterior_n.plot()\n",
"posterior_n.mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The posterior mean of $n$ is close to 137.5, which is consistent with Equation 6.134.\n",
"The MAP is 132, which is one less than the analytic result, 133.\n",
"But again, there are two values with the same probability except for floating-point errors."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"132"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_n.idxmax()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.914335439641036e-16"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_n[132] - posterior_n[133]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Under B's model, the data from the first interval updates our belief about $s$, which influences what we believe about `n2`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Going the other way\n",
"\n",
"That might not seem surprising, but there is an additional point Jaynes makes with this example, which is that it also works the other way around: Having seen `c2`, we have more information about $s$, which means we can -- and should! -- go back and reconsider what we concluded about `n1`.\n",
"\n",
"We can do that by imagining we did the experiments in the opposite order, so\n",
"\n",
"1. We'll start again with a joint prior based on a uniform distribution for $s$,\n",
"\n",
"2. Update it based on `c2`,\n",
"\n",
"3. Use the posterior distribution of $s$ to form a new joint prior,\n",
"\n",
"4. Update it based on `c1`, and\n",
"\n",
"5. Extract the marginal posterior for `n1`."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"169.94393251129674"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"joint = make_joint(prior_s, ns)\n",
"posterior = update(joint, phi, c2)\n",
"posterior_s = marginal(posterior, 0)\n",
"posterior_s.mean()"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"131.49999935386944"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"joint = make_joint(posterior_s, ns)\n",
"posterior = update(joint, phi, c1)\n",
"posterior_n2 = marginal(posterior, 1)\n",
"posterior_n2.mean()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The posterior mean is close to 131.5, which is consistent with Equation 6.133.\n",
"And the MAP is 126, which is one less than the result in Equation 6.132, again due to floating-point error."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"126"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_n2.idxmax()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.8163916471489756e-16"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"posterior_n2[126] - posterior_n2[127]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's what the new distribution of `n1` looks like compared to the original, which was based on `c1` only."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"joint = make_joint(prior_s, ns)\n",
"posterior = update(joint, phi, c1)\n",
"posterior_n = marginal(posterior, 1)"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABKjElEQVR4nO3deXxU9b34/9d7sgIhYUnCkpV9XwUMigqugAtaa4utotXWWu36a29rv97ba1t76623661Xq63FFaRakSoWV4qgIGHft0BIIJANCFv29++Pc4JDyDJJZiYzyfv5eCSznHM+n8+ZMzPvOZ/P53w+oqoYY4wxocbT3gUwxhhjGmIByhhjTEiyAGWMMSYkWYAyxhgTkixAGWOMCUkWoIwxxoQkC1AdhIh8WUTeCWD680XkUff+ZSKyy49pvy0id7n37xaRlX5MO6CvS0uJSBcR+YeInBCRv7V3eVpCRA6IyNUBTF9FZLB7/ykR+Q+vZd8QkaMickpEeovIpSKyx318c6DKFIpEZLmIfNXHdc+9puHIAlQruB/Us+6H46iI/FVE4tqQ3iMi8mJbyqSqL6nqtW1JowV5faSqw5pbz9f9UtVZqvpcW8slIpnuBzLSK+2gvS4++jzQB+itqre1dGMR+aaIZItIhYjM93vpQoSq3q+qPwcQkSjgN8C1qhqnqiXAz4A/uo8XB7Ns3j/WTGBZgGq9G1U1DpgITAb+vb0K4v2F3IptRUTa5X3Qnnm3owxgt6pWt3L7w8CjwLP+K1JgteX96eoDxALbvJ7LqPc4mOUxwaKq9tfCP+AAcLXX48eBN937N+F8cI4Dy4ERXuv9CDgEnAR2AVcBM4FKoAo4BWxy100A/gIUuNs8CkS4y+4GVgG/BUrdZXcDK73yugRYC5xwby/xWrYc+IWbxllgcAP7OAFY75b1FWAh8Ki7bDqQ38r9uiBv97mv1tu3/3XLvhO4qonX/hHgRff+QUDd/E4BU1vxuvzczf8k8A6Q6C6LBV4EStxjuxbo08j7Y4Sb1nGc98JN7vM/rfea3NvAto8Ai4Dn3TJsAyY1sN6jwPxm3qfnXhv3cab7+kQ2t7/u8juBXHefH/Z+7XF+3D4E7HOXLwJ61cvnXveYrGikfP+G8/4+DNzjbjPYXTbf3cehwGmv4/qBm2ctzvvnFBBDyz8vMcD/uOU7CjwFdPF+fwPfBwrdNL/iLrvPPX6Vbt7/aGTfFHgA2OO+tj8HBgGfAGXu6xXttf7XgL1u+ZYA/b2WXYPzOTgB/BH4F+7nxV1+D7ADOAYsAzLqlaPuNZ0NbHfLcwj4QXt/lzb7XdveBQjHv3of1DScL5Gfe32YrgGigB+6b7poYBiQV/fGcz/Eg9z7j+D1ReI+txj4E9ANSAY+Bb7uLrsbqAa+BUQCXfD6IgZ6uW/WO93lt7uPe7vLl7sfzFHu8qh6eUfjfDF9z92Pz7sfygsCVCv264K8uTBAVXvl/UX3g9mr/mtfPw/qfQF7pdeS12Wfexy7uI8fc5d9HfgH0BWIAC4C4ht4b0ThHPP/576OV+J8IQxr7DWpt/0jQDnOl0kE8EtgdQPr+StANba/I3G+gC/H+TL/jXtc6t733wVWA6nu8j8BC+rl8zzO+7dLA2WbiRMYRrvrvEwDAaqJ41r/fbCYln1efocTCHoB3d1j+0uv93c1TjVilHsszgA965etidde3fTjcd7rFcD7wECcYLoduMtd90qgGKc2Jgbnx9kKd1kiTkD7vFuW77llq/u83Izzfhvh7tu/Ax/XK0fda1oAXObe7wlMbO/v0ub+Olv1ij8tFpHjwEqcXzT/hfNl+paqvquqVTi/0Lrg/GqvwXnzjRSRKFU9oKr7GkpYRPoAs4DvquppVS3E+fU312u1w6r6v6parapn6yVxPbBHVV9wly/A+QV2o9c681V1m7u8qt72WTgfht+papWqvopzxtAQn/fLx7zB+dVal/crOGdl1zeTpi98eV3+qqq73dd0ETDefb4K6I3zYa9R1XWqWtZAHllAHM4XfaWqfgC8iRMMfbVSVZeqag3wAjCuJTvZQo3t7+dxagVWqGoF8B84Zy11vg48rKr57vJHgM/Xqz57xH3/1n9/AnzBzXurqp52t2+Vln5ecH4AfA34nqqWqupJnM+v9/pVwM/c9+BSnGDdbLtrPf+tqmWqug3YCryjqjmqegJ4G6eWAuDLwLOqut59LX8MTBWRTNyzHlV91f2s/A444pXH13EC6w533/4LGC8iGQ2UpwrncxqvqsdUdX0L9yfoLEC13s2q2kNVM1T1AfdD2B/nzAMAVa3FObtIUdW9OL86HwEKRWShiPRvJO0MnABRICLH3UD4J5xfhnXymijbeeVw5QIpLdj+kLo/tby2v0AL98uXvGkk7+bS9IUvr4v3h/8MTrABJ1AsAxaKyGER+ZXbeN9QHnnusW8sj+bUL0NsANtNGtvf/ngdJzeIlHitmwG87vX+3IHzY6WP1zrNvce8lzf4/vJRSz8vSThnwuu81v+n+3ydEj2/ndD7tfHVUa/7Zxt47P1ae39vnMJ5rVO48DhovX3JAH7vtR+lgNDw++1WnICXKyL/EpGpLdyfoLMA5V+Hcd4wgNMJAKcK8BCAqr6sqtPcdRT4b3dVrZdOHk6VQKIbBHuoaryqjvJap/42jZbDlV5XDh+2LwBS3PJ7b9+gFuyXL3nTSN6H3funcb5c6vRtQbq+vC4Ncn9J/1RVR+KcEd8AzGskj7R6nT98yiMAmnqtmlOA894FQES64pxB1skDZnm9P3uoaqyqtuQ9lub1uNH3lw9a+nkpxgkQo7zWT1Cn05MvmnuftVT9741uOK/1IS48DsL5r1seTlWm93HooqofX1Bo1bWqOgcncC/GOWMOaRag/GsRcL2IXOX+uv4+zgfnYxEZJiJXikgMThXDWZxfnOD8ssqs+1JT1QKcButfi0i8iHhEZJCIXOFjOZYCQ0XkSyISKSJfxGlTeNPH7T/Bqef+trv954ApDa3Ykv1qgWQ37ygRuQ2nfn2pu2wjMNddNgmnKqpOEU411MBG0m316yIiM0RkjIhE4LQJVPHZfnpbgxMYfuiWcTpOFeLC5vLwhVvuWJz2qQgRaersaiNwuYiki0gCTtWRr14FbhCRaSISjdMe430cnwJ+UVeVJCJJIjKnBekvAu4WkZFu8PvPFmx7npZ+Xtyz22eA34pIslv+FBG5zscsj9L4e6w1Xga+IiLj3c/RfwFrVPUA8BYwSkQ+5x7nb3P+D42ngB+LyCgAEUlwPzPnEZFoca4JTHCrCsto+P0bUixA+ZGq7gLuwGnkLMb5YrpRVStx2mkec58/gvMl/P/cTesu2CwRkbp64Xk4jezbcRryXwX6+ViOEpxf+N/HqSr4IXCDqhb7uH0l8DmcxuVjOG1rf29k9Zbuly/WAEPcNH8BfN7dJ3DaQga55fopzoe7rtxn3PVXuVUeWfX2qy2vS1+cY1CGU531L5xefedxX7ubcNpEioH/A+ap6k6f9rx5/47zI+AhnPfaWRq5xEFV38XpgbkZWIfvP1Bw200exHl9C3Be73yvVX6P0wngHRE5idNh4uIWpP82TnvKBziN/B/4um0jWvp5+ZGb72oRKQPew/c2pr/gtOUcF5HFrS6xS1Xfx3lfv4bzWg/CbQ9z35u34XzGSnA+F6u8tn0dp8ZiobsfW3Heew25Ezjgrnc/zvsnpMn5Vf3GGGNMaLAzKGOMMSHJApQxxpiQZAHKGGNMSLIAZYwxJiSF1aCJiYmJmpmZ2d7FMMYY40fr1q0rVtWk+s+HVYDKzMwkOzu7vYthjDHGj0SkwZFErIrPGGNMSPIpQInITBHZJSJ7ReShBpaLiPzBXb5ZRCZ6LXtWRApFZGsD233LTXebiPyqbbtijDGmI2k2QLlDuzyBc3XySOB2ERlZb7VZOFc4D8GZL+VJr2XzcYbWr5/uDGAOMNYdM+t/WlF+Y4wxHZQvbVBTgL2qmgMgIgtxAst2r3XmAM+7I+2uFpEeItJPVQtUdYU4w8bX9w2cKQkqANwh8o0xJqxUVVWRn59PeXl5excl5MXGxpKamkpUVEMTAVzIlwCVwvnDu+dz4ZhbDa2TgjOuVGOGApeJyC9wBhn9gapeMOeQiNyHc1ZGenpbBjw2xhj/y8/Pp3v37mRmZnL+IPzGm6pSUlJCfn4+AwYM8GkbX9qgGnrF6w/g58s69UXizOqYhTP18yJp4Oiq6tOqOklVJyUlXdAL0Rhj2lV5eTm9e/e24NQMEaF3794tOtP0JUDlc/78I6l8NjdPS9ZpKN2/q+NTnGkSEn0ojzHGhBQLTr5p6evkS4BaCwwRkQHuvDBzcYbZ97YEmOf25ssCTrhztDRlMXClW+ihOEPl+zQdhDFNqq6EfR/Ahhfh4GqoDflpb4wxDWg2QLnTHn8TZ7rrHcAiVd0mIveLyP3uakuBHJz5VZ4BHqjbXkQW4EyAN0xE8kXkXnfRs8BAt/v5QuAutbk/TFsdWAl/nAQv3AJvPAjPXgd/uhwKNrd3yYwJGBHhzjvvPPe4urqapKQkbrjhhhalk5mZSXFx0+cJTa2zYcMGRIRly5a1KN/G+DSShKou5bMZTeuee8rrvuJMbtbQtrc38nwlYTBhlgkju96GV+6EnhnwxReh7xjI/Rje/xn8dRZ8aRFkXtrepTTG77p168bWrVs5e/YsXbp04d133yUlJSXo5ViwYAHTpk1jwYIFXHedrxMUN85GkjAdQ+EOePUeJyh99X0YcSP0zITxX4KvfQjxKbDwdijZ194lNSYgZs2axVtvvQU4geL22z87NygtLeXmm29m7NixZGVlsXmzU6NQUlLCtddey4QJE/j617+OdyXWiy++yJQpUxg/fjxf//rXqalpuqpcVXn11VeZP38+77zzjl+63YfVWHzGNKimGl77GkTHwdyXoUuP85fH94MvL4Knp8Nr98K970GEvfWN//30H9vYfrjMr2mO7B/Pf944qtn15s6dy89+9jNuuOEGNm/ezD333MNHH30EwH/+538yYcIEFi9ezAcffMC8efPYuHEjP/3pT5k2bRo/+clPeOutt3j66acB2LFjB6+88gqrVq0iKiqKBx54gJdeeol58+Y1mv+qVasYMGAAgwYNYvr06SxdupTPfe5zbdp3O4My4W/NU3B0C9zwGycYNaRnJtzwWzi8AVY/EdTiGRMMY8eO5cCBAyxYsIDZs2eft2zlypXn2qiuvPJKSkpKOHHiBCtWrOCOO5yWluuvv56ePXsC8P7777Nu3TomT57M+PHjef/998nJyWky/wULFjB37lzACZYLFixo8z7Zz0gT3s4ehxW/gsFXw/BmGoRH3QKbXoF/PQ7jvgRxdl2d8S9fznQC6aabbuIHP/gBy5cvp6Sk5NzzDfU/q+vy3VDXb1Xlrrvu4pe//KVP+dbU1PDaa6+xZMkSfvGLX5y7KPfkyZN07969lXtjZ1Am3K1+EspPwJX/Ab5cY3Htz6H6LHxkQz+ajueee+7hJz/5CWPGjDnv+csvv5yXXnoJgOXLl5OYmEh8fPx5z7/99tscO3YMgKuuuopXX32VwkJnBLrS0lJycxucEQOA9957j3HjxpGXl8eBAwfIzc3l1ltvZfHixW3aHwtQJnxVnoFPn4Zhs6H/eN+2SRwCY74A65+HM6UBLZ4xwZaamsp3vvOdC55/5JFHyM7OZuzYsTz00EM899xzgNM2tWLFCiZOnMg777xzbji5kSNH8uijj3LttdcyduxYrrnmGgoKGr+0dcGCBdxyyy3nPXfrrbfy8ssvt2l/JJwuPZo0aZLahIXmnOy/wpvfhbuXtqz7eOEO+L8smP5jmH7B7DHGtMiOHTsYMWJEexcjbDT0eonIOlWdVH9dO4My4WvdfOgzBjIuadl2ySNg6ExY8yfnLMwYE5IsQJnwVLgDCjbChC/71vZU36XfhbOlsPElf5fMGOMnFqBMeNq0ECQCRn++ddtnTIW+Y2HDC/4tlzHGbyxAmfBTWwObF8GQa9rWVXzCHVCwCY5s8V/ZjDF+YwHKhJ/9K+DkYRg3t23pjLkNIqJhg1XzGROKLECZ8LPlVYhJgKGz2pZO115OF/XNr0B1hX/KZozxGwtQJrzU1sDufzrVe1GxbU9vwp1OZ4nd/2x7Wsa0k1CYbiMzM5MxY8Ywfvx4xowZwxtvvNGivBtiAcqEl0Pr4EwxDGvj2VOdQTOgWxJs/bt/0jOmHXhPtwG023QbH374IRs3buTVV1/l29/+dpvTswBlwsuupeCJdMbe8wdPhDM1x5537JooE9bae7oNb2VlZecGnm0LnwaLFZGZwO+BCODPqvpYveXiLp8NnAHuVtX17rJngRuAQlUd3UDaPwAeB5JU1aZ8N03b9U9In3rhlBptMXIOZD8Le9+DkTf5L13T+bz9kP97hfYdA7Mea3a19p5uA2DGjBmoKjk5OSxatKjNu95sgBKRCOAJ4BogH1grIktUdbvXarOAIe7fxcCT7i3AfOCPwPMNpJ3mpnuw9btgOo3S/VC0Ayb+l3/TzZgGXXvD9jcsQJmw1dx0G6+99hpw4XQbf/+7U73d2HQbAGfPniU5ObnZMnz44YckJiayb98+rrrqKqZPn05cXFyr98mXM6gpwF5VzQEQkYXAHMA7QM0Bnnenfl8tIj1EpJ+qFqjqChHJbCTt3wI/BNremmY6vt3LnNuhM/2bbkSkM1XH1tegqtw/nS9M5+TDmU4gtdd0G/UNGjSIPn36sH37dqZMmdKqNMC3NqgUIM/rcb77XEvXOY+I3AQcUtVNPpTBGNj3AfQaBL0H+T/tkXOg8pSThzFhqr2m26ivsLCQ/fv3k5GR0ab98eUMqqGBzuqHY1/W+Wxlka7Aw8C1zWYuch9wH3BuKHjTCdVUQe4qGPuFwKQ/4HKIiYfdb8Pw2c2vb0wIamq6ja985SuMHTuWrl27njfdxu23387EiRO54oorGpxuo7a2lqioKJ544olmA86MGTOIiIigqqqKxx57jD59+rRpf5qdbkNEpgKPqOp17uMfA6jqL73W+ROwXFUXuI93AdNVtcB9nAm8WddJQkTGAO/jdKgASAUOA1NU9UhjZbHpNjqxvLXwl6vhtvnOzLiBsGgeHFwD39/ZugFoTadk0220jL+n21gLDBGRASISDcwFltRbZwkwTxxZwIm64NQQVd2iqsmqmqmqmThVghObCk6mk9u/3LnNvDxweQydCaeOOOPzGWPaXbMBSlWrgW8Cy4AdwCJV3SYi94vI/e5qS4EcYC/wDPBA3fYisgD4BBgmIvkicq+f98F0BvtXOHM/desduDwGXwOIc02UMabd+XQdlKouxQlC3s895XVfgQcb2fb2hp6vt06mL+UwnVRVuVP1Nvmrgc0nLglSJjq9Ba/4YWDzMh2KqjbYG86cr6UzuNtIEib05X8KNRUw8IrA5zV0pjOc0qmiwOdlOoTY2FhKSkpa/OXb2agqJSUlxMb6fhmHT2dQxrSr/R+BeJwRJAJtyLXw4S9g3/ttn87DdAqpqank5+dTVGQ/apoTGxtLamqqz+tbgDKh7+AnznAvsfGBz6vvWGdUiZzlFqCMT6KiohgwYEB7F6NDsio+E9pqqpwqt7Ss4OTn8cCAK5wAZVU2xrQrO4MyQbGv6BRvbDjE9oIyqmuVQUlxzBrdl4syejbduHx0K1SdgfSLG1/H3wZOh21/h6JdkDw8ePkaY85jAcoE1JnKan7+5g4Wrj2IR4QhyXFERgif7CvhLyv3c8mg3vzyc2PI6N2t4QQOrnFu04IcoMA5i7IAZUy7sQBlAqbwZDnz/vIpu4+e5J5LB3D/FYNI6h4DwOmKav6Wncdv3t3N9X9YyR+/NIHpwxoYLTlvNcSnQoLvDatt1jMDeg5wAlTW/c2ubowJDGuDMgFx/EwlX3pmDQdLzzD/K1P4jxtGngtOAN1iIrn70gG8/d3LSe/Vla8+l82HOwvPT0TVOYMKZvVenUEz4MBKpw3MGNMuLEAZv6upVR58eT0HS87w7N2TuXxoUqPrpvTowitfz2J4v+48+PJ6tuSf+GzhiTw4eTh4HSS8DZwOlSfh0Prg522MASxAmQD466r9rNpbws9vHkXWwOaHJuoeG8Wzd0+mZ9dovjJ/LUfLyp0FeZ86t+1xBpV5GSCQ82Hw8zbGABagjJ/tLz7N48t2cfWIZL4wKc3n7ZK7xzL/K5M5VVHFv7262bkqP38tRHWF5FEBLHEjuvaC/uOddihjTLuwAGX8RlX50aubiYn08ItbxrR4bLIhfbrz8OwRrNhdxAurc53qtX7jnRlv28PA6U6QrDjZPvkb08lZgDJ+8872o3x6oJQfzx5Bn/jWTZt+R1YG04cl8aulW9CCTc7gre1l4HSorYbcT9qvDMZ0YhagjF/U1Cq/fmcXA5O6cdtFre8SLiI89rmxDJM8pKaifQNU6hTwREHuyvYrgzGdmAUo4xdLNh1i99FTfP+aYURGtO1t1Tchlm8NKwPg06p2HOMsuiukXAQHVrVfGYzpxCxAmTarqVV+994eRvWPZ9bovn5J87KuuRwnnv/3QRnVNbV+SbNVMi+Fwxug4lT7lcGYTsoClGmz93YcJbfkDA/OGIzH459J2yIKNlLVdzx7i07z9/WH/JJmq2RcCloDeWvarwzGdFI+BSgRmSkiu0Rkr4g81MByEZE/uMs3i8hEr2XPikihiGytt83jIrLTXf91EenR5r0x7WL+qgOk9OjCtSP7+CfBipNQuIPEYVMZl5rA79/fQ2V1O51FpU0BiYBcq+YzJtiaDVAiEgE8AcwCRgK3i8jIeqvNAoa4f/cBT3otmw/MbCDpd4HRqjoW2A38uKWFN+1vR0EZn+SUcOfUjDa3PZ1TsAlQJGUS37tmKIeOn2VRdp5/0m6pmO7O9VDWDmVM0PnyjTIF2KuqOapaCSwE5tRbZw7wvDpWAz1EpB+Aqq4ASusnqqrvqGq1+3A1EMTRQI2/zF91gNgoD3Mn+35RbrMOb3Bu+0/giqFJTMroyR8/2EtFdY3/8miJjEudOakqz7RP/sZ0Ur4EqBTA++drvvtcS9dpyj3A2w0tEJH7RCRbRLJtSuXQcrK8ijc2HeKWCSn06Brtv4QLNkF8CsQlISJ85+ohHCkrZ/GGdmqLypwGtVXORbvGmKDxJUA11Opdf6pRX9ZpOHGRh4Fq4KWGlqvq06o6SVUnJSU1PuioCb63txyhvKqW21owpJFPCjY7U6+7pg1OZFT/eP60Iofa2naY5TY9C8Rj7VDGBJkvASof8P4GSgUOt2KdC4jIXcANwJdVbX7tcPPq+nwGJnZjQloP/yVaeRqKd0O/zwKUiPD1KwaRU3Sad3cc9V9evopNgL5jrB3KmCDzJUCtBYaIyAARiQbmAkvqrbMEmOf25ssCTqhqQVOJishM4EfATapqlfth5mDJGT7dX8qtF6W2eMy9Jh3dBij0G3fe07NH9yWtVxee+tc+/+XVEhnTnCq+6or2yd+YTqjZAOV2ZPgmsAzYASxS1W0icr+I1E03uhTIAfYCzwAP1G0vIguAT4BhIpIvIve6i/4IdAfeFZGNIvKUv3bKBN5r6/MRgVsmtKSp0QcFm5xbryo+gMgID1+dNpANB4+zKe+4f/P0RealUFPhdJYwxgSFT8NEq+pSnCDk/dxTXvcVeLCRbW9v5PnBvhfThBJVZfHGQ1wyqDf9e3Txb+IFm6BLrwaneP/cxBR+9c+dPP9JLr/2Z7WiL9LcOakOroaMS4KbtzGdlI0kYVpse0EZuSVnuHFsf/8nfmSz0/7UQLVh99goPjcxlX9sPkzp6Ur/592UbomQONQJUMaYoLAAZVrs7S1HiPAI147yz7h751RXwtHtF7Q/ebtzagaV1bW8srYdLtxNu9gZ8qi2HccGNKYTsQBlWkRVWbqlgKyBvejVzY/XPgEU7XSuN6rX/uRtaJ/uTB3YmxdX51IT7C7n6VOh/DgU7wpuvsZ0UhagTIvsPnqKnOLTzBrdz/+J13WQ6De+ydXmTc3g0PGzfLCz0P9laEp6lnN70CYwNCYYLECZFlm6pQARuM7f1XvgtD9Fx0GvgU2uds3IPvRLiOX5Tw74vwxN6TUQuiXDQRvZ3JhgsABlWmTZtiNMzuxFUvcY/ydesAn6jAZP02/LyAgPcyen89GeYvJKg3gJnQikX2xnUMYEiQUo47NDx8+y88hJrhnhp2k1vNXWwJGtTXaQ8Pb5SamIwN/W5fu/LE1JnwrHc6Gs2YFSjDFtZAHK+KyuzefKEcn+T7w0B6pOnzfEUVNSenThsiFJvJqdF9zOEufaoay7uTGBZgHK+OyDHUfJ7N2VgYnd/J94IyNINOWLk9I4fKKclXuL/V+exvQdC1FdbYZdY4LAApTxyZnKalbtK+HK4X38O/ZenaPbwBMJScN83uTqkcn07BoV3MkMI6Ig5SJrhzImCCxAGZ98vLeEyupargpE9R44ASpxKET63vkiJjKCmyek8O62oxwL5sgS6VPhyBZnanpjTMBYgDI+eX9nIXExkUzO7BWYDI5uc3rwtdAXJ6dRWVPL68GczDA9C7TWJjA0JsAsQJlmqSordhdx6eDeREcG4C1z9hiU5UOfUS3edHjfeMalJgS3mi91sjOBoV0PZUxAWYAyzTpQcoZDx88ybUiAZjQ+ut25bcUZFMCtF6Wy88hJdhSU+bFQTYiNd4KptUMZE1AWoEyzVu4pApyp1wPi6DbnthVnUAA3jO1PpEdYHNRqvqmQnw011cHL05hOxgKUadbKvcWk9OhCZu+ugcng6FZnDqjurRs+qVe3aKYPS2LxxkPBuyYqPcu5buvoluDkZ0wn5FOAEpGZIrJLRPaKyEMNLBcR+YO7fLOITPRa9qyIFIrI1nrb9BKRd0Vkj3vbs+27Y/ytuqaWj/eVMG1wYmC6l4PbQWJUg3NA+eqWCakcLavgk30lfixYE9Lsgl1jAq3ZACUiEcATwCxgJHC7iIyst9osYIj7dx/wpNey+cDMBpJ+CHhfVYcA77uPTYjZcugEJ8urmTYkQNV7tbVQuL3V7U91rhqRTPeYyOD15ktIgYR0a4cyJoB8OYOaAuxV1RxVrQQWAnPqrTMHeF4dq4EeItIPQFVXAKUNpDsHeM69/xxwcyvKbwJs5R5nlIZLBvUOTAbH9kPVmVa3P9WJjYpg9ph+/HNrAWcra/xUuGakZzlnUBrkeamM6SR8CVApgHcf3nz3uZauU18fVS0AcG8DdAWoaYuVe4sZ1T+e3nEBGL0c2txBwtstE1M4XVnDO9uPtDktn6RnwamjcOxAcPIzppPxJUA11DBQ/yejL+u0iojcJyLZIpJdVFTkjySNj05XVLP+4LHA9d4DJ0CJB5KGtzmpKZm9SOnRhb+vD1I1X/pU59baoYwJCF8CVD6Q5vU4Fag/14Av69R3tK4a0L1tcHpUVX1aVSep6qSkpABdh2Ma9OmBUqpqNHDtT+D04Os1CKLb3kPQ4xFuntCfj/YUUXSywg+Fa0bScIhNsHYoYwLElwC1FhgiIgNEJBqYCyypt84SYJ7bmy8LOFFXfdeEJcBd7v27gDdaUG4TBKv2FBMd6Qnc8EbwWQ8+P7llQgq1Cks2BWG+Jo8H0i62MyhjAqTZAKWq1cA3gWXADmCRqm4TkftF5H53taVADrAXeAZ4oG57EVkAfAIME5F8EbnXXfQYcI2I7AGucR+bEPJJTgkXpfckNioiMBlUnHI6SbSxB5+3wcndGZOSELyLdtOzoHgXnGmoH5Axpi0ifVlJVZfiBCHv557yuq/Ag41se3sjz5cAV/lcUhNUJ85Wsb2gjO9cNSRwmRTucG79eAYFMGd8fx59awf7ik4xKCnOr2lfwLsdavjswOZlTCdjI0mYBmUfKEUVLh4QoO7l4LQ/gd8D1E3j+uMReCMYZ1H9J0JEtLVDGRMAFqBMg9bsLyU6wsOE9B6By+ToNojuDj3S/ZpscnwslwxKZPHGw2igr1GKinWClLVDGeN3FqBMg9bklDAuLSFw7U/glyGOGjNnfH8Olp5hQ95xv6d9gfQsOLwBqs4GPi9jOhELUOYCpyqq2Xq4LLDVe6p+78HnbebovkRHeoJTzZc+FWqr4NC6wOdlTCdiAcpcYF3uMWpqlYsHBrB7+Yl8qDgRsADVPTaKq0ck8+bmAqpqagOSxznpFzu31g5ljF9ZgDIXWJNTQqRHuCgjgAPMnxviyH9dzOubMz6FktOVrNpbHLA8AOjSE5JHQq4FKGP8yQKUucCa/aWMSU2ga7RPVyG0Tl0PvuQRActi+rAk4mMjeWNjEC7aTc+CvE+hNkgD1RrTCViAMuc5W1nD5vzjTBkQwOo9cM6gemQ406cHSExkBNeP7ceybUc4UxngmW/TL4HKk5+dGRpj2swClDnP+oPHqKpRsgLZQQLcDhKBq96rM2d8Cmcqa3h3+9HAZpReN4GhVfMZ4y8WoMx51uwvxSMwKTOA7U9V5VCyJ2AdJLxNyexF/4TYwFfz9UiD+FQLUMb4kQUoc57sA6WM6BdP99iowGVStBO0NigByuMRbhzfnxW7iyg9XRnYzGwCQ2P8ygKUOae6ppaNeceZFMjeexCUHnzebh6fQnWt8tbmAJ9FZUyFkwU2gaExfmIBypyzo+AkZypruCiQ02uAE6Aiu0CvAYHNxzWiXzzD+nRncaCr+WwCQ2P8ygKUOWddrjNlRECvfwKni3nyCPAEcBilem4a3591ucfIKz0TuEySRtgEhsb4kQUoc0527jH6JcSS0qNL4DJRdQJUENqfvM0Z3x8I8ESGHg+kZVmAMsZPLECZc9blHgv82dOpQjhTErT2pzqpPbsyObMnizccCuwI5+lZULwbTgd49ApjOgELUAaAw8fPUnCiPAjVe1uc2z4jA5tPA+aMT2FP4Sm2F5QFLpO6dqi8NYHLw5hOwqcAJSIzRWSXiOwVkYcaWC4i8gd3+WYRmdjctiIyXkRWi8hGEckWkSn+2SXTGtm5xwCYlBGEDhIQ9DMogOvH9CPSI4G9JiplIkTEWDWfMX7QbIASkQjgCWAWMBK4XUTq//ydBQxx/+4DnvRh218BP1XV8cBP3Memnaw7UErX6AhG9Ose2IyOboP4FOga4EDYgJ7dopk+LIklGw9TUxugar7IGCdI2cCxxrSZL2dQU4C9qpqjqpXAQmBOvXXmAM+rYzXQQ0T6NbOtAnUDsSUAQRjR0zRm3cFjjE/rQWREgGt9jwS/g4S3OeNTOFJWzpr9JYHLJH0qFGyEilOBy8OYTsCXb6MUIM/rcb77nC/rNLXtd4HHRSQP+B/gxw1lLiL3uVWA2UVFRT4U17TU6YpqdhScDHz7U3UlFO9ql+q9OleP6EO36Aje2BDA30OZ06C22tqhjGkjXwJUQ/Nx168faWydprb9BvA9VU0Dvgf8paHMVfVpVZ2kqpOSkpJ8KK5pqY15x6mp1cAHqOJdzhd3O55BdYmO4LpRfVm6tYDyqgBNjZF2MXgi4cDKwKRvTCfhS4DKB9K8HqdyYXVcY+s0te1dwN/d+3/DqQ407SD7wDFEYGIHG+KoMXMmpHCyvJrluwoDk0FMHPSfCLmrApO+MZ2ELwFqLTBERAaISDQwF1hSb50lwDy3N18WcEJVC5rZ9jBwhXv/SmBPG/fFtNK6g8cY1qc78YEcIBacC3QjYqD34MDm04xLB/UmMS46sL35MqfBoXVQeTpweRjTwTUboFS1GvgmsAzYASxS1W0icr+I3O+uthTIAfYCzwAPNLWtu83XgF+LyCbgv3B6/5kgq6lVNuQeC/zZEzgdJJKHQ0QAZ+r1QWSEhxvG9uf9nYWUlVcFJhNrhzKmzXz6plDVpThByPu5p7zuK/Cgr9u6z68ELmpJYY3/7T56kpMV1YEfwRycKr4h1wY+Hx/cPCGF+R8f4J9bjvCFyWnNb9BS3u1Qg670f/rGdAI2kkQnty5YF+ieKoTThe3aQcLbuNQEMnt3ZfHGQ4HJoK4dyjpKGNNqFqA6uXW5x0iMiyGtVwAHiAWn/QlCJkCJCHPGp/BJTglHTpQHJpPMS60dypg2sADVyWXnljIpoyciDV0R4Ech0oPP25zx/VGFfwRqhHNrhzKmTSxAdWKFZeXklZ5lUmaQOkh07wfdegc+Lx8NTIpjbGpC4Kr50rJAIqyaz5hWsgDVidW1PwX8Al1wzqBC6OypzpzxKWw7XMbewpP+TzwmzhmXzwKUMa1iAaoTy849Rkykh1H9EwKbUU0VFO0MmfYnbzeO64dHYHGghj7KnAaH1ls7lDGtYAGqE8vOPca41B5ERwb4bVC8G2qroO+YwObTCsndY7l0cCJvbArQRIaZ05x9z/vU/2kb08FZgOqkyqtq2HboBBcFo/3pXAeJ0DuDAqeaL6/0LOsPHvd/4tYOZUyrWYDqpDblHae6VrkoPRgdJLZARHS7D3HUmOtG9SEm0sMbgegsca4d6iP/p21MB2cBqpPKDnYHiaThEBHgsf5aqXtsFFeP7MObmwuoqqn1fwYDLof8bCgP4FTzxnRAFqA6qXW5xxiU1I2e3aIDn1mI9uDzdvP4FEpPV7JyT7H/Ex84A7TGRjc3poUsQHVCtbXKutxjgR/eCOB0MZw6ErLtT3WuGJpEQpeowFwTlTYForrCvg/9n7YxHZgFqE4op/gUJ85WBaeDRMEm57bf2MDn1QbRkR5mj+nHO9uOcrqi2r+JR8ZAxiWQYwHKmJawANUJZR+oGyA2iAGqb2gHKICbx/fnbFUNy7Yd8X/iA2c43e1PBGjUCmM6IAtQnVB27jF6dYtmQGK3wGd2ZDP0yIAuPQKfVxtNzuxFeq+uLMrO83/iA6c7tznL/Z+2MR2UBahOaF3uMSamB2GAWICCzSFfvVfH4xG+ODmN1TmlHCj288gPfUZBtySr5jOmBXwKUCIyU0R2icheEXmogeUiIn9wl28WkYm+bCsi33KXbRORX7V9d0xzik9VsL/4dHAGiC0vg9J90G9c4PPyk1snpuIR/H8WJeKcReUsh9oAdGU3pgNqNkCJSATwBDALGAncLiIj6602Cxji/t0HPNnctiIyA5gDjFXVUcD/+GOHTNM+m6AwGNc/uXNA9Q2fANU3IZYZw5J5dV0+1f6+JmrgDDhdBIXb/ZuuMR2UL2dQU4C9qpqjqpXAQpzA4m0O8Lw6VgM9RKRfM9t+A3hMVSsAVLXQD/tjmrEu9xjRER5GpwR4gFgImx589X1hchqFJytYvqvIvwmfa4eyaj5jfOFLgEoBvOs78t3nfFmnqW2HApeJyBoR+ZeITG4ocxG5T0SyRSS7qMjPXxidUPaBUsakJhAbFRH4zAo2Q1wf6N438Hn50ZXDk0mMi2HhWj9X8yWkQOJQ2PeBf9M1poPyJUA11JJef9jnxtZpattIoCeQBfwbsEgaaLVX1adVdZKqTkpKSvKhuKYx5VU1bD1UFpzhjcDpwRcG3cvri4rwcOtFKXy4q5DCMj9PBz/4ajiwyqbfMMYHvgSofCDN63EqUH/ynMbWaWrbfODvbrXgp0AtkOh70U1LbT10gsqa2uAEqKpyKNwRVh0kvH1xUho1tcpr6/183dKQa6GmAvbb4LHGNMeXALUWGCIiA0QkGpgLLKm3zhJgntubLws4oaoFzWy7GLgSQESGAtFAAAZCM3WCOkBs4XZn/Lkwa3+qMzApjimZvViUneffeaIyLoGobrDnHf+laUwH1WyAUtVq4JvAMmAHsEhVt4nI/SJyv7vaUiAH2As8AzzQ1LbuNs8CA0VkK07nibs0IDPGmTrZB44xILEbiXExgc8sjEaQaMwXJqexv/g0q3NK/ZdoZIzTWWLPu2Bvd2OaFOnLSqq6FCcIeT/3lNd9BR70dVv3+UrgjpYU1rSeqrL+4DGuHJ4cnAyPbIaYBOiZGZz8AuD6Mf34+ZvbeXF1LlMH9fZfwkOvhV1vQdFOSB7hv3SN6WBsJIlOIqf4NKWnK4Nz/RN8NoJEMEarCJAu0RF8YVIqy7Yd4ag/O0sMvsa5tWo+Y5pkAaqTWFc3QGwwRpCoqXYu0g3j6r06d2RlUKPKy2sO+i/RhBRnfqw97/ovTWM6IAtQnUR2bikJXaIYmBgX+MyKdkB1uTPVeZjL6N2NK4YmseDTg/6dbXfINXDwEyg/4b80jelgLEB1Etm5x7gooyceTxCq3A6tc247QIACmDc1g8KTFf6dhmPIdVBbbZMYGtMEC1CdQOnpSnKKTgfvAt1D66BLT+g5IDj5BdgVQ5NJ69WF5z/J9V+iqZMhNsHaoYxpggWoTmB9MAeIBTi0HlIuCusOEt4iPMIdF2fw6f5Sdh4p81Oikc5Fu7vedtrsjDEXsADVCazZX0J0hIdxaT0Cn1nlaeci3ZSLAp9XEH1hUhoxkR5e8OdZ1Igb4Wyp0xZljLmABahOYM3+Usan9QjSALGbQGs7XIDq2S2aG8f15/UNhygrr/JPooOvhshY2Pmmf9IzpoOxANXBnSyvYuuhE1w8sFdwMqzrING/Y3SQ8Hb3JZmcqaxh4ad+6nIe3Q0GXQk73rRRJYxpgAWoDm5d7jFqFS4e4MeREJpyaB30SIe4jjfy/OiUBLIG9uKvqw74r8v5iBuhLB8Ob/BPesZ0IBagOrg1+0uJ9AgTM3oEJ8ND6zpc9Z63+y4fSMGJcpZuKfBPgkNngkTAjn/4Jz1jOhALUB3cmpwSxqQm0DXap2EX2+ZUERw/2KED1PShyQxK6sbTK3L8M8p5116QOc3aoYxpgAWoDuxMZTWb808Er3rv8HrntgMHKI9H+NplA9l2uIxPckr8k+iIG6F4NxTt8k96xnQQFqA6sA0Hj1Ndq1w8IIgdJMQTtpMU+urmCSkkxkXzzIoc/yQ4/Hrn1qr5jDmPBagObE1OCR4J0gCx4ASo5JFO77QOLDYqgnlTM/lwVxG7jpxse4Lx/SF1Cmx7ve1pGdOBWIDqwFbvL2VU/wS6x0YFPrPaWidA9Z8Q+LxCwLypGXSLjuCPH+71T4JjbnNGgD+63T/pGdMBWIDqoMqratiYd5wpwareK9kDZ49BelZw8mtnPbpGc+fUTN7cfJh9RafanuCoW5zefFv+1va0jOkgfApQIjJTRHaJyF4ReaiB5SIif3CXbxaRiS3Y9gcioiKS2LZdMd425h2nsro2eO1PdcP1pHWOAAXw1csGEBPp4cnl+9qeWFwSDJoBW161i3aNcTUboEQkAngCmAWMBG4XkZH1VpsFDHH/7gOe9GVbEUkDrgH8OBucAVi1txiPwMUDg9SD7+Bq6JoIvQcFJ78QkBgXw+1T0nl9wyHySs+0PcExt8GJg5D3advTMqYD8OUMagqwV1VzVLUSWAjMqbfOHOB5dawGeohIPx+2/S3wQ8B+MvrZyr3FjE3tQUKXILQ/gROg0rM6zAjmvvr65YOIEOHJf/nhLGr49RDZBbYsantaxnQAvgSoFCDP63G++5wv6zS6rYjcBBxS1U1NZS4i94lItohkFxUV+VBcU1Zexaa841w2JEi1piePwLH9nab9yVvfhFhum5TK37LzyD/WxrOomO4wbJbTm6/GTwPSGhPGfAlQDf0krn/G09g6DT4vIl2Bh4GfNJe5qj6tqpNUdVJSUscb3y0QVu8roVbh0sFBClAHVzu36VODk1+IeXDGYESE37+3p+2JjbkNzpRAzvK2p2VMmPMlQOUDaV6PU4HDPq7T2PODgAHAJhE54D6/XkT6tqTwpmEr9xbTJSqCCek9gpPhwdVO1VTfscHJL8T079GFO7MyeG19PnsL23hd1OCrndmIN77kn8IZE8Z8CVBrgSEiMkBEooG5wJJ66ywB5rm9+bKAE6pa0Ni2qrpFVZNVNVNVM3EC2URVPeKvHevMVu4t5uKBvYiJDML8TwB5qyF1EkRGBye/EPTA9EF0iYrg1+/sbltCkdEw7nZnCo7Txf4pnDFhqtkAparVwDeBZcAOYJGqbhOR+0Xkfne1pUAOsBd4BnigqW39vhfmnMPHz5JTdJppwareKz/hTFLYSav36vSOi+Heywby9tYjbM4/3rbEJt4FtVWwaYFfymZMuPLpOihVXaqqQ1V1kKr+wn3uKVV9yr2vqvqgu3yMqmY3tW0D6Weqqv1c9INVe52XMWjtT7kfOzPoDrg8OPmFsK9dNoCeXaN4fFkbB31NHg5pF8O65+yaKNOp2UgSHcyqvcUkxkUzrE/34GS4/yOIiIHUycHJL4R1j43iwRmD+WhPMR/uKmxbYhfd7YzOkfuxX8pmTDiyANWBqCor95Zw6eBEPJ4gXY+0fwWkXwxRscHJL8TNm5rJgMRuPPrm9rbNujvyZohJgPXP+a1sxoQbC1AdyM4jJyk+VRG86r3TJXB0C2Ra9V6d6EgPD88ewb6i07y4OrcNCXWFsbfBtsVwptRv5TMmnFiA6kA+2OlUK00fGqTrxXJXOrfW/nSeq0Ykc9mQRH733h6Ona5sfUIX3Q01FbBpod/KZkw4sQDVgXy4s5DRKfEkxwepum3/CojqBikTm1+3ExER/v36kZwsr+K377Wh23nfMU5niU//BLU1/iugMWHCAlQHcex0JesPHuPKYcnByzTnX5AxFSKCNN5fGBnWtzt3ZGXw4urctnU7z3oAjh2A3f/0V9GMCRsWoDqIFXuKqFWYMTxIAerYAaeX2aCrgpNfGPrBdcNIjIvhode2UN3aDhPDb4CEdPjk//xbOGPCgAWoDuKDnYX07hbNuNQewclwz7vO7ZBrg5NfGIqPjeKRm0axvaCMv6460LpEIiLh4vuc9r6CJsdVNqbDsQDVAdTUKv/aXcQVw5KC1718z7vQM7NTzf/UGrNG9+Wq4cn85t3drZ8zasKdEN0dVv7Or2UzJtRZgOoA1h4o5fiZKq4a3ic4GVaVOx0khlzb6eZ/aikR4Wc3j0YEHl68FW3NyBBdesDke51pOIr3+r2MxoQqC1AdwD+3HiEm0sP0YUHsXl591qr3fJTSows/mjmcFbuLeGlNKyePnvogRMbAyt/6t3DGhDALUGGutlZZtu0Ilw9NoltMZHAy3fMeRMZC5rTg5NcB3JmVwWVDEvnFWzvYX3y65QnEJTuDyG5eCMdbGeSMCTMWoMLc5kMnKDhRzsxRQZpKSxV2vw2Zl0FUl+Dk2QF4PMLjnx9HdKSH772ysXW9+i79NogHVjzu/wIaE4IsQIW5t7cWEOkRrh4RpPanI1ucLuYjbgxOfh1I34RYHr15NBvzjvPHD1vRlpSQCpPugQ0vWVuU6RQsQIUxVWXZ1iNMHdSbhK5Bulh2xxLnV/zw64OTXwdz47j+3DIhhT+8v+fc1CgtctkPnOrVDx/1f+GMCTEWoMLYtsNlHCg5w+wx/YKX6fY3IONS6BakAWk7oEdvHs3ApDi+vWADR06Ut2zjuCSY+oDTo+/whsAU0JgQ4VOAEpGZIrJLRPaKyEMNLBcR+YO7fLOITGxuWxF5XER2uuu/LiI9/LJHncjiDYeIjvAwe3SQAlThTijeDSPnBCe/DqpbTCRP3TGRs1U1PPjy+pZPy3HJt6BrIrz9kE1oaDq0ZgOUiEQATwCzgJHA7SIyst5qs4Ah7t99wJM+bPsuMFpVxwK7gR+3eW86kZpa5Y1Nh5kxPCm41XuItT/5weDk7vz3rWNZl3uMR9/c3rKNYxPgqp9A3mrY+lpgCmhMCPDlDGoKsFdVc1S1ElgI1P8JPQd43p36fTXQQ0T6NbWtqr6jqtXu9quBVD/sT6fx8b5iik5WcMuElOBkqOp8GaZnQfcg9Rjs4G4c15+vThvAc5/kMn/V/pZtPOEO6DcO3v0JVLai27oxYcCXAJUC5Hk9znef82UdX7YFuAd4u6HMReQ+EckWkeyioiIfits5vL7hEPGxkUwP1ujlh9dD0U4YNzc4+XUSP549gmtG9uFnb27n/R1Hfd/QEwGzfgVlh2D5LwNXQGPakS8BqqGxbOpXfDe2TrPbisjDQDXwUkOZq+rTqjpJVSclJQVppIQQd7qimmVbjzB7TD9ioyKCk+nGl53eY6NuCU5+nUSER/j93PGM6p/AtxZsYOuhE75vnJ4FF30FPnkCDq0PXCGNaSe+BKh8IM3rcSpw2Md1mtxWRO4CbgC+rK0apKxzWrzxEKcra/jC5LTmV/aHqnLY8jen7Sk2ITh5diJdoyP5y12T6Nk1mrue/ZQ9R0/6vvE1P4VuybDkW1BTFbhCGtMOfAlQa4EhIjJARKKBucCSeussAea5vfmygBOqWtDUtiIyE/gRcJOqtnKY585HVXlx9UFG9ItnQlqP4GS6aymUn4DxXwpOfp1QcnwsL371Yjwe4ct/XsMBX4dDik2A638NR7daVZ/pcJoNUG5Hhm8Cy4AdwCJV3SYi94vI/e5qS4EcYC/wDPBAU9u62/wR6A68KyIbReQp/+1Wx7Uh7zg7Csq4IysdCdZI4uufh/gUGHBFcPLrpAYkduOlr15MVU0tX/7zGt+n5xhxg9Np4qPfwIGVgS2kMUEk4VSzNmnSJM3Ozm7vYrSr7y/axD+3FrDm4auJC8bgsIU74P+y4Mp/h8v/LfD5GbYeOsGXnllNt5hIXrh3CoOTuze/UcUp+NPlUF0O96+Err0CX1Bj/ERE1qnqpPrP20gSYaToZAX/2HyYWyamBCc4Aax+0ukccdE9wcnPMDolgYX3TaWqRrntqU/YmHe8+Y1i4uDWP8PpInj1HqitCXg5jQk0C1Bh5NlV+6mqqeXeaQODk+HpEtj8itO1vFvv4ORpABjZP57XvjGVuNhIvvTMapbvKmx+o5SJTntUzofw/s8CX0hjAswCVJgoK6/ixU9ymT26HwMSuwUn0+xnnSqji78RnPzMeTJ6d+O1+y8ho3c37pm/lj9/lNP8jLwT5zkjnq/6HWx4MSjlNCZQLECFiRc+yeVkRTXfmD4oOBmePQ6f/BGGXAfJw4OTp7lAcnwsr31jKteN6sujb+3g+3/bRHlVM9V3M/8bBl0JS74Nu98JTkGNCQALUGHgVEU1z67cz2VDEhmdEqTrkD7+Xyg/7nSOMO2qa3QkT3xpIt+7eih/X3+Im59Yxe6mrpWKjIYvPA99R8Pf7rKefSZsWYAKA08u30vJ6Up+cO2w4GR4qtDpHDHqc9BvbHDyNE3yeITvXD2E+V+ZTPGpCm7835W8sDq38Sq/mO7w5VchIQ1eus2ClAlLFqBC3OHjZ/nzR/uZM74/44J1Ye7yx5y2pxkPByc/47Ppw5J5+zuXkzWwN/+xeGvTF/XGJcPdbzpB6sXPw65/BrewxrSRBagQ9/iyXSjwb9cF6ewp71Onc8SUr0Hi4ODkaVokqXsMf717Mo/ePJot+Se47ncreOLDvVRUN9A2FZcMd78FScNg4e2w9i/BL7AxrWQBKoQt31XI6xsO8bXLBpDas2vgM6yudMZ0i0+xtqcQ5/EId2Rl8N73r+DK4ck8vmwXV//mX7y5+fCF1X5xSU6QGnw1vPX/wRsPQqWNLmZCnwWoEHXibBUPvbaFwclxfOvKIcHJ9F+POVNq3PAbpw3DhLw+8bE8ecdFPHfPFLpFR/LNlzdwy/99zPJdhecHqpg4uH0hXP5D2PAS/PlqKN7TfgU3xgcWoEKQqvLTJdsoOlXBr28bF5wpNXa9DR/92hnTbeh1gc/P+NUVQ5N469uX8atbx1JYVs7df13LTX9cxdtbCqium1LeEwFXPgx3vAqnjsBTl8HK39ko6CZk2Vh8IejZlfv52Zvb+e7VQ/ju1UMDn2HJPnh6BvTKhHuWQVSXwOdpAqayupbXN+Tzf8v3kVtyhn4Jsdw+JZ25k9NIjo91Vio7DEv/DXa+CcmjnLPm9Kz2LbjptBobi88CVIhZvquQe+av5ZqRfXjyyxfh8QR4xPJjB2D+Dc604fcth54Zgc3PBE11TS3v7SjkpTW5fLSnmEiPcNWIZG4al8KVw5PpEh0BO96Et3/ozMw75DqY8f+g//j2LrrpZCxAhYEVu4u474VsBiTG8er9U+kW6AFhS/fDczdBRRnc9Q+75qkD2198mpfX5PL6hsMUn6qgS1QEV41IZtbofkxLjyVhy7Ow6g/OxdlDZ8KU+2DgDPBYK4AJPAtQIW7plgK+u3Ajg5LjeOHeKSTGxQQ2wz3vwWv3OvfnvWG/mjuJmlplzf4S3tpcwNtbj1B6uhKPwPi0Hlw1oAtzyheTsucl5Ewx9BoEE++EUbdAz8z2LrrpwCxAhajyqhoee3sn8z8+wMT0Hvz17ikkdI0KYIYn4MNfwpqnoM8o+OIL0CtIo6ObkFJdU8vGvOOs2F3Ev/YUszn/OKrQ1VPNPb22cJsuI+P0ZgC0/0XI8NkwaAb0G+90uDDGTyxAhRhVZemWIzy+bCcHSs5wz6UD+NGsYcREBuiDX17mzIy76vfOnEGT7oFrH4XoIFxfZcLCsdOVrMs9xoa8Y2w4eJxNecfpWXWE2Z7V3BixhjGeHADKIxM4njgR7T+BbgOm0H3QxYhNkGjaoE0BSkRmAr8HIoA/q+pj9ZaLu3w2cAa4W1XXN7WtiPQCXgEygQPAF1T1WFPl6AgBquDEWd7aXMAra/PYU3iKoX3i+PfrR3L50CT/Z3amFA5+Ajv+4TSGV56EjGlw3aPQf4L/8zMdSk2tsvvoSXYUlLHr6EmOHM6jx5FPGHV2HRM9exgoBXjE+f4okV4URqdzvFsGZ+MHoPFpRPVMpWvvVOKTUunVvQvxsVFER1qblrlQqwOUiEQAu4FrgHxgLXC7qm73Wmc28C2cAHUx8HtVvbipbUXkV0Cpqj4mIg8BPVX1R02VJZQDlKqiCtW1Sll5FSfOOn9FJyvYW3iKfUWn2JR3nH1Fzrhp49N6cGdWBjdPSCGifk+9umNy7tg08LjqDFScdKb6rjjpVN2VHYITeU7PvEProXSfs3psAgy/ASZ/1ZnUzpg2OFVRTV7pGY4cLaQibx2RRzbRtWwvPc8cIKUmn3jOHxuwRoUienBM4zgpcZz2dKc8ojvlUfFURcYj0V0gqqtzeUNUFySqC57orkhUFyKiuxAdHUVkZLR7G4VERhEZEYUnIpKIyCg8kVFERETiiYwiKjKCCE8EkR4hIsLj3HoEjzh/CLg3eETc++6te98jICII7nMS4J60pk0BairwiKpe5z7+MYCq/tJrnT8By1V1gft4FzAd5+yowW3r1lHVAhHp527f5IBzbQ1Qmx+7mvTyHYj7hS/qfPWfe3zuts75z3sv83Wbc+uLOm943Df9ea+7H6tZxeMMVdR3LKReBKmTIS3LmYLBmEBTpfZUMaeKcjlVlEd5aR7Vxw8hJwuIqDhOZGUZ0VUniKkuo2vNSWK0IuBFqlXn06iI++fcdzS8TOt96s/f5kLBDmEXlKUdY+iBy37N+KvmtimNxgKUL/2YU4A8r8f5OGdJza2T0sy2fVS1AMANUsmNFPw+4D6A9PR0H4rbuDNpl7HrWLoTIc5lUK/KQQSpO9oC4Dnv4Dt3Pd6R59yziPOrLDoygphIDzFRkXSJiqBHt2iiIjxe63m/m+o/5+PjqC7OcEQxcRAT79zv3g/i+0NEADtZGNMUETzdk4jvnkT8wAu+by5UXQnVZ6Hq/L/qyjNUl5+mquIMldXVVFdVUV1VRU11FbW1VWh1NVpbjdY4t7U11VBThdbWUquKqlKrtdTWKqq1Tu2D+3zdD8Jz971vve6rs5LXOnBua/ffZ/cb19afn/XPIaSFKbb5528zJzGJSW37Xm6KLwGqodhcv8SNrePLtk1S1aeBp8E5g2rJtvVlffk/27K5McbfIqOdv9jzJ+KMdP9i26VQJlT40mKZD6R5PU4FDvu4TlPbHnWr9nBvC30vtjHGmI7OlwC1FhgiIgNEJBqYCyypt84SYJ44soATbvVdU9suAe5y798FvNHGfTHGGNOBNFvFp6rVIvJNYBlOV/FnVXWbiNzvLn8KWIrTg28vTjfzrzS1rZv0Y8AiEbkXOAjc5tc9M8YYE9bsQl1jjDHtqrFefHbVnDHGmJBkAcoYY0xIsgBljDEmJFmAMsYYE5LCqpOEiBQBuW1MJhEo9kNx2pPtQ/sL9/KD7UOosH2ADFW9YMTssApQ/iAi2Q31Fgkntg/tL9zLD7YPocL2oXFWxWeMMSYkWYAyxhgTkjpjgHq6vQvgB7YP7S/cyw+2D6HC9qERna4NyhhjTHjojGdQxhhjwoAFKGOMMSGp0wQoEZkpIrtEZK+IPNTe5fGViBwQkS0islFEst3neonIuyKyx73t2d7l9CYiz4pIoYhs9Xqu0TKLyI/d47JLRK5rn1Kfr5F9eEREDrnHYqOIzPZaFlL7ICJpIvKhiOwQkW0i8h33+bA5Dk3sQzgdh1gR+VRENrn78FP3+XA6Do3tQ+CPg7rTIHfkP5ypPvYBA4FoYBMwsr3L5WPZDwCJ9Z77FfCQe/8h4L/bu5z1ync5MBHY2lyZgZHu8YgBBrjHKSJE9+ER4AcNrBty+wD0Aya697sDu91yhs1xaGIfwuk4CBDn3o8C1gBZYXYcGtuHgB+HznIGNQXYq6o5qloJLATmtHOZ2mIO8Jx7/zng5vYryoVUdQVQWu/pxso8B1ioqhWquh9nTrEpwShnUxrZh8aE3D6oaoGqrnfvnwR2ACmE0XFoYh8aE4r7oKp6yn0Y5f4p4XUcGtuHxvhtHzpLgEoB8rwe59P0Gz2UKPCOiKwTkfvc5/qoM2Mx7m1yu5XOd42VOdyOzTdFZLNbBVhXLRPS+yAimcAEnF++YXkc6u0DhNFxEJEIEdkIFALvqmrYHYdG9gECfBw6S4CSBp4Ll/71l6rqRGAW8KCIXN7eBfKzcDo2TwKDgPFAAfBr9/mQ3QcRiQNeA76rqmVNrdrAc6G6D2F1HFS1RlXHA6nAFBEZ3cTq4bQPAT8OnSVA5QNpXo9TgcPtVJYWUdXD7m0h8DrOqfJREekH4N4Wtl8JfdZYmcPm2KjqUfeDWgs8w2fVFiG5DyIShfPF/pKq/t19OqyOQ0P7EG7HoY6qHgeWAzMJs+NQx3sfgnEcOkuAWgsMEZEBIhINzAWWtHOZmiUi3USke9194FpgK07Z73JXuwt4o31K2CKNlXkJMFdEYkRkADAE+LQdytesui8U1y04xwJCcB9ERIC/ADtU9Tdei8LmODS2D2F2HJJEpId7vwtwNbCT8DoODe5DUI5De/YOCeYfMBunF9A+4OH2Lo+PZR6I0xtmE7CtrtxAb+B9YI9726u9y1qv3AtwTvmrcH5N3dtUmYGH3eOyC5jV3uVvYh9eALYAm90PYb9Q3QdgGk61ymZgo/s3O5yOQxP7EE7HYSywwS3rVuAn7vPhdBwa24eAHwcb6sgYY0xI6ixVfMYYY8KMBShjjDEhyQKUMcaYkGQByhhjTEiyAGWMMSYkWYAyxhgTkixAGWOMCUkWoIwJISKS6c5/9Iw798477tX7xnQ6FqCMCT1DgCdUdRRwHLi1fYtjTPuwAGVM6Nmvqhvd++uAzPYrijHtxwKUMaGnwut+DRDZXgUxpj1ZgDLGGBOSLEAZY4wJSTaauTHGmJBkZ1DGGGNCkgUoY4wxIckClDHGmJBkAcoYY0xIsgBljDEmJFmAMsYYE5IsQBljjAlJ/z/LTG8v7cyqagAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"from utils import decorate\n",
"\n",
"posterior_n.plot(label='Model A')\n",
"posterior_n2.plot(label='Model B')\n",
"decorate(title='Posterior distributions of n1 under different models')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the additional information from `c2`:\n",
"\n",
"* We give higher probability to large values of $s$, so we also give higher probability to large values of `n1`, and\n",
"\n",
"* The width of the distribution is narrower, which shows that with more information about $s$, we have more information about `n1`.\n",
"\n",
"This is one of several examples Jaynes uses to distinguish between \"logical and causal dependence.\" In this example, causal dependence only goes in the forward direction: \"$s$ is the physical cause which partially determines $n$; and then $n$ in turn is the physical cause which partially determines $c$\".\n",
"\n",
"Therefore, `c1` and `c2` are causally independent: if the number of particles counted in one interval is unusually high (or low), that does not cause the number of particles during any other interval to be higher or lower.\n",
"\n",
"But if $s$ is unknown, they are not *logically* independent. For example, if `c1` is lower than expected, that implies that lower values of $s$ are more likely, which implies that lower values of `n2` are more likely, which implies that lower values of `c2` are more likely.\n",
"\n",
"And, as we've seen, it works the other way, too.\n",
"For example, if `c2` is higher than expected, that implies that higher values of $s$, `n1`, and `c1` are more likely.\n",
"\n",
"If you find the second result more surprising -- that is, if you think it's weird that `c2` changes what we believe about `n1` -- that implies that you are not (yet) distinguishing between logical and causal dependence."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Tags",
"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.8.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
![](\"https://github.com/AllenDowney/ThinkBayes2/raw/master/examples/jaynes177.png\")
\n",
"\n",
"For A, the \"different intervals are logically independent\", so the update with `c2 = 16` starts with the same prior."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"c2 = 16\n",
"posterior2 = update(joint, phi, c2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the posterior marginal distribution of `n2`."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"168.947980523708"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAp5ElEQVR4nO3deXxV9Z3/8dcnN3tCEiAJCSGQAAFkUYSAKFqV1grYilo71S7a1bFqZ9pOO7XT+XXamV87Haedzji1+tNOFzu2ajelFbWudUUIyh6WGBICSUhCNiBk//7+yI3GmOWGLOcu7+fjkUfuck7u+x5I3vd877nfY845REQk8kR5HUBERLyhAhARiVAqABGRCKUCEBGJUCoAEZEIFe11gJFIT093eXl5XscQEQkp27Ztq3POZfS/PaQKIC8vj6KiIq9jiIiEFDMrH+h2DQGJiEQoFYCISIRSAYiIRCgVgIhIhFIBiIhEKBWAiEiEUgGIiESokPocgMh46ezqZsuhenYebaKlvYuctHhWz01nxuREr6OJjBsVgES8Tbuq+M5jxRxtPP2u+96/cBrfunIR09MSPEgmMr5UABKxWju6+KdH9/BQUQWLpqfwjSvO4sKCdJJjoymtO8XGHZX85MVSLv/hC3znmiVcec50ryOLjCkVgESkjq5ubvvV6zxdXMMtl8zhS5fNI8b39lticzOT+fJl8/jQshz+7uEd/M2v36CxpZ0bzs/zLrTIGNObwBJxnHP83cM7eLq4hn+5ajF/v3bBO/749zVrahIPfO483nfWNL756B5+uXnAKVVEQpIKQCLOz14uY+OOSr56+Xw+sWrWsMvHRfu4++PLeO+CTL61cQ+vvFk3ASlFxp8KQCLK7qNNfO/xfbzvrExuuWROwOvF+KL4z+uWkp+exK0PvD7gG8YioSagAjCztWa238xKzOz2Ae43M7vTf/9OM1vW576fmlmNme3ut84UM3vKzA76v08e/dMRGVxnVzdf+c0O0hJjuOPaczCzEa0/KT6G+24opL2zm7//7Q6cc+OUVGRiDFsAZuYD7gLWAQuB681sYb/F1gEF/q+bgLv73PdzYO0AP/p24BnnXAHwjP+6yLh54LXD7Ks+wT9vWMSUpNgz+hn56Un8wxVn8XLJcR547fAYJxSZWIHsAawESpxzpc65duBBYEO/ZTYA97sem4E0M8sGcM69ANQP8HM3AL/wX/4FcNUZ5BcJSP2pdn7w5/1cODedyxdljepnfXTlTC6cm853NxVT3dQ6RglFJl4gBZADVPS5fsR/20iX6W+ac64KwP89c6CFzOwmMysys6La2toA4oq8253PHORUexf/9MGFIx766c/M+O7VS+jsdvzbE/vGKKHIxAukAAb6bek/+BnIMmfEOXevc67QOVeYkfGuU1qKDKumuZVfbTnMtctmUDBt0pj8zJlTE/ncRfn84Y2jbCtvGJOfKTLRAimAI0Bun+szgMozWKa/Y73DRP7vNQFkERmxe/5SSle349ZL547pz73lkrlkTorjX/60V28IS0gKpAC2AgVmlm9mscB1wMZ+y2wEbvAfDbQKaOod3hnCRuBG/+UbgUdHkFskILUn2njgtXKuPjeHmVPHdmK3pLhovnzZPLZXNPLsPr1+kdAzbAE45zqB24AngWLgYefcHjO72cxu9i+2CSgFSoD7gFt61zezXwOvAvPN7IiZfcZ/1/eAy8zsIHCZ/7rImHrgtXLau7pHdMz/SHxo+QxmTknkB38+QHe39gIktAQ0F5BzbhM9f+T73nZPn8sOuHWQda8f5PbjwHsDTioyQu2d3fzv5sNcMi+D2RnJ4/IYMb4ovvi+Ar788A6e2FPN+iXZ4/I4IuNBnwSWsLVpVxV1J9v45Or8cX2cDUtzmJ2exI+fL9F7ARJSVAAStn72ShmzM5K4aG76uD6OL8r47EWz2X20mVdLj4/rY4mMJRWAhKW9lc3sqGjkE6tmERU1uuP+A3HNshzSk2O594XScX8skbGiApCw9JttFcT6orhq6XCfRxwb8TE+bjw/j+f317K/+sSEPKbIaKkAJOy0d3bzyBtHuWzhNCaf4Zw/Z+Ljq2aREOPTXoCEDBWAhJ1nio/R0NLBtYUzJvRxJyfF8pEVuWzccZRjzZojSIKfCkDCzu9eP8K0lDjeUzDxU4d8anUeHV2OB7dUDL+wiMdUABJWmk538JcDtXzw7On4JuDN3/5mTU3iooJ0Htx6mM6u7gl/fJGRUAFIWHlq7zE6uhwfOGe6Zxk+vmoWVU2tmh5Cgp4KQMLKYzsryUlL4JwZqZ5leO+CTLJS4vlfnTBGgpwKQMJGU0sHLx6s4wNnZ496zv/RiPZFcd3KXF44UEv58VOe5RAZjgpAwsaTe6vp7HZccbb38/Fct2ImvijjV9oLkCCmApCw8djOKnKnJLAkx7vhn15ZqfG8d0Emv3v9CB16M1iClApAwkLDqXZeLqnjiiXTPR3+6euvCnOpO9nO8/t1KlMJTioACQtPFR/rGf4JoumYL56fQXpyLL/dps8ESHBSAUhYeLa4huzUeBbnpHgd5S0x/rmInimu4fjJNq/jiLyLCkBCXltnFy8erGXNgsygGf7pdW3hDDq7HRt3DHeKbJGJpwKQkPdaaT2n2rt471mZXkd5lwVZKSzJSeU3RUe8jiLyLioACXnP7qshPiaKC+aM74lfztS1y2ewt6qZPZVNXkcReQcVgIQ05xzP7DvG6jnpxMf4vI4zoCvPmU6Mz3h0u4aBJLioACSkHaw5SUX9adYE4fBPr8lJsVw8L4ON2yvp7tY5gyV4qAAkpD1T3DPh2poFwVsAAFcuzaG6uZUtZfVeRxF5iwpAQtqz+46xaHoK2akJXkcZ0vvOyiQx1qdhIAkqKgAJWU2nO9hW3hD0r/4BEmOjuWzhNB7fXUV7p6aGkOCgApCQ9eqbdXQ7uMiDM3+diQ1Lp9PY0sGLBzU1hAQHFYCErBcO1pEU6+PcmWleRwnIRQUZTE6M0TCQBA0VgISslw7Wcf6cqcT4QuO/cYwvivVLsnlq7zFOtXV6HUdEBSChqfz4KQ7Xt4TM8E+vDUtzON3RxdPFx7yOIqICkND04sE6AC4sCM5P/w6mcNZkpqfGaxhIgkJABWBma81sv5mVmNntA9xvZnan//6dZrZsuHXNbKmZbTaz7WZWZGYrx+YpSSR46WAd01PjmZ2e5HWUEYmKMj5wznRePFhL0+kOr+NIhBu2AMzMB9wFrAMWAteb2cJ+i60DCvxfNwF3B7DuHcC3nXNLgW/6r4sMq6vb8cqbdVxUkBF0s38GYv2SbDq6HE/t1TCQeCuQPYCVQIlzrtQ51w48CGzot8wG4H7XYzOQZmbZw6zrgN7J21MB7RNLQHYeaaS5tTPkhn96nTMjlZy0BDbtqvI6ikS4QAogB+h7SqMj/tsCWWaodb8I/LuZVQDfB74+0IOb2U3+IaKi2lodPy094/9msHpuaBaAmbF+SZaGgcRzgRTAQPvY/We0GmyZodb9PPAl51wu8CXgfwZ6cOfcvc65QudcYUZGaB3xIePj1TePszA7hSlJsV5HOWO9w0BPaxhIPBRIARwBcvtcn8G7h2sGW2aodW8Efu+//Bt6hotEhtTa0cXrhxs4f/ZUr6OMytLcNA0DiecCKYCtQIGZ5ZtZLHAdsLHfMhuBG/xHA60CmpxzVcOsWwlc7L+8Bjg4yuciEWBHRSNtnd2sCvECMDPWLc7ixYN1NLdqGEi8MWwBOOc6gduAJ4Fi4GHn3B4zu9nMbvYvtgkoBUqA+4BbhlrXv87ngB+Y2Q7gu/QcPSQypM2l9ZjBivwpXkcZtfVnZ9Pe1a1hIPFMdCALOec20fNHvu9t9/S57IBbA13Xf/tLwPKRhBXZXHqcRdNTSE2I8TrKqJ2bm8b01Hg27arimmUzvI4jEUifBJaQ0drRxbbDDazKD+3hn15mxrol2bxwQMNA4g0VgISM7RWNtIfB+H9f65f0DAM9o7mBxAMqAAkZm0uPh834f69zc9PITo3nsZ3VXkeRCKQCkJARTuP/vaKijHWLs3nhYC0nNAwkE0wFICGh5/j/xrAZ/+/rirOzaO/s5tl9NV5HkQijApCQEI7j/73OzZ3MtJQ4Ht+lYSCZWCoACQmbS48TFWbj/716h4Ge21+jM4XJhFIBSEjoGf9PDavx/77WLs6irbOb5/drwkOZOCoACXpvjf/PDr9X/71W5E0hPTmWTbs1N5BMHBWABL1dR5to7+xmRV74FoAvyrh8URbP7auhtaPL6zgSIVQAEvS2ltUDUBjGBQCwbnE2Le1d/OWAhoFkYqgAJOhtK2tgTkZSSM//H4jzZk9hcmIMj2uKaJkgKgAJat3djqLyBgpnhferf4AYXxTvX5jFM8U1tHVqGEjGnwpAgtqbtSdpOt1BYd5kr6NMiHVLsjjR1snLJXVeR5EIoAKQoLa1rAEI//H/XhfMSWdSfDSb9KEwmQAqAAlqReX1pCfHkjc10esoEyI2OorLFk7jz3uqae/s9jqOhDkVgAS1orIGls+ajJl5HWXCrF+cTXNrJ6+WHvc6ioQ5FYAErZrmVg7Xt4T18f8DubAgneS4aB0NJONOBSBBq6i8Z/x/+azIeAO4V3yMjzULMvnz3mN0dmkYSMaPCkCCVlFZA/ExUSyanup1lAm3fkkW9afa2XKo3usoEsZUABK0isrrOWdGGrHRkfff9OJ5mSTE+DQ3kIyryPvNkpDQ0t7JnsrmiBv/75UQ2zMM9MTuY3R1O6/jSJhSAUhQ2n64ka5ux/II+QDYQNYuzqLuZBvb/O+FiIw1FYAEpaLyBsxg2czILYBLF2QSFx3FJh0NJONEBSBBaWtZPfOnTQrbE8AEIjkumovnZfDE7mq6NQwk40AFIEGnq9vxxuHGiJn/Zyjrl2RT3dzKGxWNXkeRMKQCkKCzr7qZk22dETED6HDWnJVJjM94QkcDyThQAUjQ6X3TU3sAkBIfw0UFGWzaVY1zGgaSsRVQAZjZWjPbb2YlZnb7APebmd3pv3+nmS0LZF0z+4L/vj1mdsfon46Eg61lDWSlxJOTluB1lKCwbnEWRxtPs+tok9dRJMwMWwBm5gPuAtYBC4HrzWxhv8XWAQX+r5uAu4db18wuBTYAZzvnFgHfH4snJKFvW1k9hXmRNQHcUC5bOI3oKNMU0TLmAtkDWAmUOOdKnXPtwIP0/OHuawNwv+uxGUgzs+xh1v088D3nXBuAc65mDJ6PhLijjaepbGqlMMLm/xlKWmIs58+ZyuO7qzQMJGMqkALIASr6XD/ivy2QZYZadx5wkZm9ZmZ/MbMVIwku4akoQk4AP1Lrl2RTfryF4qoTXkeRMBJIAQy0H97/Zchgywy1bjQwGVgFfBV42AbY5zezm8ysyMyKamtrA4groayorIGkWB8LsiZ5HSWovH/hNKIMHtfRQDKGAimAI0Bun+szgMoAlxlq3SPA7/3DRluAbiC9/4M75+51zhU65wozMjICiCuhbGtZPctmTSbapwPU+pqaHMd5+VN5bJeGgWTsBPJbthUoMLN8M4sFrgM29ltmI3CD/2igVUCTc65qmHUfAdYAmNk8IBbQmbAjWHNrB/uPnYi4+f8DtX5JFqW1pzhYc9LrKBImhi0A51wncBvwJFAMPOyc22NmN5vZzf7FNgGlQAlwH3DLUOv61/kpMNvMdtPz5vCNTi9tItrr5Q04R8TOADqcyxdlYQaP62ggGSPRgSzknNtEzx/5vrfd0+eyA24NdF3/7e3Ax0cSVsLbtvIGfFHG0tw0r6MEpcyUeFbMmsLju6v42/cVeB1HwoAGWiVobC2rZ2F2CklxAb0uiUhrF2exr/oEb9ZqGEhGTwUgQaGjq5vtFY0a/x/G2sVZADyxW8NAMnoqAAkKeyqbae3o1vj/MKanJXDuzDQdDipjQgUgQeHtD4BpD2A46xdns/toM4ePt3gdRUKcCkCCQlFZA7lTEpiWEu91lKDXOwykvQAZLRWAeM45R1F5g+b/D1DulESW5KSySe8DyCipAMRz5cdbqDvZpuGfEVi3JIsdFY0cbTztdRQJYSoA8VxR7wlgtAcQsHWLswEdDSSjowIQzxWV1ZMSH01BZrLXUUJGfnoSC7Im8fguvQ8gZ04FIJ4rKm9g+azJREXpBDAjsX5JNkXlDVQ3tXodRUKUCkA8VX+qnZKak6zI1/DPSF1xds8w0J929p+cVyQwKgDxVO/x//oA2MjNyUhmSU4qj25XAciZUQGIp4rKG4j1RbEkJ9XrKCFpw9Lp7DrapLmB5IyoAMRTW8vqOXtGKvExPq+jhKQPnjMdM7QXIGdEBSCeOd3exe6jTRr/H4VpKfGcP3sqj24/qjOFyYipAMQz2ysa6ehyrNAHwEblqqU5lB9vYceRJq+jSIhRAYhnet8AXj5TewCjcfniLGJ9UTzyxlGvo0iIUQGIZ7aWNzB/2iRSE2O8jhLSUhNiWLMgkz/trKKzq9vrOBJCVADiia5ux+vlDazI1/DPWLjq3OnUnWzjlTePex1FQogKQDyxr7qZk22dOv5/jFwyP5NJ8dE8sl3DQBI4FYB4Yuuh3hPAqADGQnyMj3WLs3hydzWtHV1ex5EQoQIQT2wtb2B6ajw5aQleRwkbG5bmcKq9i6eLj3kdRUKECkAmnHOOorJ6vfofY6tmT2VaShx/eF3DQBIYFYBMuCMNpznW3KYPgI0xX5Rx9bkzeP5ALTUnNEOoDE8FIBNu61sTwOkIoLF27fIZdHU7fSZAAqICkAm3tayeSfHRzMuc5HWUsDM3M5lzZ6bx221HNDWEDEsFIBNua1kDhToBzLj58PJcDhw7yU5NDSHDUAHIhNIJYMbfB87JJi46it9uO+J1FAlyKgCZUFsO6QQw4y0lPoa1i7N4dPtRfSZAhhRQAZjZWjPbb2YlZnb7APebmd3pv3+nmS0bwbpfMTNnZumjeyoSCjaXHic+JopzZqR5HSWsXbt8Bs2tnfpMgAxp2AIwMx9wF7AOWAhcb2YL+y22Dijwf90E3B3IumaWC1wGHB71M5GQ8NqhegpnTSE2Wjuf4+mCOelMT43XMJAMKZDfwpVAiXOu1DnXDjwIbOi3zAbgftdjM5BmZtkBrPtD4O8BHa4QARpb2tlX3cyq2Rr+GW++KOOaZTN44UAtlY2nvY4jQSqQAsgBKvpcP+K/LZBlBl3XzK4Ejjrndgz14GZ2k5kVmVlRbW1tAHElWL12qB7nej6xKuPvIytyccBDWyuGXVYiUyAFMNCxev1fsQ+2zIC3m1ki8A3gm8M9uHPuXudcoXOuMCMjY9iwErx6x//P1vj/hMidksh7CjJ4aGuFzhMgAwqkAI4AuX2uzwD6n4F6sGUGu30OkA/sMLMy/+2vm1nWSMJLaNlcqvH/ifbR82ZS3dzKs/tqvI4iQSiQ38StQIGZ5ZtZLHAdsLHfMhuBG/xHA60CmpxzVYOt65zb5ZzLdM7lOefy6CmKZc656rF6YhJcesf/z9Px/xPqvQsymZYSx6+26DgLebdhC8A51wncBjwJFAMPO+f2mNnNZnazf7FNQClQAtwH3DLUumP+LCTovTX+P0fj/xMp2hfFRwpz+cuBWirqW7yOI0EmOpCFnHOb6Pkj3/e2e/pcdsCtga47wDJ5geSQ0PVaab1//D/V6ygR5yMrZ/Kj50p4aGsFX7l8vtdxJIhoMFYmxObS4yyfNZm4aJ/XUSJOTloCl87P5KGiCjr0ZrD0oQKQcdfY0k5xdTOr8jX845WPnjeT2hNtPL1XnwyWt6kAZNxt0fi/5y6Zn0lOWgK/eLXM6ygSRFQAMu42a/zfc74o48YLZrG5tJ49lZomWnqoAGTcvfJmHYWzpmj832MfKZxJQoyPn79c5nUUCRIqABlXNSda2Vd9ggsLNNmr11ITY7h2+Qwe3V5J3ck2r+NIEFAByLh6uaQOgAvnqgCCwSdX59He1c2vXtMHw0QFIOPspYPHmZwYw8LsFK+jCDAnI5lL5mfwy83ltHXqZDGRTgUg48Y5x0sltVwwN13n/w0in16dT+2JNh7bWeV1FPGYCkDGzZu1JznW3MZFGv4JKhcVpDM3M5mfvnyIng/xS6RSAci4efFgz/j/ahVAUDEzPr06n91Hm3n1zeNexxEPqQBk3Lx0sI68qYnkTkn0Oor0c82yHDInxXHX8yVeRxEPqQBkXHR0dbO59Lhe/Qep+Bgfn7toNi+XHGd7RaPXccQjKgAZF28cbuRUe5cO/wxi1583k9SEGH78nPYCIpUKQMbF8/tr8EUZq/UBsKCVHBfNJy/I4897j3Hg2Amv44gHVAAyLp7bX0vhrMmkxMd4HUWG8MkL8kiM9XHP8296HUU8oAKQMVfVdJriqmYuXZDpdRQZxuSkWD66ciaP7qjUGcMikApAxtzz+2sBWKMCCAmfvWg2PjN+rCOCIo4KQMbcc/tqyElLoCAz2esoEoCs1HiuX5nLb4qOUH78lNdxZAKpAGRMtXV28VJJHZfMz8BM0z+EilsvnUu0z/ivZw56HUUmkApAxtTWQw20tHdp+CfEZKbEc+P5eTzyxlFKanREUKRQAciYenZfDbHRUZyv0z+GnL++eA4JMT5++JT2AiKFCkDGjHOO5/bXcF7+FBJjo72OIyM0JSmWz1yYz2O7qth1RKeNjAQqABkzB2tOcqjuFO9flOV1FDlDn3vPbKYmxfJ/H9urmUIjgApAxsyTu6sxg8sXTvM6ipyhSfExfPGyebx2qJ6ni2u8jiPjTAUgY+aJPdUsmzmZzJR4r6PIKFy/Ipc5GUn866ZiOrq6vY4j40gFIGOior6FPZXNrNXwT8iL9kXxD+vPorTulM4dHOZUADImntxTDcDlKoCwsGZBJhfMmcp/PHWA+lPtXseRcaICkDHxxO5qzspOYeZUnfwlHJgZ37pyEafaOrnjiX1ex5FxElABmNlaM9tvZiVmdvsA95uZ3em/f6eZLRtuXTP7dzPb51/+D2aWNibPSCZcTXMr2w43aPgnzMybNolPX5jPg1sreP1wg9dxZBwMWwBm5gPuAtYBC4HrzWxhv8XWAQX+r5uAuwNY9ylgsXPubOAA8PVRPxvxxJN7j+EcrF2sAgg3f/PeArJS4vnmo7vp6tZhoeEmkD2AlUCJc67UOdcOPAhs6LfMBuB+12MzkGZm2UOt65z7s3Ou07/+ZmDGGDwf8cAft1cyJyOJedM0+Vu4SY6L5h8/cBa7jzbzi1fKvI4jYyyQAsgBKvpcP+K/LZBlAlkX4NPA4wM9uJndZGZFZlZUW1sbQFyZSBX1LWwpq+fqc3M0+VuYumJJNpfOz+Dfn9zP4eM6Z0A4CaQABvqt7r8vONgyw65rZt8AOoEHBnpw59y9zrlC51xhRkZGAHFlIm3cUQnAhqUD9bqEAzPju9csITrK+NrvduoTwmEkkAI4AuT2uT4DqAxwmSHXNbMbgQ8AH3P6XxVynHP84Y2jrMibTO4UHf0TzrJTE/iHK87i1dLj/HpLxfArSEgIpAC2AgVmlm9mscB1wMZ+y2wEbvAfDbQKaHLOVQ21rpmtBb4GXOmc035lCNpT2UxJzUmuOlev/iPBdStyWT13Kt95bK+GgsLEsAXgf6P2NuBJoBh42Dm3x8xuNrOb/YttAkqBEuA+4Jah1vWv8yNgEvCUmW03s3vG7mnJRHjkjaPE+IwrlmR7HUUmgJlxx7XnEBVl/O1Db2iaiDBgoTTyUlhY6IqKiryOIUBXt+P8f32Gc3LTuO+GQq/jyAT6445KvvDrN/ibNXP58vvnex1HAmBm25xz7/pF1SeB5Yz85UANNSfauEbDPxHng+dM50PLZvCj50rYcqje6zgyCioAOSO/eq2C9OQ43qepnyPStzcsYuaURL7w69epOdHqdRw5QyoAGbGqptM8u+8YHy6cQYxP/4UiUXJcNHd/fDlNpzu47Vd6PyBU6bdXRuzhrUfodnD9ipleRxEPnZWdwveuOZsth+r5t8c1YVwo0olbZUQ6urp5cOthLipI18yfwlXn5rC9opGfvHSIBdkpXLtcM7qEEu0ByIg8vruaqqZWPnlBntdRJEh844qzWD13Kl///U5eebPO6zgyAioACZhzjv956RCz05O4dH6m13EkSMT4ovjxx5aTNzWJm3+5jZKaE15HkgCpACRgrx9uYEdFI59anUdUlCZ+k7elJsTws0+tIDbax40/3Upl42mvI0kAVAASsHv+UkpqQgwf0jivDGDG5ER+/qkVNJ/u4GM/eU2Hh4YAFYAEpLiqmaf2HuNTq/NIjNWxAzKwxTmp/PzTK6huauUTP9lCg84nHNRUABKQHz1bQnJcNJ+6IN/rKBLkls+awk9uLOTQ8VNcd+9m7QkEMRWADOvAsRNs2l3FJy/IIzUxxus4EgJWz03nZ59cweH6Fj7y/zZzVO8JBCUVgAzrjif2kRwbzWcu1Kt/Cdzquen872dXUneyjQ/f/Qr7q3V0ULBRAciQNpce5+niGj5/6RwmJ8V6HUdCzPJZU3jwplV0djuuvfsVXjig07oGExWADKq72/HdTcVkp8bz6dV69S9nZtH0VB65dTU5kxP41M+3cv+rZTqtZJBQAcigHiqqYOeRJr56+XziY3xex5EQNj0tgd9+/gIumZfBNx/dw5ce2s6ptk6vY0U8FYAMqPZEG/+6qZjz8qdwteb8lzGQHBfNfTcU8pX3z2Pjjkquuutl9lY2ex0roqkAZED//Ke9tHZ0852rl2CmT/3K2IiKMm5bU8AvP3MeDS0dbLjrJf77mYOaTtojKgB5l0e3H+WPOyq5bc1c5mYmex1HwtDquek89aX3sHZxNj946gDX/FhHCXlBBSDvUFHfwj/+YTfLZ03mlkvmeB1HwtjkpFj++/pz+fHHlnG08TTr73yRb23cQ1NLh9fRIoYKQN5yur2LWx54HYD//MhSonW2L5kA65dk8/SXL+a6Fbnc/2oZl3z/OX65uVzDQhNAv+EC9Ez1/NXf7mB3ZRM//MhScqfoZC8ycaYkxfKdq5fwxy9cyLxpk/g/j+xmzQ+e58Eth2nvVBGMFxWA4FzP8f5/2lnF19Yu0InexTOLpqfy4E2r+MkNhUxOjOX23+/i0u8/z/2vlnFSh42OOQulD2QUFha6oqIir2OEFecc//HUAf772RJuPH8W37pykY76kaDgnOP5A7Xc+cxB3jjcSHJcNB9alsMnzs/TwQkjZGbbnHOF77pdBRC5Oru6+fYf9/LLzeVctyKX7169RCd6kaDjnGN7RSO/fLWcP+2sor2rm6W5aVy1dDofOGc66clxXkcMeioAeYf6U+18+eHtPL+/lr++eDZfu3yB/vhL0Ks72cbvth3hke2VFFc144syVs2ewpoF01izIJP89CSvIwYlFYC85fn9NXztdztpONXBP125kI+dN8vrSCIjduDYCR554yhP7T3GwZqTAOSnJ3HJ/AzOy5/C8llTyJikvQNQAQhQWnuSO57YzxN7qpmTkcSd15/LoumpXscSGbWK+hae21/DM8U1bC49Tpv/yKH89CQKZ03m7BmpnJWdwoLsFJLjIu+MdiqACNXd7dhcepxfb63gsZ2VxEZH8YU1BXz2onziojXBm4Sf9s5udh1toqisnq1lDWwrr6ehz4fLZk1N5KysFOZNS2bm1CTypiYya2oS6cmxYXsAxKgKwMzWAv8F+ICfOOe+1+9+89+/HmgBPumce32odc1sCvAQkAeUAX/lnGsYKocKIDCn2jrZcqiel0vqeHx3NUcbTzMpPprrV87kcxfN1m6xRBTnHFVNreytbKa4qpni6maKq05QfvwU3X3+/CXF+sidkkhmSjxZKXFMS4l/6ytzUhyTE2NJTYxhUlx0yL1fdsYFYGY+4ABwGXAE2Apc75zb22eZ9cAX6CmA84D/cs6dN9S6ZnYHUO+c+56Z3Q5Mds59bagskV4AzjlaO7o51d5JS1sXza0dHGtu5VhzG9XNrbxZc5L9x05wqO4UXd2O2Ogozp89lQ8tn8H7F07TlM4ifbR3dnO08TRlx09RXneK8voWKupPU3OileqmVupOtr2jIHpFGaQlxpKWEENqYgxpCTGkJMSQGBtNYqyPpFgfCbHRJMX5SIjxkRQXTUKsj7joKOKio4jxvf319nUj1n851hc15gUzWAEEMhi2EihxzpX6f9CDwAZgb59lNgD3u5422WxmaWaWTc+r+8HW3QBc4l//F8DzwJAFcKbufOYgG3dUvnUSirf+Td07vr3rfvfW/e6d1/v9pwh4vX7390kywDLv/JmdXY5T7Z3veuxeZpA7OZH5WZNYtziL8/KnUpg3WX/0RQYRGx1FfnpSz5FD8999f2dXN3Un2znW3ErNiTaaTnfQ2NJOY0sHjad7vjed7qDuZDuldadoae+ipa2Tlo6uQX9PAxVlEB0VRVQU+MzwRRk//thyLixIH90P7ieQAsgBKvpcP0LPq/zhlskZZt1pzrkqAOdclZllDvTgZnYTcBPAzJkzA4j7bpmT4pg/bZL/B77j21tjfm9fH/r+t9e3QZYf5P5+P6D/ekOuC/iiokiO85EYF/3WK4zkuGimpcSRlRpPRnKc5u4RGUPRviiyUuPJSo0f0Xq9e+ot7Z20tHdxqr2TU21dtHd209HV89Xe2U17VzcdXe6t6x1dPbe1d3bT2eXoco6u7re/slLHfug2kAIYaF+kf78Ntkwg6w7JOXcvcC/0DAGNZN1e162cyXUrz6w8RERGwsxIiPWREOtjqtdhhhHIS8YjQG6f6zOAygCXGWrdY/5hIvzfawKPLSIioxVIAWwFCsws38xigeuAjf2W2QjcYD1WAU3+4Z2h1t0I3Oi/fCPw6Cifi4iIjMCwQ0DOuU4zuw14kp5DOX/qnNtjZjf7778H2ETPEUAl9BwG+qmh1vX/6O8BD5vZZ4DDwIfH9JmJiMiQ9EEwEZEwN9hhoDpsREQkQqkAREQilApARCRCqQBERCJUSL0JbGa1QPkZrp4O1I1hnPEWSnlDKSso73gKpawQWnlHk3WWcy6j/40hVQCjYWZFA70LHqxCKW8oZQXlHU+hlBVCK+94ZNUQkIhIhFIBiIhEqEgqgHu9DjBCoZQ3lLKC8o6nUMoKoZV3zLNGzHsAIiLyTpG0ByAiIn2oAEREIlREFICZrTWz/WZW4j//cFAxszIz22Vm282syH/bFDN7yswO+r9P9jDfT82sxsx297lt0Hxm9nX/tt5vZpcHSd5vmdlR/zbe7j+Pted5zSzXzJ4zs2Iz22Nmf+u/Pei27xBZg3XbxpvZFjPb4c/7bf/tQbdth8k7ftvXORfWX/RMQ/0mMBuIBXYAC73O1S9jGZDe77Y7gNv9l28H/s3DfO8BlgG7h8sHLPRv4zgg37/tfUGQ91vAVwZY1tO8QDawzH95EnDAnynotu8QWYN12xqQ7L8cA7wGrArGbTtM3nHbvpGwB/DWSe2dc+1A74npg90G4Bf+y78ArvIqiHPuBaC+382D5dsAPOica3POHaLnHBErJyJnr0HyDsbTvM65Kufc6/7LJ4Bies6lHXTbd4isg/F62zrn3En/1Rj/lyMIt+0weQcz6ryRUACDnbA+mDjgz2a2zcxu8t82zfWcVQ3/90zP0g1ssHzBvL1vM7Od/iGi3t3+oMlrZnnAufS88gvq7dsvKwTptjUzn5ltp+eUs08554J62w6SF8Zp+0ZCAYz6xPQTYLVzbhmwDrjVzN7jdaBRCNbtfTcwB1gKVAE/8N8eFHnNLBn4HfBF51zzUIsOcNuE5h0ga9BuW+dcl3NuKT3nI19pZouHWDxY847b9o2EAgjkpPaecs5V+r/XAH+gZzfumJllA/i/13iXcECD5QvK7e2cO+b/5eoG7uPtXWXP85pZDD1/UB9wzv3ef3NQbt+Bsgbztu3lnGsEngfWEqTbtq++ecdz+0ZCAQRyUnvPmFmSmU3qvQy8H9hNT8Yb/YvdCDzqTcJBDZZvI3CdmcWZWT5QAGzxIN879P7C+11NzzYGj/OamQH/AxQ75/6jz11Bt30HyxrE2zbDzNL8lxOA9wH7CMJtO1Tecd2+E/UOt5df9Jyw/gA975J/w+s8/bLNpued/B3Ant58wFTgGeCg//sUDzP+mp5dzw56XnV8Zqh8wDf823o/sC5I8v4S2AXs9P/iZAdDXuBCenbbdwLb/V/rg3H7DpE1WLft2cAb/ly7gW/6bw+6bTtM3nHbvpoKQkQkQkXCEJCIiAxABSAiEqFUACIiEUoFICISoVQAIiIRSgUgIhKhVAAiIhFKBSAyCmaW558f/z7/HO5/9n+KUyToqQBERq8AuMs5twhoBD7kbRyRwKgAREbvkHNuu//yNiDPuygigVMBiIxeW5/LXUC0V0FERkIFICISoVQAIiIRSrOBiohEKO0BiIhEKBWAiEiEUgGIiEQoFYCISIRSAYiIRCgVgIhIhFIBiIhEqP8PCB5ZjWUGb/UAAAAASUVORK5CYII=\n",
"text/plain": [
"