Created
May 2, 2021 07:58
-
-
Save peilingjiang/b34bcbadcb352139c75b198913ecde05 to your computer and use it in GitHub Desktop.
Bayesian Modeling Problem 10.6
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"id": "retired-substitute", | |
"metadata": {}, | |
"source": [ | |
"# Problem 10.6" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"id": "chemical-syracuse", | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"import numpy as np\n", | |
"import scipy.stats as stats\n", | |
"from math import log\n", | |
"import seaborn as sns\n", | |
"import matplotlib.pyplot as plt" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"id": "essential-capture", | |
"metadata": {}, | |
"source": [ | |
"## Question 1" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"id": "egyptian-container", | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"p_same = .5\n", | |
"sig1 = 3\n", | |
"sig2 = 10\n", | |
"sigs = 10\n", | |
"\n", | |
"j1 = 1 / sig1**2\n", | |
"j2 = 1 / sig2**2\n", | |
"js = 1 / sigs**2" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"id": "caring-webster", | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"d = lambda x1, x2: log(p_same / (1 - p_same)) + .5 * log(1 + j1 * j2 / (js * (\n", | |
" j1 + j2 + js))) - .5 * j1 * j2 / (j1 + j2 + js) * (j1 * x1**2 /\n", | |
" (j1 + js) + j2 * x2**2 /\n", | |
" (j2 + js) - 2 * x1 * x2)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"id": "arctic-grace", | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"x1 = np.linspace(-30, 30, 500)\n", | |
"x2 = np.linspace(30, -30, 500)\n", | |
"\n", | |
"d_matrix = np.vectorize(d)(x1[:, np.newaxis], x2)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 5, | |
"id": "decreased-microwave", | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAD5CAYAAABPqQIFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABhgklEQVR4nO29d5glR3ku/n59zswmrXa12pxXCWUhES7OBAEO2Lo2wiYbfM0asDFcgg0IYytg4wvYj7EBsxh+XNsCjBFBQiLKYPDFMiAjUM5abc5hdjbMnNP1+6O7+lRVV+yu0zOz6u955pnTlbpntPPq/d5662tijKGNNtpoY6ZGMtUP0EYbbbRRJ1oQa6ONNmZ0tCDWRhttzOhoQayNNtqY0dGCWBtttDGjowWxNtpoY0ZHt4mbPO1ZH2OUEGikg2QkQTKSAJ0EnVkdUIdA3fx73k8JIRnpoNtNQB1Ct5ug2+0g6RA63QTdboKkQ+h2O+gUnxN08jHdboIkIXQ7CTr596wNGOkk6Hay/iSh7DofnxBhpMv7gW7elhBhJCHhetDXFdo7hKJt8L2HBH0QpSCk6NBk8TmhSSSUIqEeCH10qAdCCqJ8XP45oUkQSwGWAmk//86ve/lnBqQ9MLEvbxtc52393mCNXg9IWfZZ/er1gTQdrJmmg7G9fF1Auwbrs6ytz8D498k+0M8sPWwyHyteF2PTYg5Slq0FZJ8n+9nnfCzLv6f8Or83Sxl6kynefcdjuPqitSCiop0Ja6b5eixlSKX2tGiXxqes+DGz68GPXlyr/co1ADCmWcfwOf/Rrf0vZfdT8B+mEp+iJ3n7rWLcL1Y0AmLJSAIkSQFgNJKDlgbAMhDLgKUzkoNRt4NuDm7dboLuSIIkya9HBoCkAthIN8lAJ6ECwArwkvpkwEsSFODVTQidHLhEAOsmidCWARhv4yBGSLUAxkHLB8AI/QGAFUCmATCWlgFMuuagJqyh+0tTwEoLYHxeMdcNYOiLc+S/FZZ6/O2kwhoWABPjM4/txovWLS4ATL0XH6/en1+r7WnpWv5RlNuXrov1p6k3M5mheVkjjy0BlPi52ym3OQAs6RCSJMkYWT6m0+2gkzO3bs60OpwNdRN0u1QAmApufA4fzxkY/+oUwJVIrCsRWJcIYBzUChDyADA+pgAt8TMZQKgAKFYGN97GGZcJwHR/hSq1KPoMbE0BNROAlZhVmipjLCxMbbcEB7aDxyax+cgJXHLaKYM+Dk6aNUQWJq2nAy4DOPF+n/AdZwoOODGBJ0n8v6ZTNMLEOrM6QGfAxDgro4SQzO6CkoyVjYx2teljIjCuTg5qHIA6I53iM08HOVh1uxkIdbuDFFJKL2mQUmZzZJDqKOmiLp0ckfo4iA1SSP45oZ42TUwgp5alFDJVUsi0JwOYLoXUAZiYQopgxdNCDdti/V6pzcTYbAwsBoCJYUsjeWx6YAd+96zlBTipACamnaa1i7Y+K7EwMXxBzTRO154k9YEuNLqNoEH8aOSxefrI2ZdO/0pyICv0KQXAut1OkT7yMZx9lfWsDJw4gHGwMgFYtyuzq4SQa2CJ1CaC1UgBaHoA04GWmiZyAMtASwNguhQyJoDZAInZxlnSSxUABPZkArBsnjndFNfRAZg67r7941g02sXpo93yGhikkcXtDOmjLc2NwcaqxrAALomochHRGgD/CGAZAAZgE2Psb+LdYRCNaWJc/+Kivpo++gr4BUDlwFfWs0hiYDylFPWwgrHR4JqnjCMJGYBLBSsR+IAOMRD1JQBL1NRQSRONfRKAmVNILYCZ5tYBMMAJYAVr4sGFfB4WQJA0KAnwZHYmCuzy/AGrYozhnx7ehXdfvLZYh4vz0py+/VoXugzb1DeMMIFXrPQucprYA/AWxth/E9F8ALcT0TcYY/dEvQsaZGK69DHTsgapYTcHJpGBdbodgY3loDcySBm5SF987gyArdsVtS5Z/1J3IFX9Syfgc/1LBTZ111H8XCmFVHchm04hAQv7KoOgVciPkUaqgGMR829+fB8uX3EaOiAjUIkiv6qF6VJJMVyamDhu2BGbkcUEMcbYDgA78s9jRHQvgFUAooNYIxJdsSM50inSRw5g3ZxxdUcytsUF+9FZ3QLARmd1JAATNS4TgPH0kTOx0ZGObLnokBHABsK+nE6KAMa/dKAVCmDZZ4uNoskUEnCkkX4AVlUHy+6vpKMo62Bi8PajJ3r4/t4x/MKyBdJ83x1JCcyEfm6r0EVVzWs6RoiwT0QbieiHwtdG07pEtB7ApQD+axjP3TwTM/i/CnE+/6zaJ8T0URXwOYBJWpeij6kCvsyuTClj2A6kCFoqSGXppbw7KfZ5AZgphfTdhfRKIcupIoC4AAYUAKZGASRKGlkapxHzP/7QTvzOmctK7EpcV2w37UiWnkUTplTSBWqh9grOthKyZuRRIoSJMcY2AdjkGkdEpwC4AcCbGGOHKz+cJRrTxEz2CVHA5/qXTsB37UBycCulkAXIcUCDlDKKAMb1LdXv1VFALSGUQEtnYg0CsLQ3AKUChGQfGFQA0xlZfQGMhwvAxHZxjspaJHASdS3NX7JurggoBnCzifmbDx8DS4G182bDLNTrWZi6lhjiruR0YFRqChkzpYy9O0lEI8gA7HrG2Ofjrj6IhiwW3RyMBsDE00fVPqHTv8wMS3bg2/xffAdS9X+p+pYIaOIOpAhq6g6k+LkwroqfBdFe7EvQk0HJVwPTsTVXCtkTtS4PBsbnlMCMXw/ASqeBFaGCm0sHy0O3GymGapH4+AM78ccXrkHJUlFcy3/pqXIvFZSrCP2xxvrEMHYoY2piREQAPg7gXsbYX8VbuRyNgJhP+mjyf6np48AaoTAtYSeSj9MJ+BmzKgv45R3HsoXCtAMpAlhCA8BSdyBFACNypIWmFNIEYDYG1tOL9VYAk9iazKjYpPBZSAtLgrjIfEwApoZlN9Ik5n9r2wE8ffF8zCLZI6DbkVSfsWQH0USJmGqGxQCUqsAUC3xI+f3VjJ8B8AoAdxLRHXnbOxljt8S8CdAUE1OOBunSxwLUEgIlAzAz6V+6HUjdmUndESLVda+ClWxgHYzjAMZ1MPGzmCZykEqQmsFNBCALgJVSyFAA06WFPgCm9uVrMIF9adNIlw6mRuBuZNGWX5+Y6OOb2w/gmkvWZb8erehvZmFqqIK+K3SgEwpENutEmppTyGm+O/kfABo5X9kYExPtEwPBvnx8SGefGAAUSg78QqwXAa/Qw0z+rwFgqWxsRAEuvgNJ1LcaV00pZNngKqSQaV8BJA4cfZmBcfYVkkKq7Cs/zF0JwPi8GjuRAGTGZrBTuHYjRTGf9Rk++fAu/PaZyyQAUw9t87Fi2PrEiAUSIYBjBLUhi/vT7TiRbzQDYgb9S0wfVf+XmD5mZtXBLqMIYKpjX9W/RMG+bFyVtS5VwB8AmFvAVy0URgAzmVj7grBf2oUMBDCVfU0BgKlRMrTqdhANaaQ0XwCqHUdP4MhkH2fNnwObmK/6wtQ1+WfxHjoWJmXYgeAWW7/iMZ2F/aaimXTScXxIJ+Dr0kd9dYoBgMnpo/kIkS6d1FWh8Cmjo+5AOk2sNne9FsAqppCaLz8fWCCAAVp64K2DWfrUNFLdXdz04E685bxVKBlWA6wTPgL+dAjVaiGC1zR17DcWjYDY6Cwu0pvTR32FCfMBbpN9QmRa5t1Fu/7lK+CLVSg6MPdJACamkBoLxYxLISEAgSaN1AKYRQfz3Y387vaDePLCeZgj/OWZjK28T17PwMgkH5n8PX/0UujGqeHLmFS9yzZvOmtiTUZDmpgg3juOD4n614gBrHwrUKgAZvJ/lQEsLZ+BVAT8IADTecA0An4tAGswhQQgAVJxX4sOVoRDBxO/m3YjT0z08ZVt+3GtIOabdhbFPpe51RZ1wYLPJyLJ8FpFpB+WNtaCmCVE932I/8uVPor6l5o+lnUv2W3fET67AMzGwEo7kMLupBXA1BSy2J1kZgCrkEICqARg2TwRjDTMxSHke+tg4vf+AMCkccL1/314F377jEzM166l2ZHUhe2Y0XQKG9DFZGMtiNluUsP/JYr7pgoU8s6irsqELNarR4h0An7IDqQ9hRQATDWxigzMBmD9vBqEN4ANmJlUkYJHCeBYqV37B24qG20S8n10MGG+DuBUYX7b2HGM9/o4+9Q5WmZl0rtsgr5u/GBeqckrfXSFD/hYU0mBjUXTxKZNwemwaEYTG+1k3q+OrG+Z0sfBjmQ5fRQPcIsMS/V/6XYhTf4vXRFDtQqrj4BfAFgqsCrpXKMMYCzVOPYlwBPYl4+RtW4KCUjO+cLUatuJBGQw8tTBxDmiEG87WgRkxQ7fev5qCcBcWlixcyluEhhYGI+YWhMQBlo2xiUJ+pHTynZ30hJq7S/X8SFZrC+zL/EAtwpg4gFu8xEiwcCq1MEXP0spYmFw9ajEatqBVBhYFAArpYyOFFLsUz5L9b9ElmUBMHncAMCkEAFMDNXAqtmNFIHn37YdwNNOn4953Q50looYu47WCq5TkGW6UsiYQNamk5bQue/V9FGtNmFKH1UB31Rh1cXAjEUMPRiYuxKrRgNLU5RMrMYU0gPAdCmhKYU0ApiSYvI2ixdMCp2Qn4d5TjnFNO1GinFsoo9bc2e+CjQ+O5I+9fOnKny1LpsJNs5zzMx8spl0Mmdivu577Y6kpnyOKtjr2JiuBn5IEcPaO5A2E6vrHCTvK1JGJhc0BOqlkECxtreZFYBOyHfqZmoaCUhpJIDSbmRRlbXP8PEHd+J/nbW8EPNtxQtNWph2rCa9HMyVv8cKvkPpZFm6FFJoi/1s1GlBzBiq+16XPnY4eFmOD+nOP5rYV/kIkR7AfIsYVgKwOiZWFZAAMxiFABhfR/ju4wUrQgUwMQKEfJcOJoLKo0eOA4xh/SmzoTO2Dj7bWZhJ0B/8aEz69TQZLs1MB2Q+80KCWiZmDpP7Xlc+R9p11Bwf0tX/UrUu10s8fBmYCmBqFQrtDqTLxBoMYJpUULkuOfHFsfwzD6Xf6gUT5pSsFDxUHUwn5IuhSTuNZyNThrSX4uMP7MC7LloL4w5jgBY2FQ79UItEKDBFq2LRgpg5snSybF71dd+r9gl151GXPoo7lD4HuLX+L8PxIq2Ab6lC4e0BSwWh3vZm7lAGBgzSTr42HAwMAPhbsHUAptuJ5KFaH1SHf1p25aumVt5309b9eN7K0zCrkwSzMBPoqTEsYKsCRra00fQ52vOOzExlv5GnHumQ1n3PS+n4HB/ijnuxvr0NwDj4+fq/VABLqFe2UKhHiIYBYFwDGxaA2VJIQAYpeACYOsdTBxNd+WKIfYePT+KOfWP4uaULglmYaKtQo/TcmlDBocoegA1g6tTuEplX1EKGHfL+mk7RfDqpSR917nuX/0t3/tEl4JsYmE3AN+5AmsroiJVYAT8ACzWxugAMMKSSlhQSkFNC1UrBQyPkF+1iqFYLzQ5kKjIzzV/8h+/fgd87ZwW0Pi9JU7PTEZt4r0YJvGownSrpYxXtK94B8OkFTr7RjGPfUH1VLB+dJGUAM5lTfSwWth1Ik0lVPUJkPMSt6leqidVpocjb+evRvER5B4AB+u/KeiKAAZAArGSL4OsY2ZrF0AoMxlrsFOpuJI8f7RnDqjmjWDp7FD7HgHTivU8MI5X0TSPFc5ShQDaM8tStJmaJWaMdrXiv2id82Jd4bT4+ZK9AoQr0phpgRgMrT/906WPJ7yW48oEyg4vJwNQUEjACmMqkjOcbTV4wkw5mADhdGmnaQZycTPHZx/bgWsETVkULC0klY5yX9K04IX6uAmTZ88Y3p063NNE3mkknDexLd3g71H3vSh9dNfB1byEq+n0MrNF2ICOlkJr00Qxgwh+0hxdsMFY59+gh5IuhqxEmjvvkQzvx8g1LkRAhZfrjRaG7kFNVN6wOY7Ixr2EwsWSkE3fBhqIRYV99eYdqn+CCPRfwXf4v+eiQXf/iANahSa2Arx4v4i+yLVVh5QK+DwNL+4NxTQBYCcwcKSSg9YJpAUwIJ4CpOhgggY7JlS/2PX74OMYm+zhvwVy4dhd1aWh5zNSBl61N/CyK/DZ2pfZFZ2IJeX95rUf0i0R0PxE9RERvj/u0g2jMYmFiXyHue1f6aDvArQUvp/+LA1IqA1Kat4foX/y6rzmwDZTPQQrpYNWXehgZ2ES2tvZokMrAAP2ZyNI8DcDloUsjdQCU9lJsemAH3pG/ek3t1x7eVtiey6EvjlUjBmNy9fumlfx5VCuFCFyxU8qYmhgRdQB8CMBzAWwF8AMiupExdk+0m+TRzAFwDftSvV/+h7dl/5dX+gjP9FGtQKGCUSHEM38AE8vo9DUHtnWfqwAYnyus55tCAvA2s+rnaQAOAx2Mh00HA4Abt+zD5csXYk7uCSvGaeeUUcN2plJ73TcDX1NhAjK1b9i1xIDomtjTATzEGHsEAIjoMwCuABAdxJpJJ/M3dZdTw3oApjsypAr43O9lSh+9HPgSo8pSRW8GVmJwQ0whATeAGQ5glwBMtVLAAGBC6IT8ot0RB45N4s6D4/j5ZQNPmGqpUMO0flP6V0jK6OvtUv1j6jxdShnTYuH75RGrAGwRrrfmbdGjIYuF3joRcvZRdd9ry0c7XuBhP/8oAI+lBpjEwPgBbkCZE1H/ArQ7jdoU0uTEB7QHugdjU41eZjjUnYd2J5L3KZ91fjCmgN2H7tuO15+zUl7fEF5VWyuwq1givM8Y044jMAAyXXopXottMSIknSSijQA2Ck2bGGOb4j2NfzQDYh47i7JYbwMwu/vet/6X9wFuAcCYeCzIVgc/JoApDMuqgaEegFmtFCYvmOZgN1/bBGxq23d3HcL5C+bgtNEu1GKH0neLZcJkqzCtpYtYgBBiUhX71Gtdeqk+Z0xNLOTYUQ5YNtDaBmCNcL06b4sezaST+c4jf3WaeHRIvR6Mk48Xcf1L8nGhrz0+NNhhFB33FdLHfl+fKjYBYCmzp5C+Ir5HCgmUAay0EwkMgIuHRgcr1uZTJGF+ALR83fETPXx9+wFcsfp0rTalApjPjuQwoo7mbUslXde640nDKl4Y+djRDwCcTUQbiGgUwIsB3DiM525G2Fe8YLrvvPa9jpW5qk+E+r+kt3BbGVg6SB8Bt/4VA8AAB/tK5f78u6+NQnecyNcLBohML5XnoqyDpQoLKuYLnz98/3a89uzlJdbhkw7WsU/UAUFjccJIDEzHyABoWRkQj0FSRHRkjPWI6A8AfA1AB8AnGGN3R7uBEM1YLDrmF9eq4KWaV9XqE6bdR53+pStuKKWPQM64LAe4VQZWmqOAFyDbJdTD12KBQyMDU1NKQ3/eFzOFNFkpJC8YMBjr+VktscP6DLfvG8Py2aNYOXcWdGJ+sYaDhampZJm9DY+1VXHp82tA7tddi20iKzMBWp2IfeyIMXYLgFuiLqqJhhz7fuwrpPaXj/8rOH1UAYwDFSDPAczsCyizpQLQZg6A6T1kfH23DuZy5Z/op7hh815ce8m6QZ+nnytmuOwVdQR7XV92T/P6JiuFOq9OFQxTJO2xI8tNDGJ+uXSOXPs+pHx0Lf+XCGBpvsOnAh4Qnj7yOYY0MSaAFeHSwPh1AIDZhHxp/TxMOpgYH31gB15z1nIkRJW0MJ9zkrEjEV7KETNdNLE0wK89VrQHwG03UVz2vuwrpHih1X2f5qmgCEa+B7hN+pe2AoUGjEz6F6BlZqUvTb98zKdckTUKA8vDdyeSz7XpYLztJ3uPYH43KcpN68ZUjSaPGYWcbfRNJ8U23/ZYMVOLIjbIxKaAfUlGVQv7MgGYaKkQNTE+hn92gpEmfeTj+Rp1AQz8MT0BTBqvAbBUdtu7diJ9dDAAODHRx6ce2Y3rnjxII8vz9CxMHWdLBW12jNCQwMryirRQILO1AeZ2XV+UaNNJc4hnH1VNTCfYu9iYt/ve0/8lAZjK2ETBXgUvLzAyAVhqmaMHsCJ0ABZiowDMx4ny9XWHutWdyOI+Dh1MfJ6/f2AHXnO2nEb6poLDqf2lT2fjrW9nYKY2sd3VFyvadNISszpJJfZl2n3k7AuA3X0vHt4G4P0Cj1j2CaAMYLo5wjipX5wvpXU1UkhpfnUAK2lTFh2Mj/nx/iyN3GBJI321MDFipaLeFglFGwPKaWEV4FKByZU+RmdjnTadNIaJfZWMq3XZlyjeA/7sywRgavoYqn8B2Q4nX+MkATAeKoCZjhUBwPGJPj7z6B5tGjkTIqRShWl8SOoo9vn0R4mWiVluUvJ+6QHMpG91qGc8WmStvmoAMJYKwKLzf/HP/F9ORQNrEaa+EoiV+40ABv5oQ0ohxfWVz756kzj+w/dvx++do08jfbUwMXxTwDrszJeNmcbH1LuaYGLUCvvmEAHM9uZtl3lV3Z0E4K6+atW68n+FOv9X8dnClhwA5kwhpXVZqZ0J4JM9jicDA+T+fA3jWUbdZ4+dSABaHUwFuf/acxjLZ49i3bzZg/sH6Fs+qWQTUQIpi8gvzgHiMjDduCjRCvvm6BBDhyb86n6FsC8AlcV7II7+BZSqTFhf5KGbrwE4p4mVX0sMKhIDc3jBQnSwseM93LRlH67JTa2mml8hZyRtvjLTHFuI4n6wVUIBshhWCVf/sKIV9i1Rl31x7xcAeNW+1wGYmCaWGJZQ7DDU/8XXyL97Axgw+NdvAbAiphDAbEK+TQdjfYYP3rsNbzh3JciSRk63qAJkgDu9FNsBS6rahIivi1bYN0cV6wQvnQPA4v1SAEzydTF4vf+xGCuMA/zTRz4WlvSxtJ7wLzKAgRXhA2DStQBgxX31KaR0D2WODsDUENu+um0/LjltHpbOHi2N082xaWGxU0nqUPRD4CadDHCnjaYxprHDiJaJWaIO+yqsFN6VJwzsS0whq6SPgGBqlfu0NcBMjM0AcFYB31SJgq/pC2AOBqaOM52JlPoMad/OIydw+74jeOeFa/LH1LMw62aBZyrpG5SQkf2pfrGqIr2IA6Y0s+g3rDdl0Wpi5ghlX3bnvYV9iS/vAAbj+edinFDFwlf/4mOBcACTQIzJ/3ojAJjM2NwANri3H4AVwx0AVgBdL8UH792GdxgALDSmS9rpU4JHajekmeI8NRpJGw3R7k5aoop1IkrpHB/zKgC4vFy69FGnf0lrWABMx8C0etZwAEw9D6mO0+1EpgoA2XSw/++hXXjh2sWY1+1omZONhfmUno4VakrpYmOmNt4O2MEMMAOaus6URMvEzFGZfZnSR5v21RdZlnL2MXT3EdD2eR3gVudK1wIA6ADMWpAwLdpVAX+wRlwA46EyIl0p6LsOjGOin+LSRadox9pYlVSyJ5J1wpZC6kIHZEA5teThSgl1uplPhJDPaFLWlOez1aKZdLIu+/J9cYcP+wKAVGVeFvYFhKWPPNR1hTbrIW4XgBl2IAdtegAL1cDE69Kr1sSfQRg33uvjeuFwt2030hdYhp1K6gR+3XnKGDuNISRzKjT2yK9saywaSifN4j0w2H2sJN6Hpo8lEHKkj2Lq55M+Sm2eAMbDh4FB+cNuAMDEMFWUYH2Gv75nK9547kokAQX7fAV9ABKwVGVqocxMDJ+dRhtQ+RAdF9ANlSy1u5PmqC3eDwPATExKx8wAP/+XdM1K7VoBP/8clEIWcwwAxsMl4ufhAjCXkM/6DDds3oOfXrIAy+aM5j+uvKYYemtGAE2JHCY2BsjAOejzY17Z/LBnaTWx8GhGE0NPa51wal/DYl+u9zgKfVL5HN280joagAOs+pd0DSjsKbWbWIGhMjAXgAHAg4ePYev4BH5jzeL8kc3AVcdZ7xsuD5iOjZnmmMDMF6xmksxEI52pfoRK0Vg6Wck6AeRAMDj3yPjuZHZRkX3pwKcMQCX9SzdHamOlccb00QVgPuwLqAVgPOoA2PiJHj7+4A5ce8n6/Mcug5TLqGoS9GPpYVrQMgCZ+gw8bMws6y+3DYtcDg0YG2JiRPQ+AL8KYALAwwBezRg7WHW9Rv4/UTr3KIr3aX8AUv1e3pZrY/1BH+PmVc7KxDlc+E/7dgBL2eDtQ7b0Uad/GdNO/fxgAMvvoQJYEbo/7jwFDfKBCSEeJyqmGBiM1njaS/GBu7fijeeuQlfRU2JWVgXM4FEnTA51m8CdJFQAmiuSZDhfwwpKyPurZnwDwIWMsYsBPADgHXUWa9QnVtT8AgaMSXyjkM24KrIvtZSOD/syppAO9qUDL2kdVl4XlvQRcPu/gDD9a8gppBji3E8/ugfPWr4Qyy06mAu4hmGr0IVJ0De2C0BmY2a6iAG4vkAZNRpiYoyxrwuXtwG4ss56ze1ODkP7AhBUuJCPM/RVEu81/a7zj2bxPdWkix4ApoJHhBTSNhYAfrTvCI70+viZJafmj+4Cq3Ka6Tu+Sug0LhuQ2e5pSzN1MUwAGqoNIuC5iWgjgI1C0ybG2KYKd/0dAP9SYV4RzexONrHzCLid94a+gu35AF9xXQY5NX0soqaAP2irD2C2FFL1gonrin17xifw+cf34uqL1+aPrgcolxY2FWGzWLjsFzoAif3zhIJUzEPbIcJ+DlhG0CKibwJYrum6ijH2pXzMVQB6AK4Pe1I5GgGxSjW/fMALiM++SuMD2RcQbmAV5/AxngJ+sQ6/FipR+Drx1WublWIyTfH+u7fgqgvXll7gagIwEzDETCVDdhxtzEsFBRcjHCYz8gGoqJUnIgpujLHLbf1E9CoALwDwHCa+zrxCNARioh4laGB9odoqF+7zz0E7j8XaGvalY08R2Rcgp4/Z48cHMJ3+VazDrw1amQ3AfM9E8vire7Zh4zkrcEr+f21ngUIjS0vRRNjsFj7G11BQC3q2QAAaeqmchvwgRPSLAP4IwC8wxo7WXa8hEAtkXwBkkV8U9ZV1RKAxifdAGcBs4KSuIc4X+nTpY2X/F1APwPJwAZj0n8VhpRCD9Rk+8+huPO30+diQl5mO4aDXhQ9QBGleDiALum8DgrsXAxsGA2xuM+HvAMwC8I2czd/GGHtt1cWaAbGpTh+FNuvuIw8LAPJ+l/s+e1QLgBW/GncKWfTzZ4GSvuZjfQDM1wumzv/P3YdxtJfiWcsW5D+CiWG500g1mtDJfIoghgJa3fAFRBtgTdd00haMsbNirhcMYkT0XAC/CeBDjLE7iGijc1eizvseAThLRgN69iUwMyP7sqaQgewL8E4fy/MMDnzdOvz+prEoA5gKXtIYB4A9euQ4vrXzIN5+wWpUjaYEfave5bBN6NZpInxYVRVPW3DMpOMFQlRhYr8D4HUA3kVEiwA82TlD0L6Y6sTPGt3aVzFHw75suhgQ7v1Sx+R9ld33fG1d+ph/n64AdnCih03378DVlwyE/FAWVirZE0kPq1NiOgTQYoUv4NgAdKgWi24ziVnsqPLUY/kRgbcS0XsBPM05w8W+ADuAWQT6YPFeutalkPoxTaaPzh1IdSxfD9UATAwRfE70U7z3zsfxjovWYiT/v3TMNHKY4fJ+SWMjWSeqAIyL8ZnWHIoX7WRnYkT0NwDeBOBm3sYYezsRvcE1l6WaEjtZRz3wAvSFCYGydcKYMhrAK/9cm33lY4PZF78nVAA1jIUbwGw2CrW/30vx53c+jtefsxKndss7kb6u/FgsTFfjC3AbUavqXLEYT12R3gRWQ0l3T3YQAzAG4EYALwYAIno+gHczxn7GOdNX+9IBmEv7Mgn3RVs5LdSCo6avCvuSfmbdPPE7ECV9FNtCAEwMNf3763u24cq1i7Fm3qz8x/EDgFAW1qT5dVh2iWCrRCTQip5anuz1xBhj7yKilwL4NhFNADgC4O1ek2NpX6a+vN2r6uqwxXu+jqd9orQ2/30Hpo/iHF8A03nBWJ/h4w/uwNMXz8eFC+dBF02zMJ8IPRoENGSX8ACausBFsRjUyc7EiOg5AF4DYBzACgC/wxi732uyz4Ft3l8HvIpxnuxL7Ud19mU1rwrfs7ke6aNw7at/iW1VAeyzj+3B0tmj+Nn8TGT26/FLI2MyG91arnI4QPOCvS8bsmlYITuPKmC1u5Nh6eRVAP6EMfYfRHQRgH8hojczxv7NOTMVgcpxYBsoMyhTHzTaF1+7SfYlgFcxj1/b2Jcwx8v/pblHKIAZxfg+w81b9yFlDC9Ytaho9zW0lgBtiABi0sfUqCvY1wEIl/DuC1w6llUeE0m/65zku5OMsWcLn+8kol8CcAOAn3ZP1rAv3l5H+5LGGeaqY8Tr/HMt9sVDC3weAOZgX1I4diCB6gB2644D2HVsEq8+c9lgLQuATfUOpC+QqTEMi4LPTmEdtqUb1wr7g6gMvYyxHXmK6TFYkz6a0jvepnvbdn5tZV/i55jpI48q6SMQdnyogoAPVAewf995EA8dPo6NZw+KDtgAYlhppBg+5xp90svY4WttCPV6VQGu6EB2sgv7umCMHfMayM2uvrYJgy7mPDJkAi9NXyztqzQuIH2U2E0E/Usap9G/pHXz9u/sOoS7Dx7F685ZMXgUS2roA1jas5cVRH3fNxOZgMUH3Or4rap6vHwAS7e+eD2cs5NPMCYWFKG2CUACKG/T6jC0Lz4WzbAvef5wAezbOw/i3kP+AKaLYWth9V6xFvcPvY7nKwZwqXOiG15bELNEaLFCoEgnvdmXbg2d9sXXGIJ1YrCemX3J68VJH6VxngD29e0HsPnIcbz2bH8AG0Ya6Xswu2kNLs7h7PjABQzAizGGHx88im/tPIi/8HpaR5zsx46I6C8ZY3/satOGzlVvSh3zNiN4SXOqpY5ABfDK+4ZhndDN8fF/lcZ5AtiNW/bh4EQPv3uWWQOLCWAx/GExjapVtCRX+uYLWLr7+4JW1gncffAovr79AI6mKZ686BS8/vyV1mfzjicAE3suABWwfknTVo66qaOVddVLHUttLu2LX1exTgjXIexLbHMBmA28AOD6R3ZhVpLgFRuWFn2hAGaKOoexQ+cO26haBbRs82oBF4CHjhzHLVv3Y6zXwwWnzcPvnbcSp4x0IpfiOUmFfSJ6HYDXAziDiH4idM0H8D2vu4TuOqr9IcI9X6sigOnYV3luZAATYlj6F2MMH75/B845dQ4uX76w6AvVwEzPPdPDrxxOXODSzRXBa+uxCdy0ZR/2nZjE2afOwauftBwLRrultaKJ/CcxE/sUgK8A+AvIx4zGGGP7ve6SpnHAS/wcCl6AdedxysDLQ7wH4O3AL43pZ3Xx/89dW/H8lafhKYtOGazpAWC+aWRdQX/YTvvwF3CY/6B9QUvXZgMtSgh7jk/ipi37sO3oCayZPxu/eeYSLJk9agStaEeOsoeJt1aD4QQxxtghAIeI6NUAfgPAej6PiMAYu8Z5l55YatqROkrX/qkjUNa+gk2rJuFe+G7deczXDLFOiG2+4n322Q/AxiZ7+POfbMHvnrUMZ8yfM1i3AoA1FVXOQapzg+c5/oCHBVwAMN5Pccu2/bj/0FEsnjOCX117OtbkJcB9wCtaSnmygpgQXwRwCMDtAE4E3UXDrqZ015GPRXPal26OTbwHwtJH3XqbjxzHh+/fjj86fzUWzRoZrFsRwKqyMEqSaj6xob5JqBpoZXPdwKWuoQLXJGP41s6D+MG+McztdPDLaxfht85cWlrLB7iSaOnkSb47CWA1Y+wXK92lDvsqPtdkX7wtf57pkj6K83wYmDpHN4b1s5r4t+48gGufvB4jwqvVYgPYTIqmgQuQwYsxhtsPjOPWHQfAADxr+UJc9eR1SPL/PiHglUj9M5OJEdFbALwfwBLG2N6q64SA2PeI6CLG2J3Bd+n1qp11BLTjymDhAV661FE3V/g+E9kXYwz/9MhuMAa884I10rshpzKFrMrGqt7La9yQQYvPeezIcdy4ZR8OTfTw1CXz8eaL1mBWJ7GmnfxnMDGuoQj7FGkdr1vRGgDPA/B43bVCQOxnAbyaiB5Blk4SAMYYu9g1sfJZR824oew6CtcyqJl1LHlNkWn5a19im4t56UIFsKO9Pt5391Y8b8VpeMbi+dLYOgBme56w+l3xgSxU2K7y5iCXxgWUwetwr4+btuzDI0eOY9282Xjl2ctwWp7S10kXeTsNg4lRo0zsr5G9e/JLdRcKAbFfqnyXCOwLcKSOQDzTar6WfF2dfYltVbxfpnli30OHj+HjD+7Em85fhaWC/pX92MMBsCrB/1CrnaWs9kcWA7hM64jg1WcM39l9GN/bcwindDv4tXWL8YpTB4bimOAlPV80JtbYy3OvALCNMfZj9S3yVSIExB4H8DIAZzDGriGitQCWA9jsnNnrT33qGLDrqB03DcCrNC7v+9xje7Dt6ASuuWQdusI/bu3buSMCWB0rRFRrgLiu06QaD7T4vIfHjuFLj+/D0X4fz1xxWqFz+aSL6r0TDUjZgCuq2TWgnhgRbQSwUWjaJL66kYi+iQwf1LgKwDuRpZJRIgTEPgwgBfBsANcgq7l/A3zediQClwm81HEoMy/AsetoYV7FXH6tAogDvOR1DYCnuZfL96W2h7Cvscke/uqebfiFZQvwwrWLpf66ADZTos47G019PuA13k9x05Z9eGDsGM48ZTY2nrcC80fcRtQQkd4bvDqR/ocQwMRywDK+b5Yxdrn2FllB1Q0AOAtbDeC/iejpjLGdQc+bRwiI/Q/G2GVE9KP8IQ8Q0ajXTCOIKeCVf65qWAUaYl+2OfAHMBN42ebyvtv2HMbNW/fjD89dicUV0kfdur59Ux3DAC7duipwMcbwo4Pj+Nq2AxjtEK5YtxgvPXtQRDI2eEntKniJwBVNExu+sJ9vDBZn3ojoMQBPbWp3cpKIOgBYfvMlyJiZOySza3zRHqgAXrY0sAJ46dI8tT0GeB3vp/jwfduxcu4orr54rXX3UV1Pam9IA6sbITtvw2BclBAOTvRww+a92Hr0BJ6yeD7edvEajHbKu4dDBy5gAF4agb92NCvsR4sQEPsggC8AWEpE7wFwJYB3ec10mFqnve5lmwM9gNVJHU1r3r5vDDds3ovXnr2ieJVacY8IADZdomngAjS1uQj4/t4xfHPHQZzS7eCFZyzWuujVNUNEenGud8qoglfUA+DNgxhjbH3dNUJq7F9PRLcDeA4ye8X/ZIzd6zU5wKwKOHQvwN8y4dK9lDWL5xGu6wj3QBzta7zXx0fu344Vc0Zx3SXrnOxLXdO0dmj/sCL8XGO4eO8LXkf6fdyweS82j5/A006fj7dfslbaLImxwyjOrQ1eEZkYPQGYGBhj9wG4L/guLvY1g3UvcW5M9sX7v7n9AL6z6xBe96QVWD5bliBjsq8mAaz6+cb6Ohegr4h63+Fj+NKWvegmhBeuW4IzFsyR+qczeMVz7J/kx46I6KnItkfX5fMI8DO7ipqYj9cLKINXHeZVGpevWTyPcN0081LH8f5tR0/go/fvwM8sPRVXX7JO6p8u4GU7qB1Dq6nCuEz31gFXH8BXtu3HD/eN4bwFc/GmC1djTrejXXtYwCW1VwGuWDuTQCPC/jAiBHqvB/A2AHfCV9DnESraA82ZVYW2UOYltvuAl2mu2Hein+ITD+3ERJ/hjy9Yjbn5H1XxqNMEwMSIfVB7mOBFCeHAiUl85rE92Ht8Er+4ahGuvmx9kaIPA7yk9lCx3rE7Gbco4smfTu5hjN1Y5SZDOaRdzNeDV9GmSROnI/tKeym+uv0Avrf7MF595jJsOGW23B8AXrr1Q/ubDq+XcNRkXg8fOY7Pbd6DWZ0EL16/BCvny7/j0J1GINycKrUr4CXd1xe82t3JIBD7UyL6BwC3QijFwxj7vHPmpCZFrFMaWtsegXnp1hXmmoCrND5w1/HOA+P4zKO78dwVC3FNjdSx9LNU6G8yqjIuwB+4QMB/7jmMr28/gDPmz8EbL5TZbR3gEtttrEtuiwhcsV37TwAQezWAcwGMYJBOMgBOEKsk2PPx6lhf8NKu62ZeIX6v0vhA3WvL+Al88qGdOHP+HFx9yTp0Ku466tavOqaJaAK8UgC3bNuPH+wdw08vPRV/etn6ouSNeo8YrKv03B7mVCd45fOc4DXDzk7GjhAQexpj7EmV7mKyS/iAlzS/DGC17RK6dYW5oeAl3Ve9T9639/gkPvnwTszpdPC/z1uFeTV0L909QvubiibAa4IxfG7zHjxw+Bh+ZfXpuPapG6z3iF5FIsZOY0d5Dhd4xdKyAs5OTqcIrSd2PmPsnuC7TEV1CaHN6g8T2ofFvHjf4YkePvnwLkymDK8+Y5lUbRWIb1idDuAVW+sC9GL90V4fn350D7YfPYEr1y3By8/WV49Q11UPotc5zyi1DyttVIArrmM/4loNRgiIPQPAHUT0KELriYXYJSqCFzA9BXt+UPsfH96NI70+XrlhKZbPGZ7fK3TcsKIp8Brv9fGph3Zjz/FJvOSspdJ7BJraaZTbKqSNvppXh8rgFXN38gmQTlYrTQ2U2VdIYUIgPG3M52jXhp4lDQO8Dk308E+P7MKRyT5esWEpVs6tflRId4+qY4YZdcAL8E8bj6Uprs/B62VnLMX6CuZUdVwUwR4IThu9Na8kGapj/4kAYtXriansK/aOow28dOvCzb5KoBSQOu49Pol/fiRLG19ek3npnqXqmGFHE5rXJGP49KN78Pj4cbz8zGVWZ7267lA0Lx/wssyTUsR8nEnzog7ZNwFqBgsAsemUeDZST4xN9rMPETQvwJw21gWu0vgKu42ffnQ3RhPCyzYslUrktLuNmvaAtDFlDF94fC9+cmAcL92wFK960vLSGNO6tXcbp0LzAozMq/QMkdhYykbcgxD1llGimXpiU2CVENuGDV53HRjHF7fsxZJZo3jt2cuL4njAExu8bGN8mRclhH/feRDf2H4A/3PtYlx5xtJSv2nNqT7XaLRK8DGO3UYj81LAKxYbYywgnZyhIFa5npjOYV+0T9Pzjbb5QOaw/9bOQ/j3XQdx3oK5eNv5qzGq6C0n424j0Bx4PXD4KP7p4d34qSWn4tpL1yPpml8YOy3ONgrzfNiXFbz4urr1huQTYzj5NbHq9cQE8Bq23jVM4GJ9hmO9Pj7/+F48cPgYnrl0Af70oviFCX37pyJipo4m8Do40cOmB3bg9FldvPvSdUUBQnX96MxrKnYbAatgrzKv8lyKZrNIQ5jYNAovEKPsr/Q7yN7+HV5PrO/xGrRpvtO47egJfPaxPTjRT/HraxbjpesHaU3slHE6ghcwXADjutenH9mNx44cx8ZzVmDJvNHSGN2ajR/MHoZgz9d1MS8deEUyuzJ03IOmYXiBGGOMEdEtjLGLUKGemJg2TsXZxjop4/f2HMatOw5i+ZxRvOqMZVg4Ohy9qy5w2XY3AbOI7hvWncUIwv2dB8bx6Ud340XrluBlQt16df0Zf77Rwb50or3sExN+Fv6zxwKxk5mJ5fHfRPQ0xtgPQm8i7k4CqGxQrWqRCE0ZxyZ7RXXPZ5w+H++8cE1xrnEq00UXUPnMbaKKquk+JrPqR+7djiWzRnHdUzagk2hARbOeq3692F4VvCqVxLHtONbRvEzgFbF8Tj+sRuq0iaDdSQAvz99OMg4EFEUcgr9rGMeC7jwwjpu37suqe65ZjPVCOZxY4NUkcNWNKgAWon19a+dBfGvHQbz+3JVYccosqc+05smme1VmXokAuC0T847nV77LkDxeMfSuIyd6+NLj+/Dg2DFcsGAu3njuKswW/lHHeH9jCHBNJWj5Rl32dXiih7+9bzsuOm0ernnKemljpA6AxfB7TYlR1Vf3kkBM014zmtTEiOgNAH4fQB/AzYyxP6q6VgiI/bah/RrXxKnQu1zM68f7j+Ar2/ajQ4RfW306XrJ+yWD9AK0rBnhNR+AK9R75Ath3dh3CrTsO4A3nrhq6cF/baV8ndQTK7Mv0DHXYV0wQa4iJEdGzAFwB4BLG2AkiWuqaY4sQEBsXPs8G8AIAnm87cuxMohngOjjRwxce34vHx0/gogVz8aZzV2FW/g+pKcY1pemhhx4WQ8DXgdfxfoq/vWtrVjtNKAetrjstjgnF8nrVZV78uwpeKiuLFA36xF4H4L2MsRMAwBjbXWexkFe2fUC8JqL3A/ia1+QGWZcKJP1eiu/uOoTv7j6E+d0urlizCOvmmbWuUJ3LKd5PQ5ZVJarqX5QQ7jt0FJ98aCf+4NxVWHPqbKnPtF4Tu47DqOsVfceRp44qqA2BiYX4xIhoI4CNQtMmxtgmz+nnAPi53G96HMBbq2wY8qizHTEXwGqfgcMyp9rA65FDx/Clx/dibLKPn11yKt5+QfgO4xMNuIxMq6qAT8CnHtmNAxM9vOfSDeiOxDGtTnXqOGW6l5adCeNqRsjZyRywjKBFRN9EViBCjauQ4c4iZOW9ngbgs0R0BmOs0h9OyCvb7gTAb9IBsAQeehiAYOCqmi6OTfZw05b9eGjsGNbNnYVXnbEMC3JfV5oypIwZ11LX82m3rWUdn3qd1rKGWlLGOd6RStYBMB37OjLZxwfu2YrnrliIl61YJvWZ1hoaeOkAswnw0p1xjAFeQ9iZBOKmk4yxy019RPQ6AJ/PQev7RJQCWAxgT5V7hTCxFwifewB2McZ6PhNt4FWXcU1OpvjOroP43p7DmNvp4FdXLypE+jRlXumiDqSs2lcAaMUALNO6oUBmitgA9tDhY/iHB3fgzResxrJ5eutEqO8reuqYz1Md+l6WCVNhQjV1NIBXMdele5X6kgHICW2x3tzdoMXiiwCeBeBbRHQOgFEAe6suFqKJba56k9jglfZS3HXwKL66bT8mUoafW3oq3nHBmuJFEFWBy9ruCVzDAi3TveoCWWwA+9q2/bjz4Diuu3QDRjzSx1i2CRWMskVqpI4lfUzDvviaBvZl1L34OlVTR6GNKKIm1pzF4hMAPkFEdwGYAPDbVVNJICydfBGArzLGxojoXQAuA3AdY+y/XXNV4PItOKimi1vHT+DGLfuw98QkLlgwF68/ZwXmdjsZaDEU6WJdtuUDWE2ClS1iMjIeVQAMBHzk/u1YNnsUb7t4rXE9F4CFmlajg5c6bqrBy8DQCvCKqIk1xcQYYxMAXh5rvZB08k8YY/9KRD8L4HIA7wPwEWROfmuEaF0qcB04MYkvb92Px44cx8o5o7hy7WIsmZ0JkGq66DK2mtp0c8vz6oOW7h5RX/RgCNM9tGBVAcBOMIb3/uRxXLFmMS5bOt+41gCIIrOvutpX6K6jS7Tnc0N0L/GzZWdSAi/+PNF2J0/+Y0f5AUj8CrLt1JuJ6DqfiaHgNX6ih69vP4A7D4xjwUgHL1h9Ol6+IfPDVQEuU3ts4KqyKzllwBYJwPZN9PC+u7fgjeetwqr5evtEk+K9U/uaSuHeZlb10MdK7CuyV6zBdDJqhIDYNiL6KIDnAvhLIpoF+G1n+ADXZJri29sO4ra9hzErSfC8Fafh11YtAhE5gUtdy9pmAZoQ0BqmjaLqYW1d1FnDBmCUEB4ZO4aPPbATf3LxOsyf3ZX6dGv4AthUaF++nq9ibBW3fWz2FRnEnghnJ38T2RuP3s8YO0hEKwC8zWeiCbh6kylu23sY39l1CGDAzy9bUAj0acrAGMAY8wKtYQDWVPu9YoKZGKU0L5CBUUK4Y/8R3LRlH665dB1GRzpSn7pGFeNqDNOqNXXMx/myr1q7ji7w8mFeav8wQOxkr+zKGDtKRN8CcDYR/XzefNxrrrKz+MN9R/BvOw+i10/xjCWn4i3nrcJIkjgFeh+2ZXbc20ErBmDpns+0+9dE6MDP53lcAPbdXYfww31j+JNL12mPD0VLHyOxLy/bRPaDegn3gzbSPIedVelBTJ5bWCZMIKgDvwiRVt8gnNII2Z38XQBvRObSvwOZ2/Y/ATzbNZcxhh/vPYJv7DiA4/0Ul512Cv7wSSsxu5MUaSL/3hRwhRtUq/0HNs0LATfWZ0PTyWzmU6AMYF/dth9bxk/gzRet0a4x9PTRJN6XLBEeO4/8OkD78k4d+XNZ7RKW1LEYQxrQ0syNEDXrck5ZhKSTb0R2ROA2xtiziOhcAH/uM/E9P9qMCxfOw+vOXoHZwi89TcupYtU0sS5gNV0SWrzfMNiaDwsLBTB+jGvjeSu1a+gALDh9BMLF+1D2FdM2EVv3Uvtt4CX2RYjeDEWxEBA7zhg7TkQgolmMsfuI6Ek+E9950dqgHcW6wOXedZxe/7FYyioDma9HrC6Afe6xPUgBvOKc5VK7br7umYYi3gewr5jaV63UkT+7766jJe1EkgCUf0WIkz6dBLCViBYiOzLwDSI6AJ+3f0OfKjYJXKGgZaonFhq6iqamqANkagSXoPYAMAD4rTOXSu26+TYGZrV1RGBfpfVcqSNfz+X7UgGsTuooXHvtOurYF+/jABaJjZ30IMYY+/X845/lAv8CAF/1mttnTtDyBay6YBULoHxCvZcL1GxAVkcX89mN5KEC2Oc37wUD8JsOAAs9+xht93GY7IvP1bEvVVQ37Rr6CvdVwYsoGhObocVXgoR9AvAyAGcwxq4horUAngzg+665xQFwB3CFsiwXaDUJWD6RpiyInVWJ8sFq/zRSBbAvb92HiTTFS87SV6EYGoCp91D0smDX/TCE+1Dty5Q6qvNMzExNHTsd/otr08mAsR9G9sbvZwO4BsAYgBuQif3WMNXKz/r8gSsmaMX0gIUwJBeQmdiYXqiv94/XBmC37jiAvccn8eonrZDa1bk2AAuyT3TknyU4fVTHmMT70jgzgMU872gU7tV5pj4d+xIZXYSYZv/P946gtx0xxi4joh8BAGPsABGNuiYBfsAVClq+gNWEYVW9hwvUhsXIQliYDcBu23MY9x06hj+4YJV1rs9r0+TxDgCron+ZGFj2gw3Wqyveq8Ck7augfTmYGzpd/gsuA1hkYf+JsDs5SUQdICuMSERLkDEzZ9je+zgYo/8FusCq1rsYPf6jVRHbfZz2w04tqwAYANx9cBzf2XUIf3Rx2QfmW4HCyMB89a+q4n1V60TVnUfHGIl9dTt+c0XmlT2wHbxI/u9XJ54I6eQHAXwBwLK8NvaVAP7EZ2KIpmUDrbBihHH+g9Qxqw7XpOr/f19fIN56bAL/8ugeXC28Rq0qgGn9X+IcpWR06fcUyr5iWSdUoMr7zMK7foy3bQIoA5yaOmYPqwewdncyaHfyeiK6HcBz8qYrGGP3hdxsGKA1VZ4vX7OqDchisjFfsDSxsIOTPXzovm245ikbiuKSegCqCGAdKq9XM33U7T76VFutLN5rgUiTOqpr2OaqfSbtS00dkwHwxarsOkOzSTeIEdGNalP+/flEBMbYr7nWKIoiViqHMzViPQ8fcODPaLNH+K5V5RlsWphPGjnBGN571xa886K1xSvsbAzM5gGLoX/VSh+niH0NHbyyXzb/DzG4Z7s76cXEfgrAFgCfBvBfGICYd7jK6BR9HoDVdGUJ8X4uMHGBWaywpZKhAAYC/uInj+MN567CwjkjpXnqPYcGYENOH4fOvmz6FgB0u4OftbR2dzDXpn3pwCuiJnYyC/vLkdUQewmAlwK4GcCnGWN3+96kis+rLli5ALGOYA9UZ1Vq1E0pa9ULSwgfvm87rlizuHgfpG0XsooHzJuBacb6AFjtc48+4r10LY/x2nk0gaNO+5IAzQPARKZWM05aJsYY6yNz5n81L4T4EgDfJqKrGWN/53OTerXsh/OLrbszaUsRTYxsmEK/ej8fFnbTln1YM29WUVK6UQDTjfV134vsi6+rrhNinfBhX/w+vuzL1lcYVrl9ogL7ksbFArEoyziDiJ4M4O8BzEb25rTXM8acpnlTeAn7OXj9CjIAW4/BTqV3DPtQtm1+1fTOJz0cNjCpUed9kxzAKCH8ZP8RbB4/jjdcUH7/cXQA8xHwY6SPMdmXZUzQgW01nfQR7oUdSeoo8ziAJYp/LEL0m2Ni/wfA1YyxrxDRL+fXz6y6mI+w/48ALgRwS37ju0Jv4lsHvzRvyDYJHi6Q8xHuTYysrj5mA0htJQjon1NkYLuPT+BfN+/FtU9ZX5pTF8DkHUk7gFVx31dJH53GVSuIyWO0AGbSzcS+EN+XL/tKhL4I0aAkxgCcmn9eAGB7ncV8mNjLAYwjqyf2h0IlTwLAGGOnmibysGpf00BM9LZLDOmAtm+YWJgrjRQBbJIx/NXdW/GuiwdVWXU7keZnsDOwIqoAWB39S/c8SY3dxyrivQ7AVO0LiA9gM08TexOArxHR+wEkAH66zmI+mljt39B03HU0/cG6AM3GynyAzDXGl7m5ih6a+v/qzi147Xmrihd7mBiY2JcoDItKjKsMSs4dyCr6V0D6WMl5X4V9dbvltVXwAlA6sB0LvCgBoyTcMqCJkN1JItoIYKPQtIkxtkno/yayTUE1rkLmNf3fjLEbiOg3AXwc2WsgKwXVePGudzzwtPMZ0DxQ+YbTh2UAFq0VQVnLNUb3Fm3dOAlcPNJIlYVRQrhh8x4sHO3i8tWLpPEhh7mtKWTgESIf/5e3eVUZJ+lfPuwLKDnna1knTNoXX6sugCVdsHw9xhIk9HO1ceyGB3d7/4G+8Oylle9HRIcALGSMsbw6ziGfjM4Ujbwt0wZewy6X42NfcJlRQ7QtlWnFLHaohhFcNWnk3QfHsePYBK48Y6k018bofH1g2edIAMZDTJ19GFg+rt7u4+B7kHVCXc/XuEqcZXbldlFHS7rC+h2JfTGWgKGTgxhqR4Mvtd8O4BcAfBvAswE8WGexRkDMF6iq/BINMpHz3jpwc9kmSpaJyAAVysJ0c3X9Y70+rn9kN97z1A2l+6j3sAFYUAqpA7AAB76ahqrsS5pTh30pDMvqvK/Kvop7CZaKfNeROFDx9kS8TrJxSurIWALGOjmIJUhZEuUPucH6e68B8DdE1EX2xrSNjvHWaAjEpmZtG8Dx/2AhYBYDyIYh/usArPi5CHj/3Vvx1ovWFGci1WcZJoC5BHxv/Yuv4Wuf4M9kAjANyHlVnVDZV6h1oqb2xdlXyrL7pvl1jGhK2GeM/QeAp8RarxEQc4UvyLlYl21d09xQx/ww00NX6LQwm52CEsInH9qJF6w5HafPlo8UhQCYFNMMwLzsE1oQk+eV00cD8KnAOGzrhABgKRuR2BdPJWO99Hay11w+GTNmFBOrA3Y2QBNptAhovoysdG3ZgQy1Yni/zUiz5o/3H0EvZfgfy07N16oGYFU1MHFcEIBpUlCv2l8mjap0PeT0EcjTQL72QOcq0kfOivk4F/sS0se+AGYA0Gdx/oxP2mNHMSJmOunzN83vZ2Zf/msBFYDHZqEw7Ew67REWFqZb80i/j88+tgfXKTpYHQCrwsB89C9pjSr6F1/HBGASICklbAx9VvblW221CvsSBH8X+0oxSC9jxHR7J4VvTIt0MiR0gOgCK9OYNNWzsmG/zMM3dCzMlsqKOtgH7t6KN+c6mDpHXbcRANM9f6I3sBb3URkbXyNUwBfTR35f1fvF55Tmau4lifVyeRwjgEUQ71OWSQIZG+tGTydbELNE3d+NC1P8tC99vy+QNXVG0sXCijbLbuS/bt6L56w4DYtzHcy0DiVUqonvm0KazkFK48T2Ujro58Av1qirfwnt1rOP0lzP9FECtwGoSdYJteKEh3jP2ACsAJTYlwhwMaIFMUvUFe5tv1sV4NR76UBLbde3VXsrkRpVx9hYmE3Mf3T8OLYfPVG86NaVRurH2BmYy4BqTSHV9FGc73LgV0kfdbuPJvCzsS+xRI7IvqQD2hb2JQFaGPsygZeqjdWNVhOLEFXsEiLA6fBCx7RM7aaxscLm1Ff7bW1qH1+rB4ZND+wI0sGmLYCpAr74H0anU5nSRx2AKX1WAFPZV/YDyexLrTqhY18B1gmb9qUDsAgnAwG0u5PW8GFiLvDwShlzQNOxM18gk/tlNjaMlLIqC9O58j9y/w685pwVGEmS0rq2qhSlZ7HsQkqfHQBm1rMqAphulzD7JejTyfx7cOkcn/SxIeuECGA6gMtSTlk2qBptOmmJEM3KFLr0z7SmDsx8NbFhs7GQKLEjS3x/3xgWzxnBWQvnluYX15adSBVwJHCpsAsZZKGo4/+y6FiV00fRJgGAn1Uc3EcuZijV/OJ2Cl/2xezsq0gnDeAVk4m16aRnVHfYm8cZAYr5szLzfePvVtpSSZ/6+SoLG++n+PLW/UV9MB0LGzqAVfWA+Qr4thQQcBzerpE+ij4yMU0cEvvqo1uAEkNH6k9Z1pcqYJa2Fovhhw8TU8fZxtpYldpuYmU29tUUG7MxLBOL0rX/zV3b8IfnrQKRTuMq/yBOAOPREICpHrDs3p4AJv58RgamsrUAAIuVPuZVJ3Tivcqw1J3H4loAsJR1JTCLES2IeUadtNHWr2tXwSwEyORnHrAxURcTdyjrHEeqUl+MP8/Xtu3HpYtOwdJ5o6X1fIR8AGUAUwsaVgQw7yKG6g5k9gPKQCO1aYAo7zfW/nKlpVUATOf9qiHe27QxVcgXASzWW4oaLE8dNaZM2HexLNOYqiI94Aay0PV0NgVjv0aIL49PpLEuMX//ZA/f3zuGd1+2XvM8bgAjFVyAoQFYbQFfC2SDfq8Xd9jYFwB0Onrwyn6BQ2NfqnHVlD6q7CtlDL00jVZWutdrQSwofIAtJG3k7XXSwmGlkT5mVd8zkiIIfvDebXjrRWu064lt1p1INARgPGwAJqWGFQDMxLKaTB8tNb9UgPKxTojivcq+UsaQsniCfJtOWmIYWhdvDwWyEDYWO3w2CFwsjMeNW/fh2StPw6mjXc3cwQ/kEvKNAMZDsVuoACaNswGYuguZj7GaWI0gpgAY7+P3bRLAeNHCQvTXA1iJbXlYJ9Jc2Oefs3Gs+OJAFqtgcrs7aQkdsOiCjzGBni7V9HHgu4AsZvj4yIwivUPM5yC4Z2IS9xw6induWGdkdOKRIhJ2/LwArA4DU+fqUsgqO5BFqmk4wG3Tv7pCtVSpukTOnLJfDHSVJ5zsS3l3JEu6tZz3IexrUriOES0Tc4QPG4vJuHRrD4tx2fQwP1BLlGvLHAL+7t7t+ONL1hrXUQGstK7CrqKmkGq/eF1nB1LsB/wATJeSiuxLvOafg3cfBywuBvtSAczGvgbXpn8sYdFvQcwcJhDi4ZM+mtpCrqdDmADKl4V9acs+PGfFQpwy0rHqYMa1dV4wlAFsMD7RAxIf55tCCuuXqlD47EACheBfCcBi6l9FDTAzgDHYjw6l0uduyfvVL0yuA+YlAtgwNLFevz12ZA0bEwuySgQCl3gtfRZSypAdT99aYWqU3mpkEPRtO5p7J3q45+BRXHXpOu06OkOrcSdSdx4SCHbiG1PIUAuFyJYMYKSt/+WY4wQw5ZB2AV7ZL1LPvoSaX6AEKbra3UdVvI/NvsT+GNGmk5ZgjEF46S6A5oArNMS5Vdcwsi0VyAwAZkpP/+6+bXjbxYPdSJeQHwRgGhblBWA8mgQwNfUE7NVX1eNDtdLHwZeOfYn6lg3AdOxLBjA7++LWilYTQySrr0cwxoovNdK0zNR012Kb7to03zauavgYW03eMJegr5t7y9b9+LmlCzB/pKvoXEnpeXw+lw50qzpZ8dkBYCL42QCsuJcGwKR2DwDTAh/JfWLZZyIJrMB3FQVw8hLwk07BwjiApWwEKeuiz7olAMv6RgzmVhn4+vm4XsqkLy7g889ZeyqBXC/SP+o0Zd5fdYKIXkREdxNRSkRPVfreQUQPEdH9RPR8n/WmxCfGgUzHzkJZWZ0Y5i5laNhY2KFeHz/cN4Y/vWy9VQcT/WC2FFEr5PN11Jr4vjYKUV/z2YXUsqgAAOPj1bX5NSmfFQuEk4Elyi6mSf+yVZ7Q7E66zKsiCxvsQMqsLJVYWUQm1tzLre8C8BsAPio2EtH5AF4M4AIAKwF8k4jOYYz1bYtNqU9MZGUc0ELTSdPnpsJll0gUcCpZITRgpM790H3b8YYLVpXWVnUwkxesCB2AKWmmdwrpy8C0jEllXUqaOEwBXyknTUlXBi9t+piloa6XdoSkj9I8xpT0UAYvNX0UGVpMYZ/FfBmG7T6M3QuUSQyAKwB8hjF2AsCjRPQQgKcD+E/belPqExPbVN2ssmA/TUDNa6zGBsGDA9h39xzGhYvmYdGsEW0aCdgZWLG+aScS0wDA+Jp8HIYEYGr6mD20OX2sULhQLZ1jqjwhA5je+2VjXzLIRdqdDCiKSEQbIb/0dhNjbFPNR1gF4DbhemveZo3GdydtwKSmmcMS7EP7E5G9VAifI0Y6RncsTfGN7QeKEjtqv1X7sgn56m6iGOrPWDeF5N9tAGYqoxMLwJSD2nX9X1XOPqpg1mckpYO9NC12H1XA6ikpZV8Bt6k4dpQDlhG0iOibAJZruq5ijH0p/OnM0TgT07npeYhppim9DGVgPkbZqlHebdSng/KYRBpjslRQQvj7e3fgdeetLErsBKWRMAAYDwWIVA1sAEoVbRQGgLKBkXYHsivsJiaJbKvoqu9sNLEvktPH7KEh1b73TB917ntT4UIxZSynj3r2pWpdNvbVSxnSNJ6/K6Ymxhi7vMK0bQDWCNer8zZrNJZsuXYT1dDtYk6XKIGXzTtmYXE2Mf/OA+NYNKuLVfNmeaeR2iNFgMKOqNTvBDAeVVJI/t0XwPgzGueEpo8CgBX9HaDTMQOYhn3x3UcXgGWpZtcbwHRivgpgvXQAYL00RS9lmOgx9HoMvX4azWnf1O6kJW4E8GIimkVEGwCcDeD7rkmNgRiPEJuECGQ+lok6VorY2pkPC5P6BBaWAvjUI7vxyrOXl/qzdeTPtgPbXl4wdW31+UrzampgVgZm+eIsTASxUABzvbg26Q7OP1rc96kIaOJYZfeRA1ifkWCbyMBoMh3YJ3RfYvpY8oZxcb+XRmNiLGXeX3WCiH6diLYC+CkANxPR1wCAMXY3gM8CuAfAVwH8vmtnEmh4d9KVUiZJuU9NLX1TRtNzuMZVBTMr4/JkYRz4PvnwLrzynGXoJO400mpoVdiY13EiwZQq7mA2lkLW0b80te6t9gmVfRkOb1d136tHh0ypoapvqbuPHLx6/QF49dOMhcVkRmlzu5NfAPAFQ997ALwnZL1GfWIxhPqm5vAoGVUV8NEBlGqrkOYbLBV8zrZjEzjc6+P8hfOsaaRd1B8AWJCVoi6AiYCT/SKaBbAQ/5cCYFDSR1v1CRHAXIULbZUnRLalWilc7CtlDL1eWrTFiN5ke3bSGDbW5MOudMeWdGuZ2qqEukZl3UuzI2myYTDG8NEHduCdT16riPVJMc9bB9M8j9NKYQIwPhYOABPWaATAxDcQ2QBMrNhqE/Adx4dc6aPqxFfBybT7KGtikHYfuXifpgyT/QFocQDjbTFiph47atRiYUsnqwKZbv2mwgfYXOPFtq/sOIBnLV+IuXk65mOnABCmg4lzdCAl7ERKLn8YAAwQvgtg2AQD4+2AnYGZAEwj4BcgZQEwH/e9CcBUu4TKxCZTBdgUAOsJ3/uRhfYWxDwilmF1qkrs2HYTfWwV6hrinPF+ih/sHcPVT9lQSiPr6mAmK4Uk4idJKfW0FjQ0pZClI0QVAEzU1DpC+pf9EvQCvliBwibgK+WjXa9O09UGs7nvM0BKreBl8n7p9C8d+xIF/ZggVlewn6qYVumkK0zVMIYFaFVNrtpUUPRzacT8j9y/Ha8/Xz5aZLJTFGsYAKwERuqc0u5jGcBKaapJxM/nl77XZWB8Hc8dyMr1vxz6l764ofn4kJgaiumjDsA4YInsS0wfRR2sl4MWZ18igE1G8ne1TMwRPruJITuOU8XGeOjMqiEv2RWB8a5DR7Fs7iiWzRmVmBZfX1deR+wvrWsS8nVWikTwgqmgqAOw7AcdfB9WCpk9iBeAldpUBqYBMC0DCzg+FJo+qqVzZAG/bF7VgVU/zdlYDmp9NgC8GNGbdLoZpmVMWToZ8tlXG1PH+wKdOsY2p4rWZdK2koTAKPOEveep5jSSz9MK+aYjRTAAmLITaaxk4WJgqnOe/+I0db0KABPTRO7EV9cz6V8FCHah1b9MBQx5n8S+TK76rpF9mergq4e3Te57NZ002SfE9DHNQUvVwzj74oDWamINRIyjQj7hO8c0TtduO6StY2ElX5fIojTp6ace2Y2XnLUU3a6YdspppPN8Y8GGygBnAzDTHCeAmVJJcU0dAxPHGcZbAcy2A2kCMEXANx0RShVDa+mt3EYAc7vvVTamnn1UAUwGK7ktu0YBYJO9dndyytPJmEbV0AjI/gBAC0KhIYr5eyd62HZsAq9YdEppXaNdwqaD5f0SQ+NhslKoc6RxwjOoKaQ2HZT73Clkvj5nZS4AE0R9rx1IlYFp6n/ZBHw3Ayv7v3S7jyYAs6WPvE1kWyKAqXNiRAtiHhGSNvqkhKHpYwhw+upbZhtEIrXrWNvf378db7pwtTaNVHUwraFVfAaNlcK0E2krqWO0UQDldLD0maeoYp+SdppATbeDqNG/SK3YKo4xCPiq/iUK+KY+XfUJsfa9j/teN0bnvlfNq5MCExswMOV6GLuTzRVFjBrNpJNCBdWpFuQBP6AD9KkkKeCSzVErSYipoV4bu23vGC5cNK94+a2aRqr3BWA+2K3ueopppMZKwefq3PgQ1pPYl3qUqPRZA2CltNEDwExivUvAVwDM+/2PBgAzv7zDfXzIlT7KzMovfeRtPH0sGNtkGq2YYcvEHCECWdFWkXm5jK88QsBymMCqsrDJNMWXt+4rifl8rNaVbzK0CpoW7/O1UpSOE40o7C6fY9SzdKAkjtO+jag6gFkPcasMrIoDXzSwBgCYWcD3c9/b0sdeT2Fs+dheL0XaZ2BpGlTM0Bb9Xrs76QwOZNOBjdlCb1jVp42JgWnpylXzsf/4yG686pzlBRBzFqbbjSw9g2poBcw6WDGHBgCWjxfnSj+TKuTX1cCyX5I/gOlAqVi7KwNYImppcgrpswPpevu2S/8yAZjL/zXZl1NBFcBU+0TRz9lXL0V/sj+YH+nMY5tOWsLmuvfpqwt6urlqxmQKUzqoamG2HUlx7M7jEzg82cc5C+ZqgcrkyrcZWq1nIlUrhQJ40pyKAKbVwGxsLREBiBQwE0BJ58L3LGKoHtL2OcBt0r/EVNH06jRX6Zxez3720bX7qGNfvcmsL5bTvnXse4YtrZwqhiYSBh4+5xZt4r+aJnLdbNODO/HWi9aUgU9JI4s1NBYLk5APwA5ganTqA5g71TQwMBuA2cro+OxAMuUN27YSOuIBbsnMKqeFKoCp+hcHO9/jQ7odSV36KDK2tM/Q7/UL9pWmg+sY0TIxR0yXFNL1DIkBsEyhgpGpDwC+u+sQLl18Ck4Z6Qhj9HPK4JSgpIMJfSUdTGcD8fWC1U0hJWAzpJA2BmYCMLWMjuYdkC4LhekNROLnrP69b/0v+fyjr//LJ30s0sWcgYkA1u+lSNMBG4sRLRPziIJxBWhjVcBPmz462kz3kHcdB5/LdcYSaZzK2Hpg+Nq2/bguF/N1JXZK6aUmjbQaWk1eMF8rRVUAi5FCmjxgahkd5S1EKbrOHUhTBQqd/lWl+oR2h7In2ydSAbRUhqW678V+vvuopo+9Xl9gY5HqgEXymzUdU6aJmcbUOTPpC3YugmUtseNjSBWui2qtD+3Cq85ZjkQR30UAk9bVmE515lQtWxT9YMMGsFgppA+AKfNT2I4JmQV802ff8482hmbyf6nsyqR/SYL/ZFqwr4LBldgYi/aCj7QtihgWLjY21emnSdAvjVP8XSoLA4BdJyYxNtnH2QvmSuvpQFBrp4AGkMT5vkK+yQs21SmkqoGZzkEqZXREBmY7QmQCLdknVgYwnYAviv2qgVWsOqECkk7/0rnvxfSRpQM2xgFMZF9ZX0RNrE0n7VFXvK+aeppSRulvrvhsF+p15tai38LCNj2wA2/RiPkA5FeuldhYApsOln33BzCjG386pJC2N3ELKaRoYlVFelH/8i0hLb7/sVzbPsw+UZhRNemj6r43pY9pytCf7Be7l1k6maWPaT6vn4Ncr5eNiybIRzLNuoKIXgTgzwCcB+DpjLEf5u3PBfBeAKMAJgC8jTH2b671poyJDSNiMDeTdWJwDzND01kqbts7hosXlcV8nakVQHUdTOMFI3WeeC/+y+p2B8CT/YBhKaQ0JyCF1GpgCoDpjhEZRHqbgK8DMG6h8PV/maqv9tKB/tXn7alcfaJw2wvuex24qemjeC2mjz0B5FgObDGiwd3JuwD8BoCPKu17AfwqY2w7EV0I4GuYTm8AB8xsbBg+sPA1LFqXlaEl2t3EJMle0SU6801Hi0pppHouEnAL+dLNy+koX6N0qFu8VsvjmACs21HWMACYIsJ7A5jGxKq+Rk3VvFxVWDkbU99ApKaIIoCVjgoJAKYK+NyMKlZfreL/MqWPg3aRseWfI9UBayqdZIzdC6B06oYx9iPh8m4Ac4hoFmPshG29KRX2uS7mOz4kjGkk6dvVUBmVi6HxdjE+9dgevOysZV5ifmkN8R6WzYSiPzSN1DIp5Zdm0sDUOb4amG8tMBHAlBd5VEkhpXRS8wo1nwoUvjuQJgNriP+L2ycGYj7/PADJfs7SWJ9FE+SnmbD/QgD/7QIwYJqlk0MFL8tnsQy1jYWpWpiOhXEtbP9kDzuPT+CC0+YJa+vTSJ0r33SwO/teUwcTwUesRgHI6aUOwAo21i23mTQw3evUVA0swANmqkJRdQdS/+q0QZvO/6UeBVIPaPNdSV31iWLHUrBP6NJHDmBcvM+OGw3YV5rrZdGsEQHrENFGABuFpk2MsU1C/zcBLC9NBK5ijH3JsfYFAP4SwPN8nuUJZ3bl4WJh0rVGtHcZWz96/w78/gWrFDG/zOZKAKbeT0gxrTpYVSsF/2UUwFQDwHQamPV9kIoGpgEw3zI6pl1H8bMKVuoOpAhg6vsfVQG/tJvIdDuQZnAz+r8mZQbmAjA22Y+mZYWkkzlgbbL0X17lGYhoNbIX676SMfawz5zmjx3V3KWspZ2ZpS0DUJlZmDhHHXfnwXGsO3U2Fo4Ofr2u3cjSQyqgpBXpRSEfZZCU/WKW9E9NFQc/jGW8RwrpC2A6DYxkAHKV0dEd4ubAptuB7AtamP4IkbkChXoY22cH0ibgi/4v1T5hSh/ZZB8szVPJSLuKU33siIgWArgZwNsZY//Pd960SidDIzR1dPXpPF4mFqYCGE8jGWP4l8f24NqnyGK+Clym3UijH6xT/pz9EIJupvGCBZlZAaDb1TMwW0HDKimkepBb8ICJLnybidX0GrXQHUjT+x/VA9x6t31qLJ+js1z0BWYl2idUNtbLQa0/ma3P0hT94/0ifUwnU7BezsimgInVCSL6dQB/C2AJgJuJ6A7G2PMB/AGAswC8m4jenQ9/HmNst229KQcxLu77VmW1hWmMHcg0JW8stgppnApIAD6/ZR+uWL94UDNfuFkiMCOdK19fKkd4WAm8DHYK/lylXUgHgOX6mFdJ6WGkkBVMrD4Almp2IE0CvsnAqlagEM83qucfdQAmpo82AFPTR26fSHPWJaaP6YnegIlF0sSaYmKMsS8gSxnV9usAXBe6XqO7k7F0MXFr1hu4NJmStKZG0NftSJqqryYJYbzXxz2HjuJFZy7N+gT2pD0byUMEKol91dDB+P11qaIBmCoBmG8KqdYCs4j4ujTRpIGpFgq1CoXPDqTuDCQX8HV2CRHAeBsHsOJaEf1tAr54fEhMH3X6F08fRQCLtavI2le2TU2owOibSmbXZeBSw6SF6VLPjz20ExvPXVECPGuJHVMaGVPI5z+8iYGVxilgpQj+0yGFVH1fpiKGqj1CZFmqgG96gYfq/7Id4Ladf1SrT2Rpojl9TCdTCahYmjOw/gDQOLjFiPbYUUDE3KkU//5M/bbxOouE1r+laGFie5IQHh0/jjndDlbMnaUV2427kYqQD6B0sLvMrgIALEm8jggNWJiGgYm/uFgppOc5SBG0fKpQqA77vpROZgBmEvCL1FJzhCjU/2XTv2zHh8T0sRDw+6zQv0RQ46wsGhNr64nN7NCxLR0L42PV1PCTD+3Cuy5dJ4DWwBNm3Y0EtKAHoGBhRaiILAr5gB7AtDugBgDzOQ9pc+KHpJACgOnSRJOJVVcHX1dGx2cHspcy6QiRCcDUF9ja6n/1hNRQPsztD2Ci/0sEKhXA+A5ltDSwobOTsWPGglgQ+yJ3KgnYWZj4XV3nUK+Py1efhlmaGvfqGlHTSL42b+sIP6hVyNcxMI90M8SJb6uHb/CBuQBM58K3AZj++JD8FiLbGUjTESJVwFf1L78dSD//F9OwMv65ZWJZEGMz88HbaKONNgAgkjLVRhtttDE10YJYG220MaOjBbE22mhjRkcLYm200caMjhbE2mijjRkdLYhNURARI6J/Fq67RLSHiL48lc/VRBDRQiJ6/VQ/RxsnR7QgNnUxDuBCIpqTXz8XwLamH4KyaPrfwUIALYi1ESVaEJvauAXAr+SfXwLg07yDiF5ORN8nojuI6KNE1Mnbv0hEtxPR3Xl1TRDRPCK6mYh+TER3EdFv5e3rieguYc23EtGf5e33E9E/Intpwxrd/fJx9xHRJ4noASK6noguJ6L/R0QPEtHTbc+bz7+XiD6WP+/Xc9B+L4Az87Hv0/1iiOhb+dtvQETXEdHfRvy9t3EyBWOs/ZqCLwBHAFwM4HMAZgO4A8AzAXwZ2ausbgIwko/9MLJKlwCwKP8+BxkAnY6sHvnHhLUX5N/XA7hLaH8rsldlrQeQAnhG3q69Xz6uB+AiZP/Dux3AJwAQgCsAfNFz/pPz9s8CeLn6XIbfz88D+DaAlyErlNeZ6v9m7df0/Jqxx45OhmCM/YSI1iNjYbcIXc8B8BQAP8jLDs0BwAvD/WFeVA4A1gA4G8CdAD5ARH8J4MuMse963H4zY+w2x/2+A+BRxtidAEBEdwO4lTHGiOhOZGDkM/+OfNzt+Zz/cD0cY+w7lC32ZmTgvo6IrkIG0Fd6/HxtPEGiBbGpjxsBvB/ZH+rpeRsB+L+MsXeIA4nomQAuB/BTjLGjRPRtALMZYw8Q0WUAfhnAdUR0K2PsGmQsSJQMZgufx8WlDfdbD0B820wqXKcY/Pvxnd9HBnDOIKKLAKwAsI8xNgZgDMD/IqLP+cxv44kTrSY29fEJAFdztpPHrQCuJKKlAEBEi4hoHYAFAA7kAHYugGfk/SsBHGWM/TOA9wG4LF9nF4ClRHQ6Ec0C8ALDM5ju5xuh88cAzDd1EtEKANcjS1mPENEvBjxLG0+waJnYFAdjbCuADypt9xDRuwB8Pd85nATw+wC+CuC1RHQvgPsB8HTwIgDvI6I0H/u6fJ1JIroGwPeR7XzeZ3gG0/12ev4MQfMZY/vyzYG7AHyFMfY23kdEcwF8HsBbGGP3EtG1yF7f9VWfZ2njiRdtFYs2ZkQQ0ekA3oPMivIPjLG/mOJHamOaRAtibbTRxoyOVhNro402ZnS0INZGG23M6GhBrI022pjR0YJYG220MaOjBbE22mhjRkcLYm200caMjhbE2mijjRkdLYi10UYbMzpaEGujjTZmdPz/EIo28d8MQv8AAAAASUVORK5CYII=\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 2 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"plt.figure(figsize=(6, 4), dpi=72)\n", | |
"ax = sns.heatmap(d_matrix,\n", | |
" cmap='RdYlBu_r',\n", | |
" square=True,\n", | |
" xticklabels=False,\n", | |
" yticklabels=False)\n", | |
"ax.contour(d_matrix, levels=0, colors='black', linewidths=.5)\n", | |
"ax.set(xlabel=r'Measurement $x_1$', ylabel=r'Measurement $x_2$')\n", | |
"plt.show()" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"id": "spectacular-moscow", | |
"metadata": {}, | |
"source": [ | |
"## Question 2" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 6, | |
"id": "appreciated-paraguay", | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"s1 = np.zeros(300)\n", | |
"s2 = np.linspace(-30, 30, 300)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"id": "through-yield", | |
"metadata": {}, | |
"source": [ | |
"### Method 2: Numerical integration on a grid" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 7, | |
"id": "corporate-flight", | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEICAYAAABBBrPDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3lUlEQVR4nO3deZyNdf/H8dfHWGM0WSohJCkhMeGuHzdKltwRkaXum2SrYcxYUkqSJYxtGCKVZE+R7mSLdMtNjX2LrJGKm2zJ/vn9cY7uuTXGnGNmvnPO+Twfj+txzrnO91zX+6pxPufavl9RVYwxxpiryeI6gDHGmMzNCoUxxpgUWaEwxhiTIisUxhhjUmSFwhhjTIqsUBhjjEmR00IhInVFZLuI7BSRXsm8f7uILBORdSKyUUTqu8hpjDGhTFzdRyEiYcAOoDZwAPgWaKGqW5O0mQCsU9VxIlIGmK+qxVNaboECBbR48RSbGGOMucKaNWv+o6oFk3sva0aHSaIysFNVdwOIyAygIbA1SRsF8nqf3wgcvNZCixcvTmJiYhpHNcaY4CYi+672nstDT4WB/UleH/DOS6ov8LSIHADmA52TW5CItBeRRBFJPHz4cHpkNcaYkJXZT2a3ACapahGgPvCBiPwps6pOUNVIVY0sWDDZPSdjjDF+clkofgSKJnldxDsvqbbALABV/TeQEyiQIemMMcYAbgvFt0ApESkhItmB5sC8K9r8ADwMICL34CkUdmzJGGMykLNCoaoXgChgIbANmKWqW0Skn4g87m3WDWgnIhuA6UBrte5ujTEmQ7m86glVnY/nJHXSeX2SPN8KPJTRuYwxxvxXZj+ZbYwxxjGnexTGBDNVZdeuXWzYsIF9+/Zx+vRpRITChQtzxx13ULlyZXLmzOk6pjHXZIXCmDS2YcMGJk2axEcffcT+/fuv2i5HjhxUr16dNm3a0LhxY3LkyJGBKY1JPTv0ZEwa+eabb6hfvz4VKlRg3Lhx3H///YwfP55vvvmGo0ePcu7cOX7//Xd27tzJvHnzeP7559m5cyctW7akaNGijBo1irNnz7reDGP+TFWDaqpUqZIak5GOHDmi7du3VxHRggUL6oABA/TXX39N1WcvXryoixcv1lq1aimgd9xxhy5dujR9AxuTDCBRr/K9ansUxlyHFStWUL58ed555x1iY2PZtWsXL7/8MhEREan6fJYsWXjkkUdYsmQJCxYsICwsjFq1ahEdHW17FybTsEJhjB9UlWHDhlGjRg1y5crFN998Q1xcHOHh4X4tT0SoU6cO69evp3PnzsTHx1OzZk1+/vnnNE5ujO+sUBjjo4sXLxIVFUX37t154oknWLNmDRUrVkyTZd9www3Ex8cza9YsNmzYwAMPPMCWLVvSZNnG+MsKhTE+OH/+PM2bN2fs2LH06NGDmTNnkjdv3mt/0EdNmzbl66+/5uLFi/z1r39lzZo1ab4OY1LLCoUxqXThwgVatWrF7NmzGTZsGEOGDCFLlvT7J1ShQgX+9a9/kSdPHmrVqmXjrBhnrFAYkwqXLl2iTZs2fPjhhwwbNozY2NgMWW/JkiX517/+Rb58+ahXrx7fffddhqzXmKSsUBiTCr1792bKlCn0798/w4rEZUWLFmXx4sWEhYVRu3ZtDhw4kKHrN8YKhTHXMH78eN588006dOjAyy+/7CTDnXfeycKFCzl+/DgNGzbk9OnTTnKY0GSFwpgUrFy5kqioKOrVq8eYMWMQEWdZ7rvvPqZNm8a6deto3bo1aj3umwxihcKYqzh06BBNmzalWLFiTJs2jaxZ3XeN1qBBAwYPHsyHH35IXFyc6zgmRFihMCYZFy9epEWLFhw9epTZs2en+k7rjNC9e3eaNGnCyy+/zKpVq1zHMSHACoUxyejTpw9Lly5l7NixVKhQwXWc/yEiTJw4kSJFitCiRQuOHTvmOpIJck4LhYjUFZHtIrJTRHol8/4IEVnvnXaIyDEHMU2I+eqrrxg0aBDPPvssbdq0cR0nWREREcyYMYMDBw7w3HPP2fkKk66cFQoRCQMSgHpAGaCFiJRJ2kZVY1S1gqpWAEYDH2d4UBNSTp06RevWrSlRogSjRo1yHSdFVapUYdCgQXz00Ue88847ruOYIOZyj6IysFNVd6vqOWAG0DCF9i2A6RmSzISs7t27s3fvXiZNmkSePHlcx7mm2NhYatasSbdu3ez+CpNuXBaKwkDS4b8OeOf9iYgUA0oAS6/yfnsRSRSRxMOHD6d5UBMaFi5cyPjx4+nWrRvVqlVzHSdVsmTJwttvv82FCxfo0KGDHYIy6SJQTmY3B2ar6sXk3lTVCaoaqaqRBQsWzOBoJhicPHmStm3bUqZMGd544w3XcXxSsmRJBg4cyPz585kyZYrrOCYIuSwUPwJFk7wu4p2XnObYYSeTjvr06cPBgwd59913yZkzp+s4PouKiuLBBx8kOjqan376yXUcE2RcFopvgVIiUkJEsuMpBvOubCQidwM3Af/O4HwmRKxbt474+Hg6duxIlSpVXMfxS1hYGO+++y6nT5+me/furuOYIOOsUKjqBSAKWAhsA2ap6hYR6Scijydp2hyYoXbw1aSDixcv0rFjRwoUKMDAgQNdx7kupUuXpmfPnkybNo3ly5e7jmOCiATb929kZKRav/0mtd566y06derElClTaNWqles41+306dOUKVOG8PBw1q5dS7Zs2VxHMgFCRNaoamRy7wXKyWxj0tx//vMfXnrpJWrVqkXLli1dx0kTN9xwAyNHjmTz5s0kJCS4jmOChBUKE7Jee+01Tp48yejRo532CpvWGjZsSN26dXnttdf4+eefXccxQcAKhQlJW7duZfz48XTs2JEyZcpc+wMBRESIj4/nzJkz9O7d23UcEwSsUJiQ1K1bN/LkyUPfvn1dR0kXpUqVIioqivfee4+NGze6jmMCnBUKE3IWLFjAggUL6NOnDwUKFHAdJ9307t2biIgIevbs6TqKCXBWKExIuXDhAt26daNkyZK88MILruOkq3z58vHKK6+wcOFCFi1a5DqOCWBWKExI+eCDD9i6dSuDBw8mR44cruOkuxdeeIESJUrQo0cPLl5MtgccY67JCoUJGWfPnqVv375ERkbSuHFj13EyRI4cORg0aBAbN27kgw8+cB3HBCgrFCZkTJgwgR9++IGBAwcG1eWw19KsWTMiIyPp27cv586dcx3HBCArFCYk/Pbbb/Tv358aNWrwyCOPuI6ToUSEN954g3379vHuu++6jmMCkBUKExLi4+M5dOgQAwYMCKm9icvq1KnDgw8+SP/+/Tlz5ozrOCbAWKEwQe/XX39lyJAhNGjQgAcffNB1HCdEhP79+/Pjjz8yfvx413FMgLFCYYJeXFwcx44do3///q6jOFWzZk1q1qzJoEGD+O2331zHMQHECoUJaocOHWLkyJE0b96c++67z3Uc59544w1++eUX6zDQ+MQKhQlqw4cP5/fffw/arjp89dBDD1G3bl2GDBnCiRMnXMcxAcIKhQlaR44cISEhgaeeeorSpUu7jpNp9OvXjyNHjhAfH+86igkQVihM0Bo1ahSnTp2yHlSv8MADD9CgQQNGjBjBqVOnXMcxASDNCoWItPHjM3VFZLuI7BSRXldp00xEtorIFhGZdv1JTSg4fvw48fHxNG7cmLJly7qOk+n07t2bo0eP2hVQJlXSco/idV8ai0gYkADUA8oALUSkzBVtSgEvAQ+p6r1A17SJaoLd6NGjOX78uO1NXEXVqlV5+OGHiYuLs/sqzDX5VChEZONVpk3ALT6uuzKwU1V3q+o5YAbQ8Io27YAEVf0VQFUP+bgOE4JOnTrFiBEjeOyxx6hYsaLrOJlW7969+fnnn3nvvfdcRzGZnK97FLcAfwf+lsx0xMdlFQb2J3l9wDsvqbuAu0TkaxFZJSJ1fVyHCUHjxo3j6NGjvPrqq66jZGo1atTgL3/5C4MHD+b8+fOu45hMzNdC8U8gj6ruu2LaC3yZ5ukgK1AKqAG0AN4WkYgrG4lIexFJFJHEw4cPp0MMEyhOnz5NXFwctWvXpkqVKq7jZGoiQu/evdm3bx9Tp051HcdkYj4VClVtq6orrvJeSx/X/SNQNMnrIt55SR0A5qnqeVXdA+zAUziuXPcEVY1U1ciCBQv6GMMEk0mTJnHo0CFeeeUV11ECQv369alQoQKDBg2y8SrMVbm8PPZboJSIlBCR7EBzYN4Vbebi2ZtARArgORS1OwMzmgBy4cIFhg0bRtWqValWrZrrOAFBRHj55ZfZsWMHH330kes4JpNyVihU9QIQBSwEtgGzVHWLiPQTkce9zRYCR0RkK7AM6KGqvp4LMSHio48+Yvfu3fTs2TMke4j1V+PGjSldujRvvvkmquo6jsmExN8/DBGppapLLz+mcS6/RUZGamJiousYJoOpKpGRkZw6dYqtW7cSFhbmOlJAmThxIu3atWPJkiU8/PDDruMYB0RkjapGJvfe9exRxF3xaIwzS5cuZe3atXTv3t2KhB+efvppbrnlFoYOHeo6ismE0uLQk+3jG+eGDBnCLbfcwjPPPOM6SkDKmTMnXbp0YeHChWzcuNF1HJPJWF9PJuCtX7+eRYsW0bVrV3LmzOk6TsDq1KkTuXPnJi7ODhKY/2WFwgS8oUOHkidPHjp27Og6SkC76aabeO6555g+fTr79++/9gdMyLBCYQLa3r17mTlzJh06dCAiIsJ1nIAXExODqjJq1CjXUUwmcj2F4nL/xCfTIogx/hgxYgQiQteuXV1HCQrFihWjWbNmTJgwgePHj7uOYzIJvwuFqlZP+mhMRjty5AgTJ06kVatWFClSxHWcoNG9e3dOnjxpXZCbP9ihJxOwxo8fz+nTp+nevbvrKEGlYsWK1KpVi1GjRnHu3DnXcUwmYIXCBKRz586RkJBA7dq1bWCidNCjRw8OHjzItGk2VpjxfTyKV0XkhfQKY0xqffjhhxw8eJCYmBjXUYJSnTp1KFeuHMOGDbNuPUzKhUJE7hWRpP0PNwPeSabdcyLyUlqHMyY5qsrIkSMpXbo0derUcR0nKIkIMTExbN68mS+++MJ1HOPYtfYolgBJx5I8r6rJjZv4AfB0mqUyJgUrV64kMTGR6OhosmSxo6fppUWLFtx8882MHDnSdRTj2LX+lT0KDEjy+pyIFLqykaqeBWyILJMhRowYwU033cTf//5311GCWs6cOenUqROfffYZO3bscB3HOJRioVDVTaraKsmsYcAnIlIsaTsRuRmwA5km3e3du5c5c+bQvn17cufO7TpO0OvUqRPZs2e3G/BCnK8j3H0IJABrROSfItJfRAYCX2O9yJoMMGbMGESEF16wayoywi233ELLli2ZNGkSv/76q+s4xhGfD/Cq6vtACWAmkA04A7RUVRt016SrkydPMnHiRJo2bUrRokWv/QGTJrp27crp06d5++23XUcxjvh1JlBVT6rqB6r6oqr2U9Vv0zqYMVeaNGkSx48ft+46Mth9991HzZo1GT16NOfP26nIUOT0khERqSsi20Vkp4j0Sub91iJyWETWe6fnXOQ07l26dIlRo0ZRtWpVqlSp4jpOyImJieHAgQN8/PHHrqMYB5wVChEJw3O+ox5QBmghImWSaTpTVSt4p4kZGtJkGp999hm7du2yvQlHHnvsMe688067VDZEudyjqAzsVNXdqnoOmAE0dJjHZGIjRoygaNGiNGnSxHWUkJQlSxaio6NZtWoVq1atch3HZDC/CoWINBWRcO/zV0TkYxGp6ONiCgNJR0c54J13pSYislFEZouIncEMQRs2bGDZsmVERUWRNWtW13FCVuvWrbnxxhttryIE+btH8aqqnhSR/wMewdOtx7i0i/WHT4HiqloeWAy8n1wjEWkvIokiknj48OF0iGFcGjVqFDfccAPt2rVzHSWk5cmTh3bt2jF79mx++OEH13FMBvK3UFzyPj4GTFDVz4DsPi7jRyDpHkIR77w/qOoR713fABOBSsktSFUnqGqkqkYWLFjQxxgmMzt06BBTp06ldevW3HTTTa7jhLzOnTsDkJCQ4DiJyUj+FooDIjIeeAqYLyI5APFxGd8CpUSkhIhkB5oD85I2uKK7kMeBbX7mNQFq3LhxnDt3ji5duriOYoDbb7+dxo0bM2HCBE6dOnXtD5ig4G+h2AosBOqo6jEgH3DRlwWo6gUgyrucbcAsVd0iIv1E5HFvsy4iskVENgBdgNZ+5jUB6OzZs4wdO5b69etTunRp13GMV0xMDMeOHWPy5Mmuo5gMIv70NS8ia1W14hXzNnrPJTgVGRmpiYmJrmOYNPD+++/TunVrFi9ezCOPPOI6jvFSVapWrcqvv/7Kd999Zz34BgkRWaOqkcm95+vARZ1EZBNwt/dKpMvTHmBjWoQ1BjxfRiNGjODee+/l4Ycfdh3HJHF5rIrvv/+ezz//3HUckwF8/SkwDfgb8In38fJUSVVtPAqTZpYvX86GDRvo2rUrIr6e/jLprUmTJhQuXNgulQ0RvvYeexzYB/RU1X1JpqPpE8+EqpEjR1KgQAFatWp17cYmw2XLlo2oqCiWLFnCpk2bXMcx6cyf3mMVmJ8OWYwBYNeuXcybN4+OHTuSK1cu13HMVbRv355cuXLZWBUhwN+zUGtF5IE0TWKMV3x8PFmzZuX55593HcWkIF++fPzjH/9gypQp2I2uwc3fQlEF+LeI7PKezN4kInYy21y348eP8+6779K8eXMKFfrTqLsmk+nSpQtnz55l/PjxrqOYdORvoagDlARq4TmZ3cD7aMx1effddzl16hTR0dGuo5hUuOeee6hbty4JCQmcPXv22h8wAcnfgYv2ARH896qnCO88Y/x28eJF4uPjqVatGpUqJdtbi8mEYmJi+Pnnn5k1a5brKCad+Nt7bDQwFbjZO00Rkc5pGcyEnk8++YS9e/cSExPjOorxQe3atbnnnnsYMWIE/tzAazI/fw89tQWqqGofVe0DVAWsa09zXUaMGEHx4sV5/PHHr93YZBoiQteuXVm3bh0rVqxwHcekA38LhfC/fTtdxPdOAY35Q2JiIitWrCA6OpqwsDDXcYyPnnnmGfLnz8+IESNcRzHpwN9C8R6wWkT6isjrwCo8Y1IY45cRI0YQHh7Os88+6zqK8UOuXLno0KEDc+fOZffu3a7jmDTm78ns4UAb4CjwH6CNqo5Mw1wmhPz444/MmjWLtm3bkjdvXtdxjJ+ef/55wsLCGDNmjOsoJo35ezI7J1ADqOl9rOGdZ4zPEhISuHTpko05EeAKFy5Ms2bNmDhxIidOnHAdx6Qhfw89TQbuBeKBMUAZ4IO0CmVCx+nTpxk/fjyNGjWiRIkSruOY6xQTE8PJkyd57733XEcxacjfQlFWVduq6jLv1A5P4TDGJ5MnT+bo0aN2SWyQiIyM5KGHHiI+Pp6LF30ay8xkYtfT11PVyy9EpApgowUZn1y6dImRI0f+8eVigkPXrl3ZvXs3//znP11HMWnE30JRCVgpIntFZC/wb+AB6/PJ+GLBggVs376dmJgYG3MiiDRq1IhixYrZpbJBxN9CURcoAfzVO5XwzvOpzycRqSsi20Vkp4j0SqFdExFREUl2mD4TmEaOHMltt93Gk08+6TqKSUNZs2alc+fOLF++nHXr1rmOY9KAv4XiB6Aa8A9vH08K3HJ5IKPULEBEwoAEoB6ek+EtRKRMMu3CgWhgtZ9ZTSa0efNmFi9eTFRUFNmzZ3cdx6Sxtm3bkjt3bhurIkj4WyjGAn8BWnhfn8Tzpe+LysBOVd2tqueAGUDDZNq9AQwGzviZ1WRCI0eO/OMmLRN8IiIiaNOmDdOnT+fnn392HcdcJ7/Ho1DVF/B+eavqr4CvPwsLA/uTvD7gnfcHEakIFFXVz1JakIi0F5FEEUm0AVQyv0OHDjFlyhT+8Y9/kC9fPtdxTDrp0qUL58+fZ9y4ca6jmOvkb6E47z10pAAiUhC4lGapPMvMAgwHul2rrapOUNVIVY0sWLBgWsYw6eCtt97i7NmzdO3a1XUUk45KlSpFgwYNGDduHGfO2AGBQOZvoYgH5gA3i8gAYAUw0Mdl/AgUTfK6iHfeZeFAWeBL75VVVYF5dkI7sJ09e5axY8dSv359Spcu7TqOSWddu3bl8OHDTJ8+3XUUcx18LhTiuY7xK6AnMAj4CWikqh/6uKhvgVIiUkJEsgPNgXmX31TV46paQFWLq2pxPB0PPq6qdr9GAJs+fTq//PKL7U2EiJo1a1K+fHkbqyLA+Vwo1PN/e76qfqeqCao6RlW3+bGcC0AUsBDYBsxS1S0i0k9EbECCIKSqDB8+nLJly/LII4+4jmMywOWxKjZt2sSSJUtcxzF+up47sx+43pWr6nxVvUtVS6rqAO+8Pqo6L5m2NWxvIrAtXLiQTZs20b17d7vBLoS0bNmSQoUKMXToUNdRjJ/8vuoJ+LeI7BKRjXZHtkmNoUOHUrhwYVq0aHHtxiZo5MiRg+joaBYvXmw34AUofwtFHaAkUAvPndg+3ZFtQs/atWtZunQp0dHRdoNdCOrQoQPh4eHExcW5jmL84O/ARfuSm9I6nAkeQ4cOJTw8nPbt27uOYhyIiIigffv2zJw5k3377Ksi0Pi7R2FMqu3du5cPP/yQDh06cOONN7qOYxzp2rUrImKdBQYgKxQm3Y0YMQIRITo62nUU41CRIkVo2bIlb7/9NkePHnUdx/jACoVJV0ePHmXixIm0bNmSIkWKuI5jHOvevTunT5+2bj0CTFZ/PiQiOYAmQPGky1DVfmkTywSLcePGcfr0abp37+46iskEypUrR7169YiPj6dbt27kzJnTdSSTCv7uUXyCp6fXC8BvSSZj/nDmzBlGjx5NnTp1KFeunOs4JpPo0aMHhw4dYvLkya6jmFTya48CKKKqddM0iQk6H3zwAb/88gs9evRwHcVkIjVq1CAyMpK4uDjatm1LWFiY60jmGvzdo1gpIvYT0VzVxYsXiYuLo2LFitSqVct1HJOJiAg9evTg+++/Z968P3XCYDIhfwvF/wFrvMOY2p3Z5k9mz57Njh07eOmll6y7DvMnjRs35o477mDIkCHWWWAAEH/+J4lIseTmZ4ab7iIjIzUx0bqEcklVqVChAufOnWPLli1kyWIX15k/S0hIICoqiuXLl1O9enXXcUKeiKxR1WSHcfD7zmwgAk+3HX8DIjJDkTCZw2effcbGjRvp1auXFQlzVc8++yw333wz/fv3dx3FXINf/4pFJBqYCtzsnaaISOe0DGYCk6oyYMAAihcvTsuWLV3HMZlYrly56N69O4sXL2b16tWu45gU+Ptzry2ecbP7qGofPKPPtUu7WCZQffnll6xatYqePXuSLVs213FMJtexY0fy5cvHgAEDXEcxKfC3UAhwMcnri955JsQNGDCAW2+9lTZt2riOYgJAeHg4Xbt25dNPP2X9+vWu45ir8LdQvAesFpG+ItIXzzCl76RZKhOQVq9ezRdffGF33BqfdO7cmbx58zJw4EDXUcxV+HsyezjwLHDUO7VR1ZG+LkdE6novsd0pIr2Seb+j99Lb9SKyQkTK+JPXZIyBAweSL18+Onbs6DqKCSARERF07tyZ2bNns22bz6Mqmwzg9yUpqrpGVeO9k8/DVolIGJAA1APKAC2SKQTTVLWcqlYAhgDD/c1r0temTZuYN28e0dHR5MmTx3UcE2C6du1Krly5GDRokOsoJhk+FQoRWeF9PCkiJ5JMJ0XkhI/rrgzsVNXdqnoOmIGn/6g/qGrSZeYG7M6cTGrgwIHkyZOHzp3t4jfjuwIFCtCpUyemTZvGrl27XMcxV/CpUKjq/3kfw1U1b5IpXFXz+rjuwsD+JK8PeOf9DxF5QUR24dmj6JLcgkSkvYgkikji4cOHfYxhrtfWrVuZOXMmUVFR3HTTTa7jmADVrVs3smbNyuDBg11HMVfw9z6KP/2fTG5eWlDVBFUtCbwIvHKVNhNUNVJVIwsWLJgeMUwK+vXrR+7cua0rcXNdChUqxHPPPcekSZPYv3//tT9gMoy/5yhqJzOvno/L+BEomuR1Ee+8q5kBNPJxHSadbd68mVmzZhEdHU3+/PldxzEBrmfPnoiI3VeRyfh6jqKTiGwC7vZ2Bnh52gP42ingt0ApESkhItmB5sD/dCUpIqWSvHwM+N7HdZh09vrrrxMeHk5sbKzrKCYI3H777bRr14533nmHPXv2uI5jvHzdo5iG54TzKv7bz9PfgEqq+rQvC1LVC0AUsBDYBsxS1S0i0k9EHvc2ixKRLSKyHogF/uFjXpOONmzYwOzZs+natSv58uVzHccEiZdffpmsWbPSr58NmJlZ+Nt77GZVLZsOea6b9R6bcRo3bszSpUvZu3cvERERruOYINKtWzdGjhzJtm3buOuuu1zHCQlp3nssnrEoHriOTCbArVu3jjlz5hAbG2tFwqS5F198kVy5ctG3b1/XUQz+F4oqwL9FZJcNXBSaXnvtNSIiIoiOjnYdxQShm2++mc6dOzNjxgw2b97sOk7I87dQ1AFKArXwnKNo4H00IWDlypV8+umn9OjRgxtvvNF1HBOkevToQXh4OK+99prrKCHPBi4yPlFVXnzxRW699VbbmzDpKl++fMTExPDxxx+zZs0a13FCmg1cZHzy2WefsWLFCl577TVy587tOo4JcjExMeTPn59evf7UZ6jJQP5e9bQR+Iuq/uZ9nRv4t6qWT+N8PrOrntLPxYsXue+++/4YC9sGJjIZYeTIkcTExLBw4UIeffRR13GCVnpc9WQDF4WgKVOmsGXLFgYMGGBFwmSYTp06UaJECXr27MmlS5dcxwlJaTFw0evYwEVB78yZM7z66qs88MADPPnkk67jmBCSI0cOBg4cyIYNG5g6darrOCHpegYuaoNn0KL/4OfARSZwjB07lv379/Pmm28iYjuPJmM1a9aMSpUq8corr3DmzBnXcUKOvyezcwI1gJrexxreeSYIHT16lAEDBvDoo49Sq1Yt13FMCMqSJQtDhw7lhx9+YMyYMa7jhBx/Dz1NBu4F4oExeEao+yCtQpnM5fXXX+fYsWMMHTrUdRQTwmrWrEm9evUYMGAAR48edR0npPhbKMqqaltVXead2uEpHCbIbNu2jYSEBNq1a0f58s4vajMhbvDgwZw4ccI6DMxg/haKtSJS9fILEakC2DWpQSg2NpY8efLwxhtvuI5iDOXKlaNdu3aMGTOGrVu3uo4TMvwtFJWAlSKyV0T2Av8GHrA+n4LL/PnzWbBgAX369MFGDjSZRf/+/QkPD6dLly74cx+Y8Z2/N9wVS+l9l9152A13aeP8+fOUK1cOVWXTpk1kz57ddSRj/jBmzBg6d+7Mxx9/zBNPPOE6TlBI8xvurtbX0+XJ76Qm00hISGD79u0MGzbMioTJdDp27EjZsmWJjY3l999/dx0n6FlfT+ZPDh06xOuvv86jjz7KY4895jqOMX+SNWtW4uPj2bt3L8OGDXMdJ+j5e46iLVBFVfuoah+gKtDO14WISF0R2S4iO0XkT71+iUisiGz1jnnxxbUOeZm00a1bN3777TdGjRplN9eZTKtmzZo8+eSTDBw4kP3797uOE9Sc9fUkImFAAlAPz30YLUSkzBXN1gGR3s4GZwND/MxrUmnp0qVMmTKFF198kbvvvtt1HGNSFBcXh6pal/fpLC36euqLf309VQZ2qupuVT0HzAAaJm3gvUfjtPflKqCIn3lNKpw9e5ZOnTpRsmRJXn75ZddxjLmmYsWK0bdvX+bMmcPcuXNdxwlaPhcK8RyL+JD/9vV0FP/6eioMJN1fPOCddzVtgc+vkqm9iCSKSOLhw4d9jGEuGzx4MDt27GDs2LHkypXLdRxjUiU2Npby5csTFRXFiRMnXMcJSj4XCvVcTztfVdeqarx3WpcO2f4gIk8DkUCyfUio6gRVjVTVSLve3z/ff/89AwcOpHnz5tbnvwko2bJlY8KECRw8eJBXXnnFdZygdD13Zj9wnev+ESia5HUR77z/ISKPAL2Bx1X17HWu0yRDVXn++efJkSMHw4cPdx3HGJ9VqVKFqKgoxowZw+rVq13HCTr+FooqwCoR2eW9IsmfO7K/BUqJSAkRyQ40B+YlbSAi9wPj8RSJQ35mNdcwadIklixZwqBBgyhUqJDrOMb4pX///tx22220b9+e8+fPu44TVPwtFHWAO4BaeG64a+B9TDVVvQBEAQuBbcAsVd0iIv1E5HFvs6FAHuBDEVkvIvOusjjjpwMHDhATE0P16tXp2LGj6zjG+C1v3rwkJCSwceNG4uLiXMcJKj514eEdc6IjcCewCXjH+4WfaVgXHqmnqjRo0IBly5axadMmSpYs6TqSMdetWbNmzJ07l8TEROvx2Adp2YXH+3hOKm/Cc/+D3RIZwN5//33mz5/Pm2++aUXCBI2xY8eSP39+nnnmGc6etdOaacHXQlFGVZ9W1fHAk0C1dMhkMsCBAwfo2rUr1apVIyoqynUcY9JMgQIFmDhxIhs3bqRv376u4wQFXwvFH2eIMtshJ5N6Fy9e5JlnnuHChQu8++67ZMni76kqYzKnxx57jHbt2jFkyBC+/vpr13ECnq/fEPeJyAnvdBIof/m5iNidLgEiLi6OL7/8ktGjR3PnnXe6jmNMuhg2bBjFihXj73//O6dOnXIdJ6D5VChUNUxV83qncFXNmuR53vQKadLOmjVreOWVV2jatCmtW7d2HceYdBMeHs7777/Pnj17iImJcR0noNkxhxBy6tQpWrZsya233spbb71lPcOaoFetWjV69erFxIkTmTZtmus4AcsKRYhQVTp27Mj333/P5MmTyZcvn+tIxmSIfv368dBDD9GhQwd27NjhOk5AskIRIsaPH8/UqVPp168fNWvWdB3HmAyTNWtWpk+fTvbs2WnWrBlnzpxxHSngWKEIAWvWrCE6Opp69epZ9+EmJBUtWpTJkyezYcMGunTp4jpOwLFCEeSOHj3Kk08+yS233MIHH3xgl8KakPXYY4/Rq1cv3n77bSZMmOA6TkDJ6jqAST8XLlzgqaee4uDBgyxfvpz8+fO7jmSMU/3792fdunVERUVRtmxZHnzwQdeRAoL9vAxisbGxLFmyhPHjx1O1alXXcYxxLiwsjOnTp1OsWDGaNGnCjz/+aWQDkwwrFEHq7bffZvTo0cTGxtr9EsYkcdNNNzF37lxOnTpFo0aN+O2331xHyvSsUAShZcuW8fzzz1O3bl2GDBniOo4xmc69997LtGnTWLt2LS1btuTixYuuI2VqViiCzIYNG2jUqBGlS5dm+vTphIWFuY5kTKb0t7/9jVGjRjFv3jxiY2Ndx8nU7GR2ENm3bx/16tUjPDyczz//nIiICNeRjMnUoqKi2LNnD8OHD6d48eLW1cdVWKEIEkeOHKFu3br8/vvv/Otf/6Jo0aLX/pAxhqFDh7Jv3z5iY2O56aab7JxeMpweehKRuiKyXUR2ikivZN6vLiJrReSCiDzpImMgOHbsGI8++ih79uzhk08+oWzZsq4jGRMwsmTJwtSpU6lduzZt27Zl9uzZriNlOs4KhYiEAQl4RsorA7QQkTJXNPsBaA1Yb15XceLECerWrcumTZv4+OOPqV69uutIxgScHDlyMGfOHP7yl7/QsmVLFixY4DpSpuJyj6IysFNVd6vqOWAG0DBpA1Xdq6obgUsuAmZ2v/32G4899hhr1qxh1qxZ1K9f33UkYwJW7ty5+eyzzyhXrhyNGjVi/vz5riNlGi4LRWFgf5LXB7zzfCYi7UUkUUQSDx8+nCbhMrvLh5tWrlzJ1KlTadSoketIxgS8G2+8kUWLFlG2bFkaNWrE3LlzXUfKFILi8lhVnaCqkaoaWbBgQddx0t2hQ4eoWbMm3377LTNnzqRZs2auIxkTNPLnz8+SJUuoVKkSTz75JDNnznQdyTmXheJHIOmlOUW880wK9u/fT/Xq1dm+fTvz5s3jySftHL8xaS0iIoJFixbx4IMP0qJFC0aPHu06klMuC8W3QCkRKSEi2YHmwDyHeTK9tWvXUqVKFX766ScWLlxI3bp1XUcyJmiFh4ezYMECGjZsSJcuXejZsyeXLoXm6VJnhUJVLwBRwEJgGzBLVbeISD8ReRxARB4QkQNAU2C8iGxxlde1Tz/9lGrVqpEtWza+/vprqlWr5jqSMUHvhhtuYPbs2Tz//PMMHTqUVq1a8fvvv7uOleGc3nCnqvOB+VfM65Pk+bd4DkmFLFVl2LBh9OzZk0qVKvHpp59y6623uo5lTMgICwtjzJgx3H777fTq1YsdO3YwZ84cbr/9dtfRMkxQnMwOVidOnKBp06b06NGDJk2asHz5cisSxjggIrz44ovMmzeP77//nsjISJYvX+46VoaxQpFJbdmyhcqVKzN37lyGDRvGrFmzuOGGG1zHMiak/e1vf+Obb74hX758PPzwwwwYMCAkep61QpHJXLp0ifj4eCpVqsSxY8f44osviI2NRURcRzPGAHfffTerV6+madOmvPLKK9SuXZuDBw+6jpWurFBkIgcPHqR+/fpER0fzyCOPsGHDBv7617+6jmWMucKNN97ItGnTeOedd1i9ejXly5dnxowZqKrraOnCCkUmcOnSJSZMmECZMmX46quvGDt2LJ9++im33HKL62jGmKsQEZ599lnWrFlDyZIladGiBY0bN+ann35yHS3NWaFwbOvWrdSoUYMOHTpw//33s379ejp16mSHmowJEHfffTcrV64kLi6OBQsWUKZMGcaMGcOFCxdcR0szVigcOXLkCF26dKF8+fJs3ryZd955h6VLl3LXXXe5jmaM8VFYWBjdunVjw4YNREZG0rlzZypVqhQ0V0ZZochgp0+fZvjw4ZQqVYqEhATatWvH9u3befbZZ20vwpgAd9ddd7Fo0SI++ugjjh8/To0aNWjQoAHr1693He26WKHIIJcLxB133EG3bt2IjIxk/fr1jBs3jlDoyNCYUCEiNG7cmK1btzJw4EC+/vpr7r//fpo1a8Z3333nOp5frFCks6NHjzJkyBBKlChBt27dKFu2LF999RWLFi2iXLlyruMZY9LJDTfcwEsvvcSePXt49dVX+fzzz7n33ntp0qQJK1asCKgrpKxQpJNNmzbRvn17ihQpwosvvki5cuX46quvWLJkifXTZEwIiYiIoF+/fuzevZsXX3yRZcuWUa1aNapUqcL06dM5d+6c64jXZIUiDR07doy3336batWqUb58eT744ANatWrF+vXrrUAYE+IKFizIwIED2b9/P2PHjuX48eO0bNmS2267jejoaNatW+c64lVJIO3+pEZkZKQmJiZm2PpOnz7NokWLmDp1Kp9++ilnz57lrrvuom3btrRt25b8+fNnWBZjTOC4dOkSixYt4r333mPu3LmcO3eO++67j+bNm/PEE09QunTpDM0jImtUNTLZ96xQ+O6nn37is88+Y968eSxevJgzZ85QsGBBWrRowdNPP01kZKRdwWSMSbWjR48yffp0Jk+ezDfffAPAPffcQ+PGjalfvz6VK1cma9b07ezbCsV1+uWXX1i+fDlffvklX375Jdu2bQOgWLFiNGzYkMcff5zq1auTLVu2NF2vMSb07N+/n7lz5zJnzhyWL1/OpUuXyJs3LzVq1KB27dpUr16de++9l7CwsDRdrxUKH5w8eZL169ezZs0a1qxZQ2Ji4h+XtOXJk4dq1apRo0YN6tWrR9myZW3PwRiTbo4cOcKyZctYvHgxS5YsYffu3QDkzp2bBx54gKpVq1KlShUiIyMpXLjwdX0fWaFIhTVr1tCqVSt27Njxx2VrhQoVolKlSn8Uh4oVK6b77p8xxlzNnj17WLlyJatWrWL16tWsW7fuj65C8ufPT+vWrYmLi/Nr2SkVCqffeiJSFxgFhAETVfXNK97PAUwGKgFHgKdUdW96ZClUqBB33303rVq1omLFilSsWJFChQqlx6qMMcYvJUqUoESJErRq1QqA33//nXXr1rF27Vo2bNhAkSLpMyCosz0KEQkDdgC1gQPAt0ALVd2apM3zQHlV7SgizYEnVPWplJab0Vc9GWNMMEhpj8LlfRSVgZ2qultVzwEzgIZXtGkIvO99Pht4WOykgDHGZCiXhaIwsD/J6wPeecm2UdULwHHgTzcmiEh7EUkUkcTDhw+nU1xjjAlNQXFntqpOUNVIVY20DvaMMSZtuSwUPwJFk7wu4p2XbBsRyQrciOektjHGmAzislB8C5QSkRIikh1oDsy7os084B/e508CSzXYruc1xphMztnlsap6QUSigIV4Lo99V1W3iEg/IFFV5wHvAB+IyE7gKJ5iYowxJgM5vY9CVecD86+Y1yfJ8zNA04zOZYwx5r+C4mS2McaY9BN0XXiIyGFg33UsogDwnzSK41KwbAfYtmRGwbIdYNtyWTFVTfay0aArFNdLRBKvdndiIAmW7QDblswoWLYDbFtSww49GWOMSZEVCmOMMSmyQvFnE1wHSCPBsh1g25IZBct2gG3LNdk5CmOMMSmyPQpjjDEpskJhjDEmRVYoABF5Q0Q2ish6EVkkIrd554uIxIvITu/7FV1nvRYRGSoi33nzzhGRiCTvveTdlu0iUsdhzFQRkaYiskVELolI5BXvBdq21PVm3SkivVzn8YWIvCsih0Rkc5J5+URksYh87328yWXG1BKRoiKyTES2ev+2or3zA2p7RCSniHwjIhu82/G6d34JEVnt/Tub6e1H7/qpashPQN4kz7sAb3mf1wc+BwSoCqx2nTUV2/IokNX7fDAw2Pu8DLAByAGUAHYBYa7zXmNb7gFKA18CkUnmB9S24OnLbBdwB5Ddm72M61w+5K8OVAQ2J5k3BOjlfd7r8t9ZZp+AQkBF7/NwPKNslgm07fF+J+XxPs8GrPZ+R80CmnvnvwV0Sov12R4FoKonkrzMDVw+w98QmKweq4AIEcnUA2mr6iL1DPIEsApP9+3g2ZYZqnpWVfcAO/GMMphpqeo2Vd2ezFuBti2pGc0x01LVr/B0yplU0tEn3wcaZWQmf6nqT6q61vv8JLANzwBpAbU93u+kU96X2byTArXwjAYKabgdVii8RGSAiOwHWgGXOyZMzSh8mdmzePaIIPC3JalA25ZAy5sat6jqT97nPwO3uAzjDxEpDtyP59d4wG2PiISJyHrgELAYz17rsSQ/FNPs7yxkCoWILBGRzclMDQFUtbeqFgWmAlFu06bsWtvibdMbuIBnezKt1GyLydzUc5wjoK6zF5E8wEdA1yuOKATM9qjqRVWtgOeoQWXg7vRal9NuxjOSqj6SyqZT8XR9/hqpG4Uvw11rW0SkNdAAeNj7Rw8Bui1XkSm3JQWBljc1fhGRQqr6k/dw7CHXgVJLRLLhKRJTVfVj7+yA3R5VPSYiy4C/4Dk8ntW7V5Fmf2chs0eREhEpleRlQ+A77/N5wN+9Vz9VBY4n2T3NlESkLtATeFxVTyd5ax7QXERyiEgJoBTwjYuMaSDQtiU1ozkGmqSjT/4D+MRhllQTEcEzINo2VR2e5K2A2h4RKXj5ikYRyQXUxnO+ZRme0UAhLbfD9dn7zDDh+XWxGdgIfAoU1v9eWZCA59jfJpJceZNZJzwndvcD673TW0ne6+3dlu1APddZU7EtT+A5znoW+AVYGMDbUh/PFTa7gN6u8/iYfTrwE3De+/+jLZAf+AL4HlgC5HOdM5Xb8n94DittTPJvpH6gbQ9QHljn3Y7NQB/v/Dvw/GjaCXwI5EiL9VkXHsYYY1Jkh56MMcakyAqFMcaYFFmhMMYYkyIrFMYYY1JkhcIYY0yKrFCYgCYivb29Z17u/beKd35XEbkhSbv5SXvSvc51nrp2K5+XuVdECnifr0zr5V+xrscv92ArIo1EpEx6rs8EPrs81gQsEfkLMByooapnvV+02VX1oIjsxXPfy3/SYb2nVDVPGi9zL+mU94r1XL5r9/LrScA/VXX21T9lQp3tUZhAVgj4j6qeBVDV/3iLRBfgNmCZt2uDP36xi0hx8YzXMUlEdojIVBF5RES+9o5FUNnbvq+IdL+8Im//U8WTrlxEaojIP5O8HuPtPgURedM75sFGEYm7MriI5BfP2CdbRGQinps7L793yvtYSES+8u4pbRaRapffF5ER3s9+ISIFvfPbici34hmj4KPLe1TebX1LRFYDQ0SktTfrg8DjwFDvOkqKyNokOUolfW1ClxUKE8gWAUW9X/hjReSvAKoaDxwEaqpqzWQ+dycwDE8nancDLfHcsdsdePl6Q4lIfjx3ld+rquWB/sk0ew1Yoar3AnOA25Np0xLP3egVgPvw3EUMnq7wE72fXe5dFsDHqvqAqt6HpzuHtkmWVQR4UFVjL89Q1ZV4uq7ooaoVVHUXcFxEKnibtAHe82XbTXCyQmEClnr6468EtAcOAzMv/6K/hj2quklVLwFbgC/Ucwx2E1A8DaIdB84A74hIY+B0Mm2qA1MAVPUz4Ndk2nwLtBGRvkA59YyfAHAJmOl9PgVPkQMoKyL/EpFNeLrLvzfJsj5U1YupyD7Ru84w4ClgWio+Y4KcFQoT0NTT1fKXqvoanu7hm6TiY2eTPL+U5PUl/tuj8gX+999HzmSWk2wb7zmAyngGkGkALEhFpj9Rz4BB1fH0ADpJRP5+tabex0lAlKqWA16/IvNvqVztR0A9PLnXqOoRX3Ob4GOFwgQsESl9Rc+/FYB93ucn8Qx16a+9eIb/RDxjpZdIps0+oIy3F9sI4GFv+zzAjao6H4jBc9joSl/hObSEiNQD/jRGs4gUA35R1bfx/NK/PGZ7Fv7bQ2hLYIX3eTjwk7cb7Vap3M7/+e+kqmeAhcA47LCT8QqZ8ShMUMoDjPZ+SV/A02Nme+97E4AFInLwKucpruUjPF3Mb8EzAtqOKxuo6n4RmYWn9849eHrzBM8X7ycikhPPSerYKz+L5xf/dO/yVwI/JNOmBtBDRM4Dp4DLexS/AZVF5BU84yY85Z3/qjfrYe9jagrlDOBt7wUAT3rPU0zFc45lUSo+b0KAXR5rTIBJj8tzr1h+dzx7RK+m1zpMYLE9CmPMH0RkDlASqOU6i8k8bI/CGGNMiuxktjHGmBRZoTDGGJMiKxTGGGNSZIXCGGNMiqxQGGOMSdH/A17X0PVc29nVAAAAAElFTkSuQmCC\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"x1vec = np.linspace(-30, 30, 300)\n", | |
"x2vec = np.linspace(-30, 30, 300)\n", | |
"[X1, X2] = np.meshgrid(x1vec, x2vec)\n", | |
"\n", | |
"P = stats.norm.pdf(X1, s1, sig1) * stats.norm.pdf(X2, s2, sig2)\n", | |
"\n", | |
"P = P / np.sum(P) * 100\n", | |
"D = np.vectorize(d)(x1vec[:, np.newaxis], x2vec)\n", | |
"\n", | |
"p_report_same = [\n", | |
" np.sum([P[i][j] if D[i][j] else 0 for j in range(300)]) for i in range(300)\n", | |
"]\n", | |
"\n", | |
"fig1 = plt.figure(figsize=(6, 4), dpi=72)\n", | |
"plt.plot(s2, p_report_same, 'k')\n", | |
"plt.xlabel(r'Stimulus disparity')\n", | |
"plt.ylabel(r'Proportion reports $\\hat{C}=1$')\n", | |
"plt.show()\n", | |
"fig1.savefig('ten6_2_1.eps', format='eps')" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"id": "canadian-easter", | |
"metadata": {}, | |
"source": [ | |
"### Method 3: Monte Carlo simulation" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 8, | |
"id": "laughing-repeat", | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8fElEQVR4nO3dd3gVZfr/8fedQGgGAxJARGmhBekQQocA0hRwdREsoNhAxfazra6yuKv7ZXddC6KiYkMQhZWiIBYMUmISAlJM6AiBiHTppN6/P85JTCCE5HCSyTm5X9d1rpzzzGTmMxpyZ+Z55hlRVYwxxpjzCXA6gDHGmNLNCoUxxpgCWaEwxhhTICsUxhhjCmSFwhhjTIHKOR3A22rUqKH169d3OoYxxviU1atXH1TV0PyW+V2hqF+/PgkJCU7HMMYYnyIiu863zC49GWOMKZAVCmOMMQWyQmGMMaZAViiMMcYUyAqFMcaYAlmhMMYYUyArFMYYYwpkhcKYYrBixQo++ugj0tPTnY5izEWzQmGMF6WlpfHkk0/So0cPRo8eTbNmzZg+fTqZmZlORzPGY353Z7YxJWn27NnExsZSs2ZNLrvsMt544w1++ukn7rnnHgYOHMjEiRMZNWoUb7zxBj/88ANBQUFORzamyKxQGJOPrKws7r77burUqcP9999P7dq1z1lnz549jBo1irS0NLKysgCoUaMG8+bNY+jQoQAMGTKEd955h7Fjx/Lyyy/z5JNP5rs/VWXLli00btyYgAA70TeljKr61at9+/ZqzMX65ptvFFBAg4KC9I477tBdu3blWefuu+/W8uXL644dO/TEiRO6fft2PXbsWL7bGzZsmFauXFl37tyZ7/JZs2YpoPXr19cJEybo9u3bvX5MxhQESNDz/F61P12Myce0adOoXr06GzZs4K677mLWrFn06dOHgwcPArB582bee+89xo0bR4MGDahSpQoNGzYkODg43+29+uqrADz00EP5Lv/888+pUaMGjRs35vnnn6dx48ZMnz69eA7OmKI6XwXx1ZedUZiLdfDgQQ0KCtIHH3wwpy0mJkYrVKigXbt21dOnT+uNN96ol1xyie7bt6/Q2500aZICumDBgjztaWlpWrVqVb3zzjtVVTU5OVn79OmjAQEBOmfOnJz1MjIydPny5ZqWlnaRR2jMuSjgjMLxX+zeflmhMBfrlVdeUUDXrVuXp/3TTz9VQHv27KmATpgwoUjbTU1N1fDwcK1fv76mpqbmtEdHRyugn3/+eU7biRMntEuXLlq+fHmdN2+evvvuuxoWFqaAvvDCCxd1fMbkxwqFMYWUlZWlV199tXbs2DHf5f/85z8V0NDQ0PP2RxTkyy+/VEBnzpyZ0/bYY49pUFDQOds7cuSItmvXLqevpG3bttqmTRutU6eOnVUYryu1hQIYAGwGtgFP5bP8ZWCt+7UF+P1C27RCYS5GXFycAjp16tR8l2dlZel//vMfXbx4sUfbz8zM1CZNmminTp1y2po3b679+vXLd/0DBw7ogw8+qF999ZVmZWXpwoULzyk0xnhDqSwUQCCwHWgIBAHrgPAC1h8PvHeh7VqhMBfjnnvu0cqVK+vRo0eLbR+TJ09WQGNjY3X79u0K6CuvvFKo782v0BjjDQUVCidHPUUA21R1h6qmAbOAoQWsPxL4pESSmTLp9OnTfPLJJwwfPpyqVasW235Gjx5N1apVefXVV1m4cCEAgwcPLtT3BgQEMH78eOLi4oiNjS22jMbk5mShuALYnevzHnfbOUSkHtAA+P48y+8RkQQRSThw4IDXg5qyISEhgePHj3P99dcX636Cg4MZM2YMs2fP5oMPPqBp06aEhYUV+vtzFxpjSoKv3EcxApijqvlOmKOqb6tqB1XtEBoaWsLRjL+Ii4sDIDIystj3NX78eDIzM1mzZk2hzyayBQcHc+eddzJnzhxeffVV7r33Xnr06MH7779fTGlNWedkoUgBrsz1ua67LT8jsMtOppjFxsbSsGFDatasWez7atiwIddddx1Q+MtOuT3wwAOoKg8//DCzZ89m7969jBkzho8//tjbUY1xdK6nVUBjEWmAq0CMAG4+eyURaQZUA34s2XimrImNjaVnz54ltr9//OMf1K5dm+7duxf5exs2bMiGDRuoVKkS9erVIzU1lcGDB3P77bdTtWpVhgwZUgyJTVnl2BmFqmYADwBfAxuBz1Q1UUSeF5HcP+UjgFnuXnljisWePXtISUkpkctO2Vq2bMnUqVMpX768R9/fvHlz6tevj4hQsWJF5s2bR7t27Rg+fDjLly/3clpTljk6e6yqLgIWndX23Fmf/1aSmUzZlD2CqCQLhbcFBwfz1Vdf0bFjRx577LGcPhdjLpavdGYbU6xiY2OpUKECrVu3djrKRbnssst4+OGHiY+PZ82aNU7HMX7CCoUxuApF+/bt/eLBQqNGjaJSpUq89dZbTkcxfsIKhSnz0tLSWL16tU9fdsotJCSEkSNHMmPGDI4ePep0HOMHrFAYv3XixAnmz59/wfXWr1/PmTNn6NSpUwmkKhnjxo3j1KlTNlzWeIUVCuO3nn32WYYNG8batWsLXM8fOrLP1qFDBzp06MCbb76JDRg0F8sKhfFLv/32W841+u+++67AdePi4rj88su58sorC1zP14wdO5bExERWrFjhdBTj46xQGL/0n//8h7S0NGrXrs2SJUsKXDc2NpbIyEhEpITSlYwRI0Zw6aWXMnbsWKKjo52OY3yYFQrjd/bv38+bb77JzTffzA033MCyZctIS0s777rbtm3zq8tO2apUqcLMmTM5ceIEUVFRDBkyhK1btzody/ggKxTG77z00kucPn2aZ555hj59+nDq1Knz3nyWPZHewIEDSzJiiRk0aBCbNm3in//8J0uXLiUyMpJ169Y5Hcv4GCsUxq8cPHiQKVOmMGLECJo1a0avXr0ICAjI9/JTeno6kydPJioqipYtWzqQtmRUqlSJp556irVr11KlShX69OnDhg0bnI5lfIgVCuNXPvjgA06ePMkzzzwDQLVq1WjXrl2+hWLOnDmkpKTwyCOPlHRMRzRs2JDo6GgqVqxInz59SExMdDqS8RFWKIxfWbZsGU2aNKFFixY5bX369CE2NpYTJ07ktKkqL7/8Mk2aNGHQoEFORHVEo0aNiI6Opnz58vz5z3+2obOmUKxQGJ/0ww8/8Prrr+dpU1ViYmLo0qVLnvY+ffqQkZGRZ0bVmJgYVq1axcMPP0xAQNn6Z9C4cWMmTJjAxo0bWb9+vdNxjA8oW/9CjN94/vnnefjhh/NMUbF161YOHTpE165d86zbtWtXgoKC8lx+evnll6lWrRqjRo0qscylyfXXX09gYCCfffaZ01GMD7BCYXxOamoqMTExZGZm5vnlv3LlSoBzzigqV65Mly5dWLJkCb/88gsTJ05k7ty53HvvvVSpUqVEs5cWoaGhREVF8dlnn9nlJ3NBViiMz4mLi+PMmTMAfP311zntMTExhISE0KxZs3O+p0+fPqxdu5aGDRsyceJE+vbty6OPPlpimUuj4cOHs23btgKnOElPTy+5QKbUskJhfE50dDQiQu/evVm8eHHOX8QxMTF07tw53z6Hm2++maioKF544QV27tzJ119/TWhoaElHL1UudPnpzTffpFq1aixatCjf5abssEJhfE50dDRt27ZlxIgRJCcns2nTJo4cOUJSUtI5/RPZGjZsyJIlS3j66ae56qqrSjhx6XTZZZfRt2/ffC8/paen8+KLL3Ly5EmGDh3KJ5984lBKUxpYoTA+5cyZM8TGxtKrVy/69+8PwOLFi/nxxx+Bc/snTMGGDx/Ojh07znka3uzZs9mzZw8zZ86ka9eu3HLLLbzxxhsOpTROs0JhfMqPP/5IamoqvXv3pl69ejRv3pzFixcTExNDYGAgERERTkf0KcOGDaNcuXJ5Lj+pKi+99BLNmjXjpptu4quvvmLw4MHcf//97Nixw8G0ximOFgoRGSAim0Vkm4g8dZ51hotIkogkisjMks5oSpfo6GgCAgLo3r07AAMGDGDZsmUsWbKENm3alNlRTJ6qXr06/fr144MPPmDTpk2A6x6VNWvW8OijjxIQEEClSpV44YUXAHLO3EzZ4lihEJFAYAowEAgHRopI+FnrNAb+AnRV1RbAwyWd05Qu0dHRtG/fnksvvRRwFYrsy1F22ckz//d//wdAt27diI2N5aWXXiI0NJRbb701Z53w8HAqV65MfHy8UzGNg5w8o4gAtqnqDlVNA2YBQ89a525giqoeAVDV/SWc0ZQi2bPA9urVK6etR48eVKpUCbD+CU+1atUqZ2hxVFQUX375Jffdd1/Of1eAcuXK0b59eysUZZSTheIKYHeuz3vcbbk1AZqIyEoRiRWRAfltSETuEZEEEUk4cOBAMcU1TouJiSE9PZ3evXvntFWsWDGncJxvxJO5sEaNGrFy5UqaN29OlSpVuO+++85ZJyIigp9++um8z/Yw/strhUJE7vDWtnIpBzQGegEjgXdEJOTslVT1bVXtoKodyvrYeH8WHR1NYGAg3bp1y9P+0EMPcc899/jdo0xLWq1atYiJiWHz5s3UrFnznOURERGkpqbaFOVlkDfPKCYWcf0UIPe/7Lruttz2AAtUNV1VfwG24CocpoxRVRYsWEDnzp0JDg7Os6x///5MnTrVoWT+pUKFClxxxdkn9i7ZI8rs8lPZU6RCISLrz/PaANQq4r5XAY1FpIGIBAEjgAVnrTMP19kEIlID16UoG59XBq1bt46ff/6Zm2++2ekoZVa9evUIDQ21QlEGlSvi+rWA/sCRs9oFiCnKhlQ1Q0QeAL4GAoH3VDVRRJ4HElR1gXvZNSKSBGQCj6vqoSJmNn5g+vTplC9fnuHDhzsdpcwSESIiIli1apXTUUwJK2qh+BK4RFXXnr1ARJYWdeequghYdFbbc7neK/Co+2XKqIyMDGbOnMngwYO57LLLnI5TpnXs2JFFixZx/Pjxcy4BGv9VpEtPqnqnqq44zzK7JmCKxZIlS/jtt9+47bbbnI5S5kVERKCqrF692ukopgTZFB6m1Js+fTohISEMHjzY6ShlXseOHQHr0C5rrFCYUu3EiRPMnTuX4cOHU6FCBafjlHk1atSgYcOGVijKGI8LhYhE5f5qzMVYvHgxsbGx57R//vnnnDp1yi47lSIRERFWKMqYizmj+M9ZX43x2MiRI+natSvPPfccGRkZqCrz58/n2WefpUGDBnbXdSkSERHB7t27SU5OdjqKKSHeuPQkXtiGKcOOHj3K77//zlVXXcXf//53evbsSd++fRk2bBhVqlTh448/RsR+zEqLQYMGUaFCBW677TabzqOMsD4K47jsv0wnTZrEjBkz2LBhA+vWreP1119n3bp1NtlfKdO0aVOmTZvGsmXLuO+++855Ol625cuXc/r06RJOZ4pDUe+jMMbrsgvFVVddRWRkJNdccw0VKlSwcfql2C233MKmTZv4xz/+QYsWLXjkkUfyLP/111/p2bMnTz31FC+++KJDKY232BmFcVzuQgGukTVWJEq/iRMncsMNN/DYY4+dM1FgYmIiqspHH31EZmamQwmNt1xMoTjh/nrcG0FM2ZWcnEz58uWpXbu201FMEQQEBPDKK6+QlZXF8uXL8yxLSkoCICUlhW+//daJeMaLPC4Uqtoj91djPJWcnEzdunUJCLATXF9zxRVXEBIScs4ZRVJSEiEhIVSvXp0PPvjAmXDGa6yPwjguOTk557KT8S0iQsuWLc8pFBs3bqRFixa0bduWd955hyNHjlCtWjWHUpqLZX/CGcdZofBtLVu25Oeff84z+ikpKYnw8HDuuOMOUlNTmTVrloMJzcUq6vMonhWR+4srjCl7MjIySElJsULhw1q2bMnRo0fZvdv1ZOMDBw5w6NAhwsPDadu2LS1btuT99993OKW5GAUWChFpISIzcjUNB6bls95dIvIXb4cz/m/v3r1kZmZaofBhLVu2BMi5/JTdkd28eXNEhDvuuINVq1aRmJjoWEZzcS50RvEd8Eyuz+mqeiaf9aYDt3otlSkzzh4aa3xPixYtgD8KxcaNGwEIDw8HXPdciAjz5s1zJJ+5eBcqFNcAL+T6nCYil5+9kqqmAuneDGbKBisUvi8kJIQrr7wyzxnFJZdcQt26dQGoWbMmNWrUYM+ePU7GNBehwEKhqhtU9ZZcTS8B80WkXu71RKQmkP99/MYUILtQXHnllQ4nMRcj98injRs35lx2ylarVi1+++03p+KZi1Sk4bGqOltEKgOrRSQWWIur2PwZ+JvX0xm/l5ycTLVq1exObB/XsmVLvv32W9LT00lKSqJfv355lteuXdsKhQ8r8vBYVf0QaAB8CpQHzgA3q+qMAr/RmHzY0Fj/0LJlS9LT04mPj+fXX3+lefPmeZZbofBtHt1wp6rHcXVgG3NRkpOTqVev3oVXNKVa9sin2bNnA390ZGfLLhSqalPG+yBHb7gTkQEisllEtonIU/ksv11EDojIWvfrLidymuJjZxT+oVmzZpQrVy6nUOR3RnHmzBmOHTvmRDxzkRwrFCISCEwBBgLhwEgRCc9n1U9VtY379W6JhjTF6tixYzkPLDK+LSgoiKZNm/Lrr79SoUIFGjRokGd5rVq1ANi3b58T8cxFcvKMIgLYpqo7VDUNmAUMdTCP8aKTJ09ecJ3sO3mtUPiHq6++GnCdXQQGBuZZlj0zsPVT+CaPCoWI/FlEgt3v/yoin4tIuyJu5gpgd67Pe9xtZ7tBRNaLyBwRyXcMpYjcIyIJIpJw4MCBIsYw3rZ161ZCQkL44YcfClwve2is9VH4h+x+irMvO4EVCl/n6RnFs6p6XES6AX1xTevxpvdi5fgCqK+qrYBvgQ/zW0lV31bVDqraITQ0tBhimKKIi4sjIyODRYsWFbjerl27ADuj8BfZheLsjmywQuHrPC0UWe6vg4G3VXUhEFTEbaQAuc8Q6rrbcqjqIfdd3wDvAu09yGpKWPZcP8uWLStwveTkZMqVK2cPLPITnTp1onr16vTs2fOcZdWrV6dcuXJWKHyUp8+j2CMiU4F+wCQRqQAUdczbKqCxiDTAVSBGADfnXkFELlfVve6PQ4CNHuY1JSh78reEhAROnjxJlSpV8l0v+4FFZ1/PNr6pVq1aHDp0KN9lAQEB1KxZ0zqzfZSnZxRJwNdAf1X9HagOFOnBuKqaATzg3s5G4DNVTRSR50VkiHu1B0UkUUTWAQ8Ct3uY15SgpKQkLrvsMjIyMoiNjT3vejY0tmzJ76a7P/3pT0ycONGhRKawPC0UfVX1c1XdCuD+q798UTeiqotUtYmqNlLVF9xtz6nqAvf7v6hqC1Vtraq9VXWTh3lNCTl9+jQ7duxg1KhRBAQEFHj5yQpF2XJ2ocjMzGThwoV8+eWXDqYyhVHUBxeNE5ENQDP3SKTs1y/A+uKJaHzJ5s2bycrKonPnzrRp0+a8heKrr74iOTmZJk2alHBC45SzC8XOnTtJS0vj559/JiMjw8Fk5kKKekYxE7gOmO/+mv1qr6r2PAqT05EdHh5Ojx49iI2NJTU1Nc8669atY/jw4bRp04ZHHnnEiZjGAbVr12bfvn1kZbnGwmzevBmAM2fOsG3bNiejmQsoUqFQ1aPALuAJVd2V63W4eOIZX5OYmEi5cuVo3LgxPXr04MyZMyQkJOQsT0lJYfDgwYSEhPDll19yySWXOJjWlKTatWuTmZmZ0+G9adMfV5LXr7cLEqWZJ7PHKlDwAHlTZiUlJdG4cWOCgoLo1q0b8Mcw2WPHjnHttddy9OhRvvzyS+rUqeNkVFPCzp7GY/PmzYSEhBAYGGiFopTztDN7jYh09GoS4xcSExNzbrgKDQ0lPDycZcuWcfr0aYYMGcLPP//M7Nmzad26tcNJTUk7+6a7TZs20aJFC5o2bWqFopTztFB0An4Uke3uzuwNImL/p8u4M2fOsH379jx35vbo0YOVK1dy0003sWzZMj766CMGDBjgYErjlPwKRbNmzWjdurUVilLO00LRH2gEROHqzL7W/dWUYdkjnlq0aJHT1qNHD44fP84XX3zBlClTGDlypIMJjZNyF4ojR46wf/9+mjVrRqtWrdi1axdHjx51OKE5H48KharuAkL4Y9RTiLvNlGG5Rzxli4qKonbt2rz44ouMGzfOqWimFAgODqZSpUr89ttvOSOemjZtSqtWrQBynrltSh+PpvAQkYeAu4HP3U0fi8jbqjrZa8mMz0lMTCQwMDDPvRG1atUiJSWFgABHn5FlSgERoVatWuzbty9nxFOzZs2oVKkS4Bo2nT0AwpQuns71dCfQSVVPAojIJOBHwApFGZaUlERYWBgVKlTI025FwmTLvulu8+bNlC9fngYNGhAYGEi1atWsn6IU8/RfsJB3bqdMij4poPEzSUlJ+U4xbUy27EKxadMmwsLCKFeuHCJCq1atrFCUYp4WiveBOBH5m4hMBGJxPZPClFGpqals27YtT0e2MWfLfUbRrFmznPZWrVqxYcOGnLu2TeniaWf2f4E7gMPAQeAOVX3Fi7mMj9myZQuZmZl2RmEKVLt2bQ4ePMjWrVtp2rRpTnurVq04efIkv/zyi4PpzPl4+ijUikAvoLf7ay93mymjFi5cCEBERITDSUxpln13dkZGxjlnFGBTeZRWnl56+ghoAbwGvA6EA9O9Fcr4lqysLN5991169OhBo0aNnI5jSrHcTzPMXSiuvvpqROSChSI+Pp7XXnut2PKZ/Hk66ulqVc19jSFaRJK8Ecj4nujoaLZv324PoDEXlLtQ5L70VLlyZRo3bsz7779PrVq1GDFiBCEhIed8//jx41m1ahW33XYb1apVK4nIhoub6yky+4OIdAISCljf+LF33nmHatWqccMNNzgdxZRy2YWiVq1a5xSCV199lapVqzJu3Dguv/xynn/++TzL4+PjiY+PR1VZuXJlSUU2eF4o2gMxIrJTRHbiuoeio835VPYcOHCAzz//nFGjRlGxonVTmYJl91HkPpvINmDAANatW8eqVavo378/EyZMID4+Pmf55MmTCQ4OJigoqMAnJxrv8/TSk83qZgD46KOPSE9P5+6773Y6ivEBlSpVombNmrRs2TLf5SJChw4dmD59Ok2aNOHhhx9m5cqV7N+/n08//ZSxY8eyZs0ali9fXsLJyzZPzyiSge7AaPccTwrUyn6QkdfSmVJNVXnnnXfo2rWr3T9hCm3p0qXnXFY6W3BwMC+++CI//vgjn3zyCW+//Tbp6ek88MADdO/enYSEBE6ePFlCiY2nheINoDOQPRXocWBKUTciIgNEZLOIbBORpwpY7wYRURHp4FlcUxxWrFjB5s2b7WzCFEnz5s2pXr36BdcbPXo07du358knn+TNN9+kf//+NGnShB49epCRkUFcXFwJpDVwEc+jUNX7gTMAqnoECCrKBkQkEFdxGYhreO1IETnnbi0RCQYeAuynopSZM2cOFStW5MYbb3Q6ivFDAQEBvPrqq+zZs4e9e/cyfvx4ALp06YKI2OWnEuRpoUh3/6JXABEJBYp6730EsE1Vd6hqGjALGJrPen8HJuEuSqZ0UFXmzZvHNddcQ5UqVZyOY/xU165dGT16NFdffTUDBw4E4NJLL6VNmzbWoV2CPC0UrwFzgZoi8gKwAnixiNu4Atid6/Med1sOEWkHXKmqCwvakIjcIyIJIpJw4MCBIsYwnli3bh3JyckMHZpfbTfGe9577z1++umnPLMQd+/endjYWNLT0x1MVnYUuVCIiADLgCeAfwJ7gWGqOtubwUQkAPgv8P8utK6qvq2qHVS1Q2hoqDdjmPOYN28eIsK1117rdBTj5wICAihXLu8Aze7du3Pq1CnWrFnjUKqypcjDY1VVRWSRqrYENl3EvlOAK3N9rutuyxYMXA0sddUmagMLRGSIqtrNfQ6bP38+Xbt2pWbNmk5HMWVQ9+7dAVi2bBmdOnVyOI3/u5g7szte5L5XAY1FpIGIBAEjgAXZC1X1qKrWUNX6qlof11TmViRKgV27drF27Vq77GQcU6tWLZo0acKyZctISkpi8uTJPP/88zZNeTHx9Ia7TsAtIrILOInroUWqqq0KuwFVzRCRB4CvgUDgPVVNFJHngQRVXVDwFkxJ2b17N9988w233347gYGBzJ8/H8AKhXFUjx49ePfdd/nyyy9z2m644Qa7p6cYeFoo+ntj56q6CFh0Vttz51m3lzf2aYru8ccf59NPP2XevHnMnDmT+fPnEx4eTuPGjZ2OZsqw+++/H1UlMjKSGjVqcP3117Nx40YrFMXAo0Jhd1+XHYcPH2bu3Lm0bduWr776iq5du5KUlMQTTzzhdDRTxrVp04Z3330XgFOnTiEiJCXZJNbFwZ56bwo0c+ZM0tLSeP/991m4cCG7du0iMzOTYcOGOR3NmByVK1emQYMGViiKiaeXnkwZ8d5779GuXTtat25N69atiY2NZenSpXTseLFjGYzxrvDwcCsUxcQKhTmvtWvX8tNPPzF58uSctubNm9O8eXMHUxmTv/DwcL755hsyMjLOue/CXByP/muKSAXgBqB+7m2oasFTQhqf8v777xMUFMTNN9/sdBRjLig8PJy0tDR27NhBkyZNnI7jVzzto5iPa16mDFzDY7Nfxk+kpqby8ccfc/311xdqpk9jnBYe7ppT1C4/eZ+n52d1VdUeXuTHFixYwOHDhxkzZozTUYwplGbNmgGuQmGDLbzL0zOKGBHJ/xFVxi/MnTuXWrVq0adPH6ejGFMowcHBXHXVVXZGUQw8PaPoBtwuIr8AqXhwZ7Yp3VauXEnPnj0JDAx0OooxhWYjn4qHp4VioFdTmFJl9+7dJCcn8//+3wUn7jWmVAkPD2fp0qVkZmbaHzle5NGlJ/ed2SHAde5XiN2t7T9WrlwJQLdu3RxOYkzRhIeHc+bMGXbtsl9H3uRRoRCRh4AZQE3362MRGe/NYMY5K1asoEqVKrRqZVcSjW+xkU/Fw9PO7DtxPTf7OfckfpHA3d6LZZy0cuVKOnfubDctGZ+TfTOoFQrv8rRQCJCZ63Omu834uGPHjrF+/Xq6du3qdBRjiiwkJIQ6deqwceNGp6P4FU//ZHwfiBORue7Pw4BpXklkHBUbG0tWVpYVCuOzbOST93namf1fYAxw2P26Q1Vf8WIu45AVK1YQEBBAZGSk01GM8Uh2oVBVp6P4DY8vQqvqamC1F7OYUmDlypW0bt2a4OBgp6MY45Hw8HBOnDjB7t27ueqqq5yO4xeKdEYhIivcX4+LyLFcr+Micqx4IpqSkp6eTmxsrA2LNT6tTZs2ACQkJDgbxI8UqVCoajf312BVrZrrFayqVYsnovG2xYsXc9dddxEVFUW9evXo3r07y5YtY926dZw6dcr6J4xPa9OmDUFBQcTGxjodxW94eh/FpMK0mdJHVbnrrruYPXs2qampdOvWjZ07d9KzZ0+GDx8OYIXC+LQKFSrQrl07KxRe5Onw2H75tBV5Wg8RGSAim0Vkm4g8lc/ysSKyQUTWisgKEQn3KK3JsX37dlJSUpg0aRIrV65kxowZbNmyhRdffJGDBw/StGlT6tat63RMYy5KZGQkCQkJpKenOx3FLxS1j2KciGwAmonI+lyvX4D1RdxWIDAFV4EJB0bmUwhmqmpLVW0D/Av4b1H2Yc4VHR0NQO/evXPaKlWqxF/+8hd27drFsmXLnIpmjNdERkZy+vRpNmzY4HQUv1DUUU8zga+Bd4E7crUfV9XDRdxWBLBNVXcAiMgsXA9DyhkAraq5O8irADbe7SItXbqU2rVr5/sEsGrVqjmQyBjvyx7eHRsbS7t27RxO4/uK2pl91P2Lvaaq7sr1KmqRALgC2J3r8x53Wx4icr+IbMd1RvGgB/spk1SVhIQEjh49mqctOjqa3r17I2I30hv/ddVVV1G7dm3rp/AST/soVotIR68mOQ9VnaKqjYAngb/mt46I3CMiCSKScODAgZKIVWqlpKTwj3/8g7CwMDp27Mi4ceNylm3ZsoW9e/fSq1cv5wIaUwJEhMjISCsUXuJpoegE/Cgi2919FBtEpEh9FEAKcGWuz3XdbeczC9dUIedQ1bdVtYOqdggNDS1iDP+RmppK+/btefbZZ6lXrx6DBg3is88+Y8+ePYDrshPk7Z8wxl9FRkaydetWDh065HQUn+dpoegPNAKicD2P4lr316JYBTQWkQYiEgSMABbkXkFEGuf6OBjY6mHeMmHz5s3s27eP9957j++//54pU6agqrz++uuAqyO7Tp06hIWFOZzUmOKX3U8RFxfncBLf59iDi1Q1A3gAV+f4RuAzVU0UkedFZIh7tQdEJFFE1gKPAqM9yVtWZE+E1r59ewDq16/Pn/70J95++21OnDjB0qVLrX/ClBkdOnQgICDALj95gaMPLlLVRaraRFUbqeoL7rbnVHWB+/1DqtpCVduoam9VTfQkb1mRlJREQEAATZs2zWl75JFHOHLkCE8//TT79u2z/glTZmQ/fKsohSI9PZ2RI0fyxRdfFGMy3+PppIDZDy46CTl3Zf8ITPZWMFN0iYmJhIWFUaFChZy2zp07ExERweTJrv811j9hypLIyEhmzpxJVlYWAQEX/rt48uTJzJo1i6CgIK67rqhX0/2XPbjIjyQlJeU8CjKbiPDII48AULduXRo2bOhENGMcERkZybFjx/jpp58uuG5KSgoTJkwAYNu2bcUdzad4WiiyH1z0NxGZCMRiDy5yVGpqKlu3bqVFixbnLLvhhhto1KgRgwcPtv4JU6YMHDiQatWq8cADD5CRkVHguo8++igZGRn06dPHCsVZLubBRXfgemjRQezBRY7bunUrmZmZ55xRAJQvX56ffvqJ1157zYFkxjinZs2avPHGG8TGxvKvf/3rvOt9++23fPbZZzz99NP069eP/fv3c+yYPTkhm6ed2RWBXkBv99de7jbjkMREVz9/fmcUAMHBwQQFBZVkJGNKhREjRnDTTTcxYcKEfC9BnTx5kvvvv5+wsDAef/xxGjVqBLgm0DQunl56+ghoAbwGvI5rUr/p3gplii6/EU/GGJcpU6YQGhrKbbfdxpkzZ3LaVZVx48axbds2pk6dSsWKFXPuM7JC8QdPC8XVqnqnqka7X3fjKhzGIYmJiTRq1IiKFe3EzpizXXbZZUybNo3ExET69evHwYMHAXj33XeZPn06EyZMICoqCiDnjML6Kf7gaaFYIyKR2R9EpBNgzx10UH4jnowxfxg4cCCffPIJq1atolOnTsyaNYvx48fTr18//vrXP6aRCw4OplatWlYocvG0ULQHYkRkp4jsxHUPRUcP53wyFyktLe28I56MMX8YMWIES5cu5cSJE4wcOZIaNWowY8YMAgMD86wXFhZmhSIXT2+4G+DVFOaibN26lYyMDDujMKYQIiMjiY+P5+mnn+bhhx8mv4lEw8LC+O677xxIVzp5da6n7JcX85lCyJ7jyQqFMYVTr149ZsyYQceO+T8toVGjRqSkpHD69OkSTlY6OTrXk/GOxMRERIRmzZo5HcUYv5A98mnHjh0OJykdPO2jyJ7r6TlVfQ6IBO72XixTFElJSTRs2JBKlSo5HcUYv5BdKKyfwsXmevIDiYmJ1pFtjBdZocjL087s7Lme5ro/D8PmenJEeno6W7ZsYciQIRde2RhTKNWqVaN69epWKNyKXCjENavcbGAp0M3dfIeqXnh6RuN10dHRZGRknLdTzhjjGRsi+4ciFwpVVRFZpKotgTXFkMkUwYcffkhISAiDBg1yOooxfqVRo0b8+OOPTscoFS7mzmz7E9Zhx44dY+7cuYwYMcKm7jDGy8LCwkhOTiYtLc3pKI7ztFB0AmJFZLuIrLc7sp0xe/ZsTp8+zejR9ihxY7wtLCyMrKwsdu7c6XQUx3namd3fqymMRz788EOaNm1Kp06dnI5ijN/JPfKpSZMm7Nq1i+DgYKpXr+5wspJXpDMKEakoIg8Dj+OaxiPF7sguPnv37mX06NH89a9/ZdasWWzatCln2fbt21m+fDmjR4+2p9YZUwyyC8WcOXO47rrrqF+/PrfffruzoRxS1DOKD4F0YDkwENdzKB7ydOciMgB4FQgE3lXV/ztr+aPAXUAGcAAYU5YK0pw5c/joo48IDAwkM9N128qQIUN45ZVX+OijjxARbrvtNodTGuOfQkNDCQ4O5v3336dGjRq0adOGJUuWkJaWVuYeAlbUQhHuHu2EiEwD4j3dsYgEAlOAfsAeYJWILFDVpFyr/QR0UNVTIjIO+Bdwk6f79DVxcXHUqVOHHTt2sHnzZhYuXMgLL7xAeHg4lSpVok+fPtStW9fpmMb4JRHhjTfe4PTp09x66618/fXXXH/99cTHx9OtW7cLb8CPFLUzOz37jaoW/KTyC4sAtqnqDlVNA2YBQ3Ov4H4o0in3x1igTP1WjIuLo1OnTlSoUIFWrVrxl7/8hU2bNjF06FCOHDnCvffe63REY/zarbfeyt13302lSpXo2bMnAQEBLFmyxOlYJa6ohaK1iBxzv44DrbLfi0hRn0R+BbA71+c97rbzuRP4Kr8FInKPiCSISMKBAweKGKN0OnToENu2bSMiIiJPe926dZk1axb79+/nxhtvdCidMWVPtWrVaNeuHd9//73TUUpckQqFqgaqalX3K1hVy+V6X7W4QorIrUAH4N/nyfW2qnZQ1Q75zS3vi+LjXVf1zjeiyV+O0xhfEhUVxY8//sjJkyedjlKiPL2PwhtSgCtzfa7rbstDRPoCzwBDVDW1hLI5Li4uDhGhQ4cOTkcxxrhFRUWRnp7OypUrnY5SopwsFKuAxiLSQESCgBHAgtwriEhbYCquIrHfgYyOiYuLo0WLFgQHBzsdxRjj1q1bN8qXL1/mLj85VijcneEPAF8DG4HPVDVRRJ4XkeypUP8NXALMFpG1IrLgPJvzK6pKfHy83UhnTClTpUoVIiMjy1yHtpNnFKjqIlVtoqqNVPUFd9tzqrrA/b6vqtZS1Tbul9/Mpb148WLGjh2Lqp6zbPv27Rw+fNgKhTGlUJ8+fVizZg1HjhxxOkqJcbRQlGVTp05l6tSpbNmy5ZxlcXFxwPk7so0xzomKiiIrK4tly5axevVqRo4cybhx45yOVaw8nevJXARVzZm+eOHChTRt2jTP8ri4OKpUqWJPrTOmFOrUqROVK1fm7rvv5sCBAzkzJ9xyyy1+eyOenVE44JdffmHfvn0ALFq06JzlcXFxdOjQgcDAwJKOZoy5gKCgIK677joCAwOZNGkSu3fvpmbNmkycONHpaMXGCoUDYmJiAOjfvz/Lli3j+PHjOctSU1NZu3atXXYyphT75JNPSElJ4YknnuDyyy/niSee4LvvvmPFihU566gqqan+MaLfCoUDYmJiCA4O5sknnyQ9PZ1vv/02Z9natWtJS0uzQmFMKSYiBAT88etz7Nixec4qDh06xODBgwkNDeXDDz/Md9CKL7FC4YCYmBgiIyPp1q0bl156aZ7LT9OnTycgIIDOnTs7mNAYUxRVqlTJOat4/fXXad++PUuWLKFx48bcfvvtDB8+nEOHDjkd02NWKErY8ePH2bBhA507d6Z8+fJcc801LFq0CFVlzZo1vPnmm9x3331cfvnlTkc1xhRB9lnF+PHjUVVWrFhBfHw8kyZNYv78+bRu3ZqNGzc6HdMjVihKWHx8PFlZWXTp0gWAwYMHs3fvXtasWcN9991HjRo1+Pvf/+5wSmNMUVWpUoXXX3+d0aNHs3r1ajp27EhgYCBPPPEEcXFxZGZm0qtXL37++WenoxaZFYoSFhMTg4jk9EEMGDAAgDFjxhAXF8e///1vQkJCHExojPHUn//8Zz744ANq1KiRp71t27YsXbqUcuXK0atXL9atW+dQQs9YoShhMTExhIeH5xSDWrVq0bFjR9avX0+3bt3siXXG+KmmTZvyww8/UKlSJaKionyqz8IKRTE4fvw4v//++zntWVlZxMbG5lx2yjZs2DDKlSvHlClT7PnXxvixsLAw5syZw+HDh/niiy+cjlNoViiKwbBhw6hduzajRo0iNjY2Z2jcpk2b+P33388pFI899hhbt26lVatWTsQ1xpSgiIgIrrzySubNm+d0lEKzKTy87Pjx4/zwww+Eh4czb948pk+fTpMmTejVqxeZmZkA5wx9DQoKon79+g6kNcaUNBFh6NChTJs2jVOnTlG5cmWnI12QnVF42fLly8nMzOTll18mJSWFt956i7CwMGbNmsW0adOoWbMmTZo0cTqmMcZBQ4cO5fTp03luti3NrFB42ffff09QUBBdunQhODiYe++9l4ULF3L48GFWr17NDz/8YP0QxpRxPXv25NJLL2X+/PlORykUu/TkZd9//z1dunShUqVKedoDAwNp166dQ6mMMaVJ+fLlGTx4MF988QWZmZmlfgJQO6PwosOHD7N27Vp69+7tdBRjTCk3dOhQDh48mDNJ6PHjx5k6dSrHjh1zONm5rFB40Q8//ICqEhUV5XQUY0wpN2DAAMqXL8/8+fPZvn07nTt3ZuzYsYwcOTJn4EtpYYXCi77//nsqV65MRESE01GMMaVc1apViYqK4uOPP6Zjx47s3buX8ePHs2jRIp544gmn4+VhfRRe9P3339O9e3eCgoKcjmKM8QHDhg3j66+/5uqrr2b+/Pk0bNgQVeW///0vLVq0YMyYMU5HBBw+oxCRASKyWUS2ichT+SzvISJrRCRDRG50ImNh/fbbbyQlJdllJ2NMod1xxx1MmzaNmJgYGjZsCMDLL79M3759GTt2LGvWrHE4oYtjhUJEAoEpwEAgHBgpIuFnrZYM3A7MLNl0Rbd06VIA68g2xhRahQoVGDNmDMHBwTlt5cqV47PPPqNKlSpMmjTJwXR/cPKMIgLYpqo7VDUNmAUMzb2Cqu5U1fVAlhMBi+L777/n0ksvpW3btk5HMcb4uGrVqjFmzBj+97//kZKS4nQcRwvFFcDuXJ/3uNuKTETuEZEEEUk4cOCAV8IV1tGjR3n88cf54IMP6Nu3L+XKWbePMebi3X///WRlZfHWW285HcU/Rj2p6tuq2kFVO4SGhnp9+/Pnz+eXX37J05aVlcW0adNo0qQJL730ErfeeitvvPGG1/dtjCmbGjZsyODBg3n77bdJTU11NIuThSIFuDLX57rutlJl3bp1DBs2jBYtWjBp0iTS09PZuHEjPXv25K677iIsLIz4+Hjee+89atas6XRcY4wfeeCBB9i/fz+zZ892NIeThWIV0FhEGohIEDACWOBgnnzNmDGDcuXK0a9fP5566inCw8Np3bo1iYmJTJs2jRUrVtChQwenYxpj/FC/fv1o0qQJkydP5sCBAzz77LNcfvnlXHPNNcTHx5dcEFV17AUMArYA24Fn3G3PA0Pc7zvi6rs4CRwCEi+0zfbt26u3ZGRk6BVXXKHXXnutqqrOmzdPGzdurLfeeqvu27fPa/sxxpjzee211xTQihUrKqADBgzQ0NBQBXTIkCE6Y8YMXb58ue7atUvT09M93g+QoOf5vSrqfqiOv+jQoYMmJCR4ZVvR0dFERUXxySefMGLECK9s0xhjiuLYsWP07duXFi1a8MQTT9C8eXNOnDjBq6++yr///W+OHj2as+6wYcOYO3euR/sRkdWqmu/lESsUBbjrrrv49NNP2bdvn088XMQYU7acOXOGnTt3kpycTHJyMnXq1GHQoEEebaugQmFjOc/jzJkzzJkzh+uvv96KhDGmVKpYsSLNmjWjWbNmxbofvxgeWxwWLVrE0aNHufXWW52OYowxjrJC4bZ//34GDRrEK6+8QlJSEjNmzKBWrVo2d5MxpsyzS09uu3fvZseOHTzyyCM5bQ899JDdaW2MKfPst6Bb+/bt2bRpE7t27eLbb79l1apVeYqGMcaUVTbqyRhjTIGjnqyPwhhjTIGsUBhjjCmQFQpjjDEFskJhjDGmQFYojDHGFMgKhTHGmAJZoTDGGFMgKxTGGGMK5Hc33InIAWDXRWyiBnDQS3Gc5C/HAXYspZG/HAfYsWSrp6qh+S3wu0JxsUQk4Xx3J/oSfzkOsGMpjfzlOMCOpTDs0pMxxpgCWaEwxhhTICsU53rb6QBe4i/HAXYspZG/HAfYsVyQ9VEYY4wpkJ1RGGOMKZAVCmOMMQWyQgGIyN9FZL2IrBWRb0SkjrtdROQ1EdnmXt7O6awXIiL/FpFN7rxzRSQk17K/uI9ls4j0dzBmoYjIn0UkUUSyRKTDWct87VgGuLNuE5GnnM5TFCLynojsF5Gfc7VVF5FvRWSr+2s1JzMWlohcKSLRIpLk/tl6yN3uU8cjIhVFJF5E1rmPY6K7vYGIxLl/zj4VkSCv7FBVy/wLqJrr/YPAW+73g4CvAAEigTinsxbiWK4ByrnfTwImud+HA+uACkADYDsQ6HTeCxxLc6ApsBTokKvdp44FCHRnbAgEubOHO52rCPl7AO2An3O1/Qt4yv3+qeyfs9L+Ai4H2rnfBwNb3D9PPnU87t9Jl7jflwfi3L+jPgNGuNvfAsZ5Y392RgGo6rFcH6sA2T38Q4GP1CUWCBGRy0s8YBGo6jeqmuH+GAvUdb8fCsxS1VRV/QXYBkQ4kbGwVHWjqm7OZ5GvHUsEsE1Vd6hqGjAL1zH4BFVdBhw+q3ko8KH7/YfAsJLM5ClV3auqa9zvjwMbgSvwseNx/0464f5Y3v1SIAqY42732nFYoXATkRdEZDdwC/Ccu/kKYHeu1fa423zFGFxnROD7x5Kbrx2Lr+UtjFqqutf9/jeglpNhPCEi9YG2uP4a97njEZFAEVkL7Ae+xXXW+nuuPxS99nNWZgqFiHwnIj/n8xoKoKrPqOqVwAzgAWfTFuxCx+Je5xkgA9fxlFqFORZTuqnrOodPjbMXkUuA/wEPn3VFwWeOR1UzVbUNrqsGEUCz4tpXueLacGmjqn0LueoMYBEwAUgBrsy1rK67zVEXOhYRuR24Fujj/qEHHz2W8yiVx1IAX8tbGPtE5HJV3eu+HLvf6UCFJSLlcRWJGar6ubvZZ49HVX8XkWigM67L4+XcZxVe+zkrM2cUBRGRxrk+DgU2ud8vAEa5Rz9FAkdznZ6WSiIyAHgCGKKqp3ItWgCMEJEKItIAaAzEO5HRC3ztWFYBjd0jUoKAEbiOwZctAEa7348G5juYpdBERIBpwEZV/W+uRT51PCISmj2iUUQqAf1w9bdEAze6V/PecTjde18aXrj+uvgZWA98AVyhf4wsmILr2t8Gco28Ka0vXB27u4G17tdbuZY94z6WzcBAp7MW4liux3WdNRXYB3ztw8cyCNcIm+3AM07nKWL2T4C9QLr7/8edwGXAEmAr8B1Q3emchTyWbrguK63P9W9kkK8dD9AK+Ml9HD8Dz7nbG+L6o2kbMBuo4I392RQexhhjCmSXnowxxhTICoUxxpgCWaEwxhhTICsUxhhjCmSFwhhjTIGsUBifJiLPuGfPzJ79t5O7/WERqZxrvUW5Z9K9yH2euPBaRd7mThGp4X4f4+3tn7WvIdkz2IrIMBEJL879Gd9nw2ONzxKRzsB/gV6qmur+RRukqr+KyE5c970cLIb9nlDVS7y8zZ0UU96z9pN912725w+AL1V1zvm/y5R1dkZhfNnlwEFVTQVQ1YPuIvEgUAeIdk9tkPMXu4jUF9fzOj4QkS0iMkNE+orISvezCCLc6/9NRB7L3pF7/qn6uXcuIr1E5Mtcn193T5+CiPyf+5kH60XkP2cHF5HLxPXsk0QReRfXzZ3Zy064v14uIsvcZ0o/i0j37OUi8rL7e5eISKi7/W4RWSWuZxT8L/uMyn2sb4lIHPAvEbndnbULMAT4t3sfjURkTa4cjXN/NmWXFQrjy74BrnT/wn9DRHoCqOprwK9Ab1Xtnc/3hQEv4ZpErRlwM647dh8Dnr7YUCJyGa67yluoaivgH/msNgFYoaotgLnAVfmsczOuu9HbAK1x3UUMrqnwE9zf+4N7WwCfq2pHVW2NazqHO3Ntqy7QRVUfzW5Q1RhcU1c8rqptVHU7cFRE2rhXuQN4vyjHbvyTFQrjs9Q1H3974B7gAPBp9l/0F/CLqm5Q1SwgEViirmuwG4D6Xoh2FDgDTBORPwGn8lmnB/AxgKouBI7ks84q4A4R+RvQUl3PTwDIAj51v/8YV5EDuFpElovIBlzT5bfIta3ZqppZiOzvuvcZCNwEzCzE9xg/Z4XC+DR1TbW8VFUn4Joe/oZCfFtqrvdZuT5n8ceMyhnk/fdRMZ/t5LuOuw8gAtcDZK4FFhci0znU9cCgHrhmAP1AREadb1X31w+AB1S1JTDxrMwnC7nb/wEDceVeraqHiprb+B8rFMZniUjTs2b+bQPscr8/jutRl57aievxn4jrWekN8llnFxDunsU2BOjjXv8S4FJVXQQ8guuy0dmW4bq0hIgMBM55RrOI1AP2qeo7uP7Sz35mewB/zBB6M7DC/T4Y2OueRvuWQh5nnv9OqnoG+Bp4E7vsZNzKzPMojF+6BJjs/iWdgWvGzHvcy94GFovIr+fpp7iQ/+GaYj4R1xPQtpy9gqruFpHPcM3e+Quu2TzB9Yt3vohUxNVJ/ejZ34vrL/5P3NuPAZLzWacX8LiIpAMngOwzipNAhIj8FddzE25ytz/rznrA/bUwhXIW8I57AMCN7n6KGbj6WL4pxPebMsCGxxrjY4pjeO5Z238M1xnRs8W1D+Nb7IzCGJNDROYCjYAop7OY0sPOKIwxxhTIOrONMcYUyAqFMcaYAlmhMMYYUyArFMYYYwpkhcIYY0yB/j/TcDd4YPVgRQAAAABJRU5ErkJggg==\n", | |
"text/plain": [ | |
"<Figure size 432x288 with 1 Axes>" | |
] | |
}, | |
"metadata": { | |
"needs_background": "light" | |
}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"trials = 300\n", | |
"p_mc = []\n", | |
"\n", | |
"s1mc = 0\n", | |
"\n", | |
"for s2mc in np.linspace(-30, 30, 100):\n", | |
" x1mc = np.random.normal(s1mc, sig1, trials)\n", | |
" x2mc = np.random.normal(s2mc, sig2, trials)\n", | |
" d_mc = np.vectorize(d)(x1mc[:, np.newaxis], x2mc)\n", | |
" p_mc.append(np.mean(d_mc > 0))\n", | |
"\n", | |
"fig2 = plt.figure(figsize=(6, 4), dpi=72)\n", | |
"plt.plot(np.linspace(-30, 30, 100), p_mc, 'k')\n", | |
"plt.xlabel(r'Stimulus disparity')\n", | |
"plt.ylabel(r'Proportion reports $\\hat{C}=1$')\n", | |
"plt.show()" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 9, | |
"id": "mighty-toner", | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"fig2.savefig('ten6_2_2.eps', format='eps')" | |
] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "cognition_env", | |
"language": "python", | |
"name": "cognition_env" | |
}, | |
"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.5" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 5 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment