{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Stochaskell, version 0.1.0\n",
"Copyright (C) 2015-2020 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 FlexibleContexts, MonadComprehensions, NoImplicitPrelude, RebindableSyntax #-}\n",
"import Language.Stochaskell\n",
"stochaskell"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"pickStick :: RVec -> P Z\n",
"pickStick sticks = do\n",
" stick <- pmf (pickStickPMF sticks)\n",
" return stick\n",
"\n",
"pickStickPMF :: RVec -> RVec\n",
"pickStickPMF sticks =\n",
" let sticks' = vector [ 1 - (sticks!i) | i <- 1...infinity ]\n",
" rems = scanlE (*) 1 sticks'\n",
" probs = vector [ (sticks!i) * (rems!i) | i <- 1...infinity ]\n",
" in probs"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"dp :: R -> P R -> P (P R)\n",
"dp alpha proc = do\n",
" sticks <- joint vector [ beta 1 alpha | i <- 1...infinity ] :: P RVec\n",
" atoms <- joint vector [ proc | i <- 1...infinity ]\n",
" let randomDistribution = do\n",
" stick <- pickStick sticks\n",
" return (atoms!stick)\n",
" return randomDistribution"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"dpmm :: Z -> P RVec\n",
"dpmm n = do\n",
" let base = uniform 0 100\n",
" paramDist <- dp 5 base\n",
" params <- joint vector [ paramDist | j <- 1...n ]\n",
" values <- joint vector [ normal (params!j) 1 | j <- 1...n ]\n",
" return values"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"--- Generating Church code ---\n",
"(define (pmf probs)\n",
" (letrec ((go (lambda (j u)\n",
" (if (<= u (probs j)) j\n",
" (go (+ j 1) (- u (probs j)))))))\n",
" (go 1 (uniform 0 1))))\n",
"(define (bernoulli p) (if (flip p) 1 0))\n",
"\n",
"(define (model) (letrec (\n",
" \n",
" \n",
" (v_0_2 (mem (lambda (i_1_1) (letrec (\n",
" (v_1_0 (- 1 (x_church_0_0 i_1_1))))\n",
" v_1_0))))\n",
" (v_0_3 (mem (lambda (i_1_0) (if (= 1 i_1_0) 1 (letrec (\n",
" (i_1_11 (v_0_2 (- i_1_0 1)))\n",
" (i_1_12 (v_0_3 (- i_1_0 1)))\n",
" (v_1_0 (* i_1_12 i_1_11)))\n",
" v_1_0)))))\n",
" (v_0_4 (mem (lambda (i_1_1) (letrec (\n",
" (v_1_0 (* (x_church_0_0 i_1_1) (v_0_3 i_1_1))))\n",
" v_1_0))))\n",
" \n",
" (v_0_6 (mem (lambda (i_1_1) (letrec (\n",
" )\n",
" (x_church_0_1 (x_church_0_2 i_1_1))))))\n",
" \n",
" (x_church_0_0 (mem (lambda (i_1_1) (letrec (\n",
" )\n",
" (beta 1 5)))))\n",
" (x_church_0_1 (mem (lambda (i_1_1) (letrec (\n",
" )\n",
" (uniform 0 100)))))\n",
" (x_church_0_2 (mem (lambda (i_1_1) (letrec (\n",
" )\n",
" (pmf v_0_4)))))\n",
" (x_church_0_3 (mem (lambda (i_1_1) (letrec (\n",
" )\n",
" (gaussian (v_0_6 i_1_1) 1))))))\n",
" x_church_0_3))\n",
"(map (model) (iota 10000 1))\n",
"Now using node v0.10.26 (npm v1.4.3)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import Language.Church\n",
"values <- simChurchVec (dpmm 10000)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
":opt svg\n",
"import Language.Stochaskell.Plot"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEsCAIAAADfNCTgAAAABmJLR0QA/wD/AP+gvaeTAAAXIElEQVR4nO3dfXRU9Z3H8e8ND6GhuEeoSTSBhCIEhIiW0q1LQ4MuMRzrI5Wmy1bbuue4y/qwh6ZWxaU39tGm2oM9Hu1RofWwFqzoqq20jUasVmxtWSFBoGgxNJEyomhpgEyYufsHbRiSSTJzfzP36ft+HY/H/CZO7ifzy2fu3PnduZbjOAIAcKvA7w0AgHCjRgHACDUKAEaoUQAwQo0CgBFqFACMUKMAYIQaBQAj1CgAGKFGAcAINQoARoap0WTHY8svrK4sLyurOOfTzS8ddETEOdBq11dXzaiqrlvREst0BACiabi90YKSC+yNO/7U1fGyfcr91353yzE5sqlpWcv8h7fuaHt00ebrV7Z0ZzQCABE1cuibCybOu2iiiIiU1tRMe+dXf0729m5sLV781MwxUjCtYUnFJ5/ecviU4Ud662pGeZAGADyX6bHR3p0/fGh33WXnjUru64yVlJcWiIg1vqz0UGfX0QxG4vlLAAC+GmZv9Djn3dabr15bddfPLxlvxSX1A0odR8TKaGQolnXS7Y2Njc3NzX1fxmKxWCzW92VxcXFxcXFIb43FYsXFxUHbqpzfOnbs2MmTJwdtq7jV9a179uzp7u529/+G69bjf6GSLWdYf/39dy6YdfkDu3ocx3Gc+OYvn1Vz5+sJx3GSXfcunHrD890ZjMSH/AkZbUYktLW1+b0JXlAS01GTVElMx23S4V7Ux3evvurqjReseeiaaaNFRGTU7Pra/RvWbz8q8V3rHumoXTSnKIMRDowCiKphXtTHX7jrtif/GP/NZdPvEZHCC7+/9cHLF9j3nL+0YeaansKKK+9eWzdWCjIYAYCIspwAXIvJsgKxGR5ob2+fNWuW31uRd0piipqkSmKK26ScxeQpN0evQ0hJTFGTVElMcZs0ELuBevZGAUQPe6MAYIQaBQAj1CgAGKFGPZV64kSEKYkpapIqiSluk1KjnlIyHZXEFDVJlcQUahQAfEGNAoARahQAjFCjnlJyNoiSmKImqZKYwllMAOAL9kYBwAg1CgBGqFEAMEKNekrJMmYlMUVNUiUxheX3oaBkOiqJKWqSKokp1CgA+IIaBQAj1CgAGAlKjVqWZdu231uRd0rOBlESU9QkVRJTOIsJAHwRlL1RAAgpahQAjFCjAGCEGvWUkmXMSmKKmqRKYgrL70NByXRUElPUJFUSU6hRAPAFNQoARqhRADBCjXpKydkgSmKKmqRKYgpnMQFpWdZJXzLRkHMj/d4AIO/6Pq1Bwcc2wAe8qAcAI9QoABihRj2lZBmzkpiiJqmSmMLy+1BQMh2VxBQ1SZXEFGoUAHxBjQKAEWoUAIxQo55ScjaIkpiiJqmSmOI2KTXqKSXTUUlMUZNUSUyhRgHAF9QoABihRgHASFBq1LIsW8HnRihZxqwkpqhJqiSmuE0alE94UvJBebFYTMPReiUxRU1SJTHFbdKg7I0CQEhRowBghBoFACPUqKeUHGBSElPUJFUSU1h+HwpKpqOSmKImqZKYQo0CgC+oUQAwkkGNJnau/nzNjJKiwnNvb0uIiEi8dVn5uA+Vl5eXl086+0vPxUXEOdBq11dXzaiqrlvREnPSjgBABGVQo9aEudfc9cTjX/lI6lL9okvvf6Ozs7Nz77Y7F4wWObKpaVnL/Ie37mh7dNHm61e2dKcZgZqzQZTEFDVJlcSUPF5EpOC06pq5U8YXWoN/S+/Wja3FixtmjpHR0xqWVLQ+veXwgJFeF1sXOUqmo5KYoiapkpji9bWYjmy8YdbkyumfuGrV5oOOJPd1xkrKSwtExBpfVnqos+vogJG4y58EAIHm6pz6UXNWvPj6qopTEx1PLr946S1nta+S1FPiHUfESjMyFMs6cXtjY2Nzc3Pfl7FYLPUpori4OHVRArcG8NZ+/N6qWXndqvb29mA+Cjm89XjMoG1VXvNmx8nMsR3fOu+jTduO9RvufeWW6nnf+cORzV8+q+bO1xOO4yS77l049YbnuweMxAe/88w3I+za2tr83gQvBCqmiGPbf/sn5xMtUEnzR0lMx21SNy/qnYN7dsfiIpJ4+4X1v3jvzOmlY2bX1+7fsH77UYnvWvdIR+2iOUUDRka5r/roMHrGCw8lMUVNUiUxxW3SDF7UJ7tWL6lZ+eK7Bw46CyvWzGt6Yd0//qxxyR2vvJ90Rkz42BcfuO+icVKwwL7n/KUNM9f0FFZceffaurFpRqBmOiqJKWqSKokpbpNaTgA+6NOyArEZiCTLkr4PBLdtYaIh5ziLCQCMUKMAYIQa9dQQq4KiRElMUZNUSUzxevk9XFEyHZXEFDVJlcQUahQAfEGNAoARahQAjFCjnlKyjFlJTFGTVElMCftFRCzLsvsWSUeXkumoJKaoSaokpuTxZFBPcBYTgJAKyt4oAIQUNQoARqhRTylZxqwkpqhJqiSmsPw+FJRMRyUxRU1SJTGFGgUAX1CjAGCEGgUAI9Sop5QsY1YSU9QkVRJTwn4WkxJKpqOSmKImqZKYQo0CgC+oUQAwQo0CgBFq1FNKljEriSlqkiqJKSy/DwUl01FJTFGTVElMoUYBwBfUKAAYoUYBwAg16ikly5iVxBQ1SZXElLBfREQJJdMxyDEt68R/m1+5JshJc0hJTKFGgWGlXjVRwRUU4RFe1AOAEWoUAIxQo55SsoxZSUxRk1RJTAn78nvLsmwFB6uUTEclMUVNUiUxxW3SoLzF5Ji/aQoAfgjK3igAhBQ1CgBGqFFPKVnGrCSmqEmqJKZwEZFQUDIdlcQUNUmVxBRqFAB8QY0CgBFqFACMUKOeUrKMWUlMUZNUSUwJ+1lMSiiZjkpiipqkSmIKNQoAvqBGAcAINQoARqhRTylZxqwkpqhJqiSmsPw+FJRMRyUxRU1SJTGFGgUAX1CjAGCEGgUAI9Sop5QsY1YSU9QkVRJTWH4fCkqmo5KYoiapkpiSxxpN7Fz9+ZoZJUWF597eljg+5Bxoteurq2ZUVdetaIk5GY4AQARlUKPWhLnX3PXE41/5SN/l745salrWMv/hrTvaHl20+fqVLd0ZjQBAFGVQowWnVdfMnTK+0Pr7QO/Wja3FixtmjpHR0xqWVLQ+veVwBiO9+YwBAH5xc2w0ua8zVlJeWiAi1viy0kOdXUczGInndsPDSckyZiUxRU1SJTHFbVJ316lPvai844hYGY0MxbJO3N7Y2Njc3Nz3ZSwWSz3uW1xcnBo1XLce/3fQtioft0oKv7dqlgyivb3d/OfGYrHAPgq5ulVy9LsK/q39UmfIcpyM3v1J7Px2zefiP3h5ZfUI6X35pnNuKn1y0/IpBc5b911Yu+O2Vz/707nDjWxfNX/UYBthZboZQLYsS2w7zbhtC5MOOeHmRf2o2fW1+zes335U4rvWPdJRu2hOUQYjg3UoAIRaBi/qk12rl9SsfPHdAwedhRVr5jW98JNrFtj3nL+0YeaansKKK+9eWzdWCjIYAYAoCsSraT0v6o8fR/N7K/IuUDHz+qI+UEnzR0lMcZuUs5g8lXpUO8KUxBQ1SZXEFE4GBQBfUKMAYIQaBQAj1KinlBynVxJT1CRVElO4iEgoKJmOSmKKmqRKYgo1CgC+oEYBwIi7jyYBAs0a+oNwgJyiRj2l5GyQIMRMe+ZSzgUhqQeUxBTOYgoFJWeDKIkpapIqiSmcxQQAvqBGAcAINQoARqhRTyk5Tq8kpqhJqiSmsPw+FJRMRyUxRU1SJTGFGgUAX1CjAGCEGgUAI9Sop5QsY1YSU9QkVRJTwr783rIs25vT93ylZDoqiSlqkiqJKW6TBuWceiVXBgUQPUHZGwWAkKJGAcAINeopJcuYlcQUNUmVxBSW34eCkumoJKaoSaokplCjAOALahQAjFCjAGCEGvWUkmXMSmKKmqRKYkrYz2JSQsl0VBJT1CRVElOoUQDwBTUKAEaoUQAwQo16SskyZiUxRU1SJTGF5fehoGQ6KokpapIqiSnUKAD4ghoFACPUKAAYoUY9pWQZs5KYoiapkpjC8vtQUDIdlcQUNUmVxBRqFAB8QY0CgBFqFACMUKOeUrKMWUlMUZNUSUwJ+/J7y7Js2/Z7K/JOyXRUElPUJFUSU9wmHZnz7XDHcRy/NwEA3AjK3igAhBQ1CgBGqFFPKVnGrCSmqEmqJKaw/D4UlExHJTFFTVIlMYUaBQBfUKMAYIQaBQAj2ddovHVZ+bgPlZeXl5dPOvtLz8VFxDnQatdXV82oqq5b0RJz0o5ARM0yZiUxRU1SJTHF07OYii69/43Ozs7OvdvuXDBa5MimpmUt8x/euqPt0UWbr1/Z0p1mBCKiZjoqiSlqkiqJKT6eDNq7dWNr8eKGmWNk9LSGJRWtT285PGCk1/zHAEAguarRIxtvmDW5cvonrlq1+aAjyX2dsZLy0gIRscaXlR7q7Do6YCSe260GgMDI/pz6UXNWvPj6qopTEx1PLr946S1nta+S1PPhHUfESjMyDMs68S2NjY3Nzc19X8ZisdTFXMXFxak73tzKrYPdOqz29vagbTO3BuHWbFkGnwly7He3fuS/Tt3wTM39c24qfXLT8ikFzlv3XVi747ZXP/vTuSePbF81f9QQG2GZbEaYxGIxk0crLHyPaVky7OeF2baYTzrfk3pDSUxxmzTrF/XOwT27Y3ERSbz9wvpfvHfm9NIxs+tr929Yv/2oxHete6SjdtGcogEjQ3SoKqlPgBGmJKaoSaokprhNmvWL+mTXzxqX3PHK+0lnxISPffGB+y4aJwUL7HvOX9owc01PYcWVd6+tG5tmBAAiKusaHTHruideu+7kMeu0f/7aL7d/bcgRAD6zTn6PQseBNC8E5WObAXig75CxgmtNeIeTQT2l5Di9kpiiJqmSmBL2i4gooWQ6KokpYUhqDbvYMAPBj5kr1CiANHj9nm+8qAcAI9QoABihRj3l2TJmyzrxj/dYrR0xSmIKFxEJBS+no237dlCMv7qIURJTwl6jlmXZHAkHEEJBeadeyUeT5JUvr98BBKVGkRPs0Gel74knmE/iQzwvBnOD1aJGPaVkGXMoYubktMh8J027bd4/WYbiAc0Jd0mDcmxUCSXTUUlMUZNUSUyhRgHAF9QoABihRgHACDXqKSXLmJXEFDVJlcSUsC+/V0LJdFQSU9QkVRJTqFEA8AU1CgBGqFEAMEKNekrJMmYlMUVNUiUxheX3oaBkOiqJKWqSKokp1CgA+IIaBQAj1CgAGKFGPaVkGbOSmKImqZKYwvL7UFAyHZXEFDVJlcQUahQAfEGNAoARahQAjHAtJk/1Le7td7WyiF2hjNXaEaMkprhNSo16KvVByskl1YKJv7qIURJTwn4Wk2VZdvTqBIACQdkbdSL2shaAGkGpUQAeSz1Az26MCWrUU7FYTMNhJl9i9nvXzhs5T5phir5vc11/qYfQhj2cpmTeituk1KinBnuQIrZf4NdfnfdH1/ORdNgUw745mfNnFGp0aNSo/7LaL+jHl10wBB/v13qJGg29Yf9gIrarCwQNNRpxJru6ADIRlHWjSig5wKQkpqhJqiSmhH35vRJKpqOSmKImqZKYwsmgocA7QkD0UKNe4wAlEDG8qAcAI9Qoco9rTmTFsk7845e025CTDUu9k+Af1HL3gPKiHrnHSS/Z8vdQT79Vcalld/ymfoPZrj4O0WdCchYTAFNpm47Vx0PjRT0AGGFvFLlXXT0r9UvOQB0o+EcJ8ySSpyZTo4ETjXkWosNhJkwOjEb7N5NW8A8OsPw+CoI/z1wI/hOD6y1U8k5atsK7r02NhoItYvu9DW4MeynTk7/B7osZ/CcGky20bVvDNcQ2bbJra+2s/peQ/lbcPaB5eovJOdBq11dXzaiqrlvREsv1DkiIlqEN0GR+FybxTf5f2/7bPxl8w6Axw/zYpdHUlIMH1ES/VZl5+sVu2uRzzFR5DevuAc3P3uiRTU3LWub/ZOutVR3fq790Zcvv76sbm9ufkJPlbN7Ix7R291Sf7W5Xzrc853um0ehi1/za4/P+j67fA5165D0IDZCXvdHerRtbixc3zBwjo6c1LKlofXpLr6v7GXzv2v77N6TuAdlpn6MGu5Osxg3vZOj9uE2bBrlh0Jf/6cfT3k+2dz7Ylvd7Tff337Od1U5Bv41JebzSPHbZ/m4He9U52G9gkPGMfuiw0n7/YL+uLLcwuwc6V+N9+v3R9Y2f/AhmdOd9j/4gP8g++cs0D/RgGzPYnQw77o6Vjysb9zz+L+UPX9bxkyVFIj1PXV35o0WvP9owxO6oZaXfjNTxk2eeZdv9v9+2Twza9onnpUzufNhxkzuxrNQnz2G2vN+4yInx1ESDfX9Wdz7Ylqf+UNdb7no8ZXoP+uhnvoXZjqf+Wvr90GHHTxodZBbl/NeV2ztxcefmj8UQ8zzbOx94Pzn58x+ekwdHH2uY8On13cf/+8nPlSz+cfeQ35/1RgNAfrhovLwcGy04fWLJ/s59SZlS4LzT9edTys4YNdx252MzAMADeTk2Omp2fe3+Deu3H5X4rnWPdNQumjN0jQJAeOXl2KiI8/YzK5fe+PCensKKK7+39hsXlnLuPoCIylONAoAW7CUCgBFqFACMUKMAYIQaBQAj1CgAGKFGAcAINQoARvyt0fx+LKnvkh2PLb+wurK8rKzinE83v3TQkShH7n2t+RMfHPPJu/cmJcIxk7Hnvn757MmTJk088/w7thyLZlLn0CurPvPRaVOnTz1z9iVff/4dJ4IxEztXf75mRklR4bm3tyWODw3MmHFqF+fh58zh1uuq/ukb2444PTvvWlB17S/+6ufG5EFi74s/felPh5NO71uPf2HqrFt/1xvdyMfeuPeKRZ/51NQFqzoS0X1kE3sfvHhK3Z2vHko6ib/s3ftuMpJJE533Lqz81w1vJ5zkO0994cPzmncnohczEdv2q9/u+rX98Y82bTvmOE66SZtxaj/3RnP1saSBVTBx3kXnlX/AkpGlNTXT3nnrz8moRk6+9ePb/nf2V6+dNkJEovvIJjs2/Oi1+pv/c/YHLSkYN3HiqVYkk1pWQUHySPfRhCR6uns+cPoZpxyLXsyC06pr5k4ZX9j3qYcDH8rDGaf2s0aT+zpjJeWlBSJijS8rPdTZFfdxa/Kpd+cPH9pdd9l5o6IZ2Tnw1H+vLr1l+bmFfxuIZkyRY6/v2DPhrz+7qubsqhnzrvr+b993IpnUOuNzd66Qr84sKz+j6qajX/72Z4qdKMbsZ+BDeTTj1D4fG3VO/u9IXhLCebf15qvXVt11xyXjrUhG/kvr7d+TG1fM/+AgV3OISEwRkWSyZ9u2guueeHX7c7eOufcL33qlN4pJnYMv3P/I2Nu3dnZ2bf9OUfN/3PfHZBRjDjAwY6ap/azRvo8lFcnoY0lDqXvLd5fcuOfa9XdfVGJFM/KxN1/+3R9+fuOcyZVTFz/wxuamT17xg70l0YspIjKibNLEyguu+Pj4gpGlCy/7+Hvtr3WXRjBp/MX/WTeu7tLK0TJq4iWXnPV/v96aiOC87W/g32Zhxqn9rNHofyxpfPfqq67eeMGah66ZNlpEohl55NkrXnqr680333xz94Z/m3LeV59/7Nqp50YvpojIiKr6ug+8/Oxrh8V5/zfP/n5c1bTx50Qw6YjTyz607Zln9yckGWv95dayMyePieC87W/g32ZR5qk9e2csnWSs5baFZ334zCkzLrj55/sSvm5LHvQ88++njyyaUHbch7/42KFoR47/evn04+/URzZm8r2Xmi+bPWVyZeWsi7/+/IFkNJMmY882fWrmpEmVkyqqL/3mr96JYsxE54OLJ5eV/EPh6FNKyiqveOBPiTQZM03N540CgBHOYgIAI9QoABihRgHACDUKAEaoUQAwQo0CgBFqFACMUKMAYIQaBQAj1CgAGKFGAcDI/wMiCXybp8RuWQAAAABJRU5ErkJggg==",
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"toRenderable . plot . return . histToPlot $ defaultPlotHist\n",
" { _plot_hist_values = values, _plot_hist_bins = 100, _plot_hist_range = Just (0,100) }"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"import IHaskell.Display"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"svg <$> vizIR (dpmm 10000)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Haskell",
"language": "haskell",
"name": "haskell"
},
"language_info": {
"codemirror_mode": "ihaskell",
"file_extension": ".hs",
"name": "haskell",
"pygments_lexer": "Haskell",
"version": "8.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}