From a0aa2205700c8422e3e24350bd440ddf6d61922c Mon Sep 17 00:00:00 2001 From: Max Lapan Date: Sun, 23 Jan 2022 14:30:32 +0100 Subject: [PATCH] Ch11 done --- 11-Chapter 11. God Spiked the Integers.ipynb | 3005 +++++++++------- ...1-Chapter 11. God Spiked the Integers.html | 3093 ++++++++++------- 2 files changed, 3637 insertions(+), 2461 deletions(-) diff --git a/11-Chapter 11. God Spiked the Integers.ipynb b/11-Chapter 11. God Spiked the Integers.ipynb index 021d645..38c8e34 100644 --- a/11-Chapter 11. God Spiked the Integers.ipynb +++ b/11-Chapter 11. God Spiked the Integers.ipynb @@ -177,105 +177,105 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", "\n" ] }, @@ -798,104 +798,104 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", "\n" ] }, @@ -930,86 +930,86 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", "\n" ] }, @@ -1094,151 +1094,151 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n" + "\n", + "\n", + "\n" ] }, "execution_count": 17, @@ -1302,291 +1302,291 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" + "\n", + "\n", + "\n", + "\n", + "\n" ] }, "execution_count": 19, @@ -2104,291 +2104,291 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n", + "\n", + "\n", - "\n", - "\n", - "\n" + "\n", + "\n", + "\n" ] }, "execution_count": 31, @@ -2699,98 +2699,98 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" + "\n" ] }, "execution_count": 38, @@ -2923,99 +2923,99 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" ] }, - "execution_count": 48, + "execution_count": 46, "metadata": {}, "output_type": "execute_result" } @@ -8622,7 +8622,7 @@ }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 47, "id": "c7044997", "metadata": {}, "outputs": [ @@ -8632,97 +8632,97 @@ "\n", "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", "\n", - " \n", + " \n", " \n", " \n", "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n" ] }, - "execution_count": 49, + "execution_count": 47, "metadata": {}, "output_type": "execute_result" } @@ -8887,10 +8887,553 @@ "plot!(pop_seq, [λ2_mean λ2_mean], fillrange=λ2_ci, fillalpha=0.2, c=:black, lw=1.5)" ] }, + { + "cell_type": "markdown", + "id": "d45260e2", + "metadata": {}, + "source": [ + "Code 11.49" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "id": "9dc407b8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "┌───────┬───────────────────────────────────────────────────────┐\n", + "│\u001b[1m param \u001b[0m│\u001b[1m mean \u001b[0m\u001b[1m std \u001b[0m\u001b[1m 5.5% \u001b[0m\u001b[1m 50% \u001b[0m\u001b[1m 94.5% \u001b[0m\u001b[1m histogram \u001b[0m│\n", + "├───────┼───────────────────────────────────────────────────────┤\n", + "│ a[1] │ 0.8656 0.6512 -0.2068 0.864 1.8687 ▁▁▁▁▂▅█▇▄▁▁ │\n", + "│ a[2] │ 0.8979 0.8581 -0.473 0.8861 2.3729 ▁▂▄▇██▅▃▂▁ │\n", + "│ b₁ │ 0.2588 0.0341 0.2052 0.2591 0.313 ▁▁▁▁▃▇██▆▄▂▁ │\n", + "│ b₂ │ 0.2869 0.1039 0.1273 0.2828 0.4552 ▁▁▄▆▇██▆▄▂▂▁ │\n", + "│ g │ 1.0763 0.7044 0.3216 0.8842 2.5334 ▅█▅▃▂▁▁▁▁ │\n", + "└───────┴───────────────────────────────────────────────────────┘\n" + ] + } + ], + "source": [ + "@model function m11_11(T, P, cid)\n", + " a ~ MvNormal([1,1], 1)\n", + " b₁ ~ Exponential(1)\n", + " b₂ ~ Exponential(1)\n", + " b = [b₁, b₂]\n", + " g ~ Exponential(1)\n", + " λ = @. exp(a[cid]) * P ^ b[cid] / g\n", + " for i ∈ eachindex(T)\n", + " T[i] ~ Poisson(λ[i])\n", + " end\n", + "end\n", + "\n", + "m11_11_ch = sample(m11_11(d.total_tools, d.population, d.contact_id), NUTS(), 1000)\n", + "m11_11_df = DataFrame(m11_11_ch)\n", + "precis(m11_11_df)" + ] + }, + { + "cell_type": "markdown", + "id": "576bed01", + "metadata": {}, + "source": [ + "Code 11.50" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "id": "693c7e68", + "metadata": {}, + "outputs": [], + "source": [ + "Random.seed!(1)\n", + "num_days = 30\n", + "y = rand(Poisson(1.5), num_days);" + ] + }, + { + "cell_type": "markdown", + "id": "971d5d97", + "metadata": {}, + "source": [ + "Code 11.51" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "id": "da709fbb", + "metadata": {}, + "outputs": [], + "source": [ + "Random.seed!(2)\n", + "num_weeks = 4\n", + "y_new = rand(Poisson(0.5*7), num_weeks);" + ] + }, + { + "cell_type": "markdown", + "id": "4007a47f", + "metadata": {}, + "source": [ + "Code 11.52" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "id": "72fdfec5", + "metadata": {}, + "outputs": [], + "source": [ + "y_all = [y; y_new]\n", + "exposure = [repeat([1], 30); repeat([7], 4)]\n", + "monastery = [repeat([0], 30); repeat([1], 4)]\n", + "d = DataFrame(y=y_all, days=exposure, monastery=monastery);" + ] + }, + { + "cell_type": "markdown", + "id": "db84e175", + "metadata": {}, + "source": [ + "Code 11.53" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "id": "97b313f3", + "metadata": {}, + "outputs": [], + "source": [ + "Random.seed!(1)\n", + "d.log_days = log.(d.days)\n", + "\n", + "@model function m11_12(y, log_days, monastery)\n", + " a ~ Normal()\n", + " b ~ Normal()\n", + " λ = @. exp(log_days + a + b*monastery)\n", + " @. y ~ Poisson(λ)\n", + "end\n", + "\n", + "m11_12_chain = sample(m11_12(d.y, d.log_days, d.monastery), NUTS(), 1000)\n", + "m11_12_df = DataFrame(m11_12_chain);" + ] + }, + { + "cell_type": "markdown", + "id": "77430917", + "metadata": {}, + "source": [ + "Code 11.54\n", + "\n", + "Note that values are slightly different from the book. This is due to non-Normal distributions (you can check this yourself with `optimize`)" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "id": "4e35c982", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "┌───────┬──────────────────────────────────────────────────────┐\n", + "│\u001b[1m param \u001b[0m│\u001b[1m mean \u001b[0m\u001b[1m std \u001b[0m\u001b[1m 5.5% \u001b[0m\u001b[1m 50% \u001b[0m\u001b[1m 94.5% \u001b[0m\u001b[1m histogram \u001b[0m│\n", + "├───────┼──────────────────────────────────────────────────────┤\n", + "│ λ_old │ 1.4935 0.2284 1.1472 1.4783 1.8972 ▁▃▇█▅▂▁▁ │\n", + "│ λ_new │ 0.7648 0.1635 0.513 0.76 1.0218 ▁▁▂▄▇█▇▄▂▁▁▁ │\n", + "└───────┴──────────────────────────────────────────────────────┘\n" + ] + } + ], + "source": [ + "λ_old = exp.(m11_12_df.a)\n", + "λ_new = @. exp(m11_12_df.a + m11_12_df.b)\n", + "precis(DataFrame(λ_old=λ_old, λ_new=λ_new))" + ] + }, + { + "cell_type": "markdown", + "id": "674f5929", + "metadata": {}, + "source": [ + "# 11.3 Multinomial and categorical models" + ] + }, + { + "cell_type": "markdown", + "id": "08e24118", + "metadata": {}, + "source": [ + "Code 11.55" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "id": "3140ae27", + "metadata": {}, + "outputs": [], + "source": [ + "# simulate career choice among 500 individuals\n", + "N = 500\n", + "income = [1,2,5]\n", + "c_score = 0.5*income\n", + "p = softmax(c_score)\n", + "\n", + "# simulate choice\n", + "Random.seed!(34302)\n", + "w = Weights(p)\n", + "career = [sample(w) for _ in 1:N];" + ] + }, + { + "cell_type": "markdown", + "id": "dbbdcbdd", + "metadata": {}, + "source": [ + "Code 11.56" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "7ae66856", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "m11_13 (generic function with 2 methods)" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "@model function m11_13(career, income)\n", + " a ~ MvNormal([0, 0], 1)\n", + " b ~ TruncatedNormal(0, 0.5, 0, Inf)\n", + " p = softmax([a[1] + b*income[1], a[2] + b*income[2], 0])\n", + " career ~ Categorical(p)\n", + "end" + ] + }, + { + "cell_type": "markdown", + "id": "73ea3e9a", + "metadata": {}, + "source": [ + "Code 11.57" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "id": "2f75f21b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "┌───────┬─────────────────────────────────────────────────────────┐\n", + "│\u001b[1m param \u001b[0m│\u001b[1m mean \u001b[0m\u001b[1m std \u001b[0m\u001b[1m 5.5% \u001b[0m\u001b[1m 50% \u001b[0m\u001b[1m 94.5% \u001b[0m\u001b[1m histogram \u001b[0m│\n", + "├───────┼─────────────────────────────────────────────────────────┤\n", + "│ a[1] │ -1.8621 0.1702 -2.1531 -1.8514 -1.6067 ▁▁▂▅▇█▇▅▂▁▁ │\n", + "│ a[2] │ -1.7931 0.2471 -2.2586 -1.7554 -1.457 ▁▂▃▆█▆▁▁ │\n", + "│ b │ 0.1391 0.1134 0.0096 0.1104 0.3681 █▇▆▄▃▂▂▁▁▁▁ │\n", + "└───────┴─────────────────────────────────────────────────────────┘\n" + ] + } + ], + "source": [ + "Random.seed!(121)\n", + "m11_13_chain = sample(m11_13(career, income), NUTS(), 5000, n_chains=4)\n", + "m11_13_df = DataFrame(m11_13_chain)\n", + "precis(m11_13_df)" + ] + }, + { + "cell_type": "markdown", + "id": "665ffc5b", + "metadata": {}, + "source": [ + "Code 11.58" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "e591b764", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "┌────────┬─────────────────────────────────────────────────────┐\n", + "│\u001b[1m param \u001b[0m│\u001b[1m mean \u001b[0m\u001b[1m std \u001b[0m\u001b[1m 5.5% \u001b[0m\u001b[1m 50% \u001b[0m\u001b[1m 94.5% \u001b[0m\u001b[1m histogram \u001b[0m│\n", + "├────────┼─────────────────────────────────────────────────────┤\n", + "│ p_diff │ 0.0432 0.0394 0.0026 0.0315 0.1278 █▆▄▃▂▁▁▁▁▁▁ │\n", + "└────────┴─────────────────────────────────────────────────────┘\n" + ] + } + ], + "source": [ + "# logit scores\n", + "s1 = m11_13_df.\"a[1]\" + m11_13_df.b * income[1]\n", + "s2_orig = m11_13_df.\"a[2]\" + m11_13_df.b * income[2]\n", + "s2_new = m11_13_df.\"a[2]\" + m11_13_df.b * income[2] * 2\n", + "\n", + "# compute probabilities for original and counterfactual\n", + "p_orig = softmax.(eachrow(hcat(s1, s2_orig, zeros(length(s1)))))\n", + "p_orig = hcat(p_orig...)'\n", + "p_new = softmax.(eachrow(hcat(s1, s2_new, zeros(length(s1)))))\n", + "p_new = hcat(p_new...)'\n", + "\n", + "p_diff = p_new[:, 2] - p_orig[:, 2]\n", + "precis(DataFrame(p_diff=p_diff))" + ] + }, + { + "cell_type": "markdown", + "id": "8a5c480a", + "metadata": {}, + "source": [ + "Code 11.59" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "id": "fbb2ca9e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "┌───────┬────────────────────────────────────────────────────────────┐\n", + "│\u001b[1m param \u001b[0m│\u001b[1m mean \u001b[0m\u001b[1m std \u001b[0m\u001b[1m 5.5% \u001b[0m\u001b[1m 50% \u001b[0m\u001b[1m 94.5% \u001b[0m\u001b[1m histogram \u001b[0m│\n", + "├───────┼────────────────────────────────────────────────────────────┤\n", + "│ a[1] │ -2.7538 0.39 -3.3927 -2.7355 -2.1699 ▁▁▁▂▃▆██▇▅▂▁▁▁ │\n", + "│ a[2] │ -1.2474 0.2198 -1.6024 -1.2415 -0.896 ▁▁▁▂▄▆███▆▄▂▁ │\n", + "│ b[1] │ -1.6594 0.7035 -2.7403 -1.658 -0.5104 ▁▁▃▇█▇▄▂▁▁ │\n", + "│ b[2] │ -1.0656 0.4021 -1.7315 -1.0556 -0.4297 ▁▁▂▃▅▇██▆▄▂▁▁▁ │\n", + "└───────┴────────────────────────────────────────────────────────────┘\n" + ] + } + ], + "source": [ + "Random.seed!(1)\n", + "\n", + "N = 500\n", + "family_income = rand(Uniform(), N)\n", + "b = [-2, 0, 2]\n", + "career = [\n", + " begin\n", + " score = (0.5 * 1:3) + b * inc\n", + " p = softmax(score)\n", + " sample(Weights(p))\n", + " end\n", + " for inc in family_income\n", + "]\n", + "\n", + "@model function m11_14(career, family_income)\n", + " a ~ MvNormal([0, 0], 1.5)\n", + " b ~ MvNormal([0, 0], 1)\n", + " \n", + " for i ∈ eachindex(career)\n", + " income = family_income[i]\n", + " s = [a[1] + b[1] * income, a[2] + b[2] * income, 0]\n", + " p = softmax(s)\n", + " career[i] ~ Categorical(p)\n", + " end\n", + "end\n", + "\n", + "m11_14_chain = sample(m11_14(career, family_income), NUTS(), 1000)\n", + "m11_14_df = DataFrame(m11_14_chain)\n", + "precis(m11_14_df)" + ] + }, + { + "cell_type": "markdown", + "id": "2875f54a", + "metadata": {}, + "source": [ + "Code 11.60" + ] + }, + { + "cell_type": "code", + "execution_count": 60, + "id": "a911a11f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "

12 rows × 6 columns

iddeptgenderadmitrejectapplications
Int64String1String7Int64Int64Int64
11Amale512313825
22Afemale8919108
33Bmale353207560
44Bfemale17825
55Cmale120205325
66Cfemale202391593
77Dmale138279417
88Dfemale131244375
99Emale53138191
1010Efemale94299393
1111Fmale22351373
1212Ffemale24317341
" + ], + "text/latex": [ + "\\begin{tabular}{r|cccccc}\n", + "\t& id & dept & gender & admit & reject & applications\\\\\n", + "\t\\hline\n", + "\t& Int64 & String1 & String7 & Int64 & Int64 & Int64\\\\\n", + "\t\\hline\n", + "\t1 & 1 & A & male & 512 & 313 & 825 \\\\\n", + "\t2 & 2 & A & female & 89 & 19 & 108 \\\\\n", + "\t3 & 3 & B & male & 353 & 207 & 560 \\\\\n", + "\t4 & 4 & B & female & 17 & 8 & 25 \\\\\n", + "\t5 & 5 & C & male & 120 & 205 & 325 \\\\\n", + "\t6 & 6 & C & female & 202 & 391 & 593 \\\\\n", + "\t7 & 7 & D & male & 138 & 279 & 417 \\\\\n", + "\t8 & 8 & D & female & 131 & 244 & 375 \\\\\n", + "\t9 & 9 & E & male & 53 & 138 & 191 \\\\\n", + "\t10 & 10 & E & female & 94 & 299 & 393 \\\\\n", + "\t11 & 11 & F & male & 22 & 351 & 373 \\\\\n", + "\t12 & 12 & F & female & 24 & 317 & 341 \\\\\n", + "\\end{tabular}\n" + ], + "text/plain": [ + "\u001b[1m12×6 DataFrame\u001b[0m\n", + "\u001b[1m Row \u001b[0m│\u001b[1m id \u001b[0m\u001b[1m dept \u001b[0m\u001b[1m gender \u001b[0m\u001b[1m admit \u001b[0m\u001b[1m reject \u001b[0m\u001b[1m applications \u001b[0m\n", + "\u001b[1m \u001b[0m│\u001b[90m Int64 \u001b[0m\u001b[90m String1 \u001b[0m\u001b[90m String7 \u001b[0m\u001b[90m Int64 \u001b[0m\u001b[90m Int64 \u001b[0m\u001b[90m Int64 \u001b[0m\n", + "─────┼──────────────────────────────────────────────────────\n", + " 1 │ 1 A male 512 313 825\n", + " 2 │ 2 A female 89 19 108\n", + " 3 │ 3 B male 353 207 560\n", + " 4 │ 4 B female 17 8 25\n", + " 5 │ 5 C male 120 205 325\n", + " 6 │ 6 C female 202 391 593\n", + " 7 │ 7 D male 138 279 417\n", + " 8 │ 8 D female 131 244 375\n", + " 9 │ 9 E male 53 138 191\n", + " 10 │ 10 E female 94 299 393\n", + " 11 │ 11 F male 22 351 373\n", + " 12 │ 12 F female 24 317 341" + ] + }, + "execution_count": 60, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "d = DataFrame(CSV.File(\"data/UCBadmit.csv\", skipto=2, \n", + " header=[:id, :dept, :gender, :admit, :reject, :applications]))" + ] + }, + { + "cell_type": "markdown", + "id": "8afe40a9", + "metadata": {}, + "source": [ + "Code 11.61" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "id": "985ce663", + "metadata": {}, + "outputs": [], + "source": [ + "@model function m_binom(admit, applications)\n", + " a ~ Normal(0, 1.5)\n", + " p = logistic(a)\n", + " @. admit ~ Binomial(applications, p)\n", + "end\n", + "\n", + "m_binom_chain = sample(m_binom(d.admit, d.applications), NUTS(), 1000)\n", + "m_binom_df = DataFrame(m_binom_chain);\n", + "\n", + "@model function m_pois(admit, reject)\n", + " a ~ MvNormal([0, 0], 1.5)\n", + " λ = exp.(a)\n", + " admit ~ Poisson(λ[1])\n", + " reject ~ Poisson(λ[2])\n", + "end\n", + "\n", + "m_pois_chain = sample(m_pois(d.admit, d.reject), NUTS(), 1000)\n", + "m_pois_df = DataFrame(m_pois_chain);" + ] + }, + { + "cell_type": "markdown", + "id": "836cbd62", + "metadata": {}, + "source": [ + "Code 11.62" + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "id": "3c2956d0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.387568484041908" + ] + }, + "execution_count": 62, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m_binom_df.a .|> logistic |> mean" + ] + }, + { + "cell_type": "markdown", + "id": "382cc15f", + "metadata": {}, + "source": [ + "Code 11.63" + ] + }, + { + "cell_type": "code", + "execution_count": 63, + "id": "3952375d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.38722562179640474" + ] + }, + "execution_count": 63, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a1 = m_pois_df.\"a[1]\" |> mean\n", + "a2 = m_pois_df.\"a[2]\" |> mean\n", + "exp(a1)/(exp(a1)+exp(a2))" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "f14a40fd", + "id": "827b4d20", "metadata": {}, "outputs": [], "source": [] diff --git a/docs/11-Chapter 11. God Spiked the Integers.html b/docs/11-Chapter 11. God Spiked the Integers.html index d31e154..badd547 100644 --- a/docs/11-Chapter 11. God Spiked the Integers.html +++ b/docs/11-Chapter 11. God Spiked the Integers.html @@ -13307,105 +13307,105 @@

11.1 Binomial regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - - -11.1 Binomial regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - -11.1 Binomial regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + @@ -13996,104 +13996,104 @@

11.1 Binomial regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + @@ -14142,86 +14142,86 @@

11.1 Binomial regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + @@ -14328,151 +14328,151 @@

11.1 Binomial regression - + - - + - - + - - - - - - - - - - - - - - - - + + - - - + + - - - - + + - - - + + - - - - + + - - - + + - - - - + + - - - + + - - - - + + - - - + + - - - - + + - - - + + - - - - + + - - - + + + @@ -14558,291 +14558,291 @@

11.1 Binomial regression - + - - + - - + - - - - - - - - - - - - - - - - + + - - - - - - + + + + - - - + + - - - - - - + + + + - - - - + + - - - - - - + + + + - - - + + - - - - - - + + + + - - - - + + - - - - - - + + + + - - - + + - - - - - - + + + + - - - - + + - - - - - - + + + + - - - + + - - - - - - + + + + - - - - + + - - - - - - + + + + - - - + + - - - - - - + + + + - - - - + + - - - - - - + + + + - - - + + - - - - - - + + + + - - - - + + - - - - - - + + + + - - - + + - - - - - - + + + + + @@ -15390,291 +15390,291 @@

11.1 Binomial regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - + + - - - + + - - - + + - - - + + - - - + + - - - + + + @@ -15986,98 +15986,98 @@

11.2 Poisson regression - + - - + - + - - + - - - - - - - - - - - - - - - - - - - - - - -11.2 Poisson regression -11.2 Poisson regression - - - - - + @@ -16240,99 +16240,99 @@

11.2 Poisson regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - - - -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression - + - - + - - + - - - - - - - - - - - - - - - - - - - - - - -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression -11.2 Poisson regression
-
In [48]:
+
In [46]:
# get psis K values
@@ -21709,7 +21709,7 @@ 
11.2 Poisson regression
 
-    
Out[48]:

+    
Out[46]:

 
 
 
@@ -21717,97 +21717,97 @@ 
11.2 Poisson regression
 
 
-  
+  
     
   
 
-
 
-  
+  
     
   
 
-
 
-  
+  
     
   
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
 

 

-
In [49]:

+
In [47]:

 

     

 
scatter(d.population, d.total_tools, xlab="population", ylab="total tools", ylim=(0, 75), 
@@ -21996,7 +21996,7 @@ 
11.2 Poisson regression
 
-    
Out[49]:

+    
Out[47]:

 
 
 
@@ -22004,97 +22004,97 @@ 
11.2 Poisson regression
 
 
-  
+  
     
   
 
-
 
-  
+  
     
   
 
-
 
-  
+  
     
   
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
-11.2 Poisson regression
 

 
+

+

+

+

+
Code 11.49

+
+

+

+

+

+

+
In [48]:

+

+    

+
@model function m11_11(T, P, cid)
+    a ~ MvNormal([1,1], 1)
+    b₁ ~ Exponential(1)
+    b₂ ~ Exponential(1)
+    b = [b₁, b₂]
+    g ~ Exponential(1)
+    λ = @. exp(a[cid]) * P ^ b[cid] / g
+    for i ∈ eachindex(T)
+        T[i] ~ Poisson(λ[i])
+    end
+end
+
+m11_11_ch = sample(m11_11(d.total_tools, d.population, d.contact_id), NUTS(), 1000)
+m11_11_df = DataFrame(m11_11_ch)
+precis(m11_11_df)
+

+
+    

+

+

+
+

+

+
+
+

+
+    

+
+
+

+
┌───────┬───────────────────────────────────────────────────────┐
+│ param │   mean     std     5.5%     50%   94.5%     histogram │
+├───────┼───────────────────────────────────────────────────────┤
+│  a[1] │ 0.8656  0.6512  -0.2068   0.864  1.8687   ▁▁▁▁▂▅█▇▄▁▁ │
+│  a[2] │ 0.8979  0.8581   -0.473  0.8861  2.3729    ▁▂▄▇██▅▃▂▁ │
+│    b₁ │ 0.2588  0.0341   0.2052  0.2591   0.313  ▁▁▁▁▃▇██▆▄▂▁ │
+│    b₂ │ 0.2869  0.1039   0.1273  0.2828  0.4552  ▁▁▄▆▇██▆▄▂▂▁ │
+│     g │ 1.0763  0.7044   0.3216  0.8842  2.5334     ▅█▅▃▂▁▁▁▁ │
+└───────┴───────────────────────────────────────────────────────┘
+

+

+

+
+

+

+
+

+

+

+

+
Code 11.50

+
+

+

+

+

+

+
In [49]:

+

+    

+
Random.seed!(1)
+num_days = 30
+y = rand(Poisson(1.5), num_days);
+

+
+    

+

+

+
+

+

+

+

+
Code 11.51

+
+

+

+

+

+

+
In [50]:

+

+    

+
Random.seed!(2)
+num_weeks = 4
+y_new = rand(Poisson(0.5*7), num_weeks);
+

+
+    

+

+

+
+

+

+

+

+
Code 11.52

+
+

+

+

+

+

+
In [51]:

+

+    

+
y_all = [y; y_new]
+exposure = [repeat([1], 30); repeat([7], 4)]
+monastery = [repeat([0], 30); repeat([1], 4)]
+d = DataFrame(y=y_all, days=exposure, monastery=monastery);
+

+
+    

+

+

+
+

+

+

+

+
Code 11.53

+
+

+

+

+

+

+
In [52]:

+

+    

+
Random.seed!(1)
+d.log_days = log.(d.days)
+
+@model function m11_12(y, log_days, monastery)
+    a ~ Normal()
+    b ~ Normal()
+    λ = @. exp(log_days + a + b*monastery)
+    @. y ~ Poisson(λ)
+end
+
+m11_12_chain = sample(m11_12(d.y, d.log_days, d.monastery), NUTS(), 1000)
+m11_12_df = DataFrame(m11_12_chain);
+

+
+    

+

+

+
+

+

+

+

+
Code 11.54

+
Note that values are slightly different from the book. This is due to non-Normal distributions (you can check this yourself with optimize)

+
+

+

+

+

+

+
In [53]:

+

+    

+
λ_old = exp.(m11_12_df.a)
+λ_new = @. exp(m11_12_df.a + m11_12_df.b)
+precis(DataFrame(λ_old=λ_old, λ_new=λ_new))
+

+
+    

+

+

+
+

+

+
+
+

+
+    

+
+
+

+
┌───────┬──────────────────────────────────────────────────────┐
+│ param │   mean     std    5.5%     50%   94.5%     histogram │
+├───────┼──────────────────────────────────────────────────────┤
+│ λ_old │ 1.4935  0.2284  1.1472  1.4783  1.8972      ▁▃▇█▅▂▁▁ │
+│ λ_new │ 0.7648  0.1635   0.513    0.76  1.0218  ▁▁▂▄▇█▇▄▂▁▁▁ │
+└───────┴──────────────────────────────────────────────────────┘
+

+

+

+
+

+

+
+

+

+

+

+
11.3 Multinomial and categorical models

+

+

+

+

+

+

+
Code 11.55

+
+

+

+

+

+

+
In [54]:

+

+    

+
# simulate career choice among 500 individuals
+N = 500
+income = [1,2,5]
+c_score = 0.5*income
+p = softmax(c_score)
+
+# simulate choice
+Random.seed!(34302)
+w = Weights(p)
+career = [sample(w) for _ in 1:N];
+

+
+    

+

+

+
+

+

+

+

+
Code 11.56

+
+

+

+

+

+

+
In [55]:

+

+    

+
@model function m11_13(career, income)
+    a ~ MvNormal([0, 0], 1)
+    b ~ TruncatedNormal(0, 0.5, 0, Inf)
+    p = softmax([a[1] + b*income[1], a[2] + b*income[2], 0])
+    career ~ Categorical(p)
+end
+

+
+    

+

+

+
+

+

+
+
+

+
+    
Out[55]:

+
+
+
+
+

+
m11_13 (generic function with 2 methods)

+

+
+

+
+

+

+
+

+

+

+

+
Code 11.57

+
+

+

+

+

+

+
In [56]:

+

+    

+
Random.seed!(121)
+m11_13_chain = sample(m11_13(career, income), NUTS(), 5000, n_chains=4)
+m11_13_df = DataFrame(m11_13_chain)
+precis(m11_13_df)
+

+
+    

+

+

+
+

+

+
+
+

+
+    

+
+
+

+
┌───────┬─────────────────────────────────────────────────────────┐
+│ param │    mean     std     5.5%      50%    94.5%    histogram │
+├───────┼─────────────────────────────────────────────────────────┤
+│  a[1] │ -1.8621  0.1702  -2.1531  -1.8514  -1.6067  ▁▁▂▅▇█▇▅▂▁▁ │
+│  a[2] │ -1.7931  0.2471  -2.2586  -1.7554   -1.457     ▁▂▃▆█▆▁▁ │
+│     b │  0.1391  0.1134   0.0096   0.1104   0.3681  █▇▆▄▃▂▂▁▁▁▁ │
+└───────┴─────────────────────────────────────────────────────────┘
+

+

+

+
+

+

+
+

+

+

+

+
Code 11.58

+
+

+

+

+

+

+
In [57]:

+

+    

+
# logit scores
+s1 = m11_13_df."a[1]" + m11_13_df.b * income[1]
+s2_orig = m11_13_df."a[2]" + m11_13_df.b * income[2]
+s2_new = m11_13_df."a[2]" + m11_13_df.b * income[2] * 2
+
+# compute probabilities for original and counterfactual
+p_orig = softmax.(eachrow(hcat(s1, s2_orig, zeros(length(s1)))))
+p_orig = hcat(p_orig...)'
+p_new = softmax.(eachrow(hcat(s1, s2_new, zeros(length(s1)))))
+p_new = hcat(p_new...)'
+
+p_diff = p_new[:, 2] - p_orig[:, 2]
+precis(DataFrame(p_diff=p_diff))
+

+
+    

+

+

+
+

+

+
+
+

+
+    

+
+
+

+
┌────────┬─────────────────────────────────────────────────────┐
+│  param │   mean     std    5.5%     50%   94.5%    histogram │
+├────────┼─────────────────────────────────────────────────────┤
+│ p_diff │ 0.0432  0.0394  0.0026  0.0315  0.1278  █▆▄▃▂▁▁▁▁▁▁ │
+└────────┴─────────────────────────────────────────────────────┘
+

+

+

+
+

+

+
+

+

+

+

+
Code 11.59

+
+

+

+

+

+

+
In [58]:

+

+    

+
Random.seed!(1)
+
+N = 500
+family_income = rand(Uniform(), N)
+b = [-2, 0, 2]
+career = [
+    begin
+        score = (0.5 * 1:3) + b * inc
+        p = softmax(score)
+        sample(Weights(p))
+    end
+    for inc in family_income
+]
+
+@model function m11_14(career, family_income)
+    a ~ MvNormal([0, 0], 1.5)
+    b ~ MvNormal([0, 0], 1)
+    
+    for i ∈ eachindex(career)
+        income = family_income[i]
+        s = [a[1] + b[1] * income, a[2] + b[2] * income, 0]
+        p = softmax(s)
+        career[i] ~ Categorical(p)
+    end
+end
+
+m11_14_chain = sample(m11_14(career, family_income), NUTS(), 1000)
+m11_14_df = DataFrame(m11_14_chain)
+precis(m11_14_df)
+

+
+    

+

+

+
+

+

+
+
+

+
+    

+
+
+

+
┌───────┬────────────────────────────────────────────────────────────┐
+│ param │    mean     std     5.5%      50%    94.5%       histogram │
+├───────┼────────────────────────────────────────────────────────────┤
+│  a[1] │ -2.7538    0.39  -3.3927  -2.7355  -2.1699  ▁▁▁▂▃▆██▇▅▂▁▁▁ │
+│  a[2] │ -1.2474  0.2198  -1.6024  -1.2415   -0.896   ▁▁▁▂▄▆███▆▄▂▁ │
+│  b[1] │ -1.6594  0.7035  -2.7403   -1.658  -0.5104      ▁▁▃▇█▇▄▂▁▁ │
+│  b[2] │ -1.0656  0.4021  -1.7315  -1.0556  -0.4297  ▁▁▂▃▅▇██▆▄▂▁▁▁ │
+└───────┴────────────────────────────────────────────────────────────┘
+

+

+

+
+

+

+
+

+

+

+

+
Code 11.60

+
+

+

+

+

+

+
In [60]:

+

+    

+
d = DataFrame(CSV.File("data/UCBadmit.csv", skipto=2, 
+        header=[:id, :dept, :gender, :admit, :reject, :applications]))
+

+
+    

+

+

+
+

+

+
+
+

+
+    
Out[60]:

+
+
+
+

+
12 rows × 6 columns
iddeptgenderadmitrejectapplications
Int64String1String7Int64Int64Int64
11Amale512313825
22Afemale8919108
33Bmale353207560
44Bfemale17825
55Cmale120205325
66Cfemale202391593
77Dmale138279417
88Dfemale131244375
99Emale53138191
1010Efemale94299393
1111Fmale22351373
1212Ffemale24317341

+

+
+

+
+

+

+
+

+

+

+

+
Code 11.61

+
+

+

+

+

+

+
In [61]:

+

+    

+
@model function m_binom(admit, applications)
+    a ~ Normal(0, 1.5)
+    p = logistic(a)
+    @. admit ~ Binomial(applications, p)
+end
+
+m_binom_chain = sample(m_binom(d.admit, d.applications), NUTS(), 1000)
+m_binom_df = DataFrame(m_binom_chain);
+
+@model function m_pois(admit, reject)
+    a ~ MvNormal([0, 0], 1.5)
+    λ = exp.(a)
+    admit ~ Poisson(λ[1])
+    reject ~ Poisson(λ[2])
+end
+
+m_pois_chain = sample(m_pois(d.admit, d.reject), NUTS(), 1000)
+m_pois_df = DataFrame(m_pois_chain);
+

+
+    

+

+

+
+

+

+

+

+
Code 11.62

+
+

+

+

+

+

+
In [62]:

+

+    

+
m_binom_df.a .|> logistic |> mean
+

+
+    

+

+

+
+

+

+
+
+

+
+    
Out[62]:

+
+
+
+
+

+
0.387568484041908

+

+
+

+
+

+

+
+

+

+

+

+
Code 11.63

+
+

+

+

+

+

+
In [63]:

+

+    

+
a1 = m_pois_df."a[1]" |> mean
+a2 = m_pois_df."a[2]" |> mean
+exp(a1)/(exp(a1)+exp(a2))
+

+
+    

+

+

+
+

+

+
+
+

+
+    
Out[63]:

+
+
+
+
+

+
0.38722562179640474

+

+
+

+
+

+

+