{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Stochaskell, version 0.1.0\n",
"Copyright (C) 2015-2019 David A Roberts\n",
"This program comes with ABSOLUTELY NO WARRANTY.\n",
"This is free software, and you are welcome to redistribute it\n",
"under certain conditions; see the LICENSE for details.\n",
"\n",
"Using installation directory at \n",
" /home/jovyan/stochaskell"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"{-# LANGUAGE MonadComprehensions, RebindableSyntax, NoImplicitPrelude, FlexibleContexts, TypeFamilies #-}\n",
"import Language.Stochaskell\n",
"stochaskell"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"covPrior :: Z -> P RMat\n",
"covPrior n = do\n",
" tau <- joint vector [ truncated 0 infinity (cauchy 0 2.5) | i <- 1...n ]\n",
" corr <- corrLKJ 2 (1,n)\n",
" return (qfDiag corr tau)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"birats :: Z -> Z -> P (RVec,RVec,RMat,RMat,R,RMat)\n",
"birats n t = do\n",
" x <- uniforms (vector [ 0 | j <- 1...t ]) (vector [ 50 | j <- 1...t ])\n",
" betaMu <- normals (vector [ 0 | _ <- 1...2 ]) (vector [ 100 | _ <- 1...2 ])\n",
" betaSigma <- covPrior 2\n",
" beta <- joint (designMatrix 2) [ normal betaMu betaSigma :: P RVec | i <- 1...n ]\n",
" let beta' = tr' beta\n",
" yMu = asColumn (beta'!1) + outer (beta'!2) x\n",
" -- = matrix [ (beta'!1!i) + (beta'!2!i) * (x!j) | i <- 1...n, j <- 1...t ]\n",
" ySigma <- truncated 0 infinity (cauchy 0 2.5)\n",
" y <- normals yMu (matrix [ ySigma | i <- 1...n, j <- 1...t ])\n",
" return (x,betaMu,betaSigma,beta,ySigma,y)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
":opt svg\n",
"import Language.Stochaskell.Plot"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"rejecting OOB sample -3.0211867938117805\n",
"rejecting OOB sample -6.9251262844039\n",
"rejecting OOB sample -13.073243189832711\n",
"rejecting OOB sample -138.31103267470144\n",
"rejecting OOB sample -4.257210070571588"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAIAAADfNCTgAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeVxU1fsH8A8goLiUgICCAirihuCWe4niLpoi5gYK5iBl+Ks0LCu07JumqJgbuIO5zLixuILkglAKKCC7oCiKDMmiKNvMnN8fY0rIADPMcmfuef/R6+Xcy73P4+DTueeeRYsQAoqiKEpW2qoOgKIoSr3RMkpRFNUstIxSFEU1Cy2jFEVRzULLKEVRVLPQMkpRFNUstIxSFEU1Cy2jFEVRzULLKEVRVLPQMkpRFNUstIxSFEU1S2NlVJixf9GoXqYG+v1/ShECgDDNf5yVhYWFhYWFuZGB/tCN2cLq6M8s2hpbWFhYWHTp9/Wf1QDIP9FrJtrZ9rK1G786kk9n7VMUpbkaK6NaRoMXbw497TugxesPdHp/HfkgPz8/Pz+X52HpMG1KVx0ABtP35OTn5+c/TPZ31AMqrqz9LPLDI0npKScmxX3xY+RLBWdBURSlMo2VUe0OdqMGdzPU13rnSNVfx8MNZ82y1Xn3h2qSzkebuMzp0xJ6PebMtow+l1gjn2gpiqIYR+a+0cobx8+ZzXLpLq6iFed9+lpb9RzpHhBXQiAqyOebWphpA9AyNDd7kf+4Wl7xUhRFMUyLxk+p16urxy50cY220gagO3B1zL0Ay/bCvLCvnOd/2/tuAGovYkoI8G5j9g0trf8cXLFixcaNG9/8kc/n8/n8N380MTExMTFR06N8Pt/ExIRpUSni6Jv/MioqBR1t3bq1tbU106KiR5tz9OXLl7W/08aRJhCk/zps0NpkwdtPXkR4WI3fnS+sc2LNrW/tRvyWVRG3svco/3tCQojo8a5xNj5XqyVfvIkxaICUlBRVh6Ak7MmUsClZZmQqImQnIb0Veg9pM5Xtof551PE/+8ye3kkbAEjJ/Wx+NQBh0fXjF0u79zRraT9xdOHJ46mVqM48xs0bPWmgrky3oSiKqqUQmA4cAE6pOpL/aKyMih7vn9XVcvT6xGT/cZbWLvvyRQApvXT8mv3saaZar085u2J0t07m5hYDl2c57/11Slu0clyzY8yVOX1ses88N2TbT+NbKz4RiqI020nAAegN3ABsVR3Mf2gRVe/FJO4b9fPzW7NmjWojUbS7d+/27dtX1VEoA3syBZuSVV2mz4GVQBQQDIxQwv2kzVTWV0xypfJSrhy1u7E1G3syBZuSVVGmscBCYAyQDCjpuVbaTBnRGlV5DBRFMU8V4AeEAIHAVFUH0xBGtEYpiqL+KxVYAHQDkgEjVQfTCEYsTaKlpaXxHaMURTUNAQIAR8AbOMH8GgqGtEbpQz1FUQCAPGAhUAPEAd1UHUxTMaI1yhK1J05oNvZkCjYlq/hMecBQYBJwXbU1VNpMaRlVHvrvTSOxJ1lFZloEzAR+As4DviqvS7SMUhSlXi4CA4DOQDzg0PCpIiLa+tfWPjv7KCeyJmJEGaWvmCiKlSqA5cBS4DAQAOg3fHZOSc6YQ2NOpJ049QmzJoPSV0wURanE34A7YA/cBt5v+FQREe1N3Ls6evWK4StWDl+prcWI9t8bjCijLEGnu2gk9iQrv0wFgD+wDdgBfNzo2WlFaZ6hnno6erGLY20MbeQUQ0PoLCaKopgsHXAHjID9QKeGTxWIBP5x/r/d+O2bEd8wsBH6BiPCon2jFMUCBAgCPgQWAOcbraF3+XeH7xsefT/6ttdt3xG+jK2hYMhDPW2NUpSmewp8CvwDxAKNPJiLG6H+sf7rxqxbMmBJnQ0yGIgRZZSiKI3GA74AFgE/A42s4Z7CT/E449GhdYcEr4TO7TorJbzmYm47WfPQcdoaiT3JypRpGeAF+AFngfUN19AaUc2GGxucgp04Aznn559XYQ1Vy+H3LOkbpf/eNBJ7kpU+08tAPwDALWBgw6cmFSYN2TPkWt61BE4CZyBHtgjlRdpMGfFQT/tGKUqzVAJrgBBgDzC54VNrRDWb4zZvjtv8s+PPKi+gsmFEGaUoSoOkAAsAWyAFMGz41L8f/+0Z6tm1fdfbXrc7tW3k3T1jMeKhnqIojSAENgBjgRUAt+EaWimoXBW1auqRqcuHLA+fG66+NRQMKaMs6Rul0100EnuSbSzTB8AYIApIBNwavlRcflz/wP65Jblpn6cx8EGezmKiKEr5goEVwNfAyoYbZxWCirVX1gYnBW+btG1W71lKi0+haN8oRVHNwQc4QC4Q9e97eYluPLqxOHRxP9N+yd7JxgbGyolPCWgZpShKZucBDjAf4AJ6DZz3qubVT1d/CkkO2T55+4yeM5QWn3LQvlGKomTwAvACPgeOAOsbrqExD2McdjvkluQmeydrXg0FQ1qjLOkb5fP5LHkdwZ5MwaZka2X6F+AGDAOSgTYN/Ii4EXo4+fCOKTum205XTpzNJ+13yojWKEvQ6S4aiT3J8vl8oAZYA8wE/IHghmvopZxLvXb0EjdC1aiGQk1nMVEUxXz6+jnAQsAUSAA6NnDm86rnKyNXRmRF7J6627mHs9IiVBXaGqUoqlEECLK2Xgx4AecarqEX7l2w22UHIGNZBhtqKBhSRukrJopisIfAWOBgbm4w0NBQ+bKqMq8Ir8/OfrZ/+v7AqYFt9doqLUTVYsRDPUteMbHkRQTYlCk0P1kesAzwBn54//1nDZx3Lvvc0oilk2wmJXsnt9FrqM+U+egsJoqi5KIUWAbcBg4D/Rs6r7LUN8o3Kjdq77S9jlaOSouPORjxUE9RFMNEAv2A1sCthmtoRFaEuCc0aWkSO2soGi+jwoz9i0b1MjXQ7/9TihAAUB39mUVbYwsLCwuLLv2+/rMaAPknes1EO9tetnbjV0fySb2fNIT2jVIUY1QAqwAOcAgIBAwknVdSWeIV4fV/F/7v8MzDgVMD1f1BvjkaK6NaRoMXbw497TugdieqwfQ9Ofn5+fkPk/0d9YCKK2s/i/zwSFJ6yolJcV/8GPmynk8aRAihZZSiGEDc9swFEoGGmpa8NF6v7b0AJHknfWT5kbLCY6jGyqh2B7tRg7sZ6jewNV9N0vloE5c5fVpCr8ec2ZbR5xJfvfNJjVyDVlMsG6fNFpqSrADYAEwFfgC4QPt3zxBnyn/Jd+W5rrmyJnRuaODUwNa6rZUeqsIpZS+mivM+fa2teo50D4grIRAV5PNNLcy0AWgZmpu9yH9c+c4n1bLcRtNoyr+3xrEnU2hIsvcBR+AykAjMl3QSn8/npfHsd9t3a98t0StxiPkQZYaoTIqfxaQ7cHXMvQDL9sK8sK+c53/b+24Aar9oJwTQqueThtTeh3rFihUbN25880c+n187JRMTk9pjEehRZh5FLcyJSkFH62BIVNIcDRaJvios/PTZs/lACVBS78/+U/nPutvr8srzDk08NL7PeFXHrNijkBZpAkH6r8MGrU0W1Pm45ta3diN+y6qIW9l7lP89ISFE9HjXOBufqy/f+aRa8sWbGIMGSElJUXUISsKeTIl6J/uUEGdCBhOS0fB53FSu6UZTzyOelYJK5USmWtJ+p1I/1JOS+9n8agDCouvHL5Z272nW0n7i6MKTx1MrUZ15jJs3etJAg3c+aWh3aoqiVOAk4AD0Bm4AtpJOKigv+PjYxz9f/fns/LNf2n2pr6OvzBDVRiNlVpi/z8Xa3PQ9fb12puZWM/c+qk75fVovi46dOpl1tpvmd/GJkBAi4kd+P6531+7deo1ddaGg/k8kajwGTVFYWKjqEJSEPZkStUy2jBAOIV0JiWn4PG4q12SjiW+kb5WgiqhlpjKSNlPVzyASd4z6+fnRMU8UpXixwEJgDLAZkPiS/cmLJ0sjlj4se3jg4wP9zRoafk+BzqmnKNaoAvyAECAQmCrpJELInsQ9313+7tMBn5785KSuNu2QaxwjyihFUQqWCiwAugHJgJGkkx6UPlgSvuRF1YtrHtd6d+itzPjUGp1TT1GajQABgCPgDZyQVEMJIUEJQR/s+cCpq9ONxTdoDZUKI8ooS+bUyzIeTT2xJ1MwPdk8wBHgAnENLBV6v/T+2OCxB+8cvOZxzXeEr46WTr2nMTtTeZI2U0aUUcKOOfX0t1AjMThZHjAUmARcB7rVe4a4ETpkz5AJ3Sdc97ze07hnA5djcKZyJm2mtG+UojRPEeAFZAPnAQdJJ+WU5CwOXSwQCa57Xrc1kjh0lGoUI1qjFEXJz0VgANAZiJdUQ0VEFJQQNHTv0Ek2k655XKM1tJloa5SiNIZ4qdAw4DAgcfG6tKI0z1BPPR292MWxNoY2yoxPUzGiNcqSV0x11uzQYOzJFAxK9m/AASgAbkuqoQKRYMONDaMOjJrRa8aVRVekraGMyVThpM2UEbOYVB4DRakzAeAPbAN2AB9LOuku/65nqGf7Vu33OO/p8l4XZcan8ehDPUWptXTAHTACbgGd6j1DIBL4x/n7x/qvG7NuyYAltdelpOSCllGKUlME2AOsBr4HfCQt6pvCT/E449GhdYcEr4TO7TorOUSWoH2jFKWOngLOwH4gFlhebw2tEdVsuLHBKdiJM5Bzfv55WkMVhxGtUZb0jfL5fJZ00rMnU6gmWR7wBbAI+Bmof+mQpMIkjzMeHdt2TOAkWLSzkMtd2fO1SpspI1qjLEEngWgk5SZbBngBfsBZYH29NVTcCB0fMn7poKVn552VVw0Fm75WOouJojTVZcATmAjckrRU6N+P//YM9ezavuttr9ud2tb/xomSO0a0RmnfKEU1qBJYBbgDu4DAemtopaByVdSqqUemLh+yPHxuOK2hysSI1ihL+kYpSiYpwALAFkgBDOs9Iy4/zjPU087ELu3ztA4GHZQcH8WIMsoSLOmeB5syhWKTFQKbAH/AH3Cr94wKQcXaK2uDk4K3Tdo2q/cshUUCsOlrpbOYKEozPAAWAnrAAaD+10Q3Ht1YHLq4n2m/nVN2GhsYKzc86i3aN0pRDBQMfABMBi7WW0Nf1bxaFbVqNm/2r06/cl25tIaqFiMe6mlrlKL+xQc4QC4QBfSr94yYhzGeoZ4OZg7J3slGrSRurEQpDSPKKEVRAIDzAAeYD3ABvXcPv6p59dPVnw4nH94xZcd02+nKj4+qFyMe6lmCjl7WSHJK9gXgBXwOHAHW11tDL+Vc6rWjV25JbrJ3skpqKHu+VrXci4klfaP0t1AjySPZv4ABQAWQDIx69/DzqudeEV4eoR7bJ2/nunINW9U/7EnR2PO1quUsJto3SrFVDfALEATsBqbVe8aFexe8Irwmdp+YsSyjrV5bJcdHNQUjyihFsVIa4AaYAglAx3cPl1WVfRP5TWRO5P7p+8daj1V+fFQTMeKhnqJYhgBBwBjACzhXbw09l33ObqcdgGTvZFpDGY62RpWHTgLRSNIn+xBYBFQCMUD3dw+XVpb6RvlG5UYdmnHI0cpRHjHKB3u+VmkzZURrlCWvmOhvoUaSMlkeMBj4ELhebw2NyIqw22UHIGlpEqNqKNj0tUqbKSNao/QVE8UCpcAy4DZwAej/7uGSypJVUasu514+PPPwR5YSt0emGIgRrVGK0nSRQD+gNXCr3hrKS+P12t4LQJJ3Eq2haocRrVGK0lwVwFrgOHAIqOchnf+S//m5z9OK0kLnhg4xH6L8+Kjma6w1KszYv2hUL1MD/f4/pQgBQJR36qsJdlYW5uaWDrM2xpYQoDr6M4u2xhYWFhYWXfp9/Wc1APJP9JqJdra9bO3Gr47kN/bIzpK+UTp6WSM1mKy47ZkLJNZbQ3lpPPvd9t3ad0v0SmR+DWXP1yp1pqRhQn7ytZuZN9YMHbQ2WUAIIcKHMRGxj16JSM2T0x42fb+LryFVl71t5p+qrPVTr6KX2Q7/JbmCVGVsdrT1uljewB0aj0FTpKSkqDoEJWFPpkRisjWErCfEhJDD9f5UwYuCGcdm9N3Z9+bjmwoNT47Y87VKm2ljrVHtDnajBncz1H+zf6t25xFThlm00kILs1Gjejx78lT07g/VJJ2PNnGZ06cl9HrMmW0ZfS6xRrraTlFq7T7gCFwGEoH57x7mpfEcdjv0MOoRz4kf3Gmw8uOj5Ev2V0w1GQeDs8d/PEwXACrO+/S1tuo50j0groRAVJDPN7Uw0wagZWhu9iL/cbW84qUoRiPALmAIMA+4CJjXOVxQXvDxsY9/vvrz2fln1zut19fRV0mUlHzJ+IqJFEevWnjYdvOFaYZaIANXx9wLsGwvzAv7ynn+t73vBqD2CCZCAC3JlwKgpfX2+IoVKzZu3Pjmj3w+v3Y/hYmJSe0hXfQoM4+iFuZEpaCjtRQCS2pqHj14sL+qygpIrfOzvDTeZxGfTes8zW+En+4z3bvP7jIzowaO3r2rfjHLcBTSasqTvyD912H/9o0SQkh5wm9j+87Ym1lV98SaW9/ajfgtqyJuZe9R/veEhBDR413jbHyuVku+eBNj0ACFhYWqDkFJ2JMpeZvsSULMCPElpJ5f9sfPHzsfcbbfZZ9YkKjk8OSIEV/r82xycykJ76XQm0ibqfQP9dXZ+90Xnh97IHhxD/GaiKTkfja/GoCw6Prxi6Xde5q1tJ84uvDk8dRKVGce4+aNnjRQV+rbaKA6zTQNxp5MAZiYtAS8gJXACWA98J9fdkJIUEJQv139enfofYtzq79ZPYNG1YWKv9biRMS54+IQaOthbJRCbyV1po2UWWH+Phdrc9P39PXamZpbzdz7qCJqaccWBkbmYl09T70QpPw+rZdFx06dzDrbTfO7+ERICBHxI78f17tr9269xq66UCBs6A6Nx0BRzBVLSDdC3Ah58e6x+yX3nYKdhuwZkspPVX5kmkJECiLJlakk1JpkbCU1DQ37URXV78op7hj18/Njw9BRSoO8WSo0EHCuc4wQsidxz/fR3389/OsVw1foaOmoJET1JqpG3jGk/QZtPfT8P1jNgxZDpwsxIiyVl3KKklIq4AZ0BZKBurty3i+9vzh0caWg8prHtZ7GPVUSn3qreY7cA0jfhNaWcPgV5nX/L8U0dE698sjyBlA9aXSmBAgAHIGlwAnAuHayIiLa9ve2D/Z8MM12WoxnjIbVUGV8rZVPkbIGYd1RnIDR5zEuRiU1VNpMGdEaZQk+n8+Sdy+am6l4qdAqIA7oJv7oTbIPSh8sDltcUVMR4xlja2SryjAVQ7Ffa2ky0jfh8VlYL8CkBBh0VtSNmkDaTBnRGmXJnHpKzb1ZKvTamxoqRggJSgj6YM8H47uNv+55XSNrqAIVxeCqM6InoE1XTLuHgQGqraEyYERrlPaNUswmXir0Tr1LhRZVFk07No3/kn/V42ov414qiU8tERGenMXddRBWoMcyjORBp6WqY5IRI8ooRTFYJOAJzAQSgLpzN3lpvKVRS5cMWnLqk1O62nR4dNMIypGzDxlbYGCBvt/DfGpj8xyZjpZR5dHQ7sJ6aEqmFcAqIBQIAUbXOfa0/KlXhNeD0gfcadyxvVmx5ZwcvtZKPrJ3ImsHjIdixDEYD5VHXPInbaaMGDeq8hgo6h23ADegHxAItK9zjJfG++LcF4scFv3k+JOejp5K4lMz5TnI3Ib7ITCfir7fo20PVQckT4xojWppadHh9xRjCAB/YDOwBZhX5xj/JX9pxNLs4uyz888O7DhQJfGpmeIEZAagIBI2XnDOhr6RqgOSP0aUUdoapRjjPuAOGACJ7y5zx0vj+Zz3WWi/8Oiso3SNu0aI3yClbUBFAWx9MHg3WhioOiZFYUQZpShmCAZWAKsBnzovPUoqS3zO+yQWJIbNDaMLLTdCVIW840hdjxYGsPWB1Xxo+lxYRowbZQmNntvzH2qYaSEwDdgOXAeW16mhZ7PP2u20M9A1uLXk1rs1VA2TlVHjmdaUITMAYd3xkIdB2zAxHtbu6lhDpf1OGVFGWTL8nv57Y6pTgAPQG7gB/GfkfGllqVeEl895nyMuRwKnBhro1vNYqm7Jyq6hTF8+QMJyhFqjOAGOF/FROMyclBianKnlZFDaN0qpyHNgJRAFnABG1Dl2/t55r3CvSTaTkr2TW+u2Vkl8aqDkDjI2v57EOTkFBnU7lNmAEWWUolQhDnADhgNJQJvaB8qqyr6J/CYyJ/LQjEOOVvVsjEwBQFEM0jag9C5slmLQ79B9T9UBqQwtoxQLNbRU6KWcS0vCl0zsPjHZO7mNXpt6f57VRDXIP4O030Bq0PMrjDoF1k/fYkQZZcm4UU2Z29M4ZmcqcanQ51XPV0auvJRzaf/0/WOtmzoxidnJypOpoQEyA5Duj9ZdYPejBkzilITOYqIoSQiwDfgFWAdw6hyLeRjjEeoxzGLYjik72uq1VUl8zFVZiOxdrydx9v0RRnTI138wojVKUYpXz1KhYhWCirVX1oYkh+yeutu5B9MXWle2F9nI2o4HR2A1DxMT0LqLqgNiIlpGKTbgAcsAb+AH4D/DGGMfxXqEetib2qd4pxi2MlRVfEwkfoP0LB42XnDOhB79y5GIjhulNFspsABYC1wA1tSuoRWCilVRq1x5rhucNnBdubSGvkZEeByOS8PwlwfMnDA9F3ZrGFJDRSJcuoS5c9GTYZuzMKI1ypK+Uc3dWqMuxmQqcanQv/L/8gj1sDOxS/ZONmrVrMUyGJNss72exPk/6Buj51foPLPOBCQVZvroEY4cQWAgWrWCuzt+/12xt5M2U0aUUZbQnH9vjWFAphKXCq0UVK65siY4Kfj3yb+79HJp/p0YkGyzVRUhaweydsJ4CIbuh/Hwes9SfqaVlQgPR1AQEhMxaxZOnMCAAcq4Ly2jFPVmqdDbdZYKTSpMWnh6YQ+jHkneSR0MOqgqPgYpz0VmwOtlQMddQzumPC0nJCA4GEeOwMEBbm4IC0OrVqqOSTLaN0ppEgGwAZgK/Ahwa9fQGlHNhhsbxoeM/2zwZ1xXLq2hKE5EnDsuDAaAKXcxLJgJNfTpUwQEwN4e8+ahfXvExyMyEu7ujK6hYEhrlCV9o5SCSVwqNIWfsujMIrM2Zre9bndq20lV8TEDwdPLyAxAWSpsl2PwLrRQ/XIBQiH+/BNBQbh4ERMmwN8fY8dCS32G9jOijLKE2vegNZkqMq1/qVCBSOAf5+8f679uzDrOwLpD7uVCbb5WUTXyjiHtN2jroef/weo0tKT756+ITNPTcegQDh6EhQU4HOzfjzYMmH9LZzFRbFMILAGeAiF1lrlLLUpddGaRUSujvdP2WrSzUFV8qlfzHLkHkL4JrS3R2xfmqp9iUFaG48cRHIz79+HmhsWLYWOj6piagbZGKbV2CvgcWAicBN4ukFG7EbpkwBItNXo+lK/Kp8jejayd6DQRo8/j/b6qDUckQmwsQkLA42HECCxfjhkz0KLpRYgQXLmCffuQkID0dAUGKiVGlFGWLE1CyZXEpUJzSnI8zni0bNEynhPf5T22Tl4sTUb6ptfLgE5KgEFn1YaTn48//kBQEFq2hLs7MjPRQaqXfAUFCA7G3r3Q04O7O7ZuVVSgMmFEGaUP9ZSU6l8qlBCyJ3HP6ujVK4avWDl8pbYWIwaiKFvtSZzT7kGv7u7QyiQe+BkcjLg4uLiAy8VAqXZTFYkQHY2gIERHw8UFBw5g5EhFxdoMjCijLKEJ47SbRpGZSlwq9H7pfY8zHgKRIHZxrI2h8nramPK1infivLsOwgr0WIaRPOi0lO8dpMpUPPDz6FHY28PNDcePw0CqjUHF85Z27kSHDsp/90SH3zMXU/69KZ7CMq1/qVBxI/T76O+/Hv618huhqv9aBeXI2YeMLTCwQN/vFbcMaFMyLS7GiRPYtQvFxZg7FzdvwspKmntUVSEs7O28pYgI2Nk1I2QZSfudNvYLJ8zYv2hUL1MD/f4/pQjFH5F/otdMtLPtZWs3fnUknzTxk4bQ4fdUYwgQADgCS4ETtWvog9IHTiFOB+8cvO553XeEL7se5Cv5SFmDUGs8jcKIYxgXA3NnlSylLBQiKgqzZ8PaGlFR2LgRDx5g/XppamhaGlatQpcuCAoCh4OCAgQGqqSGyoI0TMhPvnYz88aaoYPWJgsIIYS8il5mO/yX5ApSlbHZ0dbrYnmTPpGs8Rg0RUpKiqpDUBJ5Z5pHiCMhwwm5V/tTkUgUGB/Y4bcO62PWC0QCud5RCqr5Wl/cI/E+hNeexLqR55nKuWe9maanE19fYmpKBg4kgYHk+XMpL1pWRgIDyYgRxNyc+PqS+/flEWlzSfudNvZQr93BblQHYcZlLVSLP6hJOh9t4hLepyW0e8yZbfnRucRX7Rr/pGb8KLZv10LJqP6lQgvKCzjhHP5L/lWPq72Me6kwPmUrTkBmAAoiYeMF52zoN2t5Kpk9f44zZxASgrQ0zJqFyEjpG44JCQgKwokTGDsWvr6YPBk66rejvZjUfaOigny+6SAzbQBahuZmL6IeVzbhk+rag/rYiiUdo5BbpqXAMuAOcAHoX/sAL4237NwyDwePU5+c0lX1fmpK+lrFb5DSNqCiALY+GLwbLaR6ZSMHJiYmhODGjbcDPzkcKQd+AigpAY+H7dtRU4NFi5CRIeXQJ2WQ9juV4RVT7dFJhABaTfqkIbVHR69YsWLjxo1v/sjn8/l8/ps/mpiY1M5QvY6K/8u0qJh6NBLwrKiYlJsbTIgucFd8VGQg8orwelD64MKCC+ba5plpmUyIGbXI/76vlwFdL4De0/dml3WcSmq0kZGr5Hz5fN2IiPdPnKhu1Urb07NFRgbEpzT1yiJRm5s3O4aH61+5ggkTsHkzxo7lFxXxCwtRWKigmGU+Cmk15clfkP7rsH/7RqvjVvYe5X9PSAgRPd41zsbn6ssmfFIt+eJNjIFijVeE+BBiScifdQ5wU7mmG019I32rBFUqiEv5qktJxlZy2oJcmUoKIlUSQmUl4XLJ1KnEyIhwOOT6dekv8egRWb+eWFnJ2nuqBqR+ralrP3F04cnjqZWozjzGzRs9aaBBEz6hT/RU09wC+gMFwO3ayy3zX/JnHp/509Wfzs4/u95pvZ6OnuoiVIqXD88ULr8AACAASURBVJCwHKHWKE6A40V8FA4zJyWHkJr69s25qysePkRgoDSD36urwePB2RkODsjNxenTiI8Hh4O2mrjraiNlVpi/z8Xa3PQ9fb12puZWM/c+EhIRP/L7cb27du/Wa+yqCwVCQpr0iUSNx0CxQg0h6wkxIeSPOge4qVyzTWa+kb6VgkqVRKZUxbdJrBvhGZJ4H/IyXwX3LyaBgaR/f9K5M/H1Jbm50l/izcv7ESNIYCB59Ur+UTKM6kuYuJr7+fmpOhCFKywsVHUISiJ9prmEjCRkPCH/KRzFFcULTi3ovaP3zcc35RiefMnta+VfJ1emkjNWJHU9qS6VzzWbTCgkkZHE1ZW0a0dcXUlkJBGJ6p7TSKavXhEulzg5kU6diK8vyclRXLSKJu13yohZTIQdc+r5Kp/uoixSZlr/UqFns896hXtN6THl1pJbBrrKfivddM39WkU1yD+DtN9AatDzK4w6BeWOPcjMxNGjOHAAHTrAzQ27d8NQwjagEjMVD106fhyDB8v08p5xpP1O1TtbSs29WSr0eu2lQksrS32jfKNyo464HPnQ8kMVxqdYNS+Qux/p/mjdBXY/Km4SZ73eDPxMTYWrK8LD0a+flJcoLQWXi507UVEBT09kZsLUVCGxMh4to5Sq1L9U6Pl7573CvSbZTEr2Tm6tq/r9LRSishDZu5C1A8ZDMeokjAYr8+Z12o4ffwxdqZq/4lWXgoMREQEnJ2zapGY7figAI8ooXW+UZepfKrSsquybyG8icyIPzTjkaOWowvgU6EU2srbjwRFYzcPEBLRW3nKoT54gJAR790JXFwsXytJ2bMHnY8OGt4uGbt0q8fmfZRhRRlnSN8qSjlE0kmn9S4Veyrm0JHzJxO4Tk72T2+gxYDueJmvq11p7GVDnTOgpqQBVVeHSJYSENGPFzn83nOsRFQVXV/B4StotXnWk/aeq+n2Q6F5MrFH/UqHPq56vjFx5KefS3ml7x1qPVWF8CiGexJn6P1T9gx7L0J0DHSVtFpyaipAQHDgAGxu4u2P+fLSWto8kKwv797/dcG7BAikXDWULRrRGKRaof6nQmIcxHqEewyyGJXsnt9XTrIHZrydx/g/6xuj5FTrPhJYylt4Qv/gJDASfj/nzEReHrl2lvIR4zfqgIKSmwt0dMTHo3l0hsWoK1bcEaWtU0xFgG/ALsA54u8VxhaBi7ZW1Ickhu6fudu6h+r0q5amqCFk7kLUTxkPQ51sYD1fCPeu8+HFzk2nJpNpr1svy+omlGNEapa+YNNdDYBFQBcQB3d58Gvso1iPUw97UPsU7xbCVBr2mKL+P9E3IOwbL2Rgfi7bKaMRlZeHIERw8CGNjuLlhyxYYSbt4nni/48BAFBVh3jzEx6MLW7cClI3cJwBIiwkxKAf7ZjFxCTEhxI+Qt2sqv6p55Rvp28m/0+n00yqKTs5eJ1uS/HYS56vHSrjvm0lD4kVD7tyR6Srx8YTDIe3bS5y6VAv7foGbihGtUZZgzyymZ89yTEy+enep0L/y//II9bAzsUv2TjZqpZr1huWu/P4Fk3Tevztx5kDvfUXfsbkDP/HOfsVZWTA2bvSH2PMLTGcxUSoX2b27G/AJkADoiz+qFFSuubImOCn498m/u/RyUW18ckLwOAKp6y2eP4TdCkXsxFmHeODnvn3Q0cGiRTJNGlKT/YrVDiPKKO0b1RQVwCogND//V2trjzefJhUmLTy9sIdRjyTvpA4GjFvqXGqiGuQdRdoGaOuh55fZLx362Eo7j1IK1dW4eBEhIbh8GVOmYOdOmSYNqXS/Yo3HiDJK6Jt6TXALcAP6Abdfvnws/qhGVLM5bvPmuM0/O/7MGchp+OfVgOAlcva+ngXvsB7mzgDI3bsKupt44OfBg+jeHe7uOHBA+oGfzNivWOMxooyyhOb2KwkAf2AzsAWYB8DEpAZACj9l0ZlFZm3Mbnvd7tS2k6qDbJ6qf5C1/fUYplEnYPTBmyNy/1rFAz+Dg/HgARYswI0b6Nat8Z+qKy0NwcE4cAD9+oHDwfTp0GvuWtea+wtcF53FRCnZfcAdMAD2A+bijwQigX+cv3+s/7ox69S+EfoyDxmbcT8E5lPRZzXa2Tb+IzKRz8DP589x7NjbGrx0qTRbxVMyYkRrlPaNqq19wLfAWmDpm0XeUotSF51ZZNTKKNEr0aKdhWrja5bSu0j/DY/PwnoBJqfAwFxB9xH3W+7eDQMDuLtj8+amvDZ/hwbtV6x2GFFGaWtUDT0HvIE7wGXgdXdb7UbokgFLtNR38bTXy4jcgs1STLsHvfaKuMmbKZfifstTp9C/f+M/VZc67Fes8RhRRil1Ew/MBYYBt4DXa1XklOR4nPFo2aJlPCe+y3tqOgeG4HEE0jag4glsl2MkV0HLiIgbjlwuBg2Std/yzdClixff7FfM8kU/VYiWUeXRiNHLbybI/w588vojQvYk7lkdvXrF8BUrh6/U1tJWv0yJCI9OImUttFug51ewmgetpv7TaHqyBQXgcrFvH6qrMWcOkpJkmnKZn48//sDu3TAyAoeDffuUttem+n2tslLL4fcs6RtV/9/CImAR8Az4G7AWf3S/9L7HGQ+BSBC7ONbG0Eb8oTplKl6H6e7PaGkKh19l2Mmj0WTrDPyUseFYXY3QUAQHIy4OLi44fRoODlJeornU6WttHrUso7RvVB1EA+6AC3BGvOeHiIi2/b3tl+u//PDhD8s+WKatpa3qCKVU8xy5B5D2GwwHYFgIjIfK/Q5vBh117tyMMe8ZGTh48O3wUS4XrZS0YinVRIwooxSzCYB1QBBwEBgv/iinJGdx6GKBSHDD80YPox4qDU96lU+RvRtZO9FpIsZeRrue8r28eL2k2gM/ZVmus6ICEREICkJaGtzcEBsr/bqhlJLQMko1LA+YDxgAiYAZAIFIsOPWjnXX1r3pCVV1hNJ4cQ9Zv7/eCmlSIgzkPB4rMRHbt+P0aUyciB9/hJMTtGX469G4/Yo1HiO+Hpb0japhv9IpwBvwBn4EtAGkFqV6hnq2bNEybnFcd0OJTSwmZlqciMytKIiU+1ZIJiYm4n5LccNxyRJkZ8s08JPx+xUz8WtVDDqLiZKLSsAXuAAcE690Jx4TKp4dr2ZjQsWDQMtSYbsc3ZaghTx3EyosxMGD2LEDXbpg+XKZGo51Zi9xOHToktphRGuUYpg0YA7gACSIN+9MLkz2DPU0aW0Sz4nv3K6zqsNrGvFecik/QVSNXl/D6nTTxzA1hfjJ+9QpzJyJs2dlWvFDvPId3a9Y/dEyStURDHwDbATcoKZLNL0ew7QOLU1g96MMY5gaIF4yacsWFBXh00+RmSl96ft3v+LXi36yYL9ijUfLKPVGGeAFZAPXARsAd57e8Qz17NS2U6JXonlbRc0olyfxGKb0jWjfH8MOwXiYHK/95AmCgrBrF/r1w5dfYuZM6ees19mv+OBBul+xZmDEa1YtLS2Nf78EgM/nqzqEBtwEBgKmQCxgUymoXBW1anzI+KWDlkbMi5C2hqog08pCpKxBWHcUJ8DxEj4Kl2MNjYnB7Nmws0NBAa5cQWQkXF3f1tDGk62sBI+HceMwevTry8XHg8NRuxrK7F9geZI2U0a0Rlnyiompk0DE8zs3AEHAVABx+XGeoZ52Jnapn6fKtli9UjMtz0Hmtn/HMCXAQG5dt+XlOHIE27ejuhre3hJXTW4oWc3ar5ipv8Dyp5azmCjVKQQWAlXALcD8zd7x2yZtU4Mdk0puI2OLIsYw3buHvXuxbx+GDpVp7ibdr5hlaBlls0jAA/gU+AHQiXkYszhssb2pfdLSJGMDGcY9KpF4DFPpXfT8PwzeLa8xTLU3fJs/H/HxsLSU5uerq3H+PA4fRlQUpkzBpk0YPZoOXWID6fpGhWn+46wsLCwsLCzMjQz0h27MFlZHf2bR1tjCwsLCoku/r/+sBkD+iV4z0c62l63d+NWR/CY8r7Okb5RJaoA1wBLgKLDmVU3VqqhVc07M2eC0gevKZW4NJSI8DsfFD3DrM3RxxbRs2C6XSw19/hxBQejbF76+cHLCo0cICJCmhiYkYPlydOmCjRsxciSys3H4MBwdaQ1lC6l2tX+r6sbXPT/4JVVAqi5728w/VVnr0KvoZbbDf0muIFUZmx1tvS6WN3wl2WNQN4WFhaoOQew+IUMJ+ZiQZ4SQqw+udt/W3ZXr+uzVM3ndQP6ZCqtI7iESbksujSD5YYSI5HXhjAzi40PatyeuriQyUsofTk8nfn4Ca2vSuzfx8yM5OfKKipkY8wuscNJmKuNDfdVfx8MNZ4Xb6kBY91BN0vloE5fwPi2h3WPObMuPziXWjB+lxv3q8sOM7nke8AXwLbC8rKrsm0ivSzmXgpyDxnUdJ8d7yDPTmhfI3Y/0TWjvgKH7YTxcLlcViXD2LLZtQ2oq3N2RnAyLpk+vf/YMJ0++XnrExUUnOJglu70z4xdYGaTNVLYyWnnj+DmzWRe76wBCoOK8T1/rL3XMP/TeGOAzzKAgn286yEwbgJahudmLqMfV4oXVKBUrB74GrgEXAfvz984vjVg6sfvEZO/ktnpKWvdXOpV8ZO98vQ6T4wW810cuVy0txaFD2LIFxsbw8cHcuU1+eS5eckm84ufkyfD1xaRJdNEQCjKW0VdXj13o4hptpQ1Ad+DqmHsBlu2FeWFfOc//tvfdANQevkRIE2aQ1J6gvWLFio0bN775I5/Prz2Gy8TEpPb/KOjRJh9NFApnv3jR+8mT4LKqV5tTXG89u3Xg4wNjrMcwMebyXGQG4MGRCtMZD63/qNHriEcEj+4288r5+SaBgeDxMHUq9u8vNjF5AiAzs5Gf1RKJzLOy3g8LQ1gYBg+Gm1vR778XlpcDQEaG6v+u6FEFHIW0DVIZOg5eRHhYjd+dL6zzcc2tb+1G/JZVEbey9yj/e0JCiOjxrnE2PlerG7yabDFQTSYiZCshHQj5gxASkRXReXNnTjjnRdULVQdWn+LbJNaNnDQlyX6k8h+5XLKqinC5xMmJmJsTPz9SVNTkn7x7l/j6ElNTMnAg2bqVsKZnkJKWDCWs7LSb1ZS9BeJOflFxblZhFSFEwI9eMaDzwrDn5NXlz3qIXzGl+zv2WHKevmJ6TRU99EWETCVkMCH3iiuKOeGcbgHd/rz/p6LvKkum/OvkylRyxpJkbCU1jfzONFFBAVm/nlhYkBEjCJdLamqa9mN5eWT9etKjB+nZk/j5kezshk+nL140j7SZSl3CRCW8OZbTD/FfvyoVpPw+rZdFx06dzDrbTfO7+ERICBHxI78f17tr9269xq66UFC30fpOBKwpoykpKcq94Z+EWBDiQ0gVN5XbcVNHn/M+5dXyqVANkyJTkZDkh5ELQ0hYd5IdSIQNP7o0VXw8cXMjhoaEwyFNjaW4mAQGkhEjiJER4XDI9etE1KTxAEr/WlWGZiqJ6kuYuG/Bz89P1YEonBJ/CwWE+BFiRsi5p+VPXY679NnR5+/8v5V196Zl+noMU09yfiDJPUREjf3/tgkqK8mhQ8TBgdjYkPXrSXFxE36mooKEhRFXV9KuHXF1JWFhpFq6Uk6Li+aRNlNGvGck7JhTryyPgPmAHpDIS4vxOe+w0H7hHy5/6Ovoqzqwf70ew7QRra3Q/zeYOzf/ko8fY88e7NwJe3v89BOmTm1s5LtIhNhYhISAy0WfPnB3l3XDOYqik0E1TSjgBSwtKF/y2dnPcopzwueGD+o0SNVR/ev1GKbtMB6GD0NhOLD5l4yJwbZtuHwZCxbg779hbd3YD6SmgsfDoUMwMIC7O9LTYWbW/DAoNqNlVHmkG0IhNfG2H6HAKV7a4y/ODVzksOj4rON6OnqKvGn96sm0/D4yt+J+MMydMT4ObW2aeYsXL3D0KH7/HQIBli6VuPzSW/n5OHkShw7h2TPMnYsLF2Br28wYxBT8tTIIzVQixfQtSEEcBhv6RhUpnRAHQlweliWNCx7nsNshsSBR1SH9q/gOiXUjPEMS70NePWn+9bKyiK8vMTIiU6eSyMjG3gOVlJBDh4iTEzE0JG5uTfgBipIaI1qjhPaNyo4AgcCPhPxvT6Lo+2inzwZ/tvrD1braDJg4Jl6H6Vk8bLwwPRe67zXnYuLllwICcOsWFi1CYmKDi89VVeHSJfB4OHsWw4eDw8H06dBTQcOcYgNGlFFKVnxgCVCQ//zowjP/q6ipuOpxtZdxLxUHJd5LLvV/qCyC7RcYyYNOy+Zcr6wMBw8iIACGhuBwwOWiVSsJp4pfHPF4OHYMNjZwd8f27WjXrjl3p6hG0TKqvi4ASwjmbf/7w5+uzVk9arXPEB9tLZXuCiOqRt4xpK5Hi1awXQ6r+dCSdrui/8jIwK5dCAmBkxNCQjBihORT09LA5SIkBC1bwtUVf/8NK6vm3Jqimo6WUeXhy20PhgpgFRD+5MWm+ad2VwurYzxjbI3k88JERoJy5OxD+ia0tkT/DXzdIc3JVCjEuXOvl1/icJCZiQ6StjJ58gQ8Hng85OVh5kyV7LIpv6+V6WimktAt7ZSn9moIzXATsCcoPnDn/+x3L5vYfeJ1j+uqrKFVRUhZg1BrPI3CqFMYFwNzZ5kz5fOxYQO6dsXatXBzQ14e1qypr4ZWVIDHg7Mz+vVDQgLWrEFeHgICVLJTsZy+VjVAM5WEEa1R+oqpaQSAP7DtafkPn5w4KiIPYhfH2hg2d+SQ7KqKkLYBOQdgOQcTbqJNoyM2G5KQgKAgcLlwdkZ4OPr1q+8k8Q7vwcGvXxy5uuL4cbXbX5PSPIwoo1QT5ALuBAb7bnt9G/XDiuErVg5fqbKe0OpiZG5D1g50nokpKWjVSfYrVSM0FAEBKCgAh4N792BkVN954i02jx1D585wc4O/v+TnfIpSNlpG1UIw8PXT8k8/PhatpxMVtziuu2F31QQiKEfWDqT7o9NETPgbbbrKfKWCAgQGYvdu2Nlh+XLMnPl25/e38vJw7Bj27UOLFpg9G7Gx6NatOeFTlCLQvlHlkal7ng9MJ/h9T+LCfrv2zbWbd2XRFdXUUMFLZAYgzAbFCRgfi2HBDdTQhjNNSIC7O/r2RUEBoqMRGQlX1//W0OJiBAVh5EgMH44nT7B/P9LSsGYNM2soS966gGYqmZbK+yW1tFQfA1NdAJYUvRr/8dHkNvqGe5z3dHlPFdudi6qQewgpa9FhBPr9jHYyvs6qrASXC39/VFbC0xNeXnj//XfOCA9HcDCuXcOECXBzo7t0UGpB9SWMltH6VACrCMKOpEz48sKpdWPWLRmwREv5u/WKqpF7EHd/Qvv+6Pcz2jvIdpncXAQFYf9+2NvDx+ed5ZdqL7Y0aBDc3DBzJl1siVIj9H/1DHQTWFBa2XP6sfatWjxI8Ero3K6zskMQ1SDvKFLWok1XmZdiIgSXLyMoCNHRmD8fN2++MyI+NRUhITh0CObmcHNDRgZMTeURPUUpFSPKqJaWlp+fHxu6RxvzekjTqYyx3hEXf3b8WQWNUCLCo5NI+g6trTDiKIw+kOEa4uWXAgKgr4+lS3Hw4H9HJT16hFOncOAASksxZw6uXYON6oZtUVSzMaKMsuShvrGpEbmA+4tqkSvXSFvrnwROgkW7pm+dLg/iApr8A/SN8UEQTB1luEZWFvbvx759IkdH7YAAODnVOlZairAwhITg9m24uGD7dowY0djqymqAzu3RPGo5i4klGpwaEQwMuXzf0GZbzsxey87NP6fcGkrwOBwXBiJjCwZtx7gYaWuoSISoKDg746OPAODo0Wwu998aWlWF8HDMng1LS/B44HBej3UaOVIDaijo3B5NpJazmNiNDyx5VZOz4JRJlVCY4JVo3tZcqfd/GoU7qyCqQd/v0cVV2p8WL7+0dSuMjMDhgMdDy5a4e7emnsWW9u1D27aKyICiVIsRZZTFfaMXCD79O7/7LB7/xw/XcQZylHrzp1FI+g7CCvT9EV1mAVK0DQnBlSvYuxfnzmHGDJw8WWs6e2qqyc6duHABrVrB3R137qBjR0WET1EMwYgyypK+0f+qAFZVCU8sCdMrqWx789M7ndrKPqVSakUxSP4RFU/RZxWsFkCaSaUFBQgOxt690NXFwoUICICxMQDg8WOcOIGQEPD5Oh9+iFOn4CDjACmKUi+MKKMsUavT+iYh828/1XfhVn070k+pjdB/4pD8I8pz0HsVui1u+nqg4sXng4Jw+TKmTMGuXf92fZaVITgUPB7++guTJmH9eowdKywqAjveRYDO7dFEdBYTwwkA/2rhpq8vtSwsH7Zjyo4OBspaYqM0GXfX4dlN9PkO3Tyh1dT/gz56hCNHsHMnOnQAh4N589CmDVBTgwsXEBKCyEg4OWHBAkyaRHfpoNiJEW/qWTKnHsgVkQ+znwUO2av9keUWritXSTW09C5iZuPKVHQYAedMdOc0pYZWVYHHw7hxcHBAbi7CwxEfDw4HbfJSsWoVOnfG2rUYMQLZ2eDx6E5HFJsx4qGeHa3R4BqRz6bYlncKPox022lsYKyMez5PR+qvKLiEnl9iWHAT90RKT8ehQzhwAP361doL7vFjBJzAwYMoLsbcuYiJQXcVrTJFUQyj+gdqFjzU84XEs+BF4tyTgi+H7p7Za6Yy7ll+H2nrkR+Knl+ixxdo0fjaxs+f48wZhIQgPR0LFsDLC9bWQEUFIiIQHIy4OEyeDHd3jB2rGeM9KUpeGNEa1WgXqoULDyUJ/rzvuGfsup5deir8hi/zkPo/5J+BjTecs6Db+L6Y4pXnT5x4vRXxjBlooS1CbCzWh+D4cQweDDc3qdaZZ890F7ApWZqpJLSMKk6FQLTiedUfn4bputnvPeIy4+7du4q94atHSN+EB3+g26eYmgm99xs+vaQEPB527EBVFTw8kJ4OExMgLQ3ruDh0CAYGcHdHVpYM79zZ8+8NbEqWZioJI8qoJg6/v1khcLl4r+x0utOeaXuMWtW7M4b8VBUh3R85e2E1H1PS0LKh3wDx0KXgYEREwMkJ/v4YOxZaJcU4cQLBwa+32DxzBvb2io2ZojQFI8qoZvWNCmqE6ysEG766qO9sG3JoxnTF3u15BjK24OEJdPPA1HToN/Tq//FjHD6MwEAYGsLNDVu3wtCgEpGR+CQE0dGYPBm+vpg8ub7dPCiKkogRZVSD5D6v+ji5MCckecxv4w4ZtjJU4K2KYpC5Dfzr6Lqw4X3lqqtx8SJCQhAdDRcXnDyJ/vYixMZiLQ9Hj8LeHm5u2L+frpRMUbKhZVRuKgSBQtH/+cfpD+x4JHBqPY1Q+fQrERGenEXq/1D1D3osw9CDDbyFz8zEgQM4cACdO4PDwcGDMHiUiaNH4XoY+vpwd0dyMszM5BDVf7GkB02MPcnSTCUiUqm67G3exsjc3NzcvLPdV9FVhBBR0WW/CX179OzRd9x3lwpF9X7SEHEYfn5+0kXCLIWF5UPvPNX75tL04opiRd2k+jnJ2ErOWJJLI0gel4gEkk6sqCBcLnFyIh07El9fkp1NSHExCQwkI0aQTp2Ijw9JSFBUkBTFPtKXUZv5pyprffIqepnt8F+SK0hVxmZHW6+L5fV80nAE0pZyhimvPlFc0XrXrfcu555X1D1ePSHJfuSEEbkylRTFNXBifDzx8SHGxsTJiXC5pPpFJQkLI66upF074upKwsJITY2igqQotmruQ31N0vloE5fwPi2h3WPObMuPziW+alf3k5rxo3SbeRuGqsgrnautffb43clegw631VPAYpold5CxGY8jYO2GiQlobVnvWWVlOH4cu3fjn38wbx4SEtClKAHBwfjiOLp3p2t9UpRCSV9GK8779LX+Usf8Q++NAT7DDAry+aaDzLQBaBmam72Ielz5zifVgOaV0fLq6Fc1MxOfCtu3OrViuLO8L0/w9DIyA1CcgO4cTMuBXvt6z7t+HUFBiIjA5Mnw98fobo+0jh6B0z7o6OCTTxAXB2trecdGUdR/SFlGdQeujrkXYNlemBf2lfP8b3vfDUDtwUqEAFr1fNKI2ru2rVixYuPGjW/+yOfzay/ob2JiUrv3V0VHBXfy3cze54bcGvpR+y0GAoM6g3Ul/az4tIbvW1SYX5UV3KFoP9HS/cfYTW/YbhMz83evXF6uc+HCezyeiUDQwtMTW9eWGcWEVq/ZJ0pOLnNyKv3++1cDBqjw7wq1OukZ+Q3K82jr1q2ta/2/iiFR0aPNOfry5UtrqdofsvYG1Nz61m7Eb1kVcSt7j/K/JySEiB7vGmfjc/XlO59UN3ihZsSgAmWVt+89M7uW1+rGQ560P5uSktLQ4Uo+SV1PTpuTy04kP0zSWfHxhMMh7dsTV1cSeUkkirpM5s0j77//uuuzuuG/bCVpJFPNwp5kaaaSSLdQHim5n82vBiAsun78Ymn3nmYt7SeOLjx5PLUS1ZnHuHmjJw00eOcTjXmiTyr8WiAamFDQtb9Z4fDOs+R23Rf3kLAc4T1QlgrHSxgTCfO6vQRlZQgKgr095s1D1664d+0Jd+AGJ28bLZ8v0K8fsrPB5cLZGboa85dNUWpDuod60eOzK2ZvuFUmIjpGH3ju3T2lLbQd1+wYM39OnwNV+pau2w6Pb13PJ+qvtCorr3S8Fp7eKz44u4+b3K77nyH0qfUOoRevG8LjwckJmzcKx2j9qbUnCBuj4eKCgwcxcqTcgqEoSjYKahU3nTgMJo8b/St/zZMXOhfuOZRXN2tM6NsnBZGQ5IeRi0NJWHeSsZXUvHz35NJSEhhI+vUjPXqQ9etJcVwG8fUlpqZk4EASGEjKGxlGplrsefojbEqWZioJI2YxEabOqS96+fD20wm9O+QUlm+Z0O2LZl7NxMQENS+Qux8ZW2BggZ5fofPMd3dDqt383PJrpWN5uNaeIGxLg5ubuiyWXLvDXuOxJ1maqSSqXzKZscs2X77/q/X7P/BfdrU3jWml2+xfoIoCQzgxTQAADAFJREFU3AtE1nYYD0Of1TAeWue4eOznjh2orISnJ7wGJbzPDQKXi0GDwOHg449pvydFMZPqSxgDyyj/5ZPI3InjuqaXVf7Pxmhlcy9Xewh9z6/eHUJfu/n5+bySDwt5Wrt2oqICnp5YtAimps0NgKIoRVJ9CRMPGmXOeqPn7m3v0Orr91uaW74fradj1Ywr/XcIve3yOkPoazc/F3uIvGyi3wv9dxFQDofu1UFR6oIRZVTlMYgVvizkpk5d0O/2i6rlXd7b1IR5AxKIqpB3HGkboK2Hnl/Cci60//M8Xrv5+cWM/JF5f2gF/bsIqJsbDBW5vB5FUfKm+hLGkDIalrVXX2dZnw6GJq0v6unYyXiVqiLk7EfW72jXC7Y+dYZ/3rtXFB3dQdz8XOJe5dUxrO3pYMTFwcUF3t5wcJBDGszAZ81uE2BTsjRTSRjxpl61npY/DUqc6TUgnmCeWZsgQKb91l/cQ9bvuB8Mc2c4XsJ7vWsfFDc/jx0znDABu33ShmUHY9sB9OsHNzdwuWjVSj6ZMAZ7/r2BTcnSTCVhexk9lf5HcQVn2WDdtvoXdLXHyHIJyUPoa/d+es9/Hs/ZYhMXCb8HWLAAf/1FFw2hKM3AiDKqki3tCsoLNsfN+eKDv9voORm2OgpIuY5cg6vQ1+793LM04YM7QQg4UTZoEN3siKI0DyPKqPL7Rk+kH8soWvLDhzDQPdxCW8rZ8ZKH0Ndufn7h+nTziuOtj+7FXSEWLkRGxqPCwvf69pV/MhRFqRQjyqgyPXnxZM0V988G/zWha7+2+icAiTvB1aP2EPoRx2oPoX/T/Bw/Vnhg/p8D4oOwOxqTJ2PLFjg5ic8xYcCbNOVgSQ+aGHuSpZlKovq35Ep7U08I2ZO4J+np17+NF7ZssU5H60sphjRJGEJfu/m5YlrWfK0jBscPoEMHcDiYN4/utUlRbCDdQnkKoqWlpeiO0QelD2afGG1r9M3G8ZatdW/raH3VtBpK8DQKV51xZTLadMW0HAwMENfQhAR4ecHaGlcvVh6ezsvsMm7JH6MNqkoQEYH4eHA4tIZSFEtofmtU3AiNvv/Nnmloo+elhZ+bNKRJwhD62s3P7yYkzK4MbnXmKOzt6bR3imItDS+jD0offH7Wc6FD+oyeLXR1jgCjGv8ZCUPo3/R+Tv+o9Bsrbq+ru1Fejnnz4OmJLl0UFD9FUcynsWVU3Ag9lb7qqIvO+60mamFn40Oaag+h771KPIT+TfOzulL0/ZhYl1chLcN5sk17p6OXNRJ7kqWZSqKZfaP3S++PDxkjEP0SMU+3fas9WghppIYWxSBmNiJHQacVpqRiWDDe6/2m9zMx/PHJDzakC2zmX/Nq2bvr6x07nJykXTqk9i5amo09mYJNydJMJWHEgCe5t0aXRow/NEPYsa2tFg42NKSpviH0ZWU4HoIdOyCsqF47/OK2USH6N6Lh4gIeDwMGyDdOiqI0ACPKqNxdXKAHfAF4SXwdX98Q+je9nwuHZIT2Omh15SDudYe7O44ehIFB/dehKIr1NLOMAqkSj7wzhL6sDMf3YMcO6Lx68fOA09v7heimpGPBAty4gW7dlBgzRVFqSTP7RutXcgdx7jjbB9UlmJiAj8IT8oaKez9zeQnnOnslFltNecXT/ZyDBw+wfr3cayhLuufBpkzBpmRpppJo7Jv6WuquQl9W0V788l3/ZfHP9ifGZuxoIarGokXw8ABrflEoipIXTX2oB/DOEPpRpxJu6wZ9gZM80XK76EvvB5k+uAidCQjwpzt2UBQlMw1tjf53CH1ZG2dx87N9+aO1tkdG3d2lbWIMDgdz56KtlOvjURRF/ZeG9o1GfojyXIyJSngv0usn557WVS8P8v5sMe5KqcNHnXO1I8JeT3unNZSiqGbTzNbom6lHnZ+nfm8Z8kHqAW2HfuBwMH069GTaI0Qe6CQQjcSeZGmmkjCiNSp34z4oa7E/6EblwAjBxKFDoX3rb0RGwtVVhTUUdBKIhmJPsjRTSTTzFdNN7aHoNgCLN8LRkb47oihKoTSzjCI9XdURUBTFFox4qFfS8HuKoigFYERrVOWvuZSDJd3zYFOmYFOyNFNJNPNNPUVRlNIw4qGeJdjTccGeTMGmZGmmEhFpCB+c/HJ8X0vzTp262Lv8dqNYREjVZW/zNkbm5ubm5p3tvoquIoSIii77Tejbo2ePvuO+u1QoauSa4jD8/PykikQdSfu3rb7YkylhU7I0U0mke6AWPbpxPt9yzFAL3adnOB/90PHo7V/srn3Wd/+4lMMz9P89qeLPL/p/35F3+TvbvC0Tp2euStg9vnUD12TPQz3NVCOxJ1maqSTSPdRrdx4xZZhFKy20MBs1qsezJ09F755Tk3Q+2sRlTp+W0OsxZ7Zl9LnEGqnu8S9J7ep6P5fqZEV/Lu0TgUJvqqrPmRMMezKV9LlCM1X0TVXyNyA12Rq91WmbRvdYeOaZiFRd9rYwtOhuZWk7wm1rbLGIVJ6aazzr+EtCCCGVYe5mLkfLG7yUpBik+lwuF1HrYOjfAKOCoX8DGvk3IIksZVT07PJXgx28Ip6KCCGi0vwHxdWEVDwI9baz9rr0svLUHKO3ZdTN1OXoy4YjoCiKYhipSqL040ZfJm6avfy+18ljU0y1AGi9Z24JALqWkz2nfv9/dx6TUZ1NC/MLROimTZ49ftrOvJNuY+HKlCZFURQjSDngqTp7v/vC82MPBC/uIV7kg5Tcz+ZXAxAWXT9+sbR7T7OW9hNHF548nlqJ6sxj3LzRkwY2XEYpiqLUmnQvpKove1tNDK42bd8SAPQn/J4UZHNw5uwNt8pERMfoA89Nu38c31GbFEX9OH/5kftV+pauWw7/MsGMDk6lKEpzsWUEA0VRlILQhiJFUVSz0DJKURTVLLSMUhRFNQstoxRFUc1CyyhFUVSz0DJKURTVLLSMUhRFNYsKyyj5J3rNRDvbXrZ241dH8jVy9KowY/+iUb1MDfT7/5QiFH+kgVmL8k59NcHOysLc3NJh1sbYEgKNTPO1yqs/jOxlbWXVpbOt43Jebg00OVkANWkbR7Zp+dG2hyJobKbV0Z9ZtDW2sLCwsOjS7+s/qyF1pjIsTSIfr6KX2Q7/JbmCVGVsdrT1utjwOlDqSchPvnYz88aaoYPWJgsIIZqZtfBhTETso1ciUvPktIdN3+/iazQyzdeEZU8LX4kIIRUZm0ZbuIW+1ORkiSBn18xJn0y1cQzIE2rmby8hhFRd9raZf6qy1idSZqqy1qi8liVlNO0OdqMGdzPU1/r3A43M+t1VaDUyzde025matNICIBIIRFra2pr5nYqJnhz9/oy9n1cPHQAa+ttbL2kzVVkZFRXk800tzLQBaBmam73If1ytqlCUSLOzrsk4GJw9/uNhupqdJqqjVgyw6WQ2LLhnwLpJBhqbLPkn/If9Zt9+1f/fjS00NlMAFed9+lpb9RzpHhBXQqTOVJV9o+T/27ubUAjiMAzg72KRNuIyZme16/ujkJxcOPg4oD2QHBQlxU05iFaJm/JRuBA3F8JBDkSkWA6SQmQXw+5SDg5E0X44rAOJbH/aej2/0xzmME9Tz+XfPPPxWvP1rYywTe27W+tomEof7DPHaRjHJCIKL+nfszkd6/WOgVHrI9ew92u9Q9RqKdS9i8MzKWnzLZv2c9vFyXK7brKuc/Up0KRBq9EQ+W2WlOhHs6Q8sE39tkI7PVwhaRjHfCckJq+mNHRx8cjLM6xb3dk9XWrNTzSlVk+cbfcUVY1dSSyT+keTY7VEkcbyxkrd4b7LF+A7DVqNav/lLCnP1J9WaHnGJCIi7+3R3uWTj8h7vz+zYDemGZgO7IblWKzXLlVVVdtcU3JB98Z8c2oey6S/MZr8pydg3/LernSVZiWlJGcWdyzdeIL3IH/H45ysTlSkmIjwaEkxVU04PBxTP6+2yGFR/p9sK0pS4/wD45frPh6vzTXKsl42ZJW1zdqffYzD+r1stWX4T+qZJnUfjJgzDbJeH5+Qbe5evg48KfZGAQCE4CsmAAAhqFEAACGoUQAAIahRAAAhqFEAACGoUQAAIahRAAAhqFEAACGoUQAAIahRAAAhqFEAACGvC7pytspgIrsAAAAASUVORK5CYII=",
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"(x,betaMu,betaSigma,beta,ySigma,ys) <- simulate (birats 5 5)\n",
"toRenderable $ sequence_ [plot $ line \"\" [sort $ list x `zip` y] | y <- list ys :: [[Double]]]"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"-- http://www.openbugs.net/Examples/BiRatsdata.html\n",
"biratsX = [8.0, 15.0, 22.0, 29.0, 36.0]\n",
"biratsY = [\n",
" [151, 199, 246, 283, 320],\n",
" [145, 199, 249, 293, 354],\n",
" [147, 214, 263, 312, 328],\n",
" [155, 200, 237, 272, 297],\n",
" [135, 188, 230, 280, 323],\n",
" [159, 210, 252, 298, 331],\n",
" [141, 189, 231, 275, 305],\n",
" [159, 201, 248, 297, 338],\n",
" [177, 236, 285, 350, 376],\n",
" [134, 182, 220, 260, 296],\n",
" [160, 208, 261, 313, 352],\n",
" [143, 188, 220, 273, 314],\n",
" [154, 200, 244, 289, 325],\n",
" [171, 221, 270, 326, 358],\n",
" [163, 216, 242, 281, 312],\n",
" [160, 207, 248, 288, 324],\n",
" [142, 187, 234, 280, 316],\n",
" [156, 203, 243, 283, 317],\n",
" [157, 212, 259, 307, 336],\n",
" [152, 203, 246, 286, 321],\n",
" [154, 205, 253, 298, 334],\n",
" [139, 190, 225, 267, 302],\n",
" [146, 191, 229, 272, 302],\n",
" [157, 211, 250, 285, 323],\n",
" [132, 185, 237, 286, 331],\n",
" [160, 207, 257, 303, 345],\n",
" [169, 216, 261, 295, 333],\n",
" [157, 205, 248, 289, 316],\n",
" [137, 180, 219, 258, 291],\n",
" [153, 200, 244, 286, 324]]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"let n = length biratsY\n",
" t = length biratsX\n",
" prior = birats (integer n) (integer t)\n",
" posterior = [ (bm,bs,b,sd) | (x,bm,bs,b,sd,y) <- prior, x == list biratsX, y == list biratsY ]"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"let initial = (list [0,0], eye 2, list (replicate n [1,1]), 1) :: (RVec,RMat,RMat,R)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"--- Generating Stan code ---\n",
"data {\n",
" vector[5] x_stan_0_0;\n",
" \n",
" \n",
" \n",
" \n",
" \n",
" matrix[30,5] x_stan_0_6;\n",
"}\n",
"parameters {\n",
" \n",
" vector[2] x_stan_0_1;\n",
" vector[2] x_stan_0_2;\n",
" corr_matrix[2] x_stan_0_3;\n",
" matrix[30,2] x_stan_0_4;\n",
" real x_stan_0_5;\n",
" \n",
"}\n",
"model {\n",
" \n",
" matrix[2,30] v_0_1;\n",
" matrix[30,1] v_0_2;\n",
" matrix[30,1] v_0_3;\n",
" \n",
" matrix[1,5] v_0_5;\n",
" matrix[30,5] v_0_6;\n",
" matrix[30,5] v_0_7;\n",
" \n",
" matrix[30,5] v_0_9;\n",
" \n",
" \n",
" \n",
" matrix[2,2] v_0_13;\n",
" vector[2] v_0_14;\n",
" vector[2] v_0_15;\n",
" \n",
" \n",
" \n",
" v_0_1 = x_stan_0_4';\n",
" v_0_2 = to_matrix(to_vector(v_0_1[1]));\n",
" v_0_3 = to_matrix(to_vector(v_0_1[2]));\n",
" \n",
" v_0_5 = to_matrix(to_row_vector(x_stan_0_0));\n",
" v_0_6 = to_matrix(v_0_3) * to_matrix(v_0_5);\n",
" v_0_7 = rep_matrix(to_vector(v_0_2), 5) + v_0_6;\n",
" \n",
" for (i_1_1 in 1:30) for (i_1_2 in 1:5) {\n",
" \n",
" v_0_9[i_1_1, i_1_2] = x_stan_0_5;\n",
" }\n",
" \n",
" \n",
" \n",
" v_0_13 = quad_form_diag(x_stan_0_3, x_stan_0_2);\n",
" for (i_1_1 in 1:2) {\n",
" \n",
" v_0_14[i_1_1] = 0;\n",
" }\n",
" for (i_1_1 in 1:2) {\n",
" \n",
" v_0_15[i_1_1] = 100;\n",
" }\n",
" \n",
" \n",
"\n",
" \n",
" x_stan_0_1 ~ normal(v_0_14, v_0_15);\n",
" for (i_1_1 in 1:2) {\n",
" x_stan_0_2[i_1_1] ~ cauchy(0, 2.5) T[0,];\n",
" }\n",
" x_stan_0_3 ~ lkj_corr(2);\n",
" for (i_1_1 in 1:30) {\n",
" x_stan_0_4[i_1_1] ~ multi_normal(x_stan_0_1, v_0_13);\n",
" }\n",
" x_stan_0_5 ~ cauchy(0, 2.5) T[0,];\n",
" to_vector(x_stan_0_6) ~ normal(to_vector(v_0_7), to_vector(v_0_9));\n",
"}\n",
"\n",
"make -C /home/jovyan/stochaskell/cmdstan /home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a\n",
"make[1]: Entering directory '/home/jovyan/stochaskell/cmdstan'\n",
"\n",
"--- Translating Stan model to C++ code ---\n",
"bin/stanc /home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a.stan --o=/home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a.hpp\n",
"Model name=model_3a4af8197fb3f3295b183f18be009847751d441a_model\n",
"Input file=/home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a.stan\n",
"Output file=/home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a.hpp\n",
"\n",
"--- Linking C++ model ---\n",
"g++ -I src -I stan/src -isystem stan/lib/stan_math/ -isystem stan/lib/stan_math/lib/eigen_3.3.3 -isystem stan/lib/stan_math/lib/boost_1.62.0 -isystem stan/lib/stan_math/lib/cvodes_2.9.0/include -Wall -DEIGEN_NO_DEBUG -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DFUSION_MAX_VECTOR_SIZE=12 -DNO_FPRINTF_OUTPUT -pipe -lpthread -O3 -o /home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a src/cmdstan/main.cpp -include /home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a.hpp stan/lib/stan_math/lib/cvodes_2.9.0/lib/libsundials_nvecserial.a stan/lib/stan_math/lib/cvodes_2.9.0/lib/libsundials_cvodes.a \n",
"make[1]: Leaving directory '/home/jovyan/stochaskell/cmdstan'\n",
"--- Sampling Stan model ---\n",
"/home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a method=sample num_samples=1000 num_warmup=1000 save_warmup=0 thin=1 adapt engaged=1 algorithm=hmc engine=nuts max_depth=10 metric=diag_e stepsize=1.0 stepsize_jitter=0.0 data file=/tmp/stan-7ab78c4890857200/stan.data init=/tmp/stan-7ab78c4890857200/stan.init output file=/tmp/stan-7ab78c4890857200/stan.csv\n",
"method = sample (Default)\n",
" sample\n",
" num_samples = 1000 (Default)\n",
" num_warmup = 1000 (Default)\n",
" save_warmup = 0 (Default)\n",
" thin = 1 (Default)\n",
" adapt\n",
" engaged = 1 (Default)\n",
" gamma = 0.050000000000000003 (Default)\n",
" delta = 0.80000000000000004 (Default)\n",
" kappa = 0.75 (Default)\n",
" t0 = 10 (Default)\n",
" init_buffer = 75 (Default)\n",
" term_buffer = 50 (Default)\n",
" window = 25 (Default)\n",
" algorithm = hmc (Default)\n",
" hmc\n",
" engine = nuts (Default)\n",
" nuts\n",
" max_depth = 10 (Default)\n",
" metric = diag_e (Default)\n",
" stepsize = 1 (Default)\n",
" stepsize_jitter = 0 (Default)\n",
"id = 0 (Default)\n",
"data\n",
" file = /tmp/stan-7ab78c4890857200/stan.data\n",
"init = /tmp/stan-7ab78c4890857200/stan.init\n",
"random\n",
" seed = 1067463274\n",
"output\n",
" file = /tmp/stan-7ab78c4890857200/stan.csv\n",
" diagnostic_file = (Default)\n",
" refresh = 100 (Default)\n",
"\n",
"\n",
"Gradient evaluation took 0.000156 seconds\n",
"1000 transitions using 10 leapfrog steps per transition would take 1.56 seconds.\n",
"Adjust your expectations accordingly!\n",
"\n",
"\n",
"Iteration: 1 / 2000 [ 0%] (Warmup)\n",
"Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:\n",
"Exception: lkj_corr_lpdf: y is not positive definite. (in '/home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a.stan' at line 71)\n",
"\n",
"If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,\n",
"but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.\n",
"\n",
"Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:\n",
"Exception: lkj_corr_lpdf: y is not positive definite. (in '/home/jovyan/stochaskell/cache/stan/model_3a4af8197fb3f3295b183f18be009847751d441a.stan' at line 71)\n",
"\n",
"If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,\n",
"but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.\n",
"\n",
"Iteration: 100 / 2000 [ 5%] (Warmup)\n",
"Iteration: 200 / 2000 [ 10%] (Warmup)\n",
"Iteration: 300 / 2000 [ 15%] (Warmup)\n",
"Iteration: 400 / 2000 [ 20%] (Warmup)\n",
"Iteration: 500 / 2000 [ 25%] (Warmup)\n",
"Iteration: 600 / 2000 [ 30%] (Warmup)\n",
"Iteration: 700 / 2000 [ 35%] (Warmup)\n",
"Iteration: 800 / 2000 [ 40%] (Warmup)\n",
"Iteration: 900 / 2000 [ 45%] (Warmup)\n",
"Iteration: 1000 / 2000 [ 50%] (Warmup)\n",
"Iteration: 1001 / 2000 [ 50%] (Sampling)\n",
"Iteration: 1100 / 2000 [ 55%] (Sampling)\n",
"Iteration: 1200 / 2000 [ 60%] (Sampling)\n",
"Iteration: 1300 / 2000 [ 65%] (Sampling)\n",
"Iteration: 1400 / 2000 [ 70%] (Sampling)\n",
"Iteration: 1500 / 2000 [ 75%] (Sampling)\n",
"Iteration: 1600 / 2000 [ 80%] (Sampling)\n",
"Iteration: 1700 / 2000 [ 85%] (Sampling)\n",
"Iteration: 1800 / 2000 [ 90%] (Sampling)\n",
"Iteration: 1900 / 2000 [ 95%] (Sampling)\n",
"Iteration: 2000 / 2000 [100%] (Sampling)\n",
"\n",
" Elapsed Time: 2.90063 seconds (Warm-up)\n",
" 2.12953 seconds (Sampling)\n",
" 5.03016 seconds (Total)\n",
"\n",
"# stan_version_major = 2\n",
"# stan_version_minor = 16\n",
"# stan_version_patch = 0\n",
"# model = model_3a4af8197fb3f3295b183f18be009847751d441a_model\n",
"# method = sample (Default)\n",
"# sample\n",
"# num_samples = 1000 (Default)\n",
"# num_warmup = 1000 (Default)\n",
"# save_warmup = 0 (Default)\n",
"# thin = 1 (Default)\n",
"# adapt\n",
"# engaged = 1 (Default)\n",
"# gamma = 0.050000000000000003 (Default)\n",
"# delta = 0.80000000000000004 (Default)\n",
"# kappa = 0.75 (Default)\n",
"# t0 = 10 (Default)\n",
"# init_buffer = 75 (Default)\n",
"# term_buffer = 50 (Default)\n",
"# window = 25 (Default)\n",
"# algorithm = hmc (Default)\n",
"# hmc\n",
"# engine = nuts (Default)\n",
"# nuts\n",
"# max_depth = 10 (Default)\n",
"# metric = diag_e (Default)\n",
"# stepsize = 1 (Default)\n",
"# stepsize_jitter = 0 (Default)\n",
"# id = 0 (Default)\n",
"# data\n",
"# file = /tmp/stan-7ab78c4890857200/stan.data\n",
"# init = /tmp/stan-7ab78c4890857200/stan.init\n",
"# random\n",
"# seed = 1067463274\n",
"# output\n",
"# file = /tmp/stan-7ab78c4890857200/stan.csv\n",
"# diagnostic_file = (Default)\n",
"# refresh = 100 (Default)\n",
"# Adaptation terminated\n",
"# Step size = 0.190624\n",
"# Diagonal elements of inverse mass matrix:\n",
"# 5.83378, 0.0112468, 3.3753, 0.00816748, 0.0546415, 28.991, 31.7809, 28.8753, 29.5893, 32.5971, 31.8214, 26.9231, 28.9886, 35.1412, 27.7419, 29.5949, 25.4521, 25.845, 28.6024, 30.1287, 27.5221, 28.7733, 30.1964, 24.443, 24.1183, 28.9358, 29.8461, 27.2327, 30.0987, 36.5679, 28.2423, 26.8473, 27.0516, 28.2851, 24.1571, 0.0558018, 0.0573982, 0.0475101, 0.0541492, 0.0559067, 0.0548206, 0.0474767, 0.0567981, 0.0604249, 0.0473151, 0.0543833, 0.0464371, 0.0470357, 0.0543212, 0.0508562, 0.0472773, 0.0483793, 0.0524393, 0.0462552, 0.0447209, 0.0490314, 0.0532829, 0.0480532, 0.05432, 0.0630577, 0.0462804, 0.0473987, 0.0466123, 0.0485278, 0.0443185, 0.00601311\n",
"# \n",
"# Elapsed Time: 2.90063 seconds (Warm-up)\n",
"# 2.12953 seconds (Sampling)\n",
"# 5.03016 seconds (Total)\n",
"# \n",
"\n",
"Extracting: x_stan_0_1, v_0_13, x_stan_0_4, x_stan_0_5\n",
"--- Removing temporary files ---"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"([104.639,6.33149],[[135.57807844,-5.944705082836486],[-5.944705082836486,0.6518915304039999]],[[99.4325,6.38386],[91.3923,7.00401],[101.771,6.81125],[107.228,5.47218],[78.4097,7.14034],[111.988,6.22571],[99.0285,5.81951],[95.9207,6.78053],[117.306,7.49154],[98.8741,5.32688],[94.8058,7.65448],[96.4247,6.10042],[104.098,6.25678],[117.772,6.69302],[131.797,5.01399],[126.627,5.56047],[99.0814,6.09183],[107.099,6.11271],[112.064,6.43936],[111.084,5.99643],[116.036,6.25897],[95.0701,5.95259],[106.196,5.46068],[119.844,5.66106],[83.2827,6.92215],[107.478,6.64509],[108.289,6.5509],[120.509,5.51061],[96.445,5.30583],[105.957,6.38429]],6.85353)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"samples <- hmcStanInit 1000 posterior initial\n",
"print $ last samples"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"import numpy as np\n",
"import pymc3 as pm\n",
"import sys\n",
"import theano\n",
"import theano.tensor as tt\n",
"import theano.tensor.basic\n",
"import theano.tensor.slinalg\n",
"from theano.ifelse import ifelse\n",
"\n",
"def ascolumn(a): return a.dimshuffle(0,'x')\n",
"def asrow(a): return a.dimshuffle('x',0)\n",
"def quad_form_diag(m,v):\n",
" d = theano.tensor.basic.diag(v)\n",
" return d.dot(m).dot(d)\n",
"def lkj_corr(name, eta, n):\n",
" C_triu = pm.LKJCorr(name + '_triu', eta=eta, n=n)\n",
" shape = n * (n - 1) // 2\n",
" tri_index = np.zeros([n, n], dtype='int32')\n",
" tri_index[np.triu_indices(n, k=1)] = np.arange(shape)\n",
" tri_index[np.triu_indices(n, k=1)[::-1]] = np.arange(shape)\n",
" return pm.Deterministic(name, tt.fill_diagonal(C_triu[tri_index], 1))\n",
"\n",
"theano.config.floatX = 'float32'\n",
"\n",
"with pm.Model() as model:\n",
" v_0_0 = 0.0 * np.ones((5), dtype='float32')\n",
" v_0_1 = 50.0 * np.ones((5), dtype='float32')\n",
" x_pm_0_0 = tt.as_tensor_variable(np.load('x_pm_0_0.npy').astype('float32'))\n",
" v_0_3 = 0.0 * np.ones((2), dtype='float32')\n",
" v_0_4 = 100.0 * np.ones((2), dtype='float32')\n",
" x_pm_0_1 = pm.Normal('x_pm_0_1', mu=v_0_3, sd=v_0_4, shape=(2-1+1))\n",
" x_pm_0_2 = pm.Cauchy('x_pm_0_2', alpha=0.0, beta=2.5, shape=(2-1+1))\n",
" x_pm_0_3 = lkj_corr('x_pm_0_3', eta=2.0, n=2)\n",
" v_0_8 = quad_form_diag(x_pm_0_3, x_pm_0_2)\n",
" x_pm_0_4 = pm.MvNormal('x_pm_0_4', mu=x_pm_0_1, cov=v_0_8, shape=(30-1+1, 2-1+1))\n",
" x_pm_0_5 = pm.Bound(pm.Cauchy, lower=0.0)('x_pm_0_5', alpha=0.0, beta=2.5, shape=())\n",
" v_0_11 = x_pm_0_4.T\n",
" v_0_12 = ascolumn(v_0_11[1-1])\n",
" v_0_13 = ascolumn(v_0_11[2-1])\n",
" v_0_14 = asrow(x_pm_0_0)\n",
" v_0_15 = v_0_13.dot(v_0_14)\n",
" v_0_16 = v_0_12 + v_0_15\n",
" v_0_17 = x_pm_0_5 * np.ones((30, 5), dtype='float32')\n",
" x_pm_0_6 = pm.Normal('x_pm_0_6', mu=v_0_16, sd=v_0_17, observed=np.load('x_pm_0_6.npy').astype('float32'), shape=(30-1+1, 5-1+1))\n",
"\n",
" trace = pm.sample(draws=1000,step=None,init=\"adapt_diag\",start={'x_pm_0_1': np.load('x_pm_0_1.npy').astype('float32'), 'x_pm_0_2': np.load('x_pm_0_2.npy').astype('float32'), 'x_pm_0_3': np.load('x_pm_0_3.npy').astype('float32'), 'x_pm_0_4': np.load('x_pm_0_4.npy').astype('float32'), 'x_pm_0_5': np.load('x_pm_0_5.npy').astype('float32')},tune=1000)\n",
" print(map(list, zip(trace['x_pm_0_1'].tolist(), trace['x_pm_0_2'].tolist(), trace['x_pm_0_3'].tolist(), trace['x_pm_0_4'].tolist(), trace['x_pm_0_5'].tolist())))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"Auto-assigning NUTS sampler...\n",
"Initializing NUTS using adapt_diag...\n",
"Multiprocess sampling (2 chains in 2 jobs)\n",
"NUTS: [x_pm_0_5, x_pm_0_4, x_pm_0_3_triu, x_pm_0_2, x_pm_0_1]\n",
"\r",
" 0%| | 0/2000 [00:00