update code pymc 3.11

code_3_11/Chp1/01 Thinking Probabilistically.ipynb

@@ -0,0 +1,229 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "from scipy import stats\n",
+    "import arviz as az"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "az.style.use('arviz-darkgrid')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([-1.75915361,  0.20230817,  0.79831042])"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "μ = 0.\n",
+    "σ = 1.\n",
+    "X = stats.norm(μ, σ)\n",
+    "x = X.rvs(3)\n",
+    "x"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mu_params = [-1, 0, 1]\n",
+    "sd_params = [0.5, 1, 1.5]\n",
+    "x = np.linspace(-7, 7, 200)\n",
+    "_, ax = plt.subplots(len(mu_params), len(sd_params), sharex=True, sharey=True,\n",
+    "                     figsize=(9, 7), constrained_layout=True)\n",
+    "for i in range(3):\n",
+    "    for j in range(3):\n",
+    "        mu = mu_params[i]\n",
+    "        sd = sd_params[j]\n",
+    "        y = stats.norm(mu, sd).pdf(x)\n",
+    "        ax[i,j].plot(x, y)\n",
+    "        ax[i,j].plot([], label=\"μ = {:3.2f}\\nσ = {:3.2f}\".format(mu, sd), alpha=0)\n",
+    "        ax[i,j].legend(loc=1)\n",
+    "ax[2,1].set_xlabel('x')\n",
+    "ax[1,0].set_ylabel('p(x)', rotation=0, labelpad=20)\n",
+    "ax[1,0].set_yticks([])\n",
+    "plt.savefig('B11197_01_01.png', dpi=300);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data = np.genfromtxt('../data/mauna_loa_CO2.csv', delimiter=',')\n",
+    "plt.plot(data[:,0], data[:,1])\n",
+    "plt.xlabel('year')\n",
+    "plt.ylabel('$CO_2$ (ppmv)')\n",
+    "plt.savefig('B11197_01_02.png', dpi=300)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "n_params = [1, 2, 4]  # Number of trials\n",
+    "p_params = [0.25, 0.5, 0.75]  # Probability of success\n",
+    "\n",
+    "x = np.arange(0, max(n_params)+1)\n",
+    "f,ax = plt.subplots(len(n_params), len(p_params), sharex=True, sharey=True,\n",
+    "                    figsize=(8, 7), constrained_layout=True)\n",
+    "\n",
+    "for i in range(len(n_params)):\n",
+    "    for j in range(len(p_params)):\n",
+    "        n = n_params[i]\n",
+    "        p = p_params[j]\n",
+    "\n",
+    "        y = stats.binom(n=n, p=p).pmf(x)\n",
+    "\n",
+    "        ax[i,j].vlines(x, 0, y, colors='C0', lw=5)\n",
+    "        ax[i,j].set_ylim(0, 1)\n",
+    "        ax[i,j].plot(0, 0, label=\"N = {:3.2f}\\nθ = {:3.2f}\".format(n,p), alpha=0)\n",
+    "        ax[i,j].legend()\n",
+    "\n",
+    "        ax[2,1].set_xlabel('y')\n",
+    "        ax[1,0].set_ylabel('p(y | θ, N)')\n",
+    "        ax[0,0].set_xticks(x)\n",
+    "plt.savefig('B11197_01_03.png', dpi=300)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "params = [0.5, 1, 2, 3]\n",
+    "x = np.linspace(0, 1, 100)\n",
+    "f, ax = plt.subplots(len(params), len(params), sharex=True, sharey=True,\n",
+    "                     figsize=(8, 7), constrained_layout=True)\n",
+    "for i in range(4):\n",
+    "    for j in range(4):\n",
+    "        a = params[i]\n",
+    "        b = params[j]\n",
+    "        y = stats.beta(a, b).pdf(x)\n",
+    "        ax[i,j].plot(x, y)\n",
+    "        ax[i,j].plot(0, 0, label=\"α = {:2.1f}\\nβ = {:2.1f}\".format(a, b), alpha=0)\n",
+    "        ax[i,j].legend()\n",
+    "ax[1,0].set_yticks([])\n",
+    "ax[1,0].set_xticks([0, 0.5, 1])\n",
+    "f.text(0.5, 0.05, 'θ', ha='center')\n",
+    "f.text(0.07, 0.5, 'p(θ)', va='center', rotation=0)\n",
+    "plt.savefig('B11197_01_04.png', dpi=300)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.figure(figsize=(10, 8))\n",
+    "\n",
+    "n_trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150]\n",
+    "data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]\n",
+    "theta_real = 0.35\n",
+    "\n",
+    "beta_params = [(1, 1), (20, 20), (1, 4)]\n",
+    "dist = stats.beta\n",
+    "x = np.linspace(0, 1, 200)\n",
+    "\n",
+    "for idx, N in enumerate(n_trials):\n",
+    "    if idx == 0:\n",
+    "        plt.subplot(4, 3, 2)\n",
+    "        plt.xlabel('θ')\n",
+    "    else:\n",
+    "        plt.subplot(4, 3, idx+3)\n",
+    "        plt.xticks([])\n",
+    "    y = data[idx]\n",
+    "    for (a_prior, b_prior) in beta_params:\n",
+    "        p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y)\n",
+    "        plt.fill_between(x, 0, p_theta_given_y, alpha=0.7)\n",
+    "\n",
+    "    plt.axvline(theta_real, ymax=0.3, color='k')\n",
+    "    plt.plot(0, 0, label=f'{N:4d} trials\\n{y:4d} heads', alpha=0)\n",
+    "    plt.xlim(0, 1)\n",
+    "    plt.ylim(0, 12)\n",
+    "    plt.legend()\n",
+    "    plt.yticks([])\n",
+    "plt.savefig('B11197_01_05.png', dpi=300)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    " <img src=\"B11197_01_06.png\" width=\"500\">"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "np.random.seed(1)\n",
+    "az.plot_posterior({'θ':stats.beta.rvs(5, 11, size=1000)})\n",
+    "plt.savefig('B11197_01_07.png', dpi=300)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<img src=\"B11197_01_08.png\" width=\"500\">"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}