# Project4_Bayesian_HardyWeinberg.r
# Multinomial Counts: Hardy Weinberg Equilibrium
# 1.0 Trinomial data from Example 8.5.1.A, p. 273 of Rice.
#
# x=(X1,X2,X3) # counts of multinomial cells (1,2,3)
# Erythrocyte antigen blood types of n=1029 Hong Kong population in 1937.
x<-c(342,500,187)
# Two estimators: Multinomial MLE and Binomial(X3) MLE
x.n=sum(x)
x.theta.mle=(2*x[3] + x[2])/(2*sum(x))
x.theta.binomialx3=sqrt(x[3]/sum(x))
print(x.theta.mle)
## [1] 0.4246842
print(x.theta.binomialx3)
## [1] 0.4262978
# 2.1 R functions
# function computing cell probabilities given Hardy-Weinberg theta parameter
fcn.probs.hardyweinberg<-function(theta0){
probs=c((1-theta0)**2, 2*theta0*(1-theta0), theta0**2)
return(probs)
}
# 3. Bayesian inference with uniform prior
dgrid=.001
theta.grid=seq(0,1,dgrid)
# Compute likelihood of x
# First, compute matrix with multinomial probabilities (3 columns) for each theta (row)
probs.grid=t(apply(as.matrix(theta.grid),1, fcn.probs.hardyweinberg))
plot(theta.grid, probs.grid[,1],type="l",col='red', xlab="theta",
ylab="Cell Probability",
main="Hardy Weinberg Cell Probabilities \n 1 (Red), 2 (Green), 3 (Blue)")
lines(theta.grid, probs.grid[,2],type="l",col='green')
lines(theta.grid, probs.grid[,3],type="l",col='blue')
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAgVBMVEUAAAAAADoAAGYAAP8AOjoAOpAAZmYAZpAAZrYA/wA6AAA6OgA6Ojo6ZmY6kLY6kNtmAABmAGZmOpBmZgBmZjpmZmZmkJBmtv+QOgCQkDqQkNuQtpCQ27aQ29uQ2/+2ZgC2/7a2///bkDrbtmbb////AAD/tmb/25D//7b//9v///9lZlrxAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3di3bjOpaYYdWpnIzb3R1Xkpm4ksmU5ihluUrv/4AR7yAJXgAC4N7A/60106dsWSYl4jdvIi8PAICXy9kTAABaEVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwRECl+Xy9XC5ffiz+28X1+aOXl61H/f4+elj9Q5c//mr/eav+9edP+4/++rY1cbddk3BQM8nDNC8xJtc+5fVXDe9u0+EyswuPtU/idGrvb7NH4xQEVJr0AW3GclefJqfDb1x9DgkBrV+gfcVzDqhjQhMF9DnHBFQIAirNCQG9m8XsemQO0cWMCAjodRy81V/lEVCngiYJaP0XjoAKQUClOSGgo0jexiVa//WnB7RdXzYs7W2o+ATU5aVPEtDr9M8bAT0RAZXmhICOdoJ2a3TtJv3d+G8fkQN6NVaXb6NVZ5t9Ae2+etuxTjsWIKBLEzO4bswkUiKg0pwQUPNAUb8OZq7vHAhg3IDeRyudzbSv1N4toO3eDIe/HgS0OARUmjMCaqxmDodk6jHarJweGK5RAzotZjPxy1PrGNB2/Xb/XlACWhwCKs1GQG9m3vovvHX/M/xINe7bgDad6X7CcljI+B318/xbP7jHg7hdPR1Gr/Hte/ud5jH9nsiuE/d5i2bPtjgrb+1/THsze9Lr6HevTe6ugA7Pb5m0bv+rsYa6MrP2923+2I19oPfL4H02vbMZfgx/EQ/shsEyAirNakBvxvh5M770du3H1LAR/sdf3RroqCv1E46Hk7GeWT30y//91j3EfLBxhKUr2Sygw2PMyXuZ7hiwP9vyrHz5YQ3odVqG+zRFi5PrFVBj0syzp/ofWJzZpfdt/tgDAbXM8OgcLxIaAQGVZi2g5jjsRmL9tX8MQ8Qc2H9vB1O/Str/wCRF/bZ+PQj//Hntfqfx4NH5lm2NpwH9p/EYY/L+6zC0p2dLmc82nRWjCF/+3RLQ5vsLm8Ebk+sTUGPSrG/F4swuvm/zx/oH1DbDo69R0AgIqDQrAR0Ph3Y8TFZu5uf1vIwPszf/Pdmvd+8G3b15YL+pOuxwm5zh0zzZNKDrk2d82fpse2Zl9lIt7A7cmtz9+0CHPyTDpM3m1ojifGY337fpC+MTUNsMT1/B2J8IKxABlWYy2hrNGOm3WJvBYuy0HB516/+rGz0v3YOM0TxdF+krfWtGZrfBbNR2OKBiHFqZB7T6RzsL78bkGSuUb4vPNpmV+/QnpwGw7GrsbU3unoAOEzCdNGNW7sYbtDazC+/b9LHe54HaZng4ONi8kZwxGhwBlWY5oMZoMbbJx2crNgOl7eO1/45x4OhuKVH//frHnz/9q90JOtTWHMLDdv0soM2mo1G2m/HlYeLsz7Y8K+0fg/0B3ZzcHQG1nGPaTcCwajrakbA5s7P3bfpY74BaZ/g2vD7Vr6Cf4RFQaVbWQKtRYGyLzgfiZG9n84/6J67j/5pFpx2V7S7Q9ljSj2HTftTdLq+WgL4bkzeUpxu5/WdG7c+2b1YGKwHdnNzdn0QyN7jHJ5y+mb9r2NCfzezK+zZ/rG9ArTPcTDUb7hERUGnWAtqbBbSryOiDQ8auz/7rfSLH2lPpu3OR5qftmMO2ft7pcDc/UG8c1R9dzKmfbPuz7ZuVwUpANyd3b0BHh++sk2bMrX1mZ88/vG/zx/oG1DrDwx5TzhuNhIBKs+dEevN0mtGazOTic8N6Z78zs366+XBqt9X7bb5mJ+iwC3R+QGc63M2oTAPaha/7+sKzjWdlWPsdz8pgOaDbk7szoF0l117lYdqsM2uYvW/zx3oG1D7Do69yDD4GAirNzhPpxwNx4eqdxiDt/vNm60Y3Ev/j+8U86PDHf34b/XtsMtxHa2XDIB+fNNV+feHZds7K+KWxrlxtT+6ugA6NtEzaMDHjfY7TmTUeM3vf5o8NF9C3x3Rzhn2g4RFQadYCOjl5Zr4puFyddhWzHmiWCxY1A/B/vfbfrXeC/nv/750B7Z75QEC3Z6U3Pw/083U4ZLI+ubtOYxqsTsxmQBfft9gBnXydggZHQKVZCWh3oaQ/fy7tS1ve7m23xusftK2zDWtIb8MT/WMyEGc/uDOg3ZfHm/CzZ1sLqG0TfnTKQf8z1Y9sT+7xgC5twk9nds/7FmgT3rqn06z32rX+4IOASrMcUOMwzXpArUdemsO096X1kH6YzT7X0qzT2S9hsTOg3QRZDyIN9s7K+AeMnaBDUTcn91BAbQeR+p24s5lded9mjw1yEGluco1CBENApVkOqDFGFgM6OvQ7OvenObflfy6thvR7/0bHgYZJmax0mT+1HdDRBPXn68+ebfzV7dOY2okeqjCctrM5uYcCunoa02xmV9632WO9A2qf4ecTW86qQjgEVJrFgJpbacboG4+clbPP+3tfWNdT+mBOHz76VHW/eT//JMxqQM0Jqh5jf7bds9Ibn55pfChoc3IPBXT9RPrJzK69b7MXxjugthke7SO2fgINRxFQaTYCanxYzxZQIyLdOuWL+UTLayFdMEfrVbOgvvTfmQ731YCaE/S2+GzWWVn5KOfklbiaD9qa3GMBbV/LhY9yjmd27X2bvzA7Azq/Ooptho2/L83vYB9oaARUmuVN+NnlJ+ZHaR6zm6xNzzVcHETGZ9lrXbS6PYzTE/ytn4VfDOigeYj12SazsnkxEetjrNch8vssvGH6Kk9n6936VePvgmnthdkO6M14GuPrthmendbqeJtmbCKg0iwH1BwO/6W/utt0aJuP+ufoGnA3Y/gu/N7ZoRFjSsbn4hj79LYD+m/Dz3XPb3u2lVmpd95aP5Q4+YMxPtCzPLkHAzrpn7kPYjazK+/b/IXZDmhfynFArTM8qSqfRwqOgEqz63J27/3HnedDe1gpextfhbj58aXjCLPdjOP9i49xCN6NL20H9L0f3MNOOMuzzWbFOLK1fHeSpfvCr0/u0YAuX1B5PrPL79v8sdsB7f9kvEy+bpnh8Ro665/hEVBpVj+J1Iy30YU6LYdf+4MHk8u4X+cPNUxvAHTvRun0iY2p2xvQ6Wch7c9mmZXxoeyFy2J0lZgeIlmb3MMBXbylh21ml963+WN3BNT4SziZ3tkMPx59bzkAHwUBLYjm25Fpnnbki4CWY3H1SqbnCp4xsQQUEhHQciycay3V6CLqyuKPUhDQYtwuylbiZmcxcho4pCGgRehPutHUoNmlpRXFH4UgoEW4q2zQ5FxzPXsfUAwCWoS7xn5Ozg3n1j6Qh4AWodkaVrgKZ37sBhCHgAKAJwIKAJ4IKAB4IqAA4ImAAoAnAgoAnggoAHgioADgiYACgCcCCgCeCCgAeCKgJbitXcju2l6p+NDV4/pbFlVX7txz2Y/29j3Ll9e7Ti6gbNxgzmG6HDlO/fJd2hZf1PfxczhfI3rvFCINAlqA4VaQFtUQf+n+Y2L/SB3u+XZdDYvxS1sLCeknq2G9i3AEu6beuP3lwmPXXtR381ncL7K/awqRCgEtwHVlxa26ccb0psmD3ddvGgJaPc3Wj41zaMvBMFnWSYu2Frpn6pfuozxYf1Hb/HsGdM8UIhkCmr/bWnPufSZtY333ut4Q0PpWRus/NvlNtogMk/WYXBY0bkF3TP3oVuv2Aq6/qMdu87RjCpEOAc3edS059ShuhuPkhvRXp1AZt22/b/1YW6Bq3a1p+3xPgTFZ06tBt7MTayVsc+rbyXnv4mjZzbH1or4Mj/K4T972FCIdApq5jWMvxmicjPV2sO/c3XadhGFtFakJT/PEN/ukmZFo56B/xs0DOIdsTn0zq+/dVNpe160Xtfmeb0C3pxDpENCsGZub1oDW32/XoaZj/T7u1jojoPV/r60gmt2Z/tL5ZLUhN6Yj7uX1t6a+fkQz9fWEzF/XtRfV+Lf3nZq3pxDJENCcGYeL7QG9G+ty1rE+W/EbB7VfmTQDaj6pjXl/+uZZpw+eTdZ4Q7l+gmZKb8003YxJs07o7Iv39p/NN4webU39dO7tq8+rL+psDdR8RZq/ef0EWObGYQoRGwHNWTP63pY2NZvudcN7NaBGioeWGbsHzIDWX955BpR9NcycLKOWxpR1T9+Us90v+r45ocMXm4AO3xj/odg39deL9aFrL+p9SOWegC7PDeeCykBAc1YNtee4XAxoPRS7ta+1TfjRweTRDzS+/N0Y43UCdh7msE6aOVlNTxa3WOv0/KPL+PaE9l+sf/E/jW90Bd0/9dd52qdTb98H+jY8bj2g1rlxen0RGQHN2a9v9VhdDGg9QF/MfwxjvRnIxu66wYvxgPnXmx7s27lnz6M5Wc2vXlzfMs8ofVuYUOsXZ6dG9a/Pvqnv2jbflF55UZvZbX/TdkCtE757CpECAS3AYkDNXYfWsd4O2uGAvHFovmlX9QNdSsebx3umzH6k35ysyZ5Y6xxc+imxTqj1i/fhh9oWdlOxb+rv4+ouTL3tPNDR2QWrAbVO+O4pRAoEtACLAb2aw3L5nG9zR+WtK+Vobco8v3GyDram+bHZ9rk5WeOATlY4+y+Yq6vTCbV+sQ1g86v7Mzsdpn78i5em3vKi9j+xGVD7hO+eQqRAQAuwFNB6rPbrnJaAtt+7Txs17G20Hiqv/7XjRJur+UsWJmtXQLvfZZ1Q6xdH0WyPtTlNfTcl61O/8qJuB9Q+4bunECkQ0AIsBXR8kZHllSXjyMfQh9GTNqP+xfq8i4xj58uTtSug7+YzTifU+sVm8tuWNVPf/Y6dU/9YOIq09aJedh9Esk+40xQiNgJaAK+ADpuI86NF1fg2B/3kjJ59A/xq1mRxsvYEtAuLdULtUz96TXwDaj2PyfKiDo29DxOwGdCFCXebQkRGQAvgENB6QA+Hh2qWkfw2DehoH914K9aue1LLCeGjybIdhTeiOvooqHVC7VM/fk3Mdb1dU9+yvbCrATU/ReURUNvqKE5FQAvgGtDRilKcgLbPaX3QeLIs1+KbBrSbCueAdj/oG1DbxwDWA2qsMhPQHBDQAjgHdLwOOt7CHT2pfRN++0zvNg4rV4ganZY5WVFdD+hkQu1TvxXQfRvItg+ibgR0OHi1END+ORcm3G0KERkBLYB7QMcXs7vaRvKhg0gL5y/Zfnz2YfXlgNon1PrFlYBuTf1tEju3NdD9AV2Y8B1TiHQIaAFWT2Oyj/X2QzBNEseVepiPsJ7GNPowo83SZUBtk9U+2AyJ8RXLjoTpL7Z+cSugy1NvHr63vbArL+ow7ZOAms8zFNY+4dtTiIQIaAEczgPt/2XuBjWPhN+7/zBPpB99bGn7PMXFKxHbJmvc8uGX2QJqnVDrF1cCujH1xlGt8Xq3deqnAW1+un7VpgFtXmvjatP2Cd+eQiREQAuw/lHO0YdmphdPb4bp8EEj4wM8zX9WP9B9Znv+UU775T7bE5gGk0Pxt/HXuvOr6mecXD9pspZmndDlqbcGdN/Uv4yvFrDzRW1n5aWfE/POSeZL+T76TaPPTfFRTkkIaAEWAzoaiZOxbl6Icno6eJOHWQdHAX23PKfxi9YCOg2E5Z5I3Q9NAmqdUOsXNwK6PPXTS3zY966ufBbe8gHZ+QH35QkfTSHORkAL4HU1JmMNc5qw+TVCmwvDmZ9cMlatpgka35PTEoPZZ70tGWrbN91PaJ1Q2xeXA7ox9ZOpme+ImL+o9tk1j0AZU/jlX4eHWOdmfMFRnIuAFmAxoPWaz9KlK9vVIuPqv9PcDStOb2Y1jGe1JWj5EzbWyTJmwdKS2YEW64RavrgY0I2pnzydZUN6/qLa53Z0CL+fvy8/btPHWF92doHKQEALsHxF+pvxjVku2sE/Xt8c96Q/ojK7hKf1ppTmE60EdDRZvW6PwejpLEeqVyZ0+qFKW0A3pn70dPb9kLMX1TTZRzLZ51z9c7wH2DI33FVOEAJatHp8h96bdh8N+JvHGYtRJmsn+VN/X8o60iOgRRvd/jKU0V0jn7/BfWszymTtJH/quSunIAS0bOYN2AOp+9FvYd69WhJhsnaSP/XjKcS5CGjZ6tEYdnNzlI+bX0oiTNZO8qf+vD8umCOghbsH39y8XkbnVPrtrQs/WTvJn/orK6CCENDCVWtLQY9IVMeI+xWk3989OxJ8snaSP/WjKcTZCGjp7oFXaG5hNl9DT9ZO8qc+0BQiDAJavGvQNZpqBSnI1mvYydpJ/tQHm0IEQUABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8CQ7oBcACCZ8ooI/Y0Bnv9oA8hK8UaGfcOLaT/mL+w9H+IMBoFjKAnqb1N+1oQQUQDiqAvr5Olt//uMvp2cgoADC0RTQ39/Hxfz17fnvP3+6PAUBBRCOpoDeZ7mskvru8hQEFEA4mgJ6nW+wP1dCnXaDElAA4SgK6HN18232xZvbNjwBBRCOooA+1zbnm+t3t8NIBBRAOAQUADwpCiib8ABkURRQDiIBkEVTQO2nMc3XSlcQUADhaAqo9UR6t48iEVAA4WgKaFPMsS8/nJ6BgAIIR1VAzUsxcTERAGdTFtAHl7MDIIa+gB5CQAGEQ0ABwBMBBQBPBBQAPCkP6Ofr2omgKW4BBaBcBNTu4yPMBALQ5+vXnQ/MOqBz+2eXggLFKiWgrlwCSkGBMu3uJwFdREGBMn0loAtcZpeCAiVy6CcBXUFBgfK49JOArqGgQGmc+klAV1FQoCxu/SSg6ygoUBLHfhLQDRQUKIdrPwnoFgIKlMK5n5oCarmhx8X1pkges0tBgUI495OAbmIjHiiDez81BbT64PsJAaWgQBHcN+B1BbReB3W6DfyM1+xSUCB/Pv3UFdCqoI73MZ7wm10KCuTOq5/KAlptxf/588DPe84uBQXy5tdPbQF93C6X9wM/7ju7FBTImWc/1QX0uRF/ZBXUe3YpKJCtr779VBfQg6ug/rNLQYFM+fdTX0CPOTC7BBTIkn8+CagDCgrk6EA/Ceh+bMQDGTrSTwLqgIIC2TmyAU9AnVBQIDPH+klAnVBQICsH+0lA3VBQICNH+0lAHVFQIBuH+0lAXX2QUCAPx/tJQJ1RUCALAfpJQD1QUEC/EP0koD4IKKBdkH4SUC8UFFAuSD8JqBc24gHdwvSTgPqhoIBmYTbgCagvCgroFaqfBNQXBQW0CtZPAuqNggI6hesnAfVHQQGNAvaTgB5AQQF9QvaTgB5BQQFtgvaTgB5CQQFdwvaTgB5DQQFNAveTgB5EQQE9QveTgB5FQAEtgveTgB5GQQElgveTgB7GRjygQ/h+EtDjKCigQfgNeAIaAgUF5IvRTwIaAgUFpIvSTwIaBAUFZIvTTwIaBgUFJIvUTwIaCAUF5IrVTwIaCgUFpIrWTwIaDAUFZIrXTwIaDgUFJIrYTwIaEAUF5InZTwIaEgUFpInaTwIaFAUFZInbTwIaFgEFJIncTwIaGAUFBIncTwIaGBvxgByx+0lAQ6OggBSxN+AJaHgUFJAhfj8JaHgUFBDga4J+EtAIKChwuiT9JKAxUFDgZEnySUDjoKDAqRL1k4DGQUGBE6XqJwGNhIICp0nWTwIaCwUFTpKunwQ0GgoKnCJhPwloPB8kFEgvZT8JaEQUFEguaT8JaFQUFEgrbT8JaFwUFEgpcT8JaGQUFEgndT8JaGwUFEgleT8JaHQUFEgjfT8JaHwUFEjhhH4S0AQIKBDfGf0koClQUCC6M/pJQFNgIx6I7ZR+EtAkKCgQ1ykb8AQ0EQoKxHRSPwloIhQUiOesfhLQVCgoEMtp/SSgyVBQII7z+klA06GgQAwn9pOAJkRBgfDO7CcBTYmCAqGd2k8CmhQFBcI6t596A3q/PL24/tTJAaWgQFAn91NbQG/Pav7x1+Px+/ul8e7282cHlIICAZ3dT10BbbP55cfjerl4FfT0gFJQIJjT+6kroF02//x/z//7+ahXSKv10f3ODygFBQI5v5+qAnpvsvn5+mzoW/8lp1VQAQGloEAQAvqpKqDXZrWzKmi/3nl1O5AkIaAUFAhAQj81BfT393a98/kffTVvbVR3EhFQCgocJqKfmgL661u7ud6X9FFtwzvtBJURUAoKHCSjnwT0HBQUOEJIPzUFNJ9N+AcFBY6Q0k9NAe0PIt2Hk5fMlu4hJqAUFPAmpp+qAtqexlR/iFPxaUwtCgr4kdNPVQHtP7/Zn0h/V3gifY+CAj4E9VNVQKvDSM1HOfuU6vsop4GCAu4k9VNXQIeLibQp1XcxkREKCjj6Kqqf2gI6uKm8nN0EBQWcCOun3oD6ERZQCgq4EJZPAno2CgrsJq6fBPRsFBTYSV4/CejpKCiwi8B+ag+oeWW7uYtFwonbiYICO0jsJwEVgIICm0T2M++AzokMKAUFtsjsp/aAupIZ0KqgJBRYJrSfBFQGCgqskNpPAioFBQWWiO0nARWDggJ2cvupMKDX/oC68yfhRQeUggJWgvupLaC3yUlJrg2VHFAKClhI7qeugH6+zk7rdLuesuyAUlBgRnQ/VQW0voyyUcz6oqBO95QTHlAKCkzI7qeqgN5nuaySqvWeSFYUFDAJ76eqgF7nG+zPlVCld+VcQEGBgfR+agpof194k9b7wi+ioEBHfD81BfS5tjnfXL/n8Fn4EQoKNOT3k4DKQ0GBioJ+agpoGZvwDwoKVDT0U1NAiziIVKOggIp+qgqo/TSm+VrpCh0BpaAono5+qgqo9UR6t48iKQkoBUXhlPRTVUCbYo59+eH0DFoCSkFRNC391BVQ81JMGV5MZISColxq+qktoI+ML2c3QUFRKj39VBjQQxQFlIKiUIr6SUAFo6AokaZ+ElDJKCjKo6qfBFQ0CorS6OonAZWNgqIsyvpJQIWjoCiJtn4SUOkoKMqhrp8EVDwKilLo6ycBlY+CogwK+0lAFaCgKIHGfhJQDSgo8qeynwRUBQqK3OnsJwHVgYIib0r7SUCVoKDI2Fet/SSgWlBQZEtvPwmoGhQUmdKbTwKqyAcJRYYUr34+CKgmFBT50d1PAqoKBUVmdOeTgCpDQZEV7f0koMpQUGREfT8JqDbsCEU29PeTgKpDQZGJDPpJQBWioMhBDv0koBpRUOiXRT8JqEoUFNrl0U8CqhMFhW6Z9JOAKkVBoVku/SSgWlFQ6JVNPwmoWhQUWuXTTwKqFwWFThn1k4AqRkGhUU79JKCaUVDok1U/CahqFBTa5NVPAqobBYUumfWTgCrHpUWgiPLLz1sQUOUoKNTIr58EVD8KCh3yyycBzQEFhQY59pOA5oCCQr4s+5kkoL/+x8/Qv8RXngGloBAvz36mCei3y+Ul9K/xk2lAKSiEy7SfyQL69Bb6N3nINaAUFKLl2s9E+0Bvl8Z76F/mKtuAUlAIlm0/0x1EujYJ/eOv0L/PSb4BpaAQK99+pjwK//t709A/TzyklHFAKSiEyrifiU9javeGnndIKeeAUlCIlHM/058H2q6HfvkR+vfuknVAKSgEyrqfqQPaHU06ay0074DywXhIk+HH30dSBvTWH0eqV0PPOK0p84BSUMiSez/TBXRyJtOvb6cckM89oGzGQ5Lc85n6PFBjpfNKQOOgoJAi/36m/CTS6PSl50b8GaczFRBQCgoZst98r6QK6ElH3adKCCgFhQRF9DNRQE//CGeniIBSUJyviHxyPdA8UVCcrJB+JloDHW3AX0/8MGchAaWgOFcp/TwnoOddUaSUgFJQnKmYfp4Q0JPOAG0UE1AKivOU08/YATU+umlgEz6BDz6VhFOUcfi9FTmg3eWXxs67Nn1BAaWgOEVR/Yy+CX8X1c+iAvpgOx7pFZXPUw4inamwgFJQJFZYPwlo5igoUiqtn5xInzsKinSK6ycBzR4FRSrl9TNyQKuD8O/zQ/GcB5oSBUUaBfaTgBaAgiKFEvtJQEvACaGIrqzTP3vK9oH++tbG9+Z3a88yA0pBEVuh/VQW0Gu79nr3XZctNKBsxiOuQvOpLKDXdq3T+HiT40posQGloIio2H6qCujna3M3+Wqfar3qWd0d2W0dtNyAUlBEU24/ExxEsvA8iHRtb4p8azpauTt+sr7ggFJQRFJwPzUF9Lm++dL+7/AEjpe3LzmgFBRRlNxPTQF9Ptl7+79DNG9uT1Z0QCkoIii6n5r2gRoBfem/eCegDjidCYGVevpSR1FAn5vub+3/ElBPFBRBld5PTQF9XNtwXkeb8OwDdUJBEU7p+dQV0Ht71ufna7Mt/xjWSvcioBQUwdBPVZ+Fr077rNdBhyNHV8cz6QkoBUUo9FNXQOsz6etN9rab1SeS3G6wREAfFBRh0M+HsoA2BR152f6h0cQR0AcFRQj0s6JpH2jlOu7n+/ZPjBDQGgXFUfSzpi2g5sn5TsffGwS0welMOKT405c6+gJ6CAFtUVAcQD87BLRUFBS+yGcvWUDb4z/n3c2jRkAHFBR+6OcgUUD9L4EcFgE1UFD4oJ+GNAG9HTp0vuK5XruyTmu7ElS4360fBYUzdn+OJAlotfneds74zxAI6CEUFI7o51iSgF7Ntc6b89nvK9YDOkdAxzgYDxfkcypFQEfXn3O+iHxQBHSCgmI/+jmTIqDdlZBbjpfwDIqAzlBQ7EQ+5who8SgodqGfFon2gb4s/istAmpBQbED/bRJdRR+uOrc/eCZoMPlRDw6TEBtKCg20U+rNBbwh2IAACAASURBVOeB3vvz55+b80f6eZuclOTaUAJq9cGxJKzi8NECRbc1tl0O1PWpCKgdBcUa+rlEU0CrW3qYP1s/u9spUQR0EQXFEvK5SFNA77NcVkl1+mAoAV1GQWFHP5dpupzddV7eZ6GddoMS0BUUFDb0c4WigFpvYcx94QOioJijn2sUBXRyPn7D8ax8ArqKQ0mY4PDRujMC+uu/e+0DJaDxUVCM0M8N6S+o7H0QiU34FCgoBuRzS5qAXoMchecgUgoUFB36uSnZBZVNnpekt5/GNF8rXZs4ArqJgqJBP7eluqDyS7W1/fKYXFzZjfVEere1WQK6AwVFhX7ukOiCytWaY3u85+Z/NTvLefmOH6wnoHtwKAkcPton0fVAq+3sz9e6ds+c+t9VbrIzlYuJxEFBi0c/90l4QeXuNKTbseuBcjm7FCho2cjnTgkD2p2G9PnKPZHko6Alo597pbylR3sp+s9XbumhAAUtF/3cLdFBpHrVsz3pnXsi6cCO0EKx+9NBktOY2gPv7c08HD88FBQBdUBBi0Q/XaQ6kb5K5/N/num8ex39CYSAOqGg5SGfTtJ9lPNtOIDufxrTUQTUDQUtDf10kyag7ZXj2xPhz1sBJaCuKGhZ6KejVJezuzcfGWpXRU9DQF2xI7Qg7P50puiCyiEQUGcUtBj00x0BxRYKWgby6YGAYhMFLQH99JEsoO01Qc87h75GQL1Q0PzRTy/pb+nheAG6sAioH3aE5u0ruz89pQnoLcQF6UMgoJ4oaM7op7dkt/Rot92N/zwDAfVGQbNFPf2luqXHsNZ546OcOlHQTNHPAxJdjclM5pWLiejEZnyO2Ho/JOX1QFtczk4rCpof+nkMAYUDCpoZ8nlQon2gL4v/SouAHkRBs0I/j0p1FH64gsj9zDNBCehRbMZnhH4eluY80Ht//vxzc/7MM+kJ6GEfJDQT7P4MIG5A2wuATrEPVDUKmgX6GQIBhTsKqh/5DIKAwgMF1Y5+hsHl7OCDzXjV2HwPhYDCCwVVjH4GQ0DhiYJqRT7DSRbQdncoF1TOBwXViX4GlCigxtGkUxNKQENiM14hNt+DShPQ0dF4TqTPBgVVh36GlSSgv7/316Gv7u1x3tXsCGhoFFQX8hlYkoDejWhWMT3vnh4ENDQKqgn9DC3R1ZiMHZ/PzXmuxpSRfDbj/7bs7EkLg8338BJdkf7N+OeNK9JnRWVBV2K519mz4Ix+RsAFlXGYkoL6JzGLmJLPGAgojhNd0EjJ01ZS+hkFm/AIQOBmfMq6yQ8pm++RcBAJIQgq6JkxkxpS+hkLpzEhDAEFldMuYRkln9FwIj0COXMlVFSuDCKmi9XPiE74KOeZH4YnoBGdU1AJidpw7iTSz5i4mAjCSVxQ+ek0nTS15DMqLmeHgJIVVFc7B8mnm37GxQWVEVKKzXil7RykmwE232NLcxrTyeudAwIaW9SCal3vtEkxK/QzukQn0p944H2EgMYXqaA5xbMTd5bIZwKJPsp53qnzYwQ0gfAFzS+dhlgVpZ8psAaK4IJuxmcdz074iJLPNFJ9Eum8Dx+NENAkQhW0iHj2Qs4t/UwkzVH4z1ch66AENJEABS0qnp1AM00+U0m0D3SMy9nl79hKaJHx7ByeeVY/0yGgiOPDO6FF17Nx6CWgnwkRUMTiU1Di2fN9KchnSnwSCfE4JpR6Tni8IKx+pkVAEZFDQYmnndvrQj8TI6CIaldBWfVctf/VIZ+pxQ/o52u12/PN+tjkCGhymyuh1HOHXa8Rq5/pxQ5ofTH6ypcfoX+RDwKa3mpBiedumy8V/TxB7IBe+0PvIgpKQM+wVFDq6Wh1ZZ18niFyQKtbIFUXEqnOZJLwWSQCegrbSij19LLUUFY/zxE5oNfuU/DPTXkJq6AE9BzTglLPIyyvHv08SdyAGtdhuom4oAgBPYtRUOp53OQ1JJ9niRtQ40qgn68SjsQT0NN0BaWegRgvJP08TbKAyriqMgE9T7UZTz2Dal5ONt9PRECRCpvu4dHPkxFQpNHUM+l944vwzCd/ls5DQJFAt+6Z4q7HRWlWP1m1Pw0BRWzmpjsFDarffGf3yEn0BvT5hO6H9QloctORTUKDGe/9JKFnIKCIyTKoKWgg86NHJDS56AG1CXJFegIq3eJmJQUNwXrwndXQxBQFNMSTEdBk1vbKsRJ62OLJSyQ0KQKKCLZGMQU9aPXkTxqajqYr0t8JqA57BjAFPWLr3HkSmoqmgI4visc+UKF2jl1WQr3t+uwRDU1CVUDry+N1V8UjoBI5DFsK6mnvZzc5NzQBZQGt77DUdJOAyuM4YimoB6ePvpPQ2LQFtL7JUr0ZT0ClcR+srIS6+up86RASGpW6gNaXZq6uzUxAZfEaqBTUjdeVl0hoRAoDWm/GvxBQUbwHKQndz/vCdSQ0Go0BrTfj//hPAirFoT1tFHSvIxf+ZGdoJCoD2mzGXwioCIfHJgXd5eB1k0loFEoDWm/GE1ABQgzLD9ZCN4W47jwJDU9rQOtTQgno2QINSQq6JdB9O0hoaHoD6oWAhhRyOFLQNeFue0RCwyKg8BR4KLISuijsbeNIaEgEFF7CD0MKuiD4bTdJaDiKLmcXZOIIaBBxhiAJtYhy12ISGorygH6+rj2Z7XeH+s0lizb8KOhMrLu+k9AwCCjcxD2fkIKORcpnhYSGoHwf6HpA5wjoQdFPx2Yl1BBr9bNFQo9THlBXBPSQFCOOgvYi9/NBQo8joNgr1WgjobX4+ayQ0GMIKPZJONIo6CNVPx8k9BjlB5GcJ46A+kk8yopPaLJ8VrjOiD99Ab32T/PiMXEE1Ef6AVZ4QZP280FC/SkL6G3yTK4NJaAezhlcJRc0cT5rJNSLqn2g9TXsDrWYgLo7a2AVuxKaevWzQ0E9aApodSF6s5j1+m1/n/hdCKirE0dVoQU9q58PEupBU0Dvs1xWSX13eQoC6ubk7boCE3piPh+nv98KJQtotbr45cfj5rbGOHKdb7A/n9VpNygBdXH+cCquoOf28yHhPdclUUCbo0lVQB1XGQ3P1c35JehvbtvwBHQ/GUOpqISens+KjPddizQBbY/GPwNanYTkWdDnk8x/8s5n4SORMowKKqiIfj7kvPUaJAlotavyz5+fr8+AVgX1PIuJgCYkaQwVklAp+axIevtlSxLQe33/tyag1Ua8+83gKmzCJyNsK66Igkrqp7glQK4kAb3WlWsD+uyg54EkDiKlIXDw5F9QUfmsCFwKJEoR0HbVsQ3oc6XRcxvefhqT0+osAd0mcuDkvRL6VdbqZ0vkgiBNioC2Oy+7gDruthxYT6R3ey4CukXqqMm5oDL7+ZC7MAiiKaC2j9Y3T7l/4gjoKsnbbbkmVGo9K5KXBxlO2IS/+u4DfZiXYuJiIuEJHy55FlRyPx/il4nTJTqIVK1ztgF1Pe4zfzIuZxeH/KGSX0KF57Mif7E4U6rTmF66gDp/fD0oArpIxTjJraAK+vlQsmicJM0nkeqz5+uA3pwvoBQUAbXTc0Hdj4waqiOfD7bjVyT9KOfZN/QgoAs0DZB8Cqqmnw9WQhcluphIfQbS+f0koDbqBkcWCdWUz4q6pSSNZJezaxN6aj4JqIWmtc9OBgXV1k+dC0p8mi6oHAABnVI6KpQnVF8+K0qXlagIaNH0DgnVBdXZz4fm5SWWtAG9nrwFT0DHVI8HtQlVm88H2/EzCQL6+/tL/1+nnsP0IKAj2seC0oJq7ueDhE7ED+htiGZ9W2LHT6+HRUB7OQwEhQlVns9KDktOMNEDeh0188p5oELkMQjUFTSDfj5yWXiCiB3Q26SY1Sn1hz4KfwwBbeQzAlQlNI98VvJZgA6KHNCql+MrHj8340/ciCegtZwWf0UFzaefeS1CR0QO6G2+vul9T6QQCOgjv2VfSUJzymclt8XIT+SAXuerm8+V0vOOxBPQHA8BqChobv3McklyFzegthvIPb923mEkAprnQi8+ofnls5LnwuQkbkCtF0+23FwzmdIDmu0SL/saTWJvenRYtgvUXgS0JDkv7oILmm0+H3kvUntEDyib8GLkvqwLTWjO+azkvliti74PdBZLDiKdJP8FXWRBc+9nCQvWivinMU1vgHQ/80z6cgNaxlIuLqH557NSxsJlFTmgn6/T1U1uKneKUhZxYQUto5/lLF5zkQNa5XJ82vyNm8qlV9LyLSeh+R57tyhpETPF/ix8dQEmc4u9uprIeSughQa0rIVbyBlNX4vqZ2kLWS/61Ziqq4l065zNzTnP+yBnmQEtb8mWUNCy6lkrb0F7pLge6P0yduL6Z5EBLXKxPjuhBebzUeailuCK9KObwp97QfoCA1rgMl05t6Bl9vNR4tKW5p5It6aep16MvlZaQEtcJ2idtyu02Hw+ClzguCtnzopbnE1nFbTkfpa3yBHQbHG1sTM25MvOZ6WsxY6A5op+PtKvhRZ26pJdUQseAc1TUQvxmqQFJZ+NgpY+ApqlgpbgTckSSj575Sx/BDRD5Sy+u6TZjmfrfaSUZZCA5qeUZXe/BAUlnxOF7IMnoNkpYrl1FTmh5NOiiIIS0MwUsdR6iFpQ+mlXwMJIQPNSwCLrK1pCyeei/BdHApqV7JfXI+IcTOLg0arcl0gCmpH8/94fFKGg5HND5gslAc1H5otqEIETSj635b1YEtBc5L2cBhNyO56t931yXjQJaCZyXkjDClZQ8rlXxgsnAc1DvktoBEESSj5dZLt8EtAcZPwXPorj2/FsvTvKdREloBnIdeGM6GBByaezTBdSAqpfnktmbAcSSj69ZLmcElDtMv3LHp/ndnxp93sPKMdFlYAql+NCmYpPQcnnARkurARUt/yWyKRcE0o9D8pueSWgmmX4Fz0xp+141j6Py22RJaCK5bYwnmJ3QslnEJkttARUr7yWxPPsKyj5DCWr5ZaAapXZX/JTba+Fks+Aclp0CahSOS2E5/tYTyhb72FltPASUJ3yWQKlWE4o+Qwvm+WXgKqUzfInyUJByWcMuSzBBFShjLaAZLEklHxGkslCTED1yWTRk2i6Hc/Wezx5LMYEVJs8ljuxzISSz7hyWJIJqDL0M7a+oOQztgyWZQKqCvlMoU4o+UxA//JMQDXRv7zp8PHB1nsa6pdoAqqI9oVNjTqfwW8hDxvlCzUBVUP9H2s16rXPwHeQxxLdyzUB1UL3cqZIv/FOQtNQvWQTUCU0L2SamDs/j9+8E7soTigB1UHvEqbK9NgRCU1Db0EJqApqly9VbIfeSWgaWpdwAqqA3r/PmiyduURCk1C6kBNQ+ZQuWsqsnPhJQVPQuZgrC+jn6+Vy+fNn/+/f3y9//OXw8xoDqnK5UmbrvHkSmoLGJV1XQK+Xxkv3hQICqnGpUubr9ueO2I5PQeGyriqgXT+HldDsA6pzu0aVnZ/aJKEJ6FvcNQW02n5/f/7vbSho7gHVt0Bp4/ChdxIan7oFXlNAr5cvP+r/eGazLWjmAdW2NKnjeM0QEhqfsmVeUUCfsez3fV7bguYdUGXLkj7ul1wiodHpWuoVBfTXt3oDvtEWNOuA6lqS9PG7Yh0JjU3Vcq81oFVBX7IOqLq9Qcr4X/GThEamadFXG9BqP+hbxgHVtBApdOyCySQ0LkULv6KAmvtAm39e3rMNqJ4lSKPj15snoXGpWf4VBbQ6fclYBa3OavryH5kGVM3yo1GY23WQ0Ki0jABNAa3OA30z/n2vz6nPMaBalh6Nwt3tiITGpGQMaApofQa9uRl/zzSgSpYdjcLeLI6ERqRjFKgKaF1Qcx20WifNLqCK9qBrE/5emyQ0HhUDQVdAqyNHb6Mv3LILqIrFRqU4tyomodFoGArKAnqU/IAqWGZ0inendxIajfzRQEBlkb/E6BQvnxUSGov48UBARRG/vOgUN58VEhqJ9BFBQCWRvrToFD+fFRIah/AxoTygn69rB5EuFgknzpWGfeb6pMlnhYRGIXtYEFAxZC8oSqXLZ4WExiB6YGQd0DnBAZW8lGiVNp8VEhqD4LGhPKCu5AZU8DKi1I5bxcVAQiOQOzoIqAiit1JUOimfFRIantjxQUAloJ+BnVfPGgkNTuoIIaACSF04tDo5nxUSGprQQaIvoMPN4V+2HzwlMaCsfoYlIJ+VDxoalsxxoiygt8lJSa4NFRhQmcuFWkLyWSGhYYkcKaoCWl29bsLpJCaJAZW4UOglKJ81EhqUwMGiKaDVXZDMYv76dmlvD7+buIAKXCT0kpbPCgkNSd5w0RTQ+yyX9Y3lXJ5CWkDlLRB6ScxnhYQGJG7AaArodb7B/lwJddoNKiyg4hYHvaTms0JCw5G2I1RRQGdXo6/c3LbhZQVU2LKgmOR8VkhoMMIKqiigz7XN+eb6XfFn4WUtCYpJz2eFhAYjatwQ0NOIWg70OvEzm244qykUSSNHUUAz24SXtBTopSafFRIaiKCxoyigWR1EErYnRylN9WyQ0CDkDB9NAbWfxjRfK10hJaByFgDF9OWzQkJDEDOANAXUeiK920eRhARUyruvmc58VkhoCELGkKaANsUc+/LD6RlkBFTIe6+Yql2fc+wMDUDGKFIVUPNSTIovJiLjnVdMeT4rJPQ4EeNIWUAfGVzOTsT7rlcG9WyQ0KMkjCR9AT1EQEAlvOt6ZZPPCgk9SMBYIqCJCXjP9cqpnjUSesz5o4mApnX+O65WViufPXaGHnL6eCKgSZ3+fquVZz4rHzT0gLNHFAFN6ex3W61s69kgof5OHlMENCH66SXflU8DDfV17qgioMmI+fSZLkXks0JCPZ06sAhoKvTTRyn1bJBQL2cOLQKaCPl0V8zK54DVUC/njS4Cmgb9dFZgPisk1Mdp44uAJkE/XZVZzwYJdXfWCCOgKdBPN4WufA5YDXV20hgjoAnQTyfF57NCQl2dM8oIaHz00wX17JBQN6eMMwIaHf3cj5XPEVZDnZwx0ghobPRzN/I5Q0JdnDDWCGhk9HMn6rmAhu6XfrQR0Ljo5y5fyecKErpb8vFGQKOinztQz200dKfUI46AxkQ/t1HPfUjoPonHHAGNh8uHbKKeLmjoHmlHHQGNhn5uIZ+uuHr9DknHHQGNhXyuo55+aOimlEOPgEZCP9dw3OgIGroh4eAjoHHQzxXU8zASuird8COgUdDPRdQzEBq6ItkAJKAx0M8l5DMgEros1RAkoBHQTzvqGRwNXZJoEBLQ8OinDceN4uCI0oI0pxES0ODo5xz1jImGWiUpKAENjX7OUM/oSKhNgrFIQAOjnxPUMxEaOhd/NBLQsOjnCJvuKbEpPxN9PBLQoOjn4Cv1TI+GTsQekQQ0JPrZo55nIaEjkcckAQ2IfraI57loqCHuqCSg4dDPGuueAnywLd+LOi4JaDD080E9BaGhnZgjk4CGQj+ppzg0tBZxbBLQQIrvJ/WUiYY+Yo5OAhpG4f2knpLR0Hjjk4AGUXQ/qad8xTc01ggloCEU3E/qqUXhDY00RgloAMX2k3rqUnRD44xSAnpcof2knhoV3NAo45SAHlZiP/mcu2LFniAaY6QS0KPK6yf1VK/QhkYYqwT0oML6STyzUWJDw49WAnpMUf2knpkpr6HBxysBPaScfhLPPH0UtjUfesQS0CNK6Sf1zFlZDQ08ZgnoAWX0k3iWoJyGhh21BNRfAf1k1bMgpayIBh23BNRb9v2knsUpo6EhRy4B9ZV3P4lnsQqIaMCxS0A9ZdzPr9SzcNkfVgo3egmon2z7STxRy7uhwcYvAfWSZz+JJ0YyjigB9RIuoGGeRw6222GT7dZ8oCFMQH3k1c+vxBNrsoxooI1IAuohow144old8otomFFMQN1l00/aCReZRTTIOCagrv6WRT9Z8YSXj4wqGmIkE1BHGfSTzXYckk1EA4xlAupGez5pJ8LIIqLHhzMBdaK6n8QTYemP6OEBTUBdqO0nW+2IRPn2/NEhTUAd6Own8URkmiN6cFAT0P309ZN2IhmtFT02rPUF9HrpvLj/8JHZ1dXPr8QTyamM6KGBrSygt8uYa0MPzK6eftJOnOlDW0aPDG1VAf18vUz98ZfTM/jPro5+0k6I8KGqogdO7tYU0N/fx8X89e357z9/ujyF9+zK7yfthDRqKlpGQO+zXFZJfXd5igMB9fzBFGgn5NJRUd8Brimg1/kG+3Ml1Gk3qO/sSu3nV9oJDcRv0vuugyoK6HN18232xZvbNrzn7ErcgKedUEZ0RT3HuKKAPtc255vrd7fDSH6zK6yfpBOKSc2o3ygnoNvE9PPrV9qJLHwIzKjXOFcU0LM24QX08yvpRIY+PkSF1GekKwroSQeRzu0n5UTu5GTUY6xrCqj9NKb5WukK99k9q5+UE4UR0FH30a4poNYT6d0+iuQ8u8n7OdlaJ50oy8epK6TO411TQJtijn354fQMrrObsJ+UE+id1VHXEa8qoOalmJJcTCR+P6fZJJyAYRrS6CXNO6CPtJezi9bPeTcpJ7Bo1tF4KXUb9PoCeojb7IbuJ9kEjrGVNGxK3TY7CeiyIBvwtmjSTeAoa0pDxNRp3BPQRX79XAgm2QQiWUipd1BdRj4BXbLnVVytJdEEUluN6d6mOhRUeUA/X9dOBJ2d9HTZP7t/2xNHggmItieoc/sLSkAXUEsgS/sKuvPJsg7oXPjZBVAu5QF1RUABhENAAcATAQUATwQUADzpC2jKz8IDwAplAb2lvRoTAKxQFdDP19lpnW7XUyagAALSFFDrFemd7ilHQAEEpCmg9nsizW91vIKAAghHU0DPuSsnACxQFNCz7gsPAHaKAvpc25xvrt/5LDyAsxBQAPCkKKBswgOQRVFAOYgEQBZNAbWfxjRfK11BQAGEoymg1hPp3T6KREABhKMpoE0xx778cHoGAgogHFUBNS/FxMVEAJxNWUAfXM4OgBj6AnqI7T6dAOAreKNCP2FIZ7/YAPISvFGhn/AUWW3r5zQzOc0LMyPWiTOTx8vI0iBUTvPCzIhFQA9iaRAqp3lhZsQioAexNAiV07wwM2IR0INYGoTKaV6YGbEI6EEsDULlNC/MjFgE9CCWBqFymhdmRiwCehBLg1A5zQszIxYBPYilQaic5oWZEYuAHsTSIFRO88LMiEVAD2JpECqneWFmxCKgB7E0CJXTvDAzYhHQg1gahMppXpgZsQjoQSwNQuU0L8yMWAT0IJYGoXKaF2ZGLAJ6EEuDUDnNCzMjFgE9iKVBqJzmhZkRi4AexNIgVE7zwsyIRUAPYmkQKqd5YWbEIqAAoA8BBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPKkM6K9vl6c3x28JtTLFn6/Vt778SD1J3rZe/ecMvSecnGPWZuZefUvRvGwvZS+pp+igX9/++Mv65cTjX2NAr5eGpSwr3xJqeYqbZaFiXVQE2nr1qzlSE52VmWmSo+iNWZmZ39+7pUzNO1N5TrbtxU8//hUG9Na94/OXcOVbQi1P8dBPLXOz+epfFQ3TlZnp+3m5/PnzlIlztTIz1/5bat6aytW6kJ0w/vUFtFp8q+X2Pt/uWPmWUCtTXC0M9abIVcncbL76d0WjdGVmqr9sdTj7d0i6lZnpvlTNlI4/05V6tXk+uWeMf30BvXYvzvPlmqypr3xLqOUpNjZ3bzqW7a1Xv1mlVhLQlZm59SuedyWroOtDpvkboGn3SrMJMB8TZ4x/dQF9vtHdK3ebrACsfEuolSm+D39Eq7+38pftrVe/2mv1vzXMSGVlZp7z0Q1O4z8lW3tnrv0cKBky7R/if7HsAz1l/KsL6POPy8vwn6MVgJVvCbUyxeYicNXQna1X/zlD7zcNM1LZt5QpsfbOjAKq472p9py82w4inTL+1QXUCMv0RVz5llA7p/iqYUVnY16ei/SbqkG6NDN3LWtqvbV3pv+ekrXpRzXJL/ZF7JTxrzGg3Ricvucr3xJq3xQ/N00UrE6vz8vzS8950BTQ1aWsOflHx7ysvjPV9vCL8b9aLAQ0/fhXF1Bzc3ayZrbyLaH2TbGOvVPr89J8RU1AV5eyLz9uqU+WOWT1nRnOltPUT2tATxn/CgM6vDKWRXvhW0LtmmIl+3NX56Xd7lUU0LWl7F91naC7sZRd9Z0FuhTQE8a/8oCO3vaVbwm1Z4qf/VQ4TMfz0u2EUBpQc6Kbjfd6brScbry6lA0n0mtY4+htBzTRkqY8oPmvgd6VrOaszsu1nQelATVnpg5ou0dFyenG6+9MMzfVbGmYlw5roL4KC6iafq7NS9/NDAJaJaffo6JxO2c0M8afgJuSTwU0CKivso7CX/Us1cvzMpyfpyagK2/M1Tikp+OUpvUh87bwLdk4Cu/rvusMPSXngW5McbVdpWEnW215XoZLPDQUNGfljbmpC+jKzJgramr+ulVsw+WU8a8uoKV8Eqn5kobx2VqeF4UBXXljjM/YKgnoysxkFVA+ibRHKZ+FV3OIorM8LwoDuvLGmGdF6NgHujIz4014DTPTsgWUz8LvUsjVmKrFQcc8dHa8+npWcvYsZYrOMFuamfvFPIikYmYa1m10rsa0RyHXA9W1QlDZ8errCeiOpaz+noK16a2lrD+NScfMNKwB5Xqgu8wuOz2sq2VwRfpuZu7jzV4N3Vl5Y4ZHaJiRyq6lTMWO9sfazBhX19cyMzUzoKeOf4UBnd34xBin+u+J1M7McKsaPQFde2MaigK6NjP9Hzc1yVmemeGz8GpmpmIPKPdE2mdy6z1znKq/K2c7M+YdkdQEdO2NqWkK6OrMXJOu5QSwMjNL13cXbSGg3JUTANQgoADgiYACgCcCCgCeCCgAeCKgAOCJgAKAJwIKAJ4IKAB4IqAA4ImAAoAnAgoAnggoMJy9EgAAA1dJREFUAHgioADgiYACgCcCCgCeCCgAeCKgAOCJgAKAJwIKAJ4IKAB4IqAA4ImAAoAnAgoAnggoAHgioADgiYACgCcCCgCeCCgAeCKgAOCJgAKAJwIKAJ4IKAB4IqAA4ImAAoAnAgoAnggoAHgioADgiYACgCcCCgCeCCgAeCKgAOCJgAKAJwIKAJ4IKAB4IqAA4ImAAoAnAgoAnggoAHgioADgiYACgCcCCgCeCCgAeCKgAOCJgEKDX/+n+v+/v1++/Nj5WCABAgoFrpe36n92BbR9LJAAAYV898v+gHaPBRIgoJCPgEIoAgr5CCiEIqCQj4BCKAIK6Z5JrL23Aa3//efPyfffJ4+tXOv/3nHgHvBDQCHdJKDXUSMfn6/tv//46zEO6O3SMWoLhERAId04oH+/mMEc+tl8wQjo0E826hELAYV8xj5QY+3yrf1KvYlelfPFfOyvb8Z/sAqKOAgo5DMD2qx4VgV9ab7TxvHZya6kb+3PvDQ//VxJbX8ICIyAQj4zoO3W+LOKVTmvwyGi9kHzo/DPtBJQxEFAIZ/lNKYmoM829lvnz6+8POYBrTf7CSjiIKCQby2gpiqmRkDv/TcIKOIgoJBvMaDGMfhZQM3vEVDEQUAh31pApwfYzaPwTTvZB4poCCjkW9uEn7axe+ytX+8koIiGgEK+xYBWR4jerY81HnpnEx6xEFDItxhQYz1zchrT8NBqW56AIg4CCvm60+XnAa33dNbroNfhI0n1Y69tNm9cTwTxEFDI1xxQf7cE1DhVqV1L7R5rfoOAIhICCvmaz8C/2QJqnK30Nnrso7tq0+W/zXeUAmEQUChQV/HFGtDusksvk8d2aX0zPxYPhEVAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcATAQUATwQUADwRUADwREABwBMBBQBPBBQAPBFQAPBEQAHAEwEFAE8EFAA8EVAA8ERAAcDT/wcluY9Awy3ROQAAAABJRU5ErkJggg==)
# Compute likelihood function given x at theta.grid values
# Issue: scaling of likelihood; use log scale and normalize
loglike.grid<-( log(probs.grid) %*% as.matrix(x))
plot(theta.grid, loglike.grid)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAaVBMVEUAAAAAADoAAGYAOjoAOpAAZmYAZpAAZrY6AAA6OgA6kJA6kLY6kNtmAABmZjpmtv+QOgCQOmaQkDqQkNuQ29uQ2/+2ZgC2/7a2/9u2///bkDrbtmbb/7bb////tmb/25D//7b//9v///80U5yZAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3dbYPaxpqgYZLNjr2zzjgzs2E3PRsc9///kYN4Fc1Lo4eSVE/puj6c43QDRkjcllCpWL0DELKa+wkAZCWgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQVHdAVwDFlE9U8UcsaO5XG2hL8UaVfsCSRvgHA1gsAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElAymmmqMbgkoNRk1JkblZbSmg6od0rNpm5kxNyvEbUTUCYydfxGMfeLSGWaDug174BJzZ278c39CjMzAaWouYs2o7lfemYgoBQwd7tqM/f6YCoCStzcnUph7pXEmASUweZOUlZzrzfKE1CeNXd/2jH3mqQYAeVTc/emWXOvWF4moNw3d2CWYu71TJiAcsPcRVmqudc7QwkoF+ZOCCsbaSICysHc1YhoamGuTLjyCRJQKqzNkpf9ynQvBoMJ6KLN3YZVrStk7lfl2tyvCDcJ6FIJwSBzvVwfzP0y8IGALpB3/eumfg2bf0GTEtBl8RYfwVQv6qJf5EoJ6HJ4R49v/NfYi14VAV0Eb+EZjPuiWwVVENDmedPObrRVYIXMTUCb5m1alzHWh/UzJwFtlndmvYqvG+tqJgLaIu/GJAqvKKttcgLaGm/BjEquNetwQgLaFG+73IqtP6tzIgLaDO+1ZhRalVbs+AS0Dd5gDSqwUq3kkQloft5WTXt99Vrb4xHQ5LyblqFEOa328gQ0M2+ipSmTT+u/GAFNy3tnsUol1IbwMgHNyXsGEa2AgCbkrcKRiM5LQLPxFuEjFZ2NgObincE9GjoDAc3DW4JPiei0BDQJ7wSepaHTEdAUvAUYSkSnIKD1s+0TpaEjE9Da2eh5jYaOSECrZnOnDBEdh4BWzGZOSRpanoBWywZOeRpaVr6Ark9r9MvwO6fZDGzajEZDy0kW0LcPa3RoQ5NsAzZqRqahZaQK6I+vV+vz178GPUKKDcDmzCQ09HWZAvrzj8ti/vP79r9/+3vIQyRY+zZkJiSir8kU0M1VLrukfh/yENWveVsw09PQsEwBXV8fsG93Qgd9DFr5arftMhMNjUkU0O3u5rerH74NO4avep3bapmVhA6XKKDbvc3rw/XNsNNIFa9xGyzz09CBBLQOtlVqoaEDJApow4fwtlKqoqHPShTQdk8i2UCpjoQ+JVNAbw9jut4rfaDC1WzbpFYa+qlMAb05kH7YpUjVrWObJVWT0McyBXRfzEu//DnoESpbwzZJ6qehD6QKaH8qpr3ck4nYGMlBQu9JFtD3hqazsyGSiITelC+gL6lo1doISUZDrwnoPGx/pCShl/IFtIVDeNseeUloT7KAtjEjve2O3J5vaOsbc6qAtjEj/cK3OJogoXuZAtrGjPRL3tpoiYS+5wpoEzPSL3ZLo0ESmimgDUwmstCtjHYtPKGJAtrAdHZL3MJo3aITmiig6SdUXt7WxUIsN6ECOpmFbVksylITmiigyQ/hl7RVsUSLTGiigKY+ibSgTYrFWmBCMwU08Yz0S9mcWLjFJTRTQPPOSL+MbQnel/YNiZkCmnZG+gVsR3CypISmCmjSGelb34bgg+UkNFlA3/NNZ9f4BgS3LSOh+QL6kunXWcsbDzyyhIQK6Nh/YatbDnyq/YLmC2iqQ/hmtxt4SusJTRbQZDPSt7nNwABtJzRVQLPNSN/iBgNDtZzQTAEdPCP9zKuquY0FYtotaKaADp6Rft411dqmAnGtJjRTQHNNJtLWdgIv+jShKd8biQKaajq7hjYRKKPFgiYKaKYJldvZQKCc9hIqoGNoZeuAwloraKKA5jmEb2TbgBG0ldBEAU1zEqmJDQPG0lJCMwU0yYz0+TcKGFc7Bc0U0Bwz0qffJGB8rSQ0U0BTzEiffHuAaTRS0FQBTTAjfe6tAabTREKTBfS98unsMm8KMLEGCpovoC8Zea0k3hBgBukLKqBFHz3rZgDzyL4TKqAlHzznNgAzyl3QZAHdTancGwz684+KLuVMuQHA7BInNFdAjyeQTqePagpowrUPVchb0FQBPZ+AP+6EVhTQfOseqpG1oJkC2h2/d/MxvZ0LWk9A0615qEnSndBMAV0frzvqruncF7SagCZb7VCdlAVNFNBtLE+ffa4PBa05oCP9TdCqhAVNFNCLCZUPBa0loJlWOdQqX0GzBrQr6JdqAppohUPNsiU0bUD3U4HWEdA0axtql6ygiQLa/wz0/fCd8AIKjUlV0EQB7YYv9b9V7sfX1S//p4aAJlnVkEOmgmYKaDcOtP8FHpvdizp7QHOsaEgjUUEzBXQ3gr5/GL+pIaApVjOkkiahqQK6K2h/H7TbJ60voMX/CliaLAXNFdDr74Z/mzug9a9iyChHQZMF9FXFF7f6FQxJpSiogL76gFWvXsgrQ0HTBvRt6Dca75Re3LpXLuRWf0EF9BV1r1vIrvqdUAF9RcUrFlpwr6C1vNME9AXVrlVoReUFFdAXVLpOoSVVF1RA4ypdpdCWmgsqoGF1rlBoTsUFFdCwGlcntKjeggpoVI1rE9pUbUHTBjRmzIAWe2Tgo0oLKqDhR6psTULT6iyogIYfqa4VCY2rsqACGn2gqlYjtK/Gggpo9IFqWouwBBUWVECDj1PROoSlqK6gAhp7mHrWICxIbQUV0NjDVLMCYVnqKqiAhh6lltUHiyOg8xFQyK6mggpo6FHqWHmwSBUVVEAjD1LHuoOFqqegAhp5kBrWHCxXNQUV0MBj1LDiYMlqKaiABh5j/tUGC1dJQQV0+EPMvtKAOgoqoMMfYu51Btws6BxPovgjln7Akl5f3ArWGVDJLqiADn4E/YQa1FBQAR38CAIKVaigoAI69AH0Eyoxf0EFdOgDCCjUYvaCCujA++sn1ENAJyWg0JSZCyqgw+6un1CVeQsqoMPuLqBQlXk/BhXQYXffPoCAQkUEdDqvLa4dUKjPnAUV0EH3Pq8fAYVKzFhQAR1yZzugUKH5CiqgA+4roFAlAZ1GyYCWe1bAS2bbBRXQAfcVUKjTXAUV0OfvKqBQKwGdQsGAFnxWwKvmKaiAPn/X3Rj604oq+KyAV81zEC+gT9/TDihUTEDHVyygJZ8UUMAcBRXQp+95WCECCnWaoaAC+uwd7YFC3QR0bKUCWvRJAUVMX9CmA3r9cr4Y0HdH8FCxyQsqoE8/kolAoXLldpme/guLP2LpByzplYDqJ9Ru6oIK6JP3O68OAYVqCeiYXgioI3io38S7oAL63N0cwUMK0xZUQJ+7myN4yEFAxxMPqCN4SGHSXVABfepeFwUt/aSAggR0NC8E9LAiBBQqN2FBBfSpe031iQrwuukKKqBP3Wt1uRsK1ExAxxFa3N4K0E9IYLJdUAF95k4+AoVUBHQULwbUETzkMFFBBfSZOxnEBLkI6Bgii9t7+QUUkpimoAL6xH0cwUM2AjoCAYWFmKSgAvrEfXoH8FUvHXAmoMUFFtdHoJDSFLugAvr5XczEBBkJaGnhPVABhWwm2AUV0M/vcn759RMSEdDCogG9+AOQwvi7oAL66T18BApJCWhZsYCayg5SGn0XVEA/vcf5pddPyEVAi3oloHZAIZuxd0EF9NN7CCikJaAlDV5cH4FCZiPvggroZ3dY3fwzkIKAFhQJqHNIkNe4u6AC+tkdfAQKmQloOQIKCzPqLqiAfnJ755AgNwEtJhTQG38GshDQYgIBXZ3/XPrZABMYsaAC+sntT6+5HVDISUBLiQfUR6CQ1XgFFdBPbr9yDgmSE9BCBi7uxSegAgpJjVZQAX18c+eQID8BPVmfXoQvw+8cCOjq9Mfhfx1QAwHde/vwOgxtaCSgPgKF7E7v4cJv5lQB/fH14z8kq1//GvQIAgpLNNYu6LgB/ef3q+IFqnf084/L++4e/be/Bz254EkkAYXURipopoBurnLZJfX7oCc3aHFXTsJDGwT0fX19x+1fMOhj0MEBdRIempAxoEfbfcdf/tz9qUvqoH3Gs+3u5rerH74NO4YfHlAn4aEF4+yCThLQH1976dscWzrUtr3X6d0M250dGtD38v9kATNIHND1xXH2OjKA832egI6wzw/MIW1AP5RvYPNOZjqEP7zaAgqpjbILmiigs5xEOv1BQCG3rAHd7jpeHsIPG7t5cnsY0/Ve6aMnN2Rxey+yfkJ2Y+yCTvQZaG8X9C10EXvn5kD6YXuzgwPqJDw04uITuUwB7S7BPISu+2PwLPzNYaUDH2tgQJ2Eh3b04pkqoB/mAAmOA+2sV5dGnUzEN3JCS0Y4hp9oMpFNeJ/xynTT2fVfaAGF7EY4hp9sNqbD4XfwBHwpwYAaxQQNKH8Mn2o6u9cNDejpDwIK6ZU/hs8X0MkO4XsvsH5CC5Ye0ClnpO+9wgIKLUgW0P3cS1ejj6IfhE47I71RTNCY4sfwmQI68Yz0RjFBa07v5gUGdOIZ6fu3FVBoQX/vs/6AljXxZCICCq0pfQyfKKBTT2fXf5EFFJqwKntuY5rJRIoMn596QuWyHzYDFSh8DD/RdHbBCewuzBHQyz8AySU8hB/6QeUdEx/C2wGFBpU9hk+0BzrxSSSjmKBBZY/hJ/kMdPPSFHb9h5lwRnon4aFBl/1MEdDuEqIC+6DTzkgvoNCiVcmDy4k+A70UPSk/6Yz0AgotKnoaKVVAJ52RXkChRf12Li2g79NNZ3f5Mg//q4AanT//TLIHWhHDmGDpLvr54ptbQO/f1Dh6aFA/nssL6ESH8HZAoU2X/VxUQKebkd44emjUqty7O1VAp5yR3kl4aNSHiLz2UKWe1OkRr35yY/hmaFj94Bnpr/7aAS+WgEKj+jnIGdDQB5iDZ6QXUOBKrwRJ90BDBZ10MhEBhVatin0IOtm18IdrLruYfuv/4HnTTmcnoNCqch+CThLQbS7P6dvs0hmYI3TaCZV7L6t+QlN67cwQ0PVFLPf/NXDX8X3yGekL/QsFVOf0rs6wB/qhfPvmDSzf+5wz0g/4C4D6FTuGTxTQqU8iCSg0KlVAt7uOl4fw3U5jIKBTzkjvI1BoV6+b9Qe0u3y9fxJpl9P18MH0U85I7yQ8NKzUh6BTnYU/hq6LXncW/m3gruPpzh+MNSO9gELDSh3DTzMO9HIOkO+7og4eB9qZbEZ6AYWG5QpofxqQLpzbw/jADujeNNPZCSg0rNSHoJPNxnQ4/N4fyW9i+5+vE1DgvdiHoKmms3udgALvxaalzxvQ7S7t8E8BDGMC3ot9CCqgd2+5+vAHoBmrMt8sJ6D3b/ryGTqgVqf3drKA/vga/F7429OKjj+ZyPBnCtTu8kPQ+KOUfE67R/zk90kCagcUWtZ/h8ff5okCuhs9OlVAfQQKTTuW87X9pFSfgXb7oKcr6Ef9DHTlJDy07XJHLPwgRZ/T+8gnkdbnS0BHDujw+wCJLDGgu2tC990UUCDu8kPQ8IMUfU7vow9j6ua02x3GCygQ19/tDO+CTnot/Pbo+23wLKBXuqmdvgso8IoiZ5EmCuh+CFIX0NWNb4YbqDuM/yKgwCtKfAg6TUAPQzi3AV2vChS0O4z/9f8KKBCXJqD7zy1/fO1OoK+HfgvHTfsZmg1jAqJKnEWaJKCbXev2AY19mceV3QzNU8zGZCA9tKnEWaSJvlSuO3V0COh2d/T1E0nvH76p7lku5QT2SpxFmuhrjbvUHQK63QUtcAwfFJhMZMRnA8yowIegUwT0n993542OAQ18I3wxeggcLDmgb6GvVRJQ4KDAWaQZDuHXRT4DFVDgJQXOIk10Eqnb5zwEdLs/Gvg+4itjB9RHoNC4AmeRphrG9OUY0G5M6Msj6d9HD6iT8NC81z8EneZKpN3o+V1AuxHwRUYxjRtQw0ChfVkCevFtHGXOwY8a0N6LKaDQqtfPIk00mUh34F60n2MHdPBdgGxeP4s02XR2h4QWGwIqoMBrLg416w5oaQIKvOjyPHzoAYo/pauf/POfl//9899rvxJJQGEJXj6LNM2VSBfjlhJcyimgsASXZ5ECb/UJL+Xc202GLKDA/Po7njUH9PR55abomaTBDGMCTi76GXmzT3IS6cfXQ0H3p+JLXIgUZCA9cHbxTq81oN1+Z1fQ3e5nkeuQolzKCfS8eBppomFM23T++v/m3v18N5kIcCFHQPc7nzPvfr77RBO40OtmzQHdFTQy9L0sAQV6+vudNQe0m4dp7v1PAQUuvHg154SXcr6tisyk/BIBBfr6R/C1BfRiGruz6seBAgvRP4UkoJ89OQEFel47DS+g925oFBMswG7Hs9KAVmfAlUgKCgvw2jgmAb19M1dywiL0dpME9DPPLe5rI8OARF46DS+gn9xKQaFpL52GH/8k0vfrU0m1n0QSUFiM/gkPAf3kyQko0HfxtUgC+vjJCSjQ99I4Jp+BfnIrAYWmvTSOSUA/uZWAQtNeGsckoDdvZRgTLMUr45gE9PbNDKSHhXhlHJNr4e/cLviZMpDMK+OYBPTeDfUTFuGVcx4CCizaKwNBfQYKLNorA0EFFFi0VwaCCiiwbBfjmIbetfiTKf2AJQko8MEL45gEFFi2F8YxTRHQG+fi5/qGeAEFPnhhHNNMAZ3rK+IFFPggY0DnKaiAAh9UHtD39x9fV7/8uftTF9Nv/R9MS0CBD14YST9JQLe5/Hb6j80unduQzrELKqDABxeHxkPvWvzJXP9ofRHL/X+9zXIiSUCBj07hrHIPdLu3+b33n5vdtfCbWa6IF1Dgo/guqIACCxcfCDpFQH/+8eEQvjt4nyKgt07+j/13AtnET8NP9BnoxUmkXU7XE3wGKqDA5yoP6I+vpzlAu2FM3Vn4t35TpzN0Rvpxnw1Qg8oD2uWy5/uuqDWPA7XDCgtSe0B3xTzowrk9jJ9lB3Tot3IqKCxAfCT9ZLMxHa7n3B/Jb+bZ/3xycXsvoYBC++oexlSRJwM68A5AauGR9AL6+EYCCu0L74LOcwg/GwEFroRH0k8U0N6MdrMmVECBK+H3/DQBvZgRdKbzRzsCClypO6A//9iN/ux0A5jm+j6PdwEFbqg7oJteNLuYfr++yUQMYwKu1B3Qdf+Dz5mmUt4zkB64UnVAtzud/cuO5plKec+lnMCV1e4b4SsdxnRzPtB5mEwEuFbzONCEAQUWJVhQh/DA4vXqWV1A051EApYlehbJMCZg8aoOaLKB9MDCVB3Qy0s557wYXkCBa3UHNNdkIsDCVB7QVNPZAQtTfUDrIKDAtegEGAIKEJwAQ0ABqvxWzouz7zWchxdQ4KYKL+UUUCCF2DG8gAIEzyL5DBQgOI5JQAEE9BkCCtwgoM8QUOAGAX2GgAI3COgzBBS4QUCfIaDADYYxPUNAgVtqHEhfHQEFbqrwUs7qCChwm4B+SkCB2yK7oAIKEPwQVEABgqfhBRQgOBA0X0DXpw8qvgy/s4ACtywioG8fZhYd2lABBW5ZQEB/fH11bmYBBW5pP6A//7gs5m6++9/+HvIQAgrc0n5AN1e57JL6fchDCChwS/sBXV8fsG93Qgd9DCqgwC3ND2Pa7m5+u/rh27BjeAEFbmp9IP12b/P6cH0z7DSSgAK3NX4p54QBDUwqACTX9mQi0x3Ch+a1AhYnUUAnO4kUm1kVSK7pQ/g7w5iu90ofeGZxg3P7A7kFDj0zBfTmQPphlyI9F9BhtwdaEDn0zBTQfTEv/fLnoEcQUOCm0KFnqoD2p2IabzIRAYUFCr3xkwX0fYLp7AQUFmghAX2JgAI3CejnBBS4aSEBneAQ3jAmWJwlBHSSGekNpIflaT+gU81I71JOWJzmhzFNNyO9fsLitD6Q3oz0wHgav5TTjPTAiJqeTMSM9EBdEgXUjPRAXQQUIChRQB3CA3VJFFAnkYBRNX0WfqIZ6YFlGj6OKVNAJ5qRHlikwEj6TAGdaEZ6YIki13KmCug0M9IDSxSZTSRZQN8nmM4OWKJlBPQlAgrcJqCfElDgtmUE1CE8MIIFBHSSGemBBWo+oFPNSA8sT+vDmAbPSH+VW5PMA/c0PpB+8Iz0Ago8r+1LOU0mAoxp8H5WooCazg6oS6KAmlAZqIuAAgQlCqhDeKAuiQLqJBJQl0wBNSM9UJVMATUjPVCVTAE1Iz0wqpYH0r+bkR4Y0eBLkZIF9N10dsBIhl8Mny+gLxFQ4I7AdEwCCtAJTAiaLKD//H446/4WOIP0LqDAXc0HdH0YtrQJzqcsoMA9rQd0fdjr3JzPlRnGBJTReEC7b/Tozrx3w0F3u57dyHoD6YEiGg/o+jD9/Nt5BNPGpZxAGW0HdLu/+eXw/+fdzrXZmIAi2h7GdJwPdPv/52i+mQ8UKKPpgfS9gJ6vQTKhMlBKy5dyHidUPh7K7wgoUEzLk4msD+FcXxzC+wwUmEmmgG4Ooz5/fD19OdLNr/l4QECBcjIFtBv2udsHPZ85Wg8cSS+gQDmZArobSb87ZD90s7siadAOqIACBaUK6L6gr8yoLKBAObkCejUl/fUXxT8moEA52QLa/2KkQeff9wQUKCdfQF8ioEA5AgoQlDagb4H56AUUKElAAYIEFCBIQAGCBBQgSEABggQU4KjdCZUvCShQ2tAp6QUUYG/wlyKlDWiMgAL3DP9aTgEF2Bn+xfACCrAjoJ8QUOAeAf2EgAL3COgnBBS4R0A/IaDAPQL6CQEF7jGM6RMCCtxlIP1jAgrct5hLOWMEFHhgKZOJxAgoUI6AAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAApysDp69dfG/v/QDliSgwAOr1aCCCijAwSmcTxZUQAH2etkU0BsEFLirH4inYtF0QFc3zP2cgGoJaJ+AAgMI6GMCCtwloI8JKHCXgD4moMBdAvqYgAJ3Gcb0mIAC9xlI/5CAAg+4lPMRAQUeMZnIAwIKlCOgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQlC+g69XRl+F3FlCgnGQBfVtdGtpQAQXKSRXQH19XH/3616BHEFCgnEwB/fnHZTH/+X3737/9PeQhBBQoJ1NAN1e57JL6fchDCChQTqaArq8P2Lc7oYM+BhVQoJxEAd3ubn67+uHbsGN4AQXKSRTQ7d7m9eH6ZthpJAEFyhFQgKBEAXUID9QlUUCdRALqkimgt4cxXe+VPiCgQDmZAnpzIP2wS5EEFCgnU0D3xbz0y5+DHkFAgXJSBbQ/FZPJRIC5JQvou+nsgGrkC+hLBBQoR0ABgvIF1CE8UIlkATUjPTCqY1yevHHxv770A56ZkR4Y1zkuz926+N9f+gFPzEgPjOsUzucKmimgZqQHRtXLZnMBNZkIMKp+IZ6pRaKAms4OGFfDATWhMjAuAX1IQIH7Gg6oQ3hgXA0H1EkkYFwtB9SM9MCoWh7GZEZ6YFwND6Q3Iz0wsnYv5Xw3Iz0wsnYnE9kxnR1QiXwBfYmAAuUIKEBQvoA6hAcqkSygZqQH6pEqoGakB2qSKaCDZ6S/yu2zYxMAnpApoINnpBdQYEyZAmoyEaAqiQJqOjugLokCakJloC4CChCUKKAO4YG6JAqok0hAXTIF1Iz0QFUyBbTIjPQA5ZTOXOUz0gOUU7pyVc9I/7SmjvVbWpiWlsXCVGvGhal6Orun2Roq1dKyWJhqNRzQadgaKtXSsliYagnoi2wNlWppWSxMtdoP6Gbcg3hbQ6VaWhYLU61WA9rNSN8NXNqNaOo8mMzuJbaGSrW0LBamWm0G9JDNX/7snUkaqaC2hkq1tCwWplptBvSYzd/+63BN0tvg7/R4lq2hUi0ti4WpVpMBPVzKuftmpG+nH42zC2prqFRLy2JhqtVkQNeHS+G3BT3td65HOpFka6hUS8tiYarVYkBP09lt/3Cq5sDp7J5ma6hUS8tiYarVYkBPEyr3JwYdOKHy02wNlWppWSxMtQT0RbaGSrW0LBamWi0G1CF8VEsL09KyWJhqtRjQ00mkzXnwUr+lRdkaKtXSsliYajUZ0MMwpt1FnIYxDdDSwrS0LBamWk0G9HT95mkg/cZA+me0tDAtLYuFqVaTAT3OSP/Ln6eUupTzGS0tTEvLYmGq1WZAz5OJnL7cw2QiT2hpYVpaFgtTrfpkps4AAAlLSURBVFYDevY29pz0AJNr6d8hgEkJKEDQZIfwA7/RGKB6AgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQa5EAggSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIShnQf35fbX0b+KtKPXjGP752v0o0k+pnr/52gb5P+HRe82hhNt2vEi3L51vZl6mf0Yv++f3Xv27+eOL3f8aArld7N8ry4FeVuv+M99tC5+amUqHPXv1uidJE58HC7JOTaMU8WJiffxy3sjRrprN92rde/Onf/wkD+nZc49cv4YNfVer+Mz73M8vSfPrqrxO9TR8szKmfq9Vvf8/y5IZ6sDDr06/SrJrO+uZGNsP7P19Au823224318cdD35VqQfPuNsYdoci6yRL8+mrv0n0Ln2wMN2/bLtwntZQ7R4szPFH3ULl+Ge6s9ttvn66c7z/8wV0fXxxti/Xhz31B7+q1P1n3DvcfcuxbX/26u93qZME9MHCvJ12PDdJdkEfv2X2/wZk+nhlfwhw/Z6Y4/2fLqDbFX185d4+7AA8+FWlHjzjzfkf0e7f2/q37c9e/e5Tq//IsCCdBwuzXY7jm7P3x5o9WjPr0xIkecsc/iH+lxufgc7y/k8X0O0/Ll/Of7zYAXjwq0o9eMb9TWCdoTufvfrbBfr+lmFBOs9tZUk8WjMXAc2xbrpPTr7fOok0y/s/XUB7Yfn4Ij74VaWefMbrDDs6nyzLdpP+lupNem9hNln21E4erZnT75LsTb93T/nL7U1slvd/xoAe34Mf1/mDX1XquWe8PTRJsDv9eFm2P9ouQ6aAPtzK9oN/cizLwzXTHQ9/6f1/FncCOv37P11A+4ezH/bMHvyqUs894xyfTj1elv1P0gT04Vb2y59vUw+WecnDNXMeLZepnzcDOsv7P2FAz6/MjU37zq8q9dQzTvJ57sNlORz3Jgroo63s33IN0P1kK1vnGwV6L6AzvP+TB/RitT/4VaWeecbbfiZ8m14uy/FDiKQB7T/p/cH7bmmyDDd+uJWdB9Jn2OM4+TygE21pyQPa/h7oJsluzsNlWR+WIWlA+wuzC+jhE5Ukw40fr5n90nSLlWFZjuyBRi0soGn6+WhZTt1sIKBdck6fqGQ8zrlYmN4/AW9JrgrYE9CoZZ2FX+fZqu8vy3l8XpqAPlgx694pvRxDmh6/Zb7d+VXdnIWP2jw1Qi/JONBPnnF3XJXhQ7ad+8tynuJhL0FzHqyYt3QBfbAw/R21NP+6dW69XWZ5/6cL6FKuRNr/KMP78+D+siQM6IMV07vGNklAHyxMUwF1JdIzlnItfJpTFEf3lyVhQB+smP6oiByfgT5YmMtD+AwLc3AroK6Ff8pCZmPqNoccy3D0xKufZyfnma0s0QizewuzWfVPIqVYmL2bx+hmY3rGQuYDzbVD0Hni1c8T0Ce2st3vEuxNf7aVnYYx5ViYvZsBNR/oU66mnT7vqzUwI/1xYTaXh70ZuvNgxZxvkWFBOk9tZSk+aH9/tDC92fWzLMxOP6Czvv8TBvTqi09679P834l0WJjzV9XkCeijFbOXKKCPFub0j1ua5NxfmPO18GkWpnM7oL4T6Tkfvnqv/z5N/62ch4XpfyNSmoA+WjE7mQL6cGHWk+7lFPBgYe7N7161OwH1rZwAaQgoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIBShX/+s/vfn3+sfvnzydvGbP+GX/+6eLDfV7/9HX88lk1AqcF69a37v6cCerhtjIBSkoBSgc3q+YAebxsjoJQkoFRguoBeEVBeIKBUQEDJSUCpgICSk4Ayu20Sd74fArr7717VNoffXt62s979+UZzt1lcrX7961jH9fb/Nrtbnj8D3d3ki4DyCgFldh8Cur5o5PuPr4f/3oWvH9C31dHHBB5/86/ngP7X6jKgh5v89v8FlDgBZXaXAf1fxyruQ3fq5/4HvYCe+/nxoH7T+80hoP9j9zBfTmfhT3f+nwJKnIBSgd5noL29y2+Hn+wO0Te7/vVu2x2Cn/5w0cDuB6fYHgJ6vMkhoMf7fPy0AAYRUCrQD+jhI8q3fS83p75tm3cs6bfDfb7s773dSb0Y2/l2vE+393oM6OFRDgG9vgkECCgV6Af0cDS+LVsXtvX5FNHhRtdn4bdpvQjo+vT56ds5oIf77APaO9v/JqDECSgVuDGMaR/Q/jny7U++vF8HdHfY3w9or6eHCveSug/ojZtAgIBSgUcBXX04I9QL6Plk0YeAHpN4HsZ0fNhTQD/eBAIElArcDWjvHPxVQPu/6we0t095L6A3bgIBAkoFHgX0Y976Z+H37fzwGag9UCYjoFTg0SH85eRJ59u+rXrXFH36GehVQH0GSgECSgXuBrQ7Q/T95m17N91cHsL37nM+C38Z0Bs3gQABpQJ3A9rbz/wwjOl80+O4+ZNTE09D7D8G9MZNIEBAqcBxuPx1QHefdO52FtfnS5KOQ5POF2VezCdyLOruLNPtgF7fBAIElArsT6h/vxHQi+vav/Vv2//F/l7Ha5VuXQv/IaDnm/yrk0jECSgV2F8D/+1WQHujlb5d3Pb9OGvT6n8fPtE8BfTGbEwfA3osqNmYeIWAUoNdFb/cDOgxh18+3PaY1m+ny+LPAT3OB3q4eOlWQI8PYxgTLxBQGnYMKIxDQGnO+cr39dXEI1CSgNKc4/n68x9gHAJKc/oXydsBZUwCSns2+sk0BJQWrS9P3MM4BBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQg6L8B8lLzgpV1Bu8AAAAASUVORK5CYII=)
loglike.grid.norm0<-loglike.grid - max(loglike.grid)
plot(theta.grid, loglike.grid.norm0)
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAaVBMVEUAAAAAADoAAGYAOjoAOpAAZmYAZpAAZrY6AAA6OgA6kJA6kLY6kNtmAABmZjpmtv+QOgCQOmaQZgCQkDqQkNuQ2/+2ZgC2kDq2/9u2///bkDrbtmbb/7bb////tmb/25D//7b//9v///9EXbr7AAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3dDZvaRtuYYZKmtds6tfsR2mwbHO///5HlewUIFt2MpLlH53kc7/s4u4AREpclNBpW7wCErOZ+AgBZCShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgTVHdAVQDHlE1X8EQua+9UG2lK8UaUfsKQR/sEAFktAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIChfQNfnaVC+DL+zgALlJAvo29VUUkMbKqCNmGmqMbiUKqA/v968Q37/e9AjeEtVLj4t42vmXm6SyhTQXz8ui/nvn9v//uOfIQ/hnVKPuVo5xNyvEbXLFNDNTS53Sf0+5CG8JeYzdwyLmPtFpDKZArq+PWDf7oQO+hjUO2BSc+dufHO/wswsUUC3u5vfbn74NuwY3hY/trmLNqO5X3pmkCig273N28P1zaPTSLbyqUzdqtrNvT6YioASN3WXUpp7JTGmRAF1CF+LuZOU1dzrjfISBdRJpJnN3Z92zL0mKSZTQPuHMd3ulT5g242YuzfNmnvF8rJMAe0dSD/sUiTb7CBzB2Yp5l7PhGUK6KGYl377a9Aj2FafM0NEWNk680kV0O5UTAcmEyltlm5wae6NgKclC+i76exGM2MwwppamBsTrnyC8gX0JTbKPnOH4taSl/3GdC8Ggwnoos3dhlWtK2TuV+XW3K8IvQR0qYRgkLleritzvwxcEdAF8q5/3dSvYfMvaFICuize4iOY6kVd9ItcKQFdDu/o8Y3/GnvRqyKgi+AtPINxX3SroAoC2jxv2tmNtgqskLkJaNO8TesyxvqwfuYkoM3yzqxX8XVjXc1EQFvk3ZhE4RVltU1OQFvjLZhRybVmHU5IQJvibZdbsfVndU5EQJvhvdaMQqvSih2fgLbBG6xBBVaqlTwyAc3P26ppr69ea3s8Apqcd9MylCin1V6egGbmTbQ0ZfJp/RcjoGl57yxWqYTaEF4moDl5zyCiFRDQhLxVOBHReQloNt4iXFPR2QhoLt4Z3KOhMxDQPLwl+JSITktAk/BO4FkaOh0BTcFbgKFEdAoCWj/bPlEaOjIBrZ2Nntdo6IgEtGo2d8oQ0XEIaMVs5pSkoeUJaLVs4JSnoWUJaJ1s2oxGQ8sR0BrZqBmZhpYhoPWxOTMJDX2dgNbGhsyERPQ1+QK6Pq/EL8PvXP2atwUzPQ0NSxbQt6uVOLShla922y4z0dCYVAH9+fVmFf7+96BHqHqd22qZlYQOlymgv35cFvPfP7f//cc/Qx6i4jVug2V+GjpQpoBubnK5S+r3IQ9R7eq2rVILDR0gU0DXtwfs253QQR+DVrqubaVURUOflSig293Nbzc/fBt2DF/niraBUh0JfUqigG73Nm8P1zfDTiNVuJptm9RKQz8loPOyWVI1CX0sUUAbPIS3SVI/DX0gUUDbO4lkYyQHCb0nU0D7hzHd7pU+UNPKtSGSiIT2yhTQ3oH0wy5FqmjV2ghJRkNvZQrooZiXfvtr0CNUs15tf6QkoZdSBbQ7FdNB0slEbHvkJaEdyQL63sR0drY7cnu+oa1vzPkC+pIa1ufCtziaIKEHAjr9U1ju1kZLJPQ9Y0CTH8IvdkujQRKaLKDZZ6Rf6FZGuxae0FQBTT8j/RK3MFq36IRmCmj2GemXt3WxEMtNaKaAJp+RfmFbFouy1IRmCmjuyUSWtFWxRItMaKKApp7ObkGbFIu1wIQmCmjmCZWXsjmxcItLqIBOYRnbErwv7RsSEwU07yH8ArYjOFtSQhMFNO1JpNa3IbiynIRmCmjOGekb34Cg3zISmimgKWekb3njgUeWkNBMAc04I327Ww58qv2Cpgpovhnpm91u4CmtJzRZQN9zTWfX5jYDA7Sd0HwBfcm0K6vFDQaGajmhAjrmX9bYxgIx7RY0X0DTHMK3tqlAXKsJTRbQRDPSt7WdwIs+TWjK90aqgCaakb6hTQTKaLGgmQI6eEb6+dZROxsIlNNeQjMFdPCM9LOtoVa2DiistYJmCmiayUQa2TZgBG0lNFFA00xn18SGAWNpKaGJApplQuX8GwWMq52CCmhp6TcJGF8rCU0U0ByH8Mm3B5hGIwVNFNAUJ5Fybw0wnSYSmimgCWakz7wpwMQaKGimgNY/I33iDQFmkL6gmQJa/Yz0eTcDmEf2ndBUAa18Rvqs2wDMKHdBkwX0vebp7FJuADC7xAnNF9CXjLk+Eq59qELeggpouYfOtu6hGlkLKqDFHjnZmoeaJN0JFdBSD5xqtUN1UhZUQEs9cKKVDlVKWFABLfS4eVY51CpfQQW0zMOmWeFQs2wJTRTQnuuQdiqYzi7N2obaJSuogBaQZWVDAqkKmiigfd9qXEVAk6xqyCFTQTMFdL8POmj2uhtjrIMcKxrSSFTQVAHdFXTg9EtXRlgFKVYzpJImobkCujuKH/QVHtcmCWjxvwKWJktBkwX0/W3V881yzxthcatfxZBRjoJmC+j2IP6VXdDii1v9CoakUhQ0W0Bf3AWdIKCF/wJYqgwFTRfQ15Re3LpXLuRWf0EF9LWHq3ndQnbV74QK6GsPV+2KhRbcK2gt7zQBfenRKl2r0IrKC5o2oG+hIfUjB7TkgwM7VRdUQOMqXaXQlpoLKqBhda5QaE7FBRXQsBpXJ7So3oIKaFSNaxPaVG1BBTSqvnUJ7aq0oAIaVN2ahKbVWVABDaptRULjqiyogMZUthqhfTUWNG1AYwQU0qqwoAIafJyK1iEsRXUFFdDYw9SzBmFBaiuogMYeppoVCMtSV0EFNPQotaw+WBwBnY+AQnY1FVRAQ49Sx8qDRaqooAIaeZA61h0sVD0FFdDIg9Sw5mC5qimogAYeo4YVB0tWS0EFNPAY8682WLhKCiqgwx9i9pUG1FFQAR3+EHOvM6C3oHM8ieKPWPoBS3p9cStYZ0Alu6ACOvgR9BNqUENBBXTwIwgoVKGCggro0AfQT6jE/AUV0KEPIKBQi9kLKqAD76+fUA8BnZSAQlNmLqiADru7fkJV5i2ogA67u4BCVeb9GFRAh919+wACChUR0Om8trh2QKE+cxZUQAfd+2P9CChUYsaCCuiQO9sBhQrNV1ABHXBfAYUqCeg0Sga03LMCXjLbLqiADrivgEKd5iqogD5/VwGFWgnoFAoGtOCzAl41T0EnCej6vEhfSv9tA70W0P3/vQso1Geeg/jxA/p2tVCzNjS+uHZAoWpNBvTn15ul+v3v0n/l84oFtOSTAgqYo6AjB/TXj8ti/vvn9r//+Kf03/m0VwPa/QNQlRkKOnJANze53CX1e+m/82nhxXUED7VrL6Dr2wP27U7ofB+Dlgpo0ScFFDF9QccN6HZ389vNDd5mPIZ/MaDvjuChYpMXdNyAbvc2bw/XNzOeRoou7n5N2AOFut0GdOT3qoA+eT9H8FC/qQvqEP7J+32sDgGFajUV0JlPIvX8axQPqCN4qN/Eu6CzDGO63SsdR7GAOoKHJKYt6BwD6We8FOn1gDqCh6q1FNBDMS/99lfpv/J58YA6gocUJt0FHX8ykfXVsiScTOS866mfUL22Avqefzq77ooQUKjchAU1ofJT95rqExXgddMVVECfutfqcjcUqJmAjiO0uJ0VoJ+QwGS7oBMFdHOxLMku5fQRKCTTVED3o0FbCKgjeMhhooJOEtDrr0VKGFCDmCCThgK62wGde/zSUWRxOy+/gEIS0xR0ioD+++ecX4N0IR5Qe6CQSVMBnWr6kM8IKCzEJAWd6BA+e0DPf9BPSKKVgL6/1fIRaGRxfQQKKU2xCzpJQHu/2WMWsYCaiQkSaiWg7z+/VnIQH94DFVDIZoJd0ImuRNoWtCPVONDuy6+fkEgjAU19JdL5hbcDCrmMvwvqSqRP7+EjUEiqiYCmvhJpZSo7yGr0XdCJBtJX0s9YQM8vvX5CLo0EtI5z8K8F1A4oZDP2LqgrkT69h4BCWg0E9P0t72QiPgKFzEbeBZ1qGFPWK5G6L7iAQjoNBLSeCe0iAXUOCfIadxd0opNIlxKNA/URKOQmoOUIKCzMqLugAvrJ7Z1DgtyyB7QioYD2/BnIQkCLCQR09fHn0s8GmMCIBZ0koOv5jtmvRAJ6fM3tgEJOyQP660clg5heCaiPQCGr8QpqMpFPbr9yDgmSyx3QvHugF5+ACigkNVpBJ/kMdLNKeimnc0jQgtwB3X0lUh37oIGArs5/HOH5ABNIHdC8A+mdQ4ImnN/Dhd/MAvr45gIKDRhrF1RAH9/84wBeQCGvkQrqSqSHt3YSHpogoEfr80sQGFw6OKBOwkMTBPT99hvmhzZ0eECdhIcWjLMLOllAjx+EvvTx58+v1y/C0McbGtD34i84MIfUAe2cR4on9NePy7vvH3TY+NLAHqiAQgsSB/TiPPxvfwUfeXOTy11SB13jFA6oz0Aht1F2Qaf6Vs5T6DaDdxo/rG/3XofOUxI8C28PFNJLG9DuruPgncaz7T2/3fxw4FfOD1rczousn5DdGLug00yofPXJZWxyu+0db8u7GfaZ6uCAOgkPjbj4RC5RQK92HQfuNJ7NEFAn4aEdnXhmCuhV+QY272yGQ3jfyAntGOEYPlFAJz+J1H2hBRSyG+EYPtEh/J1hTLd7pY+eXCygRjFBA8ofwyc6idQ/kH7Y3uzQgJ7/IKCQXvlj+ETDmHrmxRs8Kn/I4nZeYP2EFuQMaKGB9O/dqZgOxpxMpPMKCyi0IGdAL3cdX5xOebLp7IxigsYUP4bPNJlIAQMDahQTNOX8bs4V0DLT2b1uaEAj9wNq1d37TBXQYqY8hA/dD6hV6WP4ZAGdckb67ossoNCEVdlzG6kCOu2M9GU/bAYqUPgYfqKAXqYv+EHoxDPSn7spoNCKlIfwV0fewYBOOyO9HVBoUNlj+KmuRCoR0GknEzGKCRpU9hh+omvhSwxfmng6OyfhoUGX/cwQ0G35Sgz/nHhCZQGFFq1KHlxONB/ooDnn7hBQ4GVFTyMlCqhDeOBl3XamCGhv+QKmPonUfZmH/CVAvT4+/0yyBxqfg/7StDPSG8YETSo5kGmSgJY7hp9wRvrzSyug0JKSA5mmGUj/82uRfdApZ6S3AwptKjmQaaKTSJfig5qmm5HeOHpo1KrcuztZQN8nm87OSXho1FWPXnuoUk/q/Ig3Pykb0JcIKCxet50JAloRAYXF6+x3ZtgDLcshPPCSch+CJgvodDPSCyi0qtyHoKkCOuWM9J2XVT+hKZ12Liigg2ekv8ntkNeq1Gk6oDbnd3W+PdDtbuRUM9KXCmjo2QK1KnYMnymgk04mYgcUmrXEgE46nZ2PQKFdnW5mC2jYpBMqOwkPDSv1IaiAfn5DAYXGlDqGTxTQiQ/hA3cCclhgQCc9iSSg0LBSH4JmCuiUM9ILKLSs0IegmQI65Yz0Agot6x6/v1DQcQPaM4X8anD0Hj7cWDPSG8YELSv0IWiqgE47I/3q6g9AM1ZlvlkuWUDfJ5vOzvd5QMvO7+2KA3qyOR9r75L64Or1sUUCOuLTAWbS3f2sPKA/v3bOlW+Gfm5ZUuRa+DGfDzCP7js8/jafJKDri6PtdeTY+9r+jPzwb5sffhJJQaFFp3K+tp800ZfKdY/aB159eWGzOowF3QQ/TX12cVdOwkPbLk/MhB+k6HN6HzWgh93O1W9/fUxNP9KM9IaBQuPSBHTbvctD+EGXr1/c8+CP/30s5/rxjPQ9T05AgZ3LD0HDD1L0Ob3f+Qy0swv6Fhp/tHOcSXT/xXLHbm4GntMXUGCvu9sZ3gWd6iz86Vh798foWfi344OsPx7iauf28ycnoMBOkbNI04wDvfw24uA40PN0dtsInw/cR5rOTkChdSU+BJ1oMpHNx/MMjwI9n4vq7naONKGygELrEgX0fFXnC1dxfpzMX08QUMOYoG0lziIlms5uG9DjwPl19xB+rK/0MJAemlbiLFKigG6P3G8/7xx4WZNLOYGDEmeREgV0dyrq+vzTZqwZ6T9e2CEPD+RR4EPQ8aez+347qV3wg9Dd43y5/sFIM9IDrVtYQA/n8j92QtfDx0QJKHBU4CxSqoB25hXdeRvcTwEFTgqcRcr0GejOuhPft+HXhA4/iTT0bwCSKHAWKVtAXzRkGJOCQtte/xB0mslEXvoSpJIGD6RXUGhWjoD2DuB81VvomtAnF7fzYgootOr1s0gTTahc4Es8rowc0MF3AbJ5/SySPdDPbiag0KiLQ81qAzp42uNnCCjwosvz8KEHKP6Uen7WncGzEAEFXvTyWaSJPgMtNJC+Q0CBF12eRQq81QX0s5sJKLSqu+MpoE88OcOYgJOLfkbe7GmvRBo1oAbSwyJcvNOXFNAYl3ICHS+eRhLQu7fUT2iegA4hiEBHp5sC+ikBBTq6+52VBvTmLPxqhEs7nyOgQMdFP4f3YaaArkaYXuQZAgp0dY/gEwV0noIKKNDVPYVUZ0D318IfR23uYvqt+4NpCSjQ9dpp+EkCus3lx7e3b/bpHGWO0M8NHgc66rMB5rbf8aw7oOuLWB7+622WE0kDrkRSUFiA18YxTfQZaHc60M3+WvhNkSvihxp6LbyCQts6u0kC+pnnFve1kWFAIi+dhp/oKz0uD+F3B+91B3ToPYCkXjoNP9FnoBcnkfY5XVf8GaiAwmJ0T3jUGdCfX89zgO6GMe3Owr91mzodAQUuXHwtUpUB3eWy4/u+qBWPAxVQWIyXxjFNNJnIrphHu3BuD+Nn2QEVUODSS+OYJpuN6Xg95+FIfjPP/qeAAldeGsdkOrveWxnGBEvxyjgmAe2/mYH0sBCvjGMS0Du3cyknLMMr45gE9N4N9RMW4ZVzHtMH9OfXWa5BOlBE4NIrA0EFFFi0VwaCCiiwaK8MBPUZKLBsF+OYht61+JMp/YAlCShw5YVxTAIKLNsL45gEFFi2F8YxTXot/G9/vb/NMQvoBwEFrtQf0MNUIruAri6+3mNqAgpcqT6gx6mYtgFdr2YtqIACV14YST9JQH/9WK3++Ofn190cduvVjMNABRS4tuoaetfiT+b2R5v99MmHgM71ZR4HAgpcO4ezzj3QwxfIHQO63R2d70SSgALX4rugE32t8W6f8xjQ7S6oSzmBesQHgk4R0H//3J83OgV0nm+EPxBQ4Fr8NLyAAgtXd0CvDuHXCT4DNZ8yLEfdAd0mc7fPeQzodn/0S+m/82nPfyeSgsJSVB7QzWrXzENAd2NC5xtJP/BbORUUFiA+kn6aK5H2o+f3AX1braY7gl/1eO5uH38c8ekBVah7GNP5Us6D6U4hhQPa/2egTeGR9BNNJrI7cJ+8nz0EFLgR3gWdbDq7Y0JnzaeAAj3CI+lNqPz4RgIK7Qu/5ycZSP+/Lv/713+veyC9gMKy1B3QPy/GLVV/JZKAwrJUHtDOyM/dZ6HVB9QwJliS2gN6uIjzfT+mvvqAGkgPy1J1QN9/fj0W9HAqvv6v9HApJyzJav+N8NUOY9rsvw/psPs569dymkwEuFX5ONBtOn//P3Pvfr77TBPoFSzoVONAN4dnNu+3wgso0KdTzyoDui/o+VTSbAQUuBU9izTdlUhv8+9/CijQp/6A7go630zKRwIK3KozoBfT2H2ofRwosCwC+tyTE1DghoA+9+QEFLhRZ0CrI6DAregEGAIKEJwAQ0ABqv5WzmoIKNCrwks5dyeRvt+eSnISCahL7BheQAGCZ5EEFCA4jslnoAAC+gwBBXoI6DMEFOghoM8QUKBHlQF1LTyQgYA+9eQEFLhV7TAmAQWqV+NA+uoIKNCrwks5qyOgQD8B/ZSAAv0iu6ACChD8EFRAAYKn4acIaM+5+Lm+IV5AgT6xgaAzBXSur4gXUKBPsoDOU1ABBfrUG9D3959fV7/9tf/TLqbfuj+YloACfSoO6DaX387/sdmncxvSOXZBBRToU3FA1xexPPzX2ywnkgQU6FNvQLd7m987/7nZXwu/iV4Rv37lU1QBBfpUPYypVEDfXjwTJaBAr2oH0v/6cXUIvzt4jwT059dXJ3YSUKBftZdyrlcXJ5H2OV0P/wx0G+KLYu6HRw17lCcXNzCpAJBcrZOJ7PYcj+HbRW93Fv6t29QnbW5yuUvq9zu37n9yTy1u5F8iYHmmGQd6+dHl931Rh48DXd8esA8dDfXU4oY+CwGyq/QQ/uLTy104t/uSw3dAt7ubt/cZOBrqmcWNnY0Dkgscek42G9Pxes7DPuQmch3S1cn8g4Hnop4L6LDbAy2IHHomms5OQIHRhA49EwV0ukP4YbcHGhB6489zCB8z0UkkAYUFqjmgnRnt4gntH8Y06GSUgAK9Kg7oxYyg4XnsegfSD+uxgAK96g3ovnyHE0CbwRcPdfTMzDywxoYxAb3qDWj32HvwxUMX1lf9HGUyEQPpYXnqDej66sD7pamUJ5jOzqWcsDjVDmO6Gn80z1TKByYTAfrVOpC+dz7QeWgicEell3KWDagZ6YFR1DmZSMFDeDPSA/VIdRLJjPRATTINY5puRnqAJ2QaSD/ZjPQAz5jhUs7wKaSpZqQHlqnKs/DvZSYTmWo6O2CZho9jSjSd3VQTKgOLFBhJn2hCZQEFxhO5ljNRQB3CA+OJzCaSKKBOIgHjqS+gPRN4vvBB6EQz0gNL1HpAJ5qRHlii1gM60Yz0wBLVF9DiJpmRHligBQT03XR2wCgaH8ZUgoACd7Q9kL4EAQXuqfhSzmIcwgPjGNrPbAE1Iz1Qj1QBNSM9UJNMAR08I33fGNTRnh2wOJkCOnhGegEFxpQpoCYTAaqSKKCmswPqkiigJlQG6iKgAEcND6R3CA+MavDZ5kQBdRIJGNPwi+EzBdSM9MB4AtMxZQqoGemB8QQmBM0UUDPSA+NpPaBmpAdG035A301nB4xjEQF9iYACdwjoZwQUuGM5AT2cTho0hn5PQIE72h7GtI3mqZhrJ5GA0poeSP8R0M65+IE7oQIK3NXypZzngG5Wx+uPdgfyw/ZBBRS4r+HJRE4B3f7v+fKjt4Ej6QUUKCdhQDed699//TCZCDCXhAG92Os0nR0wm4QBXXcDakJlYDYpA9pppoACs0kY0Dd7oEAVEgZ00/0u+LXPQIG5JAxo98z7xoz0wGxyBXS1H0J/ruZuinoz0gNzyRfQvX0291d03n7T8SMCCpSTKKDvnYbuA/o2uJ8CChSUK6AHP7+eTscPnpNeQIFyMgb0BQIKlCOgAEFpAzp0HqYDAQXKEVCAk3YnVL4koEBpQ6ekF1CAg8FfiiSgAHvDv5ZTQAH2hn8xvIAC7AnoJwQUuEdAPyGgwD0LCmiMgAL3COgnBBS4R0A/IaDAPQsaxhQjoMBdyxlIHyOgwH2LuZQzRkCBB5YymUiMgALlCChAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECSjA2ero2VsX//tLP2BJAgo8sFoNKqiAAhydw/lkQQUU4KCTTQHtIaDAXd1APBULAQU4ENDHBBS4S0AfE1DgLgF9TECBuwT0MQEF7hLQrlWPuZ8TUC3DmLoEFBjCQPqHBBR4wKWcjwgo8IjJRB4QUKAcAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCMoX0PXq5MvwOwsoUE6ygL6tLg1tqIAC5aQK6M+vq2u//z3oEQQUKCdTQH/9uCzmv39u//uPf4Y8hIAC5WQK6OYml7ukfh/yEAIKlJMpoOvbA/btTuigj0EFFCgnUUC3u5vfbn74NuwYXkCBchIFdLu3eXu4vhl2GklAgXIEFCAoUUAdwgN1SRRQJ5GAumQKaP8wptu90gcEFCgnU0B7B9IPuxRJQIFyMgX0UMxLv/016BEEFCgnVUC7UzGZTAQYwSkuT964+F9f+gGvmM4OGM3H7tlzty7+95d+wJIEFHjgHM7nCiqgAEedbDYaUIfwwEi6hXimFskCakZ6YDxNB9SM9MCYWg6oGemBUbUcUDPSA6NqOaAmEwFG1XBATWcHjKvhYUwmVAZG1u5AegEFxtbspZwO4YHRNTuZiJNIQFUyBdSM9EBVMgXUjPRAVTIF1Iz0QFVSBdSM9EBNkgX03RwEDocAAAsASURBVHR2QDXyBfQlAgqUI6AAQfkC6hAeqESygJqRHqhHqoCakR6oSaaADp6R/ia3z17gCvCETAEdPCO9gAJjyhRQk4kAVUkUUNPZAXVJFNAyEyoDlFOucMdGlX7AEwEFalOucMdGlX7AkxKH8E9r6li/pYVpaVksTLVmXJiqTyI9zdZQqZaWxcJUq8mAFpiR/mm2hkq1tCwWplpNBrTAjPRPszVUqqVlsTDVajKgBWakf5qtoVItLYuFqVabAX19Rvqn2Roq1dKyWJhqtRrQ9xens3uaraFSLS2LhalWwwGdhq2hUi0ti4WploC+yNZQqZaWxcJUq8WA/vvnSIPm+9gaKtXSsliYajUa0NFOut+wNVSqpWWxMNVqNaBjjZu/YWuoVEvLYmGq1WhA//gfo42cv2JrqFRLy2JhqtVqQP95G3n40omtoVItLYuFqVazAT18sdz4x/G2hkq1tCwWplrtBvQ4kn7svVBbQ6VaWhYLU62WA3r6duMH3yb3OltDpVpaFgtTraYDup/YbtSpRGwN1WppWSxMtRoP6HkvdLT57GwNlWppWSxMtZoP6Pt5VpFphjUBTGDKSzl3UywLKNCMNq6FB5iBgAIECShAUEvn4gAmJaAAQQIKEDRRQN+mm1wZYCICChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBrkQCCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUICglAH998/V1reBv6rUg2f88+vuV4kugv3s1d8u0PcJn85rHi3MZverRMvy+Vb2Zepn9KJ///z9794fT/z+zxjQ9eqgpywPflWp+8/4sC3s9G4qFfrs1d8tUZroPFiYQ3ISrZgHC/Prx2krS7NmdrZPu+/Fn/79nzCgb6c1fvsSPvhVpe4/449+ZlmaT1/9daK36YOFOfdztfrjn1me3FAPFmZ9/lWaVbOz7t3IZnj/5wvobvPdbbeb2+OOB7+q1INnvNsY9oci6yRL8+mrv0n0Ln2wMLt/2fbhPK+h2j1YmNOPdguV45/pnf1u8+3TneP9ny+g69OLs325rvbUH/yqUvefcedw9y3Htv3Zq3/YpU4S0AcL83be8dwk2QV9/JY5/BuQ6eOVwyHA7Xtijvd/uoBuV/TplXu72gF48KtKPXjGm49/RHf/3ta/bX/26u8+tfqfGRZk58HCbJfj9Obs/LFmj9bM+rwESd4yx3+I/1PPZ6CzvP/TBXT7j8uXjz9e7AA8+FWlHjzj7iawztCdz1797QJ9f8uwIDvPbWVJPFozFwHNsW52n5x87zuJNMv7P11AO2G5fhEf/KpSTz7jdYYdnU+WZbtJf0v1Jr23MJsse2pnj9bM+XdJ9qbfd0/5S/8mNsv7P2NAT+/B63X+4FeVeu4Zbw9NEuxOP16W7Y+2y5ApoA+3ssPgnxzL8nDN7I6Hv3T+N4s7AZ3+/Z8uoN3D2as9swe/qtRzzzjHp1OPl+XwkzQBfbiV/fbX29SDZV7ycM18jJbL1M/egM7y/k8Y0I9XpmfTvvOrSj31jJN8nvtwWY7HvYkC+mgr+2+5Buh+spWt840CvRfQGd7/yQN6sdof/KpSzzzjbT8Tvk0vl+X0IUTSgHaf9OHgfb80WYYbP9zKPgbSZ9jjOPs8oBNtackD2v4e6CbJbs7DZVkflyFpQLsLsw/o8ROVJMONH6+Zw9LsFivDspzYA41aWEDT9PPRspy72UBAd8k5f6KS8TjnYmE6/wS8Jbkq4EBAo5Z1Fn6dZ6u+vywf4/PSBPTBill3TunlGNL0+C3z7c6v6uYsfNTmqRF6ScaBfvKMd8dVGT5k27u/LB9TPBwkaM6DFfOWLqAPFqa7o5bmX7edvrfLLO//dAFdypVIhx9leH8e3V+WhAF9sGI619gmCeiDhWkqoK5EesZSroVPc4ri5P6yJAzogxXTHRWR4zPQBwtzeQifYWGO+gLqWvinLGQ2pt3mkGMZTp549fPs5DyzlSUaYXZvYTar7kmkFAtz0HuMbjamZyxkPtBcOwQ7T7z6eQL6xFa2/12CvenPtrLzMKYcC3PQG1DzgT7lZtrpj321BmakPy3M5vKwN0N3HqyYj1tkWJCdp7ayFB+0vz9amM7s+lkWZq8b0Fnf/wkDevPFJ533af7vRDouzMdX1eQJ6KMVc5AooI8W5vyPW5rk3F+Yj2vh0yzMTn9AfSfSc66+eq/7Pk3/rZzHhel+I1KagD5aMXuZAvpwYdaT7uUU8GBh7s3vXrU7AfWtnABpCChAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgFKFf//X7v//+rH67a8nbxuz/Rt+//viwf5c/fFP/PFYNgGlBuvVt93/PBXQ421jBJSSBJQKbFbPB/R02xgBpSQBpQLTBfSGgPICAaUCAkpOAkoFBJScBJTZbZO49/0Y0P1/d6q2Of728rY76/2fe5q7zeJq9fvfpzqut/+z2d/y4zPQ/U2+CCivEFBmdxXQ9UUj339+Pf73PnzdgL6tTq4TePrNf/kI6P9dXQb0eJM//p+AEiegzO4yoP/5VMVD6M79PPygE9CPfl4f1G86vzkG9D/sH+bL+Sz8+c7/UUCJE1Aq0PkMtLN3+e34k/0h+mbfv85td4fg5z9cNHD3g3NsjwE93eQY0NN9rj8tgEEElAp0A3r8iPLt0MvNuW/b5p1K+u14ny+He293Ui/Gdr6d7rPbez0F9Pgox4De3gQCBJQKdAN6PBrflm0XtvXHKaLjjW7Pwm/TehHQ9fnz07ePgB7vcwho52z/m4ASJ6BUoGcY0yGg3XPk2598eb8N6P6wvxvQTk+PFe4k9RDQnptAgIBSgUcBXV2dEeoE9ONk0VVAT0n8GMZ0ethzQK9vAgECSgXuBrRzDv4moN3fdQPa2ae8F9Cem0CAgFKBRwG9zlv3LPyhnVefgdoDZTICSgUeHcJfTp70cdu3Veeaok8/A70JqM9AKUBAqcDdgO7OEH3vvW3nppvLQ/jOfT7Owl8GtOcmECCgVOBuQDv7mVfDmD5ueho3f3Zu4nmI/XVAe24CAQJKBU7D5W8Duv+kc7+zuP64JOk0NOnjosyL+URORd2fZeoP6O1NIEBAqcDhhPr3noBeXNf+rXvb7i8O9zpdq9R3LfxVQD9u8l+cRCJOQKnA4Rr4b30B7YxW+nZx2/fTrE2r/3r8RPMc0J7ZmK4Deiqo2Zh4hYBSg30Vv/QG9JTDL1e3PaX12/my+I+AnuYDPV681BfQ08MYxsQLBJSGnQIK4xBQmvNx5fv6ZuIRKElAac7pfP3HH2AcAkpzuhfJ2wFlTAJKezb6yTQElBatL0/cwzgEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCDo/wO14pp8e1hYGgAAAABJRU5ErkJggg==)
#
# Convert from Log scale to Original Scale of Likelihood
like.grid.norm0 =exp(loglike.grid.norm0)
like.grid.norm0<-(1/dgrid)*like.grid.norm0/sum(like.grid.norm0)
plot(theta.grid,like.grid.norm0, type="l")
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAaVBMVEUAAAAAADoAAGYAOjoAOpAAZmYAZpAAZrY6AAA6OgA6kJA6kLY6kNtmAABmZjpmtv+QOgCQOmaQkDqQkNuQ29uQ2/+2ZgC2kDq2/9u2///bkDrbtmbb/7bb////tmb/25D//7b//9v///8dmM+sAAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3d64Ibx3aYUVBhQtqhQjk2kQi2QInv/5BG49q4DAbYAHrXrlnrxzkkB4Smp7o+9h2zXwCEzLK/AYCqBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYIEFCBIQAGCBBQgSEABggQUIEhAAYLaDugM4Gmen6inv+MTZf+0gb48vVHPfsNnesE/GMCHJaAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgdMLYMj0BpQ+v+IRZeIeA0gcBJYGA0gcBJYGA0gcBJYGA0ofV0Bpdpiag9EFASSCgdGEYWaPL1ASULggoGQSULggoGQSULggoGQSULggoGQSULggoGQSULggoGQSULqxH1vAyMQGlCwJKBgGlCwJKBgGlB5uBNbxMTEDpgYCSQkDpgYCSQkDpgYCSQkDpgYCSQkDpgYCSQkDpgYCSQkDpgYCSQkDpgYCSQkDpgYCSQkDpwXZgjS/TElB6IKCkEFB6IKCkEFB6IKCkEFA6sBtX48u0BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQOCCg5BJQO7MfVADMpAaUDAkoOAaUDAkoOAaUDAkoOAaUDAkoOAaUDAkoOAaW+w7AaYCYloNQnoCQRUOoTUJIIKPUJKEkElPoElCQCSn0CShIBpT4BJYmAUp+AkkRAqU9ASSKg1CegJBFQ6hNQkggo9QkoSQSU+gSUJAJKfQJKEgGlvtGwGmGmJKDUJ6AkEVDqE1CSCCj1CShJBJT6BJQkAkp9AkoSAaU+ASWJgFKfgJJEQClvPKpGmCkJKOUJKFkElPIElCwCSnkCShYBpTwBJYuAUp6AkkVAKU9AySKglCegZBFQyhNQsggo5QkoWQSU8gSULAJKeQJKFgGlPAElS72Azmc7X+7/y6ZXjwSULMUCupgdu7ehplePBJQspQL68+vs1G9/3vUOplePBJQslQL6zx/Hxfz799XvP/91z1uYXj06GlVDzIQqBXR5lsshqd/veQuzq0cCSpZKAZ2f77CvNkLvOgxqdvVIQMlSKKCrzc1vZ3+4uG8f3uzqkYCSpVBAV1ub57vry/tOI5ldPRJQsggo5QkoWQoF1C48lwkoWQoF1EkkLhNQslQK6OXLmM63Sq8wu3okoGSpFNCLF9LfdyuS2dWh40E1xEyoUkA3xTz26cdd72B2dUhASVMqoONHMXmYCFsCSppiAf3lcXacElDS1AvoHc52+GdmV4cElDQCSnUCSpqqAV2fkb/vEqaB2dUhASVNsYAut08AXcaep2x29UhASVMqoJvNztmnH4dH03siPQJKmlIB3Z2A//z/tuWceyI9AkqeSgFdbXcO2Vx/sNy2m0tPpEdASVMpoIvtDvv8cAPSaqfew0Q+OgElTaGA7h9nt9oS3e+4e5wdAkqaQgHdP1B5vNnpgcoIKGkqBnS1Dy+g7AkoaWoFdHvh/Hy8Cy+gH52AkqZQQFd77ufHO+dOIn14AkqaQgEdrl86vWZp6Yn0nAyqMWY6lQI6PE/5y+kfeCL9hyegpKkU0M0d8IeN0Pnxb29hcnVIQElTKqBDQUef4bG4u58mV48ElDS1Arra6hztsi/ufya9ydUhASVNtYA+yOTqkICSRkCpTkBJI6BUJ6CkEVCqE1DSCCjVCShpBJTiTsfUGDMdAaU4ASWPgFKcgJJHQClOQMkjoBQnoOQRUIoTUPIIKMUJKHkElOIElDwCSnECSh4BpTgBJY+AUpyAkkdAKU5AySOgFCeg5BFQihNQ8ggoxQkoeQSU4gSUPAJKcQJKHgGlOAElj4BSnICSR0Ap7mxMDTKTEVCKE1DyCCjFCSh5BJTiBJQ8AkpxAkoeAaU4ASWPgFKcgJJHQClOQMkjoBQnoOQRUIoTUPIIKMUJKHkElNrOh9QgMxkBpTYBJZGAUpuAkkhAqU1ASSSg1CagJBJQahNQEgkotQkoiQSU2gSURAJKbQJKIgGlNgElkYBSm4CSSECpTUBJJKDUJqAkElBqE1ASCSi1CSiJBJTaBJREAkptAkoiAaU2ASWRgFKbgJJIQKntwpAaZaYioNQmoCQSUGoTUBJNEtD5bOfLs/9rdzK1uiOgJHp9QBezY6kNNbW6I6AkenVAf36dnfrtz2f/J29nanVHQEn04oD+88dxMf/+ffX7z389+795M1OrOwJKohcHdHmWyyGp35/937yZqdUdASXRiwM6P99hX22E5h0GNbW6I6Akem1AV5ub385esEjchze1uiOgJHptQFdbm+e768vE00imVncElEQCSmmXRtQoMxW78JQmoGRyEonSBJRMKZcxnW+VTsXU6o2AkinjQvrEW5FMrd4IKJlefSvnupjHPv149n/ydqZWbwSUTK9/mMj8pJ8eJsITCSiZPM6O0gSUTB6oTGkCSiYBpTQBJZOAUpqAkmmigC4beaSyqdUbASXTJAFdXw0qoLyAgJJpkoCefiySgPIsAkqmKQI6bIBmX7+0ZWr1RkDJNEVA//4982OQjphavRFQMk0U0LzHhxwztXojoGSaaBdeQHkNASXTRCeRGjkEamp1R0DJNElAL36yRwpTqzcCSqZpLqT/+bWRnXhTqzcXR9QwM5GJ7kRaFdR1oLyAgJLJnUiUJqBkcicSpQkomdyJRGkCSqaJLqRvpJ9mVncElEzuRKI0ASWTO5EoTUDJNNFJJA8T4TUElExTXcbkTiReQkDJNM2F9M080M7M6o2Akmmik0iuA+UlLg+oYWYiAkplAkoqAaUyASWVz4WnMgEllYBSmYCSapKAzvP22U+YWZ0RUFJNdCdSGxcxmVndEVBSeZgIlQkoqWyBUpmAkmqSY6DLmVs5eQkBJdVkHyrXxjaomdUZASWVC+mpTEBJJaBUJqCkElAqE1BSuROJygSUVAJKZQJKKgGlMgEl1WQB3R4ITb4r3szqjICSaqKAjs4jpSbUzOqMgJJqss9EOvj049n/yduZWZ0RUFJN9amcu5s5l6tfJd6UZGZ15o0BNc5MY6p74ffRHGKad2O8idUZASXVNA9UHh/4TH24nYnVGQEl1USPs/s2+u3isX34+f5gaqDDJlZnBJRUE93KOd5pXz5wIn5xclPovQ01sTojoKQqFdCfX2en7nwrE6szAkqqSrvw65P5x0dT7z2lb2J1RkBJVekk0vIsl3ef0jexOiOgpKp0GdP8fIf93hqbWJ0RUFIVupD+5EjAxp3HA0yszggoqRJu5YyeQjo5F7Vx5xkpE6szAkqqQg8TEVDOCCipCj3Ozi48p94aT+PMNCo9UNlJJE4IKLkqBfTyZUznW6VXmFh9EVByVQroxQvp7zsmYGL1RUDJNVFAj2/CfOA8/Kk7n85sYvVFQMk1TUBPngESP5M0P+mnh4l8bAJKrqnuRHpSQH95nB0jAkquie6Fz/koubMd/pmJ1RcBJddET2PK+ShOAe2dgJJroueB3nWt0QuZWH0RUHIVC+j6bP7oYtB7N25NrL4IKLkSHqj8gN0JpP3pIwH92ASUXJOcRHrwY+T2Difgd+8noB+bgJJrkoA+aR9+2H8fnse0OBRUQD82ASXXNBfSr9L3hG3Q+e6+o+Gezs37CejHJqDkmugk0jMupF/Fcn/sc74tqIB+bAJKrkIBPXqg8ragAvqxCSi5qgZ0KOgXAf3oBJRchR5ndxzQzaNABfRjE1ByFQro+Bjor+3nIwvox/bmeBpoJlEooMPlS+NPlfv5dfbpPwT0QxNQclUK6HAd6Ph60uX9B1TNq74IKLkqBXTzXObRbvxSQD84ASXX9AFdbUeGH243FHS8DTpskwroByag5KoV0PPnkiwE9CMTUHIVC+ijzKu+CCi5Sh0DfZx51RcBJZeAUpiAkktAKUxAySWgFCag5BJQChNQcgkodb09nAaaSbw2oGcPsnvkcXZP+ebMq54IKMkElLoElGQCSl0CSrJpjoEuZ7vPgxuS+v3CKyZiXnVFQEk2SUB/fh3dwb7ctTSDedUVASXZJAGdHz1K/vh30zKvuiKgJJvoQ+XGe+1Lx0B5DgElmYBSl4CSbIqAnnwa3Hz9ie45zKuuCCjJJjoGOtoEXcwcA+U5BJRkU52F3+22D790Fp7nEFCSTXMd6OLoOnrXgfIcAkqyiR4msjzkM3H707zqjICSbLKnMW3v6sw7Ab9mXnVFQEnmcXbUJaAkE1DqElCSCSh1CSjJXv84u+/nD7VzJxJPcWU4jTRTEFDqElCSCSh1CSjJHAOlLgElmYBSl4CSbJqHiSRfPn9gWnVFQEk20ePs8h5gd8y06oqAkmyiByrnPcDumGnVFQElmS1Q6hJQkk1yDHSZ+gi7MdOqKwJKsmnOwv/82sg2qGnVFQEl2UTHQF1IzwsIKMkElLKujaaRZgoCSlkCSjZ3IlGWgJJNQClLQMkmoJQloGQTUMoSULIJKGUJKNlSzsLP0i6rN616IqBkSwroLOnxIqZVTwSUbGkBzSmoadUTASXbZPfCf/qx/tUQ02/jP5iWadUTASXbJAFd5fLb/jfLdTqTnhFqWvVEQMk2zUd6HMVy87tFyokk06onAkq2iY6Bjh8HulzfC79MuSPetOqJgJJNQClLQMk20Ud6HO/CDzvvAsqjBJRsEx0DPTqJtM7p3DFQHiSgZJvqLPxue3O4jGk4C78YN3U6plVPro6moWYC01wHuji6hP77uqiuA+VBAkq2iR4mMhRzawjnajc+ZQPUrOqKgJJtsqcxbe/n3OzJL3O2P82qvggo2TzOjrIElGwCSlkCSjYBpSwBJZuAUpaAkm36gP786nPheQoBJZuAUpaAkk1AKUtAyeYYKGUJKNkElKquD6ahZgICSlUCSjoBpSoBJd2k98J/+vFrkfEU0AOzqiMCSrqJArp5lMgQ0NnRx3tMzazqiICSbpqAbh/FtArofJZaULOqIwJKukkC+s8fs9nnv35+HZ5hN58lXgZqVvVEQEk3SUCX68cnbwKa9WEeG2ZVRwSUdBN9qNxw6mgb0NXmaN6JJLOqIwJKuok+1njY5twGdLUJ6lZOnkBASTdFQP/+fX3eaBfQnE+E3zCrOiKgpBNQqhJQ0iXsws8dA+UZBJR0E51EGrY5twFdbY9+efZ/82ZmVUcElHRTXcb0ZRfQ4ZrQvCvpzaqOCCjpprkTaX31/Dqgi9kscQ/erOqJgJJu0ls5NxJvRDKrevLOYBprXm+ih4kMO+4N9NOk6omAkm6yx9ltE5qaT5OqKwJKOg9UpioBJd0kF9L/+/Hv//k3F9LzOAEl3TR3Ih1dt+ROJJ5CQEk34a2cG8OxUAHlCQSUdBMFdHMT56/1NfUCylMIKOkmOYn08+u2oJtT8T7Sg2cQUNJNcxZ+uf48pM3mZ+rHcppUHRFQ0k10GdMqnb/9/+zNz18mVVcElHRTXQe63NyGlPup8CZVT94bS2PN6012If1y87HGyUyqfggo+aa7E2mRv/1pUvVEQMk34a2ci1nik5S3TKp+CCj5XhvQo8fYtfBAO5OqHwJKPgGlKAEln4BSlICSz+PsKEpAySegFCWg5BNQihJQ8gkoRQko+V5/Eun7+akkJ5F4nICST0ApSkDJJ6AUJaDkcwyUogSUfAJKUe+OpcHm5QSUogSUfAJKUQJKPvfCU5SAkk9AKUpAySegFCWg5HMMlKIElHwCSlECSj4BpSgBJZ+AUpSAkk9AKUpAySeg1PT+UBpsXk5AqemGoTTavFq9gM73l5N+uf8vm1LdEFAaUCygi5NL8u9tqCnVDQGlAaUC+vProzc1mVLdEFAaUCmg//xxXMz1jaKf/7rnLUypbggoDagU0OVZLoekfr/nLUypbggoDagU0Pn5DvtqI/Suw6CmVDcElAYUCuhqc/Pb2R8u7tuHN6W6IaA0oFBAV1ub57vry/tOI5lS3RBQGiCg1CSgNKBQQO3CMyKgNKBQQJ1EYkRAaUClgF6+jOl8q/QKU6obtwyl4ebFKgX04oX0992KZEZ1Q0BpQKWAXvqMuk8/7noHM6obAkoDSgV0/CgmDxP54ASUBhQL6C+Ps2NDQGlAvYDe4dJn0md/TzyJgNIAAaUmAaUBXQf0nBnVDQGlAQJKTQJKAwSUmgSUBggoJd00koabFxNQShJQWlAooBfuQ7r7Y+XMqF4IKC0QUEoSUFpQKKCXPtVYQD8qAaUFlQK63ga96+l1Z8yoXggoLSgV0KGgdz5+6YQZ1QsBpQW1Ajrsxd/1ER6nzKheCCgtKBbQX4vZhU+Wu50Z1QsBpQXVArraiX9kE9SM6oWA0oJqAX1wE9SM6oWA0oJyAX2MGdULAaUFAkpJt42k8ea1BJSSBJQWCCglCSgtEFBKElBaIKCUJKC0QEApSUBpgYBSkoDSAgGlJAGlBQJKSQJKCwSUkgSUFggoFd04kMab1xJQKhJQmiCgVCSgNEFAqUhAaYKAUpGA0gQBpSIBpQkCSkUCShMElIoElCYIKBUJKE0QUCoSUJogoFR060AacF5KQKlIQGmCgFKRgNIEAaUiAaUJAkpFAkoTBJSKBJQmCCgVCShNEFAqElCaIKBUJKA0QUAp6OZxNOC8lIBSkIDSBgGlIAGlDQJKQQJKGwSUggSUNggoBQkobRBQChJQ2iCgFCSgtEFAKUhAaYOAUtDt42jEeSUBpSABpQ0CSkECShsElIIElDYIKAUJKG0QUAoSUNogoBQkoLRBQClIQGmDgFLPHcNoxHklAaUeAaURAko9AkojBJR6BJRGCCj1CCiNEFDqEVAaIaDUI6A0QkCpR0BphIBSzz3DaMh5IQGlHgGlEQJKPQJKIwSUegSURggo9QgojRBQ6hFQGiGg1COgNEJAKeeuUTTkvJCAUo6A0goBpRwBpRUCSjkCSisElHIElFYIKOUIKK0QUMoRUFohoJRz3ygac15HQClHQGmFgFKOgNIKAaUcAaUVAko5AkorBJRyBJRWCCjV3DmIxpzXEVCqEVCaIaBUI6A0Q0CpRkBphoBSjYDSDAGlGgGlGQJKNfcOokHnZQSUagSUZggo1QgozRBQqhFQmiGgVCOgNENAqUZAaYaAUszdY2jQeRkBpRgBpR0CSjECSjsElGIElHYIKMUIKO0QUGoJDKFR51UElFoElIYIKLUIKA0RUGoRUBoioJQSGUGjzqsIKKUIKC0RUEoRUFoioJQioLREQCklNIKGnRcRUEoRUFoioJQioLREQClFQGmJgFJJbAANOy8ioFQioDRFQKlEQGmKgFKJgNIUAaWS4AAad15DQKlEQGmKgFKJgNIUAaWQ6PgZd15DQClEQGmLgFKIgNIWAaWQ8PgZeF5CQClEQGmLgFKIgNIWAaWO+PAZeF5CQKlDQGmMgFLHA8Nn5HkFAaUOAaUxAkoZj4yekecVBJQyBJTWCChlCCitEVDKeGj0DD0vIKBU8djgGXpeQECpQkBpjoBSxYODZ+x5PgGliEfHztjzfAJKEQJKe+oFdD7b+XL/XzaJ6np47Aw+T1csoIvZsXsbag6V9fjQGXyerlRAf36dnfrtz7vewRyqavaEoTP6PFulgP7zx3Ex//599fvPf93zFqZQVc8YuWdEGMYqBXR5lsshqd/veQszqKjnDJyC8mSVAjo/32FfbYTedRjUBCpp9qzyPe2NYK1QQFebm9/O/nBx3z686VPQ7InZe+Z7QaWArrY2z3fXl9dOI52dcrpn8lz6yyQJrTDGlbibV6cnrpubd3z2G+5MG9Cpx4u3xdYXQ8sDbl6Znrx22oUHPopCAXUSCWhLpYBevozpfKv0CgEFnqdSQC9eSH/frUgCCjxPpYBuinns04+73kFAgecpFdDxo5g2PEwEyFMsoL88zg5oRr2APkRAgecRUIAgAQUIElCAIAEFCBJQgCABBQgSUIAgAQUIElCAIAEFCBJQgCABBQgSUICgDxdQgOd5eqOe/YbPlP3DBvry9EY9+w1TdLWv39PC9LQsFqZZiQvTx4/R2tConpbFwjRLQB9kbWhUT8tiYZoloA+yNjSqp2WxMM0S0AdZGxrV07JYmGYJ6IOsDY3qaVksTLME9EHWhkb1tCwWplkC+iBrQ6N6WhYL0ywBfZC1oVE9LYuFaZaAPsja0KielsXCNEtAH2RtaFRPy2JhmiWgD7I2NKqnZbEwzRLQB1kbGtXTsliYZgnog6wNjeppWSxMswT0QdaGRvW0LBamWQL6IGtDo3paFgvTLAF9kLWhUT0ti4VploA+yNrQqJ6WxcI0S0AfZG1oVE/LYmGaJaAA9QgoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgQJKECQgAIECShAkIACBAkoQJCAAgSVDOjfv89Wvt35pUZd+Y5/fh2+9OnH1N9S2Hs//dUCfZ/w23nMtYVZDl8qtCzvr2Vfpv6OHvT377/9efGPJ57/FQM6n21cKMuVLzXq7e94sy4MLq4qDXrvpz8sUZnoXFmYTXIKDcyVhfnnj91aVmZkBqtv+9IPf/r5XzCgi92In/8Ir3ypUW9/x4d+Vlmad3/680LT9MrC7Ps5m33+K+Wbu9eVhZnvv1RmaAbziytZwvyvF9Bh9R3W2+X5fseVLzXqync8rAzrXZF5kaV596e/LDRLryzM8C/bOpz7EWrdlYXZ/dGwUDX+mR6sN5vPv92M+V8voPPdD2f14zrZUr/ypUa9/R2PdncXNdbt9376m03qIgG9sjCL/Ybnssgm6PUps/k3oNLhlc0uwPmcyJj/5QK6GujdT25xsgFw5UuNuvIdLw//iA7/3ra/br/30x+OWv3fCgsyuLIwq+XYTc7RL1t2bWTm+yUoMmW2/xD/rwvHQFPmf7mArv5x+XL45dEGwJUvNerKdzxeBeYVuvPeT3+1QN8XFRZkcNtaVsS1kTkKaI2xGY6cfL90Eill/pcL6Cgspz/EK19q1I3f8bzChs47y7Japb+VmqRvLcyyypba3rWR2X+tyNb0r+Fb/nJ5FUuZ/xUDupuDp2N+5UuNuu07Xu2aFNicvr4sqz9aLUOlgF5dyzYX/9RYlqsjM+wPfxn9fxVvBHT6+V8uoOPd2ZMtsytfatRt33GNo1PXl2XzJ2UCenUt+/RjMfXFMg+5OjKHq+Uq9fNiQFPmf8GAHn4yF1btN77UqJu+4yLHc68uy3a/t1BAr61l/6fWBbrvrGXzeleBvhXQhPlfPKBHw37lS4265Tte9bPgND1elt1BiKIBHX/Tm5339dJUudz46lp2uJC+whbH3vsBnWhNKx7Q/rdAl0U2c64uy3y7DEUDOl6YdUC3R1SKXG58fWQ2SzMsVoVl2bEFGvXBAlqmn9eWZd/NDgI6JGd/RKXifs7Rwoz+CVgUuStgQ0CjPtZZ+HmdtfrtZTlcn1cmoFcGZj46pVfjkqbrU+bbG19qm7PwUcubrtArch3oO9/xsF9V4SDb2tvLcnjEw0aB5lwZmEW5gF5ZmPGGWpl/3QaXpkvK/C8X0I9yJ9LmjyrMz623l6VgQK8MzOge2yIBvbIwXQXUnUi3+Cj3wpc5RbHz9rIUDOiVgRlfFVHjGOiVhTneha+wMFuXAupe+Jt8kKcxDatDjWXYueGnX2cj55a1rNAVZm8tzHI2PolUYmE2Lu6jexrTLT7I80BrbRAMbvjp1wnoDWvZ+msFtqbfW8v2lzHVWJiNiwH1PNCbnD12+rCt1sET6XcLszze7a3QnSsDc3hFhQUZ3LSWlTjQ/uvawoyerl9lYdbGAU2d/wUDevbBJ6N5Wv8zkbYLc/iomjoBvTYwG4UCem1h9v+4lUnO2wtzuBe+zMIMLgfUZyLd5uSj98bztPyncm4XZvyJSGUCem1g1ioF9OrCzCfdynmCKwvz1vPdm/ZGQH0qJ0AZAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoDTh738f/vefP2afftz42nKu2UgAAAOWSURBVJjVf+G3P4/e7PfZ57/i78fHJqC0YD77NvzfTQHdvjZGQHkmAaUBy9ntAd29NkZAeSYBpQHTBfSMgPIAAaUBAkpNAkoDBJSaBJR0qySufd8GdP37UdWW268ev3YwX//6QnNXWZzNfvtzV8f56v+W61cejoGuX/JFQHmEgJLuJKDzo0b++vl1+/t1+MYBXcx2ThO4+8q/HgL6n7PjgG5f8vm/BJQ4ASXdcUD/ZVfFTej2/dz8wSigh36e7tQvR1/ZBvR/rN/my/4s/P4v/08BJU5AacDoGOho6/Lb9k/Wu+jLdf9Grx12wfe/OGrg8Af72G4DunvJNqC7v3N6tADuIqA0YBzQ7SHKxaaXy33fVs3blfTb9u982fzt1Ubq0bWdi93fGbZedwHdvss2oOcvgQABpQHjgG73xldlG8I2P5wi2r7o/Cz8Kq1HAZ3vj58uDgHd/p1NQEdn+xcCSpyA0oALlzFtAjo+R776ky+/zgO63u0fB3TU022FR0ndBPTCSyBAQGnAtYDOTs4IjQJ6OFl0EtBdEg+XMe3edh/Q05dAgIDSgDcDOjoHfxbQ8dfGAR1tU74V0AsvgQABpQHXAnqat/FZ+E07T46B2gJlMgJKA67twh8/POnw2sVsdE/Ru8dAzwLqGChPIKA04M2ADmeIvl987eily+Nd+NHfOZyFPw7ohZdAgIDSgDcDOtrOPLmM6fDS3XXze/sm7i+xPw3ohZdAgIDSgN3l8ucBXR/pXG8szg+3JO0uTTrclHn0PJFdUddnmS4H9PwlECCgNGBzQv37hYAe3df+bfza8Rc2f2t3r9Kle+FPAnp4yb86iUScgNKAzT3w3y4FdHS10rej1/7aPbVp9r+3RzT3Ab3wNKbTgO4K6mlMPEJAacG6il8uBnSXwy8nr92l9dv+tvhDQHfPA93evHQpoLu3cRkTDxBQOrYLKLyGgNKdw53v87MHj8AzCSjd2Z2vP/wCXkNA6c74JnkboLySgNKfpX4yDQGlR/PjE/fwGgIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKECSgAEECChAkoABBAgoQJKAAQQIKEPTffybmBYIGaqAAAAAASUVORK5CYII=)
length(like.grid.norm0)
## [1] 1001
length(theta.grid)
## [1] 1001
# For uniform prior, the posterior is the normalized likelihiood
plot(theta.grid, like.grid.norm0[,1],type="l",
xlab="theta",
ylab="Density",
xlim=c(.3,.6),
main="Figure 8.10 Posterior Density\nHardy-Weinberg Model theta" )
# Compute 95% posterior predictive interval (numerically)
index.quantile.025=which(cumsum(like.grid.norm0*dgrid) >.025)[1]
posterior.llimit=theta.grid[index.quantile.025]
index.quantile.975=which(cumsum(like.grid.norm0*dgrid) >.975)[1]
posterior.ulimit=theta.grid[index.quantile.975]
abline(v=c(posterior.llimit, posterior.ulimit), col='blue')
abline(v=x.theta.mle, col='green')
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABUAAAAPACAMAAADDuCPrAAAAhFBMVEUAAAAAADoAAGYAAP8AOjoAOpAAZmYAZrYA/wA6AAA6OgA6Ojo6ZmY6kJA6kLY6kNtmAABmOpBmZgBmZjpmZmZmkJBmtrZmtv+QOgCQZgCQkDqQkNuQtpCQ27aQ29uQ2/+2ZgC2kDq2/7a2///bkDrb/9vb////tmb/25D//7b//9v///9KAko/AAAACXBIWXMAAB2HAAAdhwGP5fFlAAAgAElEQVR4nO3dbYPixpqY4R6fOHO8u9l2djdpJ9nEbEyM7fn//y+8CCi9gqqRHlXNdX04x9MNoiTB3RIS4u0bAFneogcAUCoBBcgkoACZBBQgk4ACZBJQgEwCCpBJQAEyCShAJgEFyCSgAJkEFCCTgAJkElCATAIKkElAATIJKEAmAQXIJKAAmQQUIJOAAmQSUIBMAgqQSUABMgkoQCYBBcgkoACZBBQgk4ACZBJQgEwCCpBJQAEyCShAJgEFyCSgAJkEFCCTgAJkElCATAIKkElAATIJKEAmAQXIJKCl+/Pnt44vv15/fP6v5R0uj/vj7xO32R9v8MNv6U925zu9D924O08f+UMbnP6wzyyyQ3vAk4tiru645swTCxPQ0oUH9PBMN87DSQO6n6pjf57yEvrHT8OBnhjjawLa/WPxKe1xzZsnFiagpYsO6D596NGC7jpR2U3GsT9POQX965e3kS3cYS8N6Gc2myfGNXeeWJiAli44oMcNotTX4VvtO1tlh954WwYCmjEr50jHBfRlBU3HNXeeWJiAlm7NNzsH7K+tOG8bjey57t7aAb3c9uu1lL3qtuZpP5nmCWvG5pDO3yG7+Q8J6MYIaOliA3pJ4UfnP9tu25O3gJ4Lc97dHx59712/0TRPCQvoNfrzm/+QgG6MgJZuAwFtwrEfDOilq+2AJhkYLEJnnnZ5e8SBAb089gJrRUA3RkBLNxLQzo8vW3HvzX+cN432t83A1tbj/nLDffJKbTYhh163DwOavp15Dcz5Ts3YDkObap3BH9oBbZLc3iS9vhfb2Y1O79mbjc6sdpdk/3EGlk06xOSWt+U98tCH9DFbx966M5KMK52n//5zOsHLZF522IqnCWjpngnorWJffn0yoLt7epIEDuyT3quZtrQzjrf3VmDOY2j+cf7v7sH7qYAmR626fx/S8vQCOjAbnVkdeuOg/Ti9ZXPTDehlgTbzNfDQl4Def3EtYX9GRgL6sUseoLVEWZOAlu6JgCYv4C///lRA//H+Em4dZu+fpnTb9rlMoreVevr98V7TAe2+9CcC2jpr6vpw7VMBzlPrBnRoNjqzOnToqv043WVz1w1oet7W0EOfb/9PP3UfYWBGxgJ6uelHMrAF3nLlEQEt3eOAJu9CNh4G9P6q7pxR1H+Rtl7z/V//+fM5Df2j1OnWWXf4Q++BXv7ZPV3ofWgGv37rBXRwNoZmNX1nofs43TskegHd38Y8+NDDJ94PzchYQG9nMnybOH7H0gS0dL1zJt/vP06jc3qFXm/7ZEBPd9+9XX9z/6+O3f1VPTrKVmDSR34ioPfxp+8lHm5DTCZ+mY+BcyYHZ6Mzq8mjDj9O9w5j83cf9MfYQydnOjV/gT5GZ2TkPND9/SHtwYcR0NI9DGh6gmazifNEQL9+60xmdD8x2QQd/yjndEC75W0FdJfM1i5p1/0c0uTo1WlGWnfs/TlJZqNzhmlyo+HHmTgldTygww99SJbXvbXDMzIS0GTJDR6JYw0CWrqHAW29WXb5x+OAXlOYvjR7n2e/3+Jq9MSdjIC2JRug7+k0T484ctZlEpvh2WjParLIRh6ne4ex+bsvlo+xh07f1r0/2vCMjH0SaXe78W5gGbIKAS3dw4C2XtvJO2eTAb2+GtMXbHr20c3tfJ2pTyJ9PqDJmxGdzzO9pw1vvTGZjH14NtqzOrbIksNjnTuMzd/tB6ebDj/0Pcrp9IdnZCygt8c83+ClF4DiWQJauocHkdJDNulmy1RAr9PrH3/q5CM59p4e1ejJeA80lZ5pei/FdV5ao7xX7B6bkdlozWo6lJHH6d5hbP5uPzg+yshDp7e/L8PhGRkL6G2lnf+IObs+hICW7mFA2zXovAM4EtD2BmtL+4V6fu2mR65HNkEzjsLfDQy+/c/22T/990BHZqM1q71FNvA4nTuMzd/tPseJjTx06/b3gQ7OyOjFRK7j2g8+A1iDgJZukYBeb/4woA9bmN5uKKDPnAd6Mxa2zjh7sZkK6H1D88mADu4rDwf0+IOpgF6nlFRxaEZGA9osu/N97MHHENDSzQzobjCgyTuRAwGd2Dt8+HZm41OfRBp8tNa8fEvfP7zepBvQ3mxMB3RsF/6pgN4+iTTy0GMBHZqR0YA2f/bOv7cHH0NAS/fce6APDiKNBfThxSsevp3ZaA2id8b6g8/CD09l4PTxXdOdy02SsQ/PxmhAxx7n+YDePws//NATAe3NyPj1QC8L7zC+2FmagJbuYUDHTmNKX/PJWTUDe/wTu4fp0eT0vwdvl3s1pvbPB04vOs7W9dbJD5NJD8/GaEAnT2N6KqD7B/cZC+jQjIwH9HJa1L9MryKWJKClexjQsRPpk2g2Px4IaHpVocPIR91b70OOvJLbgZl3PdDELmn07aGTMbTeHmh1aWg2RgM6eSL9MwE93Bfz8EMPBnRkRsYDev8QmD34IAJausefhR/5KOfllX37DONwQJuX6P3T5d183H7fRPh+MLk1pnZghtv3cJ5uR6k7H7FMzm+/TDB5D7QZ7+BsjAd0+HGeDWizKNLmdR96eAt0eEZ6Ab2P4HrU3h58FAEtXfbFRPo/Hgpo5zuPJq/3ebvjo4C2r3U0fSmptvZFkvpXqksHebvt+8hsjAd0+HEeBXTwPsMPPRzQ4RkZHFf7oL09+CgCWrrHAU1fl+c3zPpXBPryr2NV6aShv6vY6kN6aH0qoM98K+fwVtUTl7Pr/bT9IZ/7LSYCOnE5u2cDOnhYvfXZ0v57oIMzko6rNU/3YdqDjyKgpXsioPeC/vDb+cXa7DPfXtlffr1fxqIXiXSzaOgcpeT3rev7TgY0CdTYtyiN7JYOXlC5tTV9m+Cu+cHXkdmYCuj4BZWfC2g6rwMP/dR5oOmV9DuXVmm9U2oPPo6Alu6ZgCZJ2yWvvuRgyURAbwEYfZk2r/rOfv90QK8pGNx2mgroyFd6XMvSulc7L73ZmAzo6Fd6PBPQwdMK0olPnMbUm5GhsN/+uRsdEmsQ0O/Mbuj1TbGsz1gCWr3jhlSyWeUFV5XpbXUWJ6DVa12l3QuuLo8+6MDCBLR+vbMLfflDLfbj7yOzCgGtX/dERC+4KtzOY/D3MJCAfgc6Z4Xb46vCwd/DDRDQ70HrJBtfPlaHg35ugIB+H7qfYKF4l1VqdyKWgAJkElCATAIKkElAATIJKEAmAQXIJKAAmQQUIJOAAmQSUIBMAgqQSUCr0P0Sov6XEj1v9+QVRy4Xam5/vdL90mqTV/p9fFnn/ZKXPRm4MMD1glXPL7MHi2liFg/Nwx7mfpb94EIGmyOgVYgI6CU612I238DW/rKzsalsJaD3B9itFtDjg+cF9H5HtkNAqxAS0ENanGuUmtf45dL3Q9+CfPvtBgJ6216+ffHw0gE9/53JCWh6R7ZDQKsQEtBWJK/7wOn3lY8OYCsB7X8J/NIB3SV/ZGYFNL0j2yGgVQgJaOtN0Os+cLNRN/A98LOsE9BbkW7X7BdQZhHQKoQEND1Q1N0Hfn4iE5NeNqD/mDzCrvm3gDKPgFYhJqDJZmZno+6ycfqJ1/sKAf0vP9+2kU+1+/KvAspsAlqFBwHdp3G7/eD9+n/3u5ze0WzKcNmmvN5j8KBQ8ijnKf3brSntfDSbp/eXf/LrQ/Oby21uObkG9PLFP63H7U1tdGbem//oZu7y23/+5Tbl06P8+B+dgDbnFXTeh+gupidm8Sb9aqqPe0A754MNTK5zx9tMa2o4Aa3CZEDTL+V8T370vru/JG+74D/8di3DLi3aeYLd9zST7czTjb/875+vN0pvfpv0PRK9gN5vkw7w68DB8YGpjc/Ml18nAnq+x/ttCl9bJxWk3wWdLNeBxfR4Fm+GA7rrzPrA5HoBHVilhBDQKkwFtP2lxvdiXN70u3QuzcU/NC/c27bW7Q69HdZbRM4v+UsLzo+a3Lz1rfRNj7sB/afkNskA//M9JN2zpdKpdWcm6c+Xf58I6OE6ifMfgo9WQIcW2vBiejyLN4MB/R+tnw1PrhvQ4dERQECrMBHQ1svxulnY2YJp9lYTX9sH2S//3Tut8/Yu3uFy09te9O72wk5qdp9cN6DTA0x+PDi1Z2amv7TePk7/d57s+S3QX9OAdkf1niyF3pQfzWJ7cTU+Bma9czzuPrnOHQdXKSEEtAqdl9TF5eW7u77GLi/M5C3L+632t/+6NuLr9UbJMaL+C/XW6f3lpX3dYU56u7u+6pP/Ggjo6R/NTKT7qMkG5fvo1Dozc+jecyygf13fBD2cb58ENHnIQ7IohxfTw1lM7W6zcpv17roZnFzrjoOrlBACWoXxgCYv5GSffN8qyyUHTR93t98kB44OQx263+I8geP9/2zeBL33Nu3Ifb++F9DLtldyxGif/Pg+vOGpjc9M07mRgJ5n9b2ZwNf0g1W7pEuXefz6YDFNzWKqH9DLPN7XzfDk0jsOr1JCCGgVJrZAT6/7ZK+y36fOy/Dyj/um1f2/Bl6nzcu6eQu0OZb0a3KCTlrea14HAvqRDPD9NsBrf25pG57aczPTW1of5ztej4K/JwFtn4Bw+/nIlB/O4sACu0+4md79eNzw5NqnMQ2tUkIIaBWmAnrTC+j1hdf62FDy1uft57dAdjWn0l/PReqfUZS+7s9T7m4Wp+88Jkf1Wxdzug18eGrPzUxvaV3ecThO4fIWaDKS9seobqMamfLDWUz1Atqb9eHJjZwHKqDRBLQKz5xIv3u7h6a1gde59Nzu1pzbW5nnyQ0d7G321ffXCV/eBL2/Bdo/7NLdTU2r1A3oNXzXn49MrT0z7dPTd+l0WqM+3vU8jI/mLdCkZ51L8V2nMbyYHs9ibyWkAb0G+fqLkckNBzRdpYQQ0Co8eSL97dXWuhJdpwxJu67/uX9rv+n2dn01XzLxf365Tu3yZuj//bn177b73XoBTSrRPm2q+fnI1J6cmc7Sum7Rfr3uyd+XWedO6VLoT/nxLKbSDrZSf/3FyOR6Ae2tUkIIaBWmAto5W+Ye0HYLhprTbGCeX9TpwY7b6/rycv9vP91+f3qdX06+vJ9hOVmXwYpkBPTxzAwsrXM6/19zobiCAjq0SgkhoFWYCOiueZX9+Hv3PdBbC0b3ept98fMd3+8TTl7X+/Y/L5P6x9u/Rz4T/2RArz9u78L3pjYV0GRmBpbW6cEvuf/oBbS/Cz+8mB7PYurJgPbfLknvOLhKCSGgVRgPaO/sxrGADh53uRwSPtwn3g3obVuo9/v+6Yt3Twb0OqTBg0h3z85Mf+mcJ/xvb8mnsUYPIn08dRBpcBZTDwM6dtWQ5MfDq5QQAlqF8YAmL7zRgLZehq0zfy7n0fxLeuO25MPh53/f9kDTMxl7930yoK0h3c7X702t/dNnT2NKz4c/3/u+zEZOYxpZTA9nMfU4oCNfJpXccXiVEkJAqzAa0HSPMIlS+1U6ce75rinM2Aeubwnq3qGZ+O2ySN/OwUi/8OOJgKZDujWuN7WnZ6a/tNJrGnXf9pg8kT6Z8sNZTD0O6PDkkjuOrFJCCGgVHgQ0+WDgUECTzyt2P/143SUf/bjgNZjd4xvtoN4/0H1+0GcDmg7pfXRqgzPz8KOc9531t+QUrNa5Bt2Pco4spkez2FtczVhHZn1wcskdR1YpIQS0CuO78L2rcvSP0XxLNjSv2qdgTlz3N/ks+9k1LdcXdfcU/8HPwo8G9K53EsB9ap2ZefJiIueHb6bXO/W9+/AfU4vp0SymbhN+H531wcmldxxepYQQ0CqMBzS9ts9/ul3brRvQ9Fb/1LoC3P7Ba7R5uXcOuSRjaZ9y0/s091RA/+1+v+v0h6Y2MTPnt28nAtqMtr/NOHzBuJHF9GAWBxbXVEAHJ5fecXiVEkJAq/DU5ew+Oh8fT7cq75tt7+1rEF/uPn7Bn97bjO1Ppn9rV+cj+dHjgH7cUnK/ENTA1HozM3LZ4/bSST5iPnT5j+ELKo8spulZbNk1N+ueFnWf9cHJte44uEoJIaBVmPwk0iVDrVNkBg713g6WdC7ivuvftGXXfp03jzb0BRX38T0b0O7HNIenNjAzt/l/FNDkgh7d5A1/pcfIYpqaxYEHv18aZSigA5P71v5TNrBKCSGgTCv7FVr26Nk8AWXS6JbUVh03HJPhCiiLElAmjZzXvV2XPe/0JM6S8k9hBJQp+7fiNuGSo1iXtxIdY2ExAsqY26k8ZRWod3HpovJPWQSUMYdCC9Q50byk9x8ojYAy5lBmPztnojvLnAUJKGMu+8JFbsClH/iB5QgoQCYBBcgkoACZBBQgk4ACZBJQgEwCCpBJQAEyCShAJgEFyCSgAJkEFCCTgAJkElCATAIKkElAATIJKEAmAQXIJKAAmQQUIJOAAmQSUIBMAgqQSUABMgkoQCYBBcgkoACZBBQgk4ACZBJQgEwCCpBJQAEybTugbwAv8/pEvXyKLxS9tLn4299eNaW///1FE3rdkPiuvLxRr57gKy3wB4MMf/vbq6b097+/aEKvGxLfEQElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElgIBSBwElwEsD+pp1KqBkEFACvKxWl2/lfMWUBJQMAkqAV9XqTUAJJaAEeGVAX7RWBZQMAkqAF9XquDovAX3BahVQMggoAV4b0JesVgElg4AS4DW1Oq3Ny2lML1ivAkoGASWAgFIHASXAS2p1XpkCSiQBJcCrA/qCFSugZBBQAggodRBQAryiVpd1KaBEElACvDygn1+zAkoGASXAC2rVrEoBJZKAEkBAqYOAEuD1Af30qhVQMggoAT5fq+uaFFAiCSgBBJQ6CCgBPl2r24q8f6XHJ9etgJJBQAkgoNRBQAkgoNRBQAmwREA/uXIFlAwCSgABpQ4CSoDP1uq+HgWUSAJKgEUC+rm1K6BkEFACCCh1EFACCCh1EFACCCh1EFACfLJWyWpMA/qp1SugZBBQAggodRBQAggodRBQAggodRBQAnyuVulaFFAiCSgBFgroZ9avgJJBQAkgoNRBQAkgoNRBQAkgoNRBQAnwqVq1VqKAEklACSCg1EFACbBUQD+xggWUDAJKAAGlDgJKAAGlDgJKgM/Uqr0OBZRIAkoAAaUOAkoAAaUOAkqAxQKav4YFlAwCSgABpQ4CSoBP1KqzCgWUSAJKAAGlDgJKAAGlDgJKgOUCmr2KBZQMAkqA/Fp116CAEklACSCg1EFACSCg1EFACSCg1EFACbBgQHPXsYCSQUAJIKDUQUAJIKDUQUAJIKDUQUAJkF2r3goUUCIJKAGWDGjmShZQMggoAQSUOggoAQSUOggoAQSUOggoAQSUOggoARYNaN5aFlAyCCgBcmvVX38CSiQBJYCAUgcBJYCAUgcBJYCAUgcBJYCAUgcBJYCAUgcBJYCAUgcBJUBmrQZW31BAs1azgJJBQAkgoNRBQAkgoNRBQAkgoNRBQAkgoNRBQAkgoNRBQAkgoNRBQAkgoNRBQAmwcEBz1rOAkkFACZBXq6G1J6BEElACCCh1EFACCCh1EFACCCh1KC+gu7err/PvLKDbIKDUobCA7t/a5jZUQLdBQKlDUQH946e3rh9+mzUFAd0GAaUOJQX0r1/axfzz5+O/f/x9ziQEdBuWDmjGihZQMpQU0EMvl6ekfsyZhIBuQ1atBleegBKppIDu+jvsx43QWW+DCug2CCh1KCigx83N994P9/P24QV0GwSUOhQU0OPWZn93/TDvMJKAboOAUgcBJYCAUoeCAmoXvh4CSh0KCqiDSPUQUOpQUkCHT2Pqb5VOENBtWDyg89e0gJKhpIAOnkg/76NIAroNObUaXncCSqSSAnopZtuXX2dNQUC3QUCpQ1EBTS/F5GIiJRNQ6lBYQL+5nF0VBJQ6lBfQGXo7/G8Cug0CSh0ElAACSh1KDej5iPy8U5hOBHQbBJQ6FBbQQ3MF0EPe9ZQFdCMyajWy6gSUSEUF9LLZ+fbl1/ul6V2RvkjLB3T2qhZQMhQV0OsB+B//oynnzhXpyySg1KGkgB63O0/ZPH+xXNPNgyvSF0lAqUNJAd03O+y7+weQjjv1LiZSIAGlDgUF9HY5u+OW6G3H3eXsiiSg1KGggN4uqJxudrqgcpEElDqUGNDjPryAlm1+rcbWnIASqayANifO79JdeAEtkIBSh4ICetxz77/fuXMQqUQCSh0KCujp/KXuOUsHV6Qv0goBnbuuBZQMJQX0dD3lr90fuCJ9iQSUOpQU0Msn4O8bobv2P58hoNsgoNShqICeCpp8h8d+dj8FdCMElDqUFdDjVmeyy76ff016Ad2G2bUaXXECSqTSAvpJAroNAkodBJQAAkodBJQAawR05soWUDIIKAEElDoIKAEElDoIKAHm1mp8vQkokQSUAAJKHQSUAAJKHQSUAAJKHQSUAKsEdN7aFlAyCCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBZtZqYrUJKJEElADrBHTW6hZQMggoAQSUOggoAQSUOggoAQSUOggoAQSUOggoAQSUOggoAQSUOggoAQSUOggoAVYK6Jz1LaBkEFACzKvV1FoTUCIJKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAHWCuiMFS6gZBBQAsyq1eRKE1AiCSgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgB5tRqep0JKJEElACrBfT5NS6gZBBQAggodRBQAggodRBQAggodRBQAggodRBQAggodRBQAsyo1YNVJqBEElACCCh1EFACCCh1EFACCCh1EFACrBfQp1e5gJJBQAkgoNRBQAkgoNRBQAkgoNRBQAnwfK0erTEBJZKAEkBAqYOAEkBAqYOAEmDFgD67zgWUDAJKAAGlDgJKAAGlDgJKAAGlDgJKgKdr9XCFCSiRBJQAAkodBJQAAkodBJQAAkodBJQAawb0yZUuoGQQUAIIKHUQUAIIKHUQUAIIKHUQUAIIKHUQUAIIKHUQUAIIKHUQUAI8W6vH60tAiSSgBFg1oM+tdQElg4ASQECpg4ASQECpg4ASQECpg4ASQECpg4ASQECpg4ASQECpg4ASQECpg4ASYN2APrXaBZQMAkqAJ2v1xOoSUCIJKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAEElDoIKAFWDugzExJQMggoAZ6r1TNrS0CJJKAEEFDqIKAEEFDqIKAEEFDqIKAEEFDqIKAEEFDqIKAEEFDqIKAEEFDqIKAEEFDqIKAEWDugT0xKQMkgoAQQUOogoAR4qlZPrSwBJZKAEkBAqYOAEkBAqcMaAf3zv/7+6gfJJaDbIKDUYZWA/vz29vXVD5NHQLdBQKnDWgE9en/1I2UQ0G0QUOqwznug+7eLj1c/2FwCug0CSh1WO4i0uyT0h99e/XizCOg2CCh1WPEo/F+/XBr6Y+AhJQHdBgGlDuuextS8Gxp3SElAt2H1gD6emICSYfXzQJvt0C+/vvpxnyKg2yCg1GHlgF6PJkVthQroNjxTq+fWlYASac2A7m/Hkc6boZmnNe0+k2AB3QYBpQ6rBbRzJtOfP2cdkE+2YLMaKqDbIKDUYd3zQJONzl1GQP/46a1r5kQEdBsElDqs+Emk1ulLx5342acznXf8k2KepztvKgK6DQJKHdYK6CuOuh96uTwlddanmwR0G9YP6MPJCSgZ1gnoaz7CObDXf2zzrLdBBXQbBJQ6FHQ90OPmZv/A/X7ePryAboOAUod1tkBbO/C7zA9zHrc2+1uyh3mHkQR0GwSUOoQENO+KIgJaDwGlDusHNPMMULvwNXmiVk+uKgEl0sIB7Z74/qnrMTmIVA0BpQ4LB/R6+aW2zA9xDp/GNGtiAroNAkodlt6FP7yun8Mn0s97O0BAtyEgoI8mKKBkCDiIlG9ge3bmlAV0GwSUOhQV0PRSTC4mUjIBpQ4FnUjf2GXX85uAboWAUofyAjrD0AGs6DFxIqDUYdmAnt60/Oi/dbnWV3MK6FY9rtWza0pAiVR1QAcGJ6CbIKDUobCAni+pnJwM+tcvPspZIgGlDmW9B3o9gHQ7fCSgZYoI6INJCigZigro/QD8dSNUQMskoNShpICe9t9P12Pa3wsqoGUSUOqwbkBPCcy8ksi30wZoc0L+6TOdl8kIaJkElDqsFNA/fz4F7/KtmrkfSzrG8vbe564pqICWSUCpwzoBPZyPvJ8vBpJ/EL51QeWmoAJapoe1enpFCSiRVgnoacvzuN15/L9j7fa5l2NqX5F+dz4WL6BlElDqsEpA95cTj5p07nI+xP6tG9DLpUAFtEwCSh3WCOixcs3u9vntz5lfY5RO5mv7n28fAlqmkIBOT1RAybDO5ezOm47XL0PKDehpCzb9Vrk/fnr78n8EtEgCSh1WDOih+QRRdkBPb6Wmb58e5h+SEtBtEFDqsGJAd80G5Mwv0kzsO1cBPQhooQSUOqz0Huj7/S3QuV+kmToVNN0GPW2TCmiBBJQ6rHIU/vwRov3lzM3r5zHz9L4bfi+gJXpUq+fXk4ASabXzQE8+LtcDyf8w56cJ6DYIKHVY55NI+9sllPaRl1MW0K0QUOqw0mfhD2+3Q/D5++8vIKDbIKDUoaTL2b2AgG6DgFIHASVATEAnJyugZBBQAggodVgpoNfj8L6Vk5MHtZqxmgSUSGsehRdQGgJKHVYJ6MH3wtMioNRhpU8ihZ78mRDQbRBQ6rDSZ+E30k8B3QgBpQ4rXY0p70s8Xk9At0FAqYOAEkBAqcN6l7PbBAHdBgGlDit9qVzgBZhaBHQbggI6NWEBJcMqAT3HmS0AABZ4SURBVN3OPryAboOAUod1TqT/46eNbIMK6DZM12rOWhJQIq10EMmJ9KQElDoIKAEElDoIKAEElDq4nB0BBJQ6CCgBBJQ6CCgBBJQ6rBbQ0xuhp2+Hjz2dSUC3QUCpw0oBvRxHOgX0LfRrOQV0G6ICOjFpASXDOgFtjsMfA7p7Cy2ogG6DgFKHVQL61y9vbz/+/sdPx4AGX1xZQLdhslazVpKAEmmtr/R4P32e8xTQ00583AfjBXQbBJQ6rPSVHqdDR01Aj5ujcQeSBHQbBJQ6rHg90Cagx01Qn0T63gkodVjpo5yn40bXgB4E9LsnoNRBQAkgoNQhYBd+5z3Q756AUoeVDiKdtjmbgB63R7+++jGfJqDbEBbQ8YkLKBnWOo3p6zWgp3NC486kF9BtEFDqsM4nkc5nz58Dun97i/x2DwHdBgGlDqt+lDP6csoCuhVTtZq3jgSUSCtdTOS0476BfgroRggodVjtcnZNQkPzKaBbIaDUwQWVCSCg1EFACSCg1GG1y9mdxZ0A2hDQbRBQ6rBCQPetrzSOvB69gG6FgFKHxQOaHH+/CP1SJAHdhriAjk5eQMmwdEAv/bxudu6jCyqg2zBRq5mrSECJtHRAu9+BdIh9J1RAt0FAqcPCAf3jp7fLJZgmfrImAd0GAaUOCwd04AuQfCcSAkollg3oX7/0P3rkcnYIKJVYPKC9Q0a+VA4BpRLLBnRwa3PnKz2+ewJKHQSUAOO1mruGBJRIAkoAAaUOAkqAwICOPYKAkkFACSCg1EFACSCg1EFACSCg1GHxgA4R0O+dgFIHASXAaK1mryABJZKAEkBAqYPvRCJAZEBHHkNAySCgBBBQ6iCgBBBQ6iCgBBBQ6iCgBBBQ6iCgBBir1fz1I6BEElACCCh1EFAChAZ0+FEElAwCSgABpQ4CSgABpQ4CSgABpQ4CSgABpQ4CSgABpQ4CSgABpQ4CSoCRWmWsnpyADj6OgJJBQAkgoNRBQAkgoNRBQAkgoNRBQAkgoNRBQAkgoNRBQAkgoNRBQAkQHNChBxJQMggoAQSUOggoAYZrlbN2BJRIAkoAAaUOAkoAAaUOAkoAAaUOAkoAAaUOAkoAAaUOAkqA6IAOPJSAkkFACSCg1EFACSCg1EFACTBYq6yVI6BEElACCCh1EFACCCh1EFACCCh1EFACCCh1EFACCCh1EFAChAe0/2ACSgYBJYCAUgcBJYCAUgcBJcBQrfLWjYASSUAJIKDUQUAJIKDUQUAJIKDUQUAJIKDUQUAJIKDUQUAJEB/Q3sMJKBkElAACSh0ElAACSh0ElAADtcpcNQJKJAElgIBSBwElgIBSBwElwAYC2n1AASWDgBJAQKmDgBJAQKmDgBJAQKmDgBJAQKmDgBKgX6vcNSOgRBJQAmwhoJ2HFFAylBfQ3dvV1/l3FtBtEFDqUFhA929tcxsqoNsgoNShqID+8dNb1w+/zZqCgG6DgFKHkgL61y/tYv758/HfP/4+ZxICug0CSh1KCuihl8tTUj/mTEJAt6FXq+wVI6BEKimgu/4O+3EjdNbboAK6DZsIaPtBBZQMBQX0uLn53vvhft4+vIBug4BSh4ICetza7O+uH+YdRhLQbRBQ6iCgBBBQ6lBQQO3C10NAqUNBAXUQqR7dWuWvFwElUkkBHT6Nqb9VOkFAt2EbAW09rICSoaSADp5IP++jSAK6DQJKHUoK6KWYbV9+nTUFAd0GAaUORQU0vRSTi4mUTECpQ2EB/eZydlUQUOpQXkBn6O3wvwnoNnRq9YnVIqBEElACCCh1qDqgfQK6DRsJaPrAAkoGASWAgFIHASWAgFIHASVAu1afWSsCSiQBJYCAUoeCAjrwOaTZXysnoNsgoNRBQAkgoNShoIAOfauxgJZpKwFNHlpAyVBSQM/boLOuXtcjoNsgoNShqICeCjrz8ksdAroNAkodygroaS9+1ld4dAnoNggodSgsoN/2bwPfLPc8Ad2GVq0+tVIElEilBfS4E/+ZTVAB3QYBpQ6lBfSTm6ACug0CSh2KC+jnCOg2CCh1EFACbCag9wcXUDIIKAEElDoIKAEElDoIKAHSWn1unQgokQSUAAJKHQSUAAJKHQSUANsJ6O3hBZQMAkoAAaUOAkoAAaUOAkoAAaUOAkqApFafXCUCSiQBJYCAUgcBJYCAUgcBJcCGAnodgICSQUAJIKDUQUAJcK/VZ9eIgBJJQAkgoNRBQAkgoNRBQAkgoNRBQAkgoNRBQAlwq9WnV8jnA9qMQUDJIKAEEFDqIKAEEFDqIKAEEFDqIKAE2FRAL4MQUDIIKAGutfr8+hBQIgkoAQSUOggoAQSUOggoAQSUOggoAZpavWB1vCKg53EIKBkElAACSh0ElAACSh0ElAACSh0ElAAbC+hpIAJKBgElgIBSBwElwKVWr1gbAkokASWAgFIHASWAgFIHASXA1gJ6HIqAkkFACSCg1EFACSCg1EFACXCu1UtWhoASSUAJIKDUQUAJsLmAfhNQcggoAQSUOggoAQSUOggoAU61es26EFAiCSgBBJQ6CCgBthdQH+Ukh4ASQECpg4AS4FirF60KASWSgBJAQKmDgBJAQKmDgBJAQKmDgBLgb3971Zp4YUA9N5hNQAkgoNRBQAkgoNRBQAmwzYB6cjCXgBLgdZt7AkokASWAgFIHASWAgFIHASXA6066fOV5oJ4dzCWgBBBQ6iCgBBBQ6iCgBNhmQD09mEtACSCg1EFAWd8Lvz9DQIkkoKxPQKmEgLK+rQbU84OZBJT1CSiVEFBW9yagVEJAWZ2AUgsBZXXbDagnCPMIKKsTUGohoKztuBIElDoIKGsTUKohoKxtywH1DGEWAWVtAko1BJSVndaBgFIHAWVlAko9BJSVbTugniLMIaCsTECph4CyrvMqEFDqIKCsa+sB9RxhBgFlXQJKRQSUdQkoFRFQVnVZAwJKHQSUVW0/oJ4kPE9AWZWAUhMBZU3NChBQ6iCgrElAqYqAsqLr8t90QD1LeJqAsiIBpS4Cynpui3/bAfU04VkCynoElMoIKKu5L/2NB9TzhCcJKKsRUGojoKymnIB6ovAcAWUtycIXUOogoKylpIB6pvAUAWUtAkp1BJSVpMt++wH1VOEZAspKBJT6CCjraC16AaUOAso6Sguo5wpPEFBW0V7yAkodBJRVlBhQTxYeEVDW0FnwJQTUk4XHBJQ1lBhQzxYeElBW0F3upQTU04VpAsryeiUqI6CeLjwioCyvt9gLCajnCw8IKIvrL/VyAuoJwxQBZXHlBtQThmkCytIGFnoxAfWMYZKAsjQBpVoCysKGlnk5AfWUYYqAsqzB4zBFBdRzhlECyrIGF3lBAVVQJggoixpe4iUF1JOGcQLKkkY234oKqGcNowSUBY3t/pYWUE8bhgkoyxktT1kBVVDGCCiLGe9OYQFVUEYIKEuZqE5pAVVQhgkoC5lqTnEBVVAGCSjLmCxOeQFVUIYIKIuY7k2BAVVQBggoC3h7UJsSA/pwpvgOlRfQ3dvV1/l39gpYxcPUFBnQ82x5ApEqLKD7t7a5DfX8X94zmSkzoDZC6SoqoH/89Nb1w2+zpuDpv7S3pzbTSg2ohNJWUkD/+qVdzD9/Pv77x9/nTMKTf2FP9qXcgDZ/ITyPOCspoIdeLk9J/ZgzCU/8Bc0oS8EBvSXUU4myArrr77AfN0JnvQ3qWb+E9D2V5+5RdEBPZs8xdSoooMfNzffeD/fz9uE93V+q84b083csPqDfPjHzVKSggB63Nvu764epw0i9Q04sY/a6rCGgjehlT7AXPf/uT6hXT/BKQLcnd11WFNBG9Jogyouef/dn0qsneGUXvh71BZTvU0EBdRCpHgJKHUoK6PBpTP2t0gkCug0CSh1KCujgifTzPookoNsgoNShpIBeitn25ddZUxDQbRBQ6lBUQNNLMV24mEiZBJQ6FBbQby5nVwUBpQ7lBfRTBHQbBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6CCgBBJQ6fHcBZQv+9rdXTenvf3/RhF43JL4rL2/Uqyf4StELmwsBpRYvb9SrJxiioH39coZazkgNdQnljDRyqOUspCnW9QLKGamhLqGckQroZ1nXCyhnpIa6hHJGKqCfZV0voJyRGuoSyhmpgH6Wdb2AckZqqEsoZ6QC+lnW9QLKGamhLqGckQroZ1nXCyhnpIa6hHJGKqCfZV0voJyRGuoSyhmpgH6Wdb2AckZqqEsoZ6QC+lnW9QLKGamhLqGckQroZ1nXCyhnpIa6hHJGKqCfZV0voJyRGuoSyhmpgH6Wdb2AckZqqEsoZ6QC+lnW9QLKGamhLqGckQroZ1nXCyhnpIa6hHJGKqCfZV0voJyRGuoSyhmpgH6Wdb2AckZqqEsoZ6QC+lnW9QLKGamhLqGckQroZ1nXCyhnpIa6hHJGKqAABRJQgEwCCpBJQAEyCShAJgEFyCSgAJkEFCCTgAJkElCATAIKkElAATIJKEAmAQXIJKAAmQQUIJOAAmQSUIBMAgqQSUABMgkoQCYBBcgkoACZBBQgk4ACZBJQgEwCCpCppID++fPb0fvQr/74qfuriRuvYMZQ//rl7eZjtQHePFpOx+F+PH3jZc0YavBCHR/q0MBCl+qMkW52oR4duqNaa6EWFNBds+6+/Nr9zX3Nfjy+8RrmDPWypqOelo+W02l0m1+oF+lQYxfqxFAHBha6VOeMdLMLtdkoOfrht8c3fq1yArq/rbzbUmoM/GUcv/Hmhnp46/1oCyNt7JJRbXWhNtKhhi7UqaH2Bxa6VGeNdLML9dbPt7cff3904xcrJqCnhXRaOqe1+LX9q9PSOm2rn/5CXhbgxI23NtTTjyL2h84eLqdD8mLZ7EK9SIcaulAnh9obWOhSnTXSzS7U24vp+upac6EWE9DddVEcF057u/y4Vdf84LgkL/81fuM1zBrq8dYR23PfmseeXk6XnbaP5268rFlDDV2ok0PtDSx0qc4a6WYX6v62MXJo/mvFhVpKQI+vj+va23feGz4upOtfmeZv5MSNVzBrqKemNmt/fY+W03FsP/zP5B2wbS7Us9ZQQxfq5FB7AwtdqrNGutmFet8quf7nmgu1lIAm6Tn+59iK3F1eQs/deCmzhnpa2evvDzcejfQU+dt+27YXamuooQt1cqi9gYUu1Vkj3exCTX71+MYvV0pAk78kp42N4V2J61+ep268mFlDPe11fJx3PwMOwj4Y6fHJ935/42vTC7U91NCFOjnU3sBCl+qskW52oR56G5lrLtSCApqc4ze8CvfXN8GeufFyZg31+F9f/qF9BHE90yO97LKlAd3sQu0MNXShTg61N7DQpTprpJtdqOdfXU5v+Xh445crJaC7+zI5/vfAMtklZ1c8vPGiZg31dsJaxPNyeqSXn+xLWKidoYYu1Mmh9gYWulRnjXSzC/X0z+tpS5etzTUXajkBvS+HJwL64MaLmjXU01/Oy01Op1ysfBBhcqTNntG+gIXaHWroQp0aan9goUt11kg3u1CP//zXW9nPBV1zoZYZ0P75aNcz1Js/69M3Xtasod5PZzq9j7fye2BTIz0O7Hpu3VBAt7RQe0MNXahTQ+0PLHSpzhrpVhfq5eV0fgJcT/tcc6GWGdCRPyqnt7j7C3BTG0tX16G27xj6Auruwl1eIwVsgfaG2v5lbOtHF9RlYFvdAm3frLsIt7RQzwFtNoib0z5tgfY9t0yav5Lbfa3fJX/Qr1Y/EXBipP1jR5tdqEOHuVq/3dIbI3eXgZUQ0P4i3NRC3SVvya7/V6mUgD55YO2yZrd8wDi9WedJ2D8fY2HjI72fSbf5o/ADQ02tvlCfXVCH+KfqrJFO/2RpE0PdJYNZf6GWEtDDc6d2XW725I0XMmuo0z9Z2PhI7xdjuNjwQh0Y6sgdV1LOU7WcZ+rUUPe9gK65UEsJ6JMfLrgsum1/aKYxeP7vuu8sjY90oEpbXagPArr+FTCeXFCXgW32k0iJLVxWZGKoh+Rowvqv/1ICOvHx1vTfl/dANvux7f5Q2x+aWHkXbnykA1Xa6kIdGGroQn1y/Qd8bLtn1kg3u1DTcwLWf/2XEtDTohm5wEqyAK//OX7jNcwZavKTw/rXM3tiObXOTt/kQr1JP7YfuFCfW//XgYUu1Tkj3exCvf8q4vVfTEDHL/F3Oo/hdh7L6tcD/ORQWz9Z/QX0xHK6B3SrC/WmdSJ93EJ9dv1fBrbR64H2R7rZhXr71f0Mf9cDHXDfVWv+EN5OBEouSJ0ckG3feLNDTX+y/tfijI80ucVGr0g/MdTYhfrk+n8fufFWR7rZhZq+i+OK9ON2zSK5vmruC/D+ZS1fx2683aHen5cRl/seH2kjPV6w1YXaSIYau1AnhjowsNClOmekm12o9y8bSU8IXWehFhTQ7hftpS+gy7pNj7dt6gskp4d6aLd/XRMjPWsdcN3sQj1rDTV0oU4NtT+wLX0r5/RIN7tQL8FMtzZ9KyfAxgkoQCYBBcgkoACZBBQgk4ACZBJQgEwCCpBJQAEyCShAJgEFyCSgAJkEFCCTgAJkElCATAIKkElAATIJKEAmAQXIJKAAmQQUIJOAAmQSUIBMAgqQSUABMgkoQCYBBcgkoACZBBQgk4ACZBJQgEwCCpBJQAEyCShAJgEFyCSgAJkEFCCTgAJkElCATAIKkElAATIJKEAmAQXIJKAAmQQUIJOAAmQSUIBMAgqQSUABMgkoQCYBBcgkoACZBBQgk4ACZBJQSvDn/zr971+/vH359cnbwgoElALs3t5P//dUQJvbwgoElO07vD0f0OttYQUCyvYJKBsloGyfgLJRAsr2CSgbJaBs3TGJZx9NQM///vH3zu8/Orc92Z3/+4kD95BHQNm6TkB3rUZ+++On5t8//PatHdD921VSW3glAWXr2gH9h7c0mPd+Xn6QBPTeTzv1LEVA2b7kPdBk6/K9+cl5F/1Uzq/pbf/8OfkPm6AsQ0DZvjSglw3PU0G/Xn7TxPHYyWtJ35v7fL3c+7iR2twJXkxA2b40oM3e+LGKp3Lu7oeImhv1j8If0yqgLENA2b6B05guAT228bZ3fvzJ12/9gJ53+wWUZQgo2zcV0NQppklAD7dfCCjLEFC2bzSgyTH4XkDT3wkoyxBQtm8qoN0D7OlR+Es7vQfKYgSU7Zvahe+28Xrb/W27U0BZjICyfaMBPR0h+hi8bXLTg114liKgbN9oQJPtzM5pTPebnvblBZRlCCjbdz1dvh/Q8zud523Q3f0jSefb7pps7l1PhOUIKNt3OaD+MRDQ5FSlZiv1etv0FwLKQgSU7bt8Bv59KKDJ2Urvrds217I7+uf+G6XwGgJKAc5V/DoY0Otll752bntN63v6sXh4LQEFyCSgAJkEFCCTgAJkElCATAIKkElAATIJKEAmAQXIJKAAmQQUIJOAAmQSUIBMAgqQSUABMgkoQCYBBcgkoACZBBQgk4ACZBJQgEwCCpBJQAEyCShAJgEFyCSgAJkEFCCTgAJkElCATAIKkElAATIJKEAmAQXIJKAAmQQUIJOAAmQSUIBMAgqQSUABMgkoQCYBBcgkoACZBBQgk4ACZBJQgEwCCpBJQAEyCShAJgEFyPT/AVK4gnWfd1vyAAAAAElFTkSuQmCC)
print(x.theta.mle)
## [1] 0.4246842
print(c(posterior.llimit, posterior.ulimit))
## [1] 0.403 0.446