{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using Turing, Distributions, Plots, Random"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n=100 # the number of observations\n",
"\n",
"# Generate a data set for testing\n",
"dist = Poisson(3)\n",
"x = [rand(dist) for i in 1:n]\n",
"\n",
"time = plot(1:n, x)\n",
"hist = histogram(x, normed=true)\n",
"plot(time, hist)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Analytic Solution"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Gamma{Float64}(α=3.0, θ=1.0)"
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Our prior belief about the probability of heads in a coin toss.\n",
"prior_belief = Gamma(3, 1)"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Import StatsPlots for animating purposes.\n",
"using StatsPlots\n",
"\n",
"# Animate posterior distribution as number of observations increases.\n",
"default(display=:inline)\n",
"for i in 1:n\n",
"\n",
" # Count the number of heads and tails.\n",
" count = sum(x[1:i-1])\n",
" \n",
" # Update our prior belief in closed form (this is possible because we use a conjugate prior).\n",
" updated_belief = Gamma(prior_belief.α + count, 1/((1/prior_belief.θ) + i))\n",
"\n",
" # Plotting\n",
" plot(updated_belief, \n",
" size = (500, 250), \n",
" title = \"Updated belief after $i observations\",\n",
" xlabel = \"λ\", \n",
" ylabel = \"p(λ)\", \n",
" legend = nothing,\n",
" xlim = (0,10),\n",
" fill=0, α=0.3, w=3)\n",
" vline!([3.0])\n",
"end\n",
"default(display=:none)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Turing Solution"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"poisson_model (generic function with 3 methods)"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@model poisson_model(x, n) = begin\n",
" \n",
" λ ~ Gamma(1,2)\n",
" \n",
" for i in 1:n\n",
" x[i] ~ Poisson(λ)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Found initial step size\n",
"│ ϵ = 0.2\n",
"└ @ Turing.Inference /Users/john/.julia/packages/Turing/xGrY5/src/inference/hmc.jl:547\n",
"\u001b[32mSampling 100%|███████████████████████████████| Time: 0:00:00\u001b[39m\n"
]
}
],
"source": [
"model = poisson_model(x, n)\n",
"chain = sample(model, NUTS(0.65), 1000);"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{ChainDataFrame,1}\n",
"\n",
"Summary Statistics\n",
". Omitted printing of 1 columns\n",
"│ Row │ parameters │ mean │ std │ naive_se │ mcse │ ess │\n",
"│ │ \u001b[90mSymbol\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mAny\u001b[39m │\n",
"├─────┼────────────┼─────────┼──────────┼────────────┼────────────┼─────────┤\n",
"│ 1 │ λ │ 3.23552 │ 0.172132 │ 0.00769796 │ 0.00206103 │ 362.874 │\n",
"\n",
"Quantiles\n",
"\n",
"│ Row │ parameters │ 2.5% │ 25.0% │ 50.0% │ 75.0% │ 97.5% │\n",
"│ │ \u001b[90mSymbol\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │\n",
"├─────┼────────────┼─────────┼─────────┼─────────┼─────────┼─────────┤\n",
"│ 1 │ λ │ 2.88478 │ 3.11605 │ 3.24064 │ 3.34866 │ 3.58392 │\n"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"describe(chain)"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(chain[:λ], seriestype = :histogram)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.2.0",
"language": "julia",
"name": "julia-1.2"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.2.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}