{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using Distributions, Plots\n",
"default(leg=:none)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate some Poisson distributed test data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n=60 # the number of samples\n",
"λ=6.0\n",
"\n",
"# Generate a data set for testing\n",
"dist = Poisson(λ)\n",
"x = rand(dist, n)\n",
"\n",
"time = plot(x)\n",
"hist = histogram(x, normed=true)\n",
"plot(time, hist)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can calculate the posterior distribution of $\\lambda$ and the predictive posterior distribution of the next sample based on the data seen up to time 't'."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = [sum(x[1:i]) for i in 1:n]\n",
"\n",
"alpha(t) = X[Int(t)]\n",
"theta(t) = 1/(t+1)\n",
"\n",
"# Posterior distribution of λ\n",
"posterior(t) = Gamma(alpha(t), theta(t))\n",
"\n",
"# Predictive posterior of next sample\n",
"predictive(t) = NegativeBinomial(alpha(t), 1-theta(t))\n",
"\n",
"T=Int(n*0.1) # observed % of test data\n",
"plot(0:1:15, [x -> Distributions.pdf(posterior(T), x), x -> Distributions.pdf(predictive(T), x)], seriestype=:bar, layout=2, leg=:none)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See how the posterior distribution behaves over time and settles close to the known value."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"λes = [mean(posterior(i)) for i in 1:n]\n",
"λlower = [quantile(posterior(i), 0.05)-λes[i] for i in 1:n]\n",
"λupper = [λes[i]-quantile(posterior(i), 0.95) for i in 1:n]\n",
"plot(λes, ribbon=(λlower, λupper), leg=:none)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Construct a whole ensemble of trajectories to use as a background guide."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"ensemble = []\n",
"for j in 1:1000\n",
"samples = Array{Real}(undef, n)\n",
"for i in 1:n\n",
" samples[i] = 12*rand() + ((i>1) ? samples[i-1] : 0)\n",
"end\n",
"push!(ensemble, samples)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to throttle the total to a specific value by only allowing a fraction of the samples through. To do that we can calculate the overshoot factor with the data up to the current time, t, and aapply it as a correcting factor to the samples."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"T = 300 # Desired total samples allowed through\n",
"y = []\n",
"for t in 1:n\n",
" λest = quantile(posterior(t), 0.9)\n",
" S = sum(x[1:t]) + (n-t)*λest # Estimated total\n",
" factor = minimum([1.0, T/S])\n",
" push!(y, x[t]*factor)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Evolve the posterior distribution over time and see how the it narrows to the actual values. Display with the ensemble and the trottles samples."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to \n",
"│ fn = /Users/john/Projects/personal/notebooks/tmp.gif\n",
"└ @ Plots /Users/john/.julia/packages/Plots/WwFyB/src/animation.jl:98\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
"Plots.AnimatedGif(\"/Users/john/Projects/personal/notebooks/tmp.gif\")"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Xpred = Array{Real}(undef, n)\n",
"lower = Array{Real}(undef, n)\n",
"upper = Array{Real}(undef, n)\n",
"@gif for t in 5:n \n",
" for i in 1:n\n",
" Xpred[i] = i