optimiz-rs/examples/notebooks/02_mcmc_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "dc5d5825",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:08.719062Z",
     "iopub.status.busy": "2026-02-16T16:15:08.718752Z",
     "iopub.status.idle": "2026-02-16T16:15:10.156553Z",
     "shell.execute_reply": "2026-02-16T16:15:10.155128Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OptimizR MCMC Module Loaded!\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "from optimizr import mcmc_sample\n",
    "\n",
    "np.random.seed(42)\n",
    "print(\"OptimizR MCMC Module Loaded!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddfbd617",
   "metadata": {},
   "source": [
    "## Example 1: Inferring Parameters of a Normal Distribution\n",
    "\n",
    "Given observed data $\\{x_1, \\ldots, x_n\\}$, infer $\\mu$ and $\\sigma$.\n",
    "\n",
    "### Likelihood\n",
    "$$L(\\mu, \\sigma | \\mathbf{x}) = \\prod_{i=1}^n \\frac{1}{\\sqrt{2\\pi\\sigma^2}} \\exp\\left(-\\frac{(x_i - \\mu)^2}{2\\sigma^2}\\right)$$\n",
    "\n",
    "### Log-Likelihood\n",
    "$$\\log L(\\mu, \\sigma | \\mathbf{x}) = -\\frac{n}{2}\\log(2\\pi) - n\\log(\\sigma) - \\frac{1}{2\\sigma^2}\\sum_{i=1}^n (x_i - \\mu)^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f929c18e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:10.160126Z",
     "iopub.status.busy": "2026-02-16T16:15:10.159742Z",
     "iopub.status.idle": "2026-02-16T16:15:10.443829Z",
     "shell.execute_reply": "2026-02-16T16:15:10.442725Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True parameters: μ=5.0, σ=2.0\n",
      "Sample mean: 4.792\n",
      "Sample std: 1.807\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1sAAAHbCAYAAAA9E4D0AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAjvlJREFUeJzs3XmcjeX/x/HXmd2YGYxZDbJlzRJjGUskGS0iJbSQFm36WkqSkC2F/FKRUkqJUKFSClFikDV7dhlmMcOs5sxy7t8fJ4dhMGPmzDnG+/l4nMfMfd3Xfd2f+z73WT7nvu/rMhmGYSAiIiIiIiJFysXRAYiIiIiIiJRESrZERERERETsQMmWiIiIiIiIHSjZEhERERERsQMlWyIiIiIiInagZEtERERERMQOlGyJiIiIiIjYgZItERERERERO1CyJSIiIiIiYgdKtkRErqBdu3aYTCZMJhOPP/64o8O5bp3bhyaTic8//9zR4ZQ4zrh/q1SpYovpjTfesJWvXr06V7xHjhxxWIznHDlyJFdMq1evdnRIIlJCKNkSkRvC9u3bef7556lfvz5ly5bFw8OD4OBg2rdvz+TJk0lKSnJ0iHKRN954I9cXYBcXFzw9PSlfvjx16tSha9euzJgxg5SUlCJdryMSlwvXaTKZcHNzw8fHh4oVK9KyZUv+97//sWnTJrvHceE+r1Klit3XVxyUSImII7k5OgAREXvKzs7mpZde4r333rtkXlxcHHFxcaxatYq3336br776io4dOzogSskPwzDIzMwkMTGRxMRE9u7dy5IlSxg+fDiffvopXbt2dXSIRSYnJ4e0tDTS0tKIjo4mKiqK999/n27duvHJJ59Qrly5XPUnTZpk+79p06bFHW6ehg8fbvsRo2XLlg6O5sr8/f1z7cPq1as7MBoRKUmUbIlIifbiiy8yY8YM23SFChV46KGHCAgIYMeOHXzzzTfk5ORw6tQpOnfuzG+//UarVq0cGHHRSE5Oxs/Pz9FhFKnXXnuNMmXKEB8fz++//85ff/0FQGJiIt26dWPu3Ln07NnTwVEWXnh4OD169CA9PZ39+/fzww8/2JKW7777jiNHjrBmzRq8vb1ty7z88suOCvcSKSkp+Pr68vTTTzs6lHzz8/Nzqn0oIiWIISJSQq1du9YAbI/GjRsbSUlJueqsXLnScHFxsdWpV6+ekZOTY5vftm1b27w+ffoYBw4cMLp37274+/sbpUqVMlq1amUsX778knUfOXLE6Nevn1GjRg3Dy8vL8PT0NCpUqGC0bNnSGDRokLF79+5Llvnjjz+MHj16GJUqVTI8PDwMX19fo0WLFsYHH3xgZGZmXlL/wm377LPPjMWLFxsRERFG6dKljTJlyhivv/66bX6VKlUuWX7Pnj252vjzzz9t83JycowvvvjCuPPOO43AwEDD3d3dCAgIMO6++25j6dKlee7vrKwsY8KECUaNGjUMDw8Po1q1asbYsWONzMzMS2LNj1GjRuVa7vDhw7nmL1682PD09LTN9/HxMeLi4mzzDx06ZAwYMMBo3bq1UbFiRcPb29vw8PAwKlSoYNx7773G999/n6u9C5/rvB433XSTre6nn35qdO/e3ahdu7ZRvnx5w83NzfD19TUaNmxovPLKK0Z8fHy+tvGcC9fTp0+fXPNOnz5tdOrUKVedoUOHXnb5i/fvZ599ZrRt29YWZ9myZY2aNWsaDz30kDFt2jTDMAxj1apVV9z2C9u98Hm56aabjFOnThnPP/+8ERYWZri4uBj/93//ZxiGYdx00022eqNGjbLFc/G6Dh06ZHzwwQdGvXr1bK+TQYMGGcnJybm24+L1Xujw4cO52ly1atUlMeT1aNu27RWXv9A333xj3H333UZwcLDh7u5ulC1b1oiIiDAmT55spKWlXfE5/eyzz4xff/3VaNeunVG6dGnDx8fH6NSpk7Fz585LlhORkkXJloiUWH369Mn1hWfFihV51uvVq1eueqtXr7bNu/ALeEREhOHv73/JFzYXFxdjwYIFtmViY2ONwMDAK37J+/DDD3PF8Nprr12xfps2bYzU1NRcy1w8/8LpMmXKGAcOHMhVtm7dulzLjxgxwjavZs2atvL09HSjQ4cOV4xn8ODBl+zHnj175ln3nnvusUuyZRiGMWnSpFx13nzzTdu8H3744aoJxOjRo/N8rq+WbDVp0uSKdcPCwozo6Oh8badhXDnZMgzDSElJMYKDg211fHx8DLPZnOfyF+7fi/fhxY/g4GDDMK492QoICDBq166dq15Bk62Lj49zj6ZNmxpnz57Nc1uKM9nKzs42HnrooSu2U6dOHePEiROXfU5btWplmEymS5YrX758rh8IRKTk0WWEIlJirVmzxvZ/uXLluOOOO/Ks16NHD+bNm5drubZt215SLyoqigoVKjB06FBSUlL49NNPMZvNWCwW+vXrR8eOHSlTpgzffvst8fHxtvX27duX8uXLc+LECfbu3ZsrLoCvv/6aN9980zYdGRlJq1atiI2NZfbs2aSmprJmzRoGDRrExx9/fNltDQgIoGfPnpQvX55du3ZRvXp1brvtNv744w8A5s6dS0REhG2ZC7e5b9++tv8HDRrEihUrAPDw8KBnz57cfPPN7Nixg4ULF2IYBlOmTKFJkyY8/PDDAHzzzTd8/fXXtjZq1KjBQw89RHR0NF9++WWeMReFJ554gldeeQXDMABYtWoVw4YNA8DNzY1GjRoRHh5OYGAgfn5+pKWlsXbtWlatWgXA2LFjefLJJwkLC+O5557j3nvvZciQIbb2e/ToQXh4OABlypSxlQcFBdG5c2eqV6+Ov78/rq6uREdHM3/+fBISEoiOjmbcuHFMnz69SLbTx8eHnj17MnXqVABSU1PZtGnTVe+F+vDDD23/d+jQgXbt2pGWlsa///7Ln3/+ydmzZwHrPUqTJk3i119/Zfny5YD12H3ttddsy+d1L9ipU6c4deoUHTp0oFWrVsTHxxMcHFygbVu6dCldunShYcOG/Pzzz7bLQ//66y8mTpzIyJEjC9TehYYPH86RI0dyvb6effZZ2z1ZlSpVumobb775JgsWLLBNt2jRgo4dO7Jnzx4WLlwIwJ49e3jkkUf47bff8mxj7dq11K5dm27durFt2zZ++uknABISEvj000959dVXr3kbRcTJOTrbExGxl1KlStl+QW7UqNFl623dujXXr83PP/+8bd6FZzvc3d1znV356quvci03c+ZMwzAMY8qUKbayZ5555pL1paamGjExMbbpW2+91Va/d+/eueouWLDANs/Nzc1ISEiwzbtw3X5+fsbRo0cvWdfnn3+e6yxGdna2YRiGsXHjRlu5q6ur7SxMQkKC4ebmZps3a9asXO09//zztnm33nqrrTwyMtJWXqZMmVxxjh8//rJnXq4kP2e2DMMwgoKCbHXq1q17yfx9+/YZX3/9tfH+++8bkydPNiZNmmR4e3vblvniiy9y1c9vrGlpacaKFSuMjz/+2JgyZYoxadIko0uXLrZlq1Wrlq/tvHideZ3ZMgzDmD59eq56F55NvVzMfn5+tvKTJ09e0ubBgwdzTV/p7FFedQBj4MCBedbL75mtp59+2jYvMzPTqFevnm1exYoV8xXblc5M5ecSwcvVycnJyXU2OyIiwvYaMgzDeOWVV3Itt3XrVtu8C8srVaqU67LIC1/z3bp1y3P/iUjJoK7fRUTyqU2bNrm6w+7Rowfu7u626c2bNwPQqlUrTCYTAB999BFNmjThscceY9y4cSxbtgw3Nzfbr//p6els27bN1sYXX3yRq5vqhx56yDYvOzubjRs35hlb7969qVy58iXlDz74IL6+vgDExsbafnm/8KxWZGQkFSpUAGDDhg1kZ2fb5j3xxBO54rnwTM22bdtIT08HyNUteadOnfD397dNP/roo3nGXFSM/85qXezIkSO0atWKWrVq0bNnT1588UVefvllhgwZYosb4Pjx4wVe55QpUwgODqZDhw7069ePwYMHM2TIEJYsWVKodq/kctt5JW3atLH9f8stt3DPPfcwcOBAZs6cyYEDB6hWrVqh43r
      "text/plain": [
       "<Figure size 1000x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate synthetic data\n",
    "true_mu = 5.0\n",
    "true_sigma = 2.0\n",
    "n_obs = 100\n",
    "\n",
    "observed_data = np.random.normal(true_mu, true_sigma, n_obs)\n",
    "\n",
    "print(f\"True parameters: μ={true_mu}, σ={true_sigma}\")\n",
    "print(f\"Sample mean: {observed_data.mean():.3f}\")\n",
    "print(f\"Sample std: {observed_data.std():.3f}\")\n",
    "\n",
    "# Plot data\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.hist(observed_data, bins=20, density=True, alpha=0.6, color='skyblue', edgecolor='black')\n",
    "x_range = np.linspace(observed_data.min(), observed_data.max(), 100)\n",
    "plt.plot(x_range, stats.norm.pdf(x_range, true_mu, true_sigma), \n",
    "         'r-', linewidth=2, label=f'True: N({true_mu}, {true_sigma}²)')\n",
    "plt.xlabel('Value', fontsize=12)\n",
    "plt.ylabel('Density', fontsize=12)\n",
    "plt.title('Observed Data Distribution', fontsize=14, fontweight='bold')\n",
    "plt.legend()\n",
    "plt.grid(alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5626ba97",
   "metadata": {},
   "source": [
    "### Define Log-Likelihood Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "37ac93a1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:10.447063Z",
     "iopub.status.busy": "2026-02-16T16:15:10.446784Z",
     "iopub.status.idle": "2026-02-16T16:15:10.453064Z",
     "shell.execute_reply": "2026-02-16T16:15:10.451933Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log-likelihood at true params: -202.57\n"
     ]
    }
   ],
   "source": [
    "def log_likelihood_normal(params, data):\n",
    "    \"\"\"\n",
    "    Log-likelihood for Normal(μ, σ²) given data.\n",
    "    \n",
    "    Args:\n",
    "        params: [μ, σ]\n",
    "        data: observed data points\n",
    "    \"\"\"\n",
    "    mu, sigma = params\n",
    "    \n",
    "    # Ensure sigma is positive\n",
    "    if sigma <= 0:\n",
    "        return -np.inf\n",
    "    \n",
    "    n = len(data)\n",
    "    residuals = (data - mu) / sigma\n",
    "    \n",
    "    log_lik = -0.5 * n * np.log(2 * np.pi)\n",
    "    log_lik -= n * np.log(sigma)\n",
    "    log_lik -= 0.5 * np.sum(residuals**2)\n",
    "    \n",
    "    return log_lik\n",
    "\n",
    "# Test the function\n",
    "test_params = [5.0, 2.0]\n",
    "print(f\"Log-likelihood at true params: {log_likelihood_normal(test_params, observed_data):.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a36ddef7",
   "metadata": {},
   "source": [
    "### Run MCMC Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "60e64783",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:10.456697Z",
     "iopub.status.busy": "2026-02-16T16:15:10.456293Z",
     "iopub.status.idle": "2026-02-16T16:15:10.801889Z",
     "shell.execute_reply": "2026-02-16T16:15:10.800458Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running MCMC sampling...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Samples generated successfully\n",
      "Generated 20000 samples after burn-in\n",
      "\n",
      "Posterior estimates:\n",
      "μ: 4.792 ± 0.187\n",
      "σ: 1.840 ± 0.136\n",
      "\n",
      "True values: μ=5.0, σ=2.0\n"
     ]
    }
   ],
   "source": [
    "# MCMC parameters\n",
    "initial_params = [0.0, 1.0]  # Start far from true values\n",
    "param_bounds = [(-10, 10), (0.1, 10)]  # μ ∈ [-10, 10], σ ∈ [0.1, 10]\n",
    "proposal_std = 0.3  # Proposal step sizes\n",
    "n_samples = 20000\n",
    "burn_in = 2000\n",
    "\n",
    "print(\"Running MCMC sampling...\")\n",
    "samples = mcmc_sample(\n",
    "    log_likelihood_fn=lambda params: log_likelihood_normal(params, observed_data),\n",
    "    initial_params=initial_params,\n",
    "    param_bounds=param_bounds,\n",
    "    proposal_std=proposal_std,\n",
    "    n_samples=n_samples,\n",
    "    burn_in=burn_in\n",
    ")\n",
    "\n",
    "print(f\"\\nSamples generated successfully\")\n",
    "print(f\"Generated {len(samples)} samples after burn-in\")\n",
    "print(f\"\\nPosterior estimates:\")\n",
    "print(f\"μ: {samples[:, 0].mean():.3f} ± {samples[:, 0].std():.3f}\")\n",
    "print(f\"σ: {samples[:, 1].mean():.3f} ± {samples[:, 1].std():.3f}\")\n",
    "print(f\"\\nTrue values: μ={true_mu}, σ={true_sigma}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2306bc9b",
   "metadata": {},
   "source": [
    "## Visualize MCMC Results\n",
    "\n",
    "### Trace Plots - Check Convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "856b7ca8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:10.805389Z",
     "iopub.status.busy": "2026-02-16T16:15:10.805001Z",
     "iopub.status.idle": "2026-02-16T16:15:12.016617Z",
     "shell.execute_reply": "2026-02-16T16:15:12.015169Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABWwAAAMWCAYAAABlYM+1AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsnXd8HNXV93+zK+1KsiTLRXIB0xKw6T08piROXlLgCQESwhMHcPyEhOSJSQgEQkwJNoQYYiCQkEBoMQSIAVNNMcFgGVfcG+5FlmX1stqVVltn3j9md3ZmdvrO7OxK5/v52NqdveXMmXvv3Dlz7rkMx3EcCIIgCIIgCIIgCIIgCIIgCNfxuC0AQRAEQRAEQRAEQRAEQRAEwUMGW4IgCIIgCIIgCIIgCIIgiAKBDLYEQRAEQRAEQRAEQRAEQRAFAhlsCYIgCIIgCIIgCIIgCIIgCgQy2BIEQRAEQRAEQRAEQRAEQRQIZLAlCIIgCIIgCIIgCIIgCIIoEMhgSxAEQRAEQRAEQRAEQRAEUSCQwZYgCIIgCIIgCIIgCIIgCKJAIIMtQRAEQRAEQRAEQRAEQRBEgUAGW4IgiCEMwzDCv4aGBrfFIQiCIAiCIAhD0DyWIIjBDBlsCYIYdBxzzDGSCZzev/r6erdFtgWlc/N4PKipqcF//dd/4dFHH0UsFnOk7oaGBsyaNQuzZs3Co48+6kgdBEEQBEEQgx2ax9I8liAIAgAYjuM4t4UgCIKwk2OOOQYHDx40nH7JkiWYMmWKcwLlCYZhdNOcf/75+OSTT+D3+7PyHDhwAMccc4yluuvr6/HVr34VAHD00UeTlwNBEARBEIQFaB6rDs1jCYIYSpS4LQBBEITdLFiwAJFIRPj+3HPP4Z///CcAYOzYsXjttdck6U899VTFcuLxODiOg8/nc05Yh7jjjjtwySWXoL+/Hy+99BL+9a9/AQBWrlyJJ598EjfddJPLEhIEQRAEQRByaB5L81iCIAiAQiIQBDEIOeecc3DhhRcK/4466ijhN7/fL/ntyCOPRE1NjbD0qqWlBdOnT0ddXR38fj+2b9+Orq4u/PznP8d5552HcePGoaysDOXl5fjiF7+In/70p9i/f7+iHG+//Ta+/e1vY+zYsfD5fBg9ejQuuOACPP/885J03d3duPvuu3H66aejsrIS5eXlOPnkkzFr1iz09fVZ0sHxxx+PCy+8EN/85jcxb948icfBkiVLDJWxefNmTJs2DUcffTT8fj+qq6vxpS99CQ899BCi0aiQ7phjjhG8EgDg4MGDijHFZs2aJRybPn26IRmmT58u5Jk1a5bkN/GSwcGyHJAgCIIgiKENzWMHzzwWADZs2IBrrrkGxxxzDMrKylRDW9TU1BgukyCIoQF52BIEQYj4yle+gj179kiOdXR04B//+EdW2n379mHfvn14/fXXsW7dOhx33HEAAI7j8OMf/xjz5s2TpO/q6sLKlStRW1uLH/3oRwCAvXv34qtf/Sqampokabdv347Zs2fj9ddfx9KlSzFy5EjL5+TxeFBdXS18NxL/a/78+Zg2bRri8bgk39q1a7F27VrMnz8fS5YsQVVVlWW5CIIgCIIgCPugeSxPocxjFy9ejEsvvVQiB0EQhFHIw5YgCEJEY2Mj7r33Xnz44Yd46qmnMHr0aIwYMQL33nsvXnnlFSxatAj19fVYuHAhrr32WgBAT08PHn74YaGMp59+WjLJveqqq7BgwQK88847uOuuu1BbWyv8du211wqT3K9+9at48803sXDhQnzlK18BAGzbtg2//vWvLZ9POBzGU089hS1btgjHzjzzTM08ra2tuP7664XJ5SWXXIKFCxfi73//O4YPHw4AWL9+PX73u98B4Jfu/eUvfxHyjx07FsuWLRP+jRs3zrL8BEEQBEEQhDFoHls489hkMokbbrhBkOP444/H/Pnz8eqrr+Kkk04S0lVUVKC+vh6LFi2yVA9BEIMYjiAIYpBzzz33cAA4ANzRRx8t+e3AgQPCbwC4v/zlL4plvPvuu9x///d/c2PHjuVKSkokeQBwZ511lpD2nHPOEY5feeWVqnJt3bpVSFdaWsp9+OGH3LJly7hly5ZxCxYskPwWCoV0z1Muk9K/MWPGcC0tLYp5Dhw4wHEcxz322GPCsdraWm5gYEBI//jjjwu/VVdXc4lEguM4jluyZImqjnPhRz/6kVDuPffcI/nt6KOPFn5bsmSJbXUSBEEQBEEUCjSPLc557GeffSaRb926dcJv27Ztk/y2bNkyW+okCGJwQSERCIIgRHzve9/LOvbcc8/h+uuv18zX09MjfN6+fbvw+bvf/a5qHnG6eDyOb37zm4rp4vE4du3ahbPPPltTBi1KSkrw7W9/Gw8//DDGjh2rmXbnzp3C53POOQdlZWXC9wsvvFD4HAwG0dzcjAkTJliWiyAIgiAIgrAHmscWzjxWHBu4vLxccv4nn3wyampqEAgEAAAHDhyQyEYQBAFQDFuCIAgJSsueHnjgAeHzt771LfziF7/AiBEjsG7dOtx8880AAJZlHZXL7KYN6d11PR4PKisr8cUvfhEVFRUOSecOHMe5LQJBEARBEETBQPPYwqGkxLipxeOhSJUEQWRDBluCIAgRDMNkHWtsbBQ+z507F6eccgoAoL6+XrGMk046CevWrQMAvPnmm0KMsDQcx4FhGJx44onCsfLycrS0tAixtcT09/dj2LBhps4jvbuuFSZNmiR8Xr9+PSKRiOCdsGLFCuG36upq4cFAPNF0atIvvg6RSATt7e3CdzLeEgRBEAQx1KF5bOHMY08++WTh88DAADZu3CjE392+fbvgXStPSxAEkYZe5RAEQeiQ3jUXAP7whz/go48+wgMPPID7779fMf1PfvIT4fMbb7yBH/zgB3jzzTfx7rvvYvbs2bjhhhsAAKeeeirOPfdcAPxE7mtf+xqef/55fPLJJ3jllVfwwAMP4Bvf+AYuv/xyB88um6uvvlrwYmhvb8dVV12Fd999F08++STuvPNOId21114reA+MGjVKON7c3IwXXngBy5cvx/r164Xjs2bNAsMwYBgG06dPNy3X/Pnz8fjjj+PDDz/E//7v/yISiQi/LVq0CJ9//rnpMgmCIAiCIAYzNI91Zx574oknYvLkycL3qVOn4tVXX8WCBQtw9dVXC8fPOussnHHGGVZPlyCIwYy7IXQJgiCcx8xmDUo8+eSTihsfTJkyRbHcZDLJXXfddaobJlx++eVC2t27d3NHHnmk5gYLX/nKVwydpzjPP//5T9N50ps1cBzH/fvf/+ZKS0tVZTr77LO53t5eIX0ikVA8jy984QtCGvF1+NGPfmRIPvGmY/J/I0aMkHyXb0pGEARBEARR7NA81lieQpzH7t69mxs/frzmJmrbt283VBZBEEMP8rAlCILQ4Wc/+xmeeOIJTJo0CWVlZTj++OPx6KOP4ve//71ieo/HgxdeeAELFizApZdeirq6OpSUlGDkyJE4//zzccUVVwhpjz/+eGzZsgW///3vceaZZ6KyshJ+vx9HHXUUvvzlL+P+++/Hk08+maczzfCDH/wAa9aswbXXXosJEyagtLQUlZWVOPvss/GnP/0Jy5cvR3V1tZDe6/XizTffxJe//GXHYoxdeumlwjWYMmUK1q5dixtvvBEVFRWoq6vDCSec4Ei9BEEQBEEQxQrNY92bxx5//PHYtGkTfvvb3+LEE09EWVkZysrKMHHiRNxyyy3YtGmTJLQEQRCEGIbjKPAfQRAEUZhMnz4dzz//PADgnnvuwaxZs9wViCAIgiAIgiAIgiAchjxsCYIgCIIgCIIgCIIgCIIgCgQy2BIEQRAEQRAEQRAEQRAEQRQIZLAlCIIgCIIgCIIgCIIgCIIoECiGLUEQBEEQBEEQBEEQBEEQRIFAHrYEQRAEQRAEQRAEQRAEQRAFAhlsCYIgCIIgCIIgCIIgCIIgCoQStwVwA5Zl0dzcjKqqKjAM47Y4BEEQBEEQhAYcxyEUCmH8+PHweIa2vwHNYwmCIAiCIIoHq/PYIWmwbW5uxoQJE9wWgyAIgiAIgjDBoUOHcOS
      "text/plain": [
       "<Figure size 1400x800 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(14, 8))\n",
    "\n",
    "# Trace plots\n",
    "axes[0, 0].plot(samples[:, 0], linewidth=0.5, alpha=0.7)\n",
    "axes[0, 0].axhline(true_mu, color='red', linestyle='--', linewidth=2, label='True μ')\n",
    "axes[0, 0].set_xlabel('Sample', fontsize=11)\n",
    "axes[0, 0].set_ylabel('μ', fontsize=11)\n",
    "axes[0, 0].set_title('Trace Plot: μ', fontsize=13, fontweight='bold')\n",
    "axes[0, 0].legend()\n",
    "axes[0, 0].grid(alpha=0.3)\n",
    "\n",
    "axes[0, 1].plot(samples[:, 1], linewidth=0.5, alpha=0.7, color='orange')\n",
    "axes[0, 1].axhline(true_sigma, color='red', linestyle='--', linewidth=2, label='True σ')\n",
    "axes[0, 1].set_xlabel('Sample', fontsize=11)\n",
    "axes[0, 1].set_ylabel('σ', fontsize=11)\n",
    "axes[0, 1].set_title('Trace Plot: σ', fontsize=13, fontweight='bold')\n",
    "axes[0, 1].legend()\n",
    "axes[0, 1].grid(alpha=0.3)\n",
    "\n",
    "# Posterior distributions\n",
    "axes[1, 0].hist(samples[:, 0], bins=50, density=True, alpha=0.6, color='skyblue', edgecolor='black')\n",
    "axes[1, 0].axvline(true_mu, color='red', linestyle='--', linewidth=2, label='True μ')\n",
    "axes[1, 0].axvline(samples[:, 0].mean(), color='green', linestyle='-', linewidth=2, label='Posterior mean')\n",
    "axes[1, 0].set_xlabel('μ', fontsize=11)\n",
    "axes[1, 0].set_ylabel('Density', fontsize=11)\n",
    "axes[1, 0].set_title('Posterior Distribution: μ', fontsize=13, fontweight='bold')\n",
    "axes[1, 0].legend()\n",
    "axes[1, 0].grid(alpha=0.3)\n",
    "\n",
    "axes[1, 1].hist(samples[:, 1], bins=50, density=True, alpha=0.6, color='orange', edgecolor='black')\n",
    "axes[1, 1].axvline(true_sigma, color='red', linestyle='--', linewidth=2, label='True σ')\n",
    "axes[1, 1].axvline(samples[:, 1].mean(), color='green', linestyle='-', linewidth=2, label='Posterior mean')\n",
    "axes[1, 1].set_xlabel('σ', fontsize=11)\n",
    "axes[1, 1].set_ylabel('Density', fontsize=11)\n",
    "axes[1, 1].set_title('Posterior Distribution: σ', fontsize=13, fontweight='bold')\n",
    "axes[1, 1].legend()\n",
    "axes[1, 1].grid(alpha=0.3)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b20c718",
   "metadata": {},
   "source": [
    "### Joint Posterior Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ad2eebc7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:12.020057Z",
     "iopub.status.busy": "2026-02-16T16:15:12.019756Z",
     "iopub.status.idle": "2026-02-16T16:15:12.343522Z",
     "shell.execute_reply": "2026-02-16T16:15:12.342502Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzQAAALCCAYAAAAIzpbHAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAmGxJREFUeJzs3XecFEXex/Fv9+xszoQlCwJKEkRBBQyoKMGcAyrmuxNUxHzmiBhRUe85E+qJWTGcCQNGUEHxVAwgIKDksMvmmel+/thl2IVddmpmdwL7efOaFzs9v66qru7pmZqqrrZc13UFAAAAAAnIjnUBAAAAACBcNGgAAAAAJCwaNAAAAAASFg0aAAAAAAmLBg0AAACAhEWDBgAAAEDCokEDAAAAIGHRoAEAAACQsGjQAAAAAEhYNGiAHdzQoUNlWZYsy9KZZ54Z6+KgAWeeeWZwfw0dOjTWxWkUm7fHsixNnTo11sWRJHXu3DlYphtvvDG4fObMmbXKu2TJkpiVcbMlS5bUKtPMmTNjXSQAiCs0aIA4s/UXqnj5AlhTfV8GQ1Vz+2o+UlNTtdNOO+nkk0/Wp59+2vgFD8GO2KCI1Nb7KSkpSZmZmerQoYMGDx6siy66SHPmzGnyctx4443BMnTu3LnJ84sGGisAELmkWBcAQNP6xz/+ocMPP1yS1KdPnxiXZvsqKiq0dOlSLV26VC+88IJuvfVWXXPNNbEuVlSdfPLJwf3UsWPHGJemboFAQCUlJSopKdGff/6pWbNm6cEHH9Sxxx6rxx57THl5ebXi77rrruDfAwcOjHZx63TNNdeosLBQkjR48OAYl2b78vPza9Vh165dY1gaAIg/NGiAHdxJJ50U6yJs14ABA3TSSSfJcRwtWLBAzzzzjCoqKiRJ1113nUaNGqX+/fvHuJRNr6ioSNnZ2RoxYoRGjBgRtXw3bdqkrKyskOM376/S0lItWLBAb775ZrBh8Oqrr2rJkiX67LPPlJ6eHlznsssua/Ryh2vz9p533nmxLkrIsrOz46oOASDuuADiyscff+xKCj6efPLJWq/PmTPHPf30093OnTu7KSkpbkZGhtu7d293woQJ7rJly7ZJ74ADDgimNWbMmODyxYsX18rn448/dp977jl3r732ctPS0tzc3Fz3+OOPd5cuXRpcZ8yYMbXWqesRiprxNcvkuq776KOP1nr9uuuuq/X68uXL3csuu8zt06ePm5GR4aakpLg77bSTO3r0aPerr77aJi+fz+fed9997j777OPm5OS4Ho/Hzc/Pd3v16uWefvrp7nPPPee6rus++eSTDW7bxx9/HEw3EAi4Tz/9tHvIIYe4rVq1cr1er9uyZUt31KhR7n//+99tyrH1fl2wYIF71113uT169HCTk5Pdo446aps6PuCAA7ZJx3T7b7jhhmB6O+20k7t27Vr3ggsucNu3b+/atu3ed999299Z7vb314YNG9wRI0bUirnyyivrXX/r4/nJJ590DzjgALdFixZuUlKSm5ub6+6yyy7uiSee6D700EN11l1dj83phrq9O+20UzDuhhtuqHc/LVq0yJ0yZYrbu3dvNyUlxW3Xrp17ySWXuEVFRdut55rqeq9tXYa6Hpv3f33r1/Tyyy+7o0aNcgsKClyv1+vm5ua6gwYNcu+++263pKRku/v0ySefdN9//3136NChbkZGhpuZmemOGDHC/fHHH7dZDwDiEQ0aIM5sr0Fz3333ubZt1/sFKCcnZ5svO6E2aPbdd9860+zevbtbVlbmum50GjQ//vhjrdfPO++84GuffPKJm5eXV2/etm2799xzT630Girz3nvv7bquWYOmtLTUHTZs2HZjJ0yYsN39ut9++9V6HkqDJpztr/lFu2XLlm6PHj1qrRNpg8Z1XXfTpk1uQUFBMCYzM9OtqKioc/2ax3PNstX1KCgoqLPuQm3QbG97Q23QHHbYYXXmN3DgwOD7Yut8o9mg8fv97oknnrjddHr27On+9ddf9e7TIUOGuJZlbbNeixYt3NWrVzd4fABArDHkDEgQn376qSZMmCDXdSVJnTp10imnnKLi4mI9+eSTKi0tVWFhoY477jgtXLhwm+sYGvL5559r4MCBGj58uD7++GN98cUXkqQFCxZo+vTpOvnkk4PXd9x+++3asGGDJOmQQw7RoYce2mjbOWvWrFrP27RpI0nauHGjjj322GC+aWlpOuuss5Sdna3nnntOf/zxhxzH0WWXXaY999xTBxxwgIqLi/Wf//wnmNZxxx2nPfbYQ4WFhfrjjz/0ySefBF8bOHCg7rrrLr3wwgvBC9x33nln/eMf/wjGbL524ZJLLtEHH3wgSUpOTtbJJ5+s7t2764cfftBLL70k13V17733as8999Spp55a53Z+9tln6t27t4444gi5riuPx7Pdegln+7e2du1arV27VsOGDdOQIUO0Zs0aFRQUbDffUGRmZurkk0/W/fffL0kqLi7WnDlzGrw25ZFHHgn+PWzYMA0dOlQlJSVatmyZPv/8c5WVlUmqqve77rpL77//vmbMmCFJysvL0z//+c/g+nVdm9MY2/vf//5XRx11lPr166d33nlH33zzjSTpm2++0Z133qnrr7/eKL2arrnmGi1ZskS33357cNnf//734HEWyjVUt99+u1588cXg83322UeHHnqofv75Z7300kuSpJ9//lmjR4/WRx99VGcaX3zxhXr06KFjjz1W8+bN09tvvy1JWrdunR5//HFdddVVYW8jAERFjBtUALZSXw/NUUcdFVyWlZXlrlq1KrjO22+/Xe+v7qH20Oy1115uZWWl67quW1lZ6bZu3bre3ob6ft0OVc18BwwY4N51113upEmT3HPPPddNSUkJvmZZlvvtt9+6rlvVO1VzvbfffjuY3qpVq9zMzMxtejvWr18fXJadnV2r18B1XddxHHfRokW1ljU05GvdunVuUlJSMOaJJ56o9foFF1wQfK1///7B5Vvv13322afWL/wN5R/O9rvutr0g48ePr3unbEfN9evqoXFd13344Ydrxb344ot1rl+zhyY7Ozu4fMWKFduk+fvvv9d6vr1eENPtDbWHpmYPYWVlpdu7d+/gax06dAipbNvrYQllOFl9MYFAwM3Pzw8uHzRokOv3+4PrXXHFFbXW++6774Kv1VzesWPHWkPo+vfvH3zt2GOPrbP+ACCeMG0zkCBq9lyMGDFCrVu3Dj4fOXKkWrVqVWdsqM4991x5vV5JktfrVZcuXYKvbe4VaApz5szR5ZdfriuvvFKPPfZYcEIAqWqa3s0TAtTcplatWmnkyJHB561bt671fHNsXl6eevfuLanqovsuXbro6KOP1uWXX66nn35af/31V63tDMVXX30lv98ffH722WfXmnb34YcfDr42b948lZaW1pnOZZddptTU1JDzDWf763LttdeGnKcJt7rn0MR+++0X/LtPnz467LDDNH78eD366KNauHChdt5554jLFen2nn766cG/vV6vTjzxxODz5cuXa9WqVRGlH4lff/1V69evDz4/7bTTavX0jRkzplZ8fcfF6aefXmtiiF122SX4d1O+9wGgsTDkDEgQNb+41DVspqCgQGvWrJEU3peQre/rkZKSEvzbcRzj9MKRnJysgoIC7bPPPrrgggtq3QcmlO3frOb2T5s2Taeccormz5+vv/76S6+//nrwNdu2dfHFF+vee+8NuYw1y9EQ13W1bt26WjN+bdajR4+Q09k6X5Ptr6lly5Zq0aKFUb6h+u2332o9b9++fYPrPPLIIzrxxBM1e/ZsrVu3LjjUabMTTzxRzz33nGw7vN/eGmN7a/5wIG1b9xs3btxm2daNu5qN9Ma09bG4dTm2fl7fcREP730AiAQNGiBB5Ofna/Xq1ZJU56/CNZeZXj8jKdg7s5llWcZphGPMmDEh3Tw0Pz8/+LfJ9vft21c//fSTfvjhB3377bdasGCBvv32W73zzjtyHEf33XefjjjiCB144IEhlbdmOaSq62natWtXb3xOTk6dyzMyMkLKr658w93/pnmGqqSkRC+88ELweVZWlgYMGNDgeh07dtS
      "text/plain": [
       "<Figure size 1000x800 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 8))\n",
    "\n",
    "# 2D histogram\n",
    "plt.hist2d(samples[:, 0], samples[:, 1], bins=50, cmap='Blues')\n",
    "plt.colorbar(label='Sample Density')\n",
    "\n",
    "# Mark true values\n",
    "plt.scatter([true_mu], [true_sigma], c='red', s=200, marker='*', \n",
    "           edgecolors='black', linewidths=2, label='True values', zorder=5)\n",
    "\n",
    "# Mark posterior mean\n",
    "plt.scatter([samples[:, 0].mean()], [samples[:, 1].mean()], \n",
    "           c='green', s=200, marker='o', edgecolors='black', \n",
    "           linewidths=2, label='Posterior mean', zorder=5)\n",
    "\n",
    "plt.xlabel('μ', fontsize=12)\n",
    "plt.ylabel('σ', fontsize=12)\n",
    "plt.title('Joint Posterior Distribution', fontsize=14, fontweight='bold')\n",
    "plt.legend(fontsize=11)\n",
    "plt.grid(alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "049df073",
   "metadata": {},
   "source": [
    "## Example 2: Logistic Regression with MCMC\n",
    "\n",
    "Bayesian inference for binary classification.\n",
    "\n",
    "### Model\n",
    "$$P(y=1 | \\mathbf{x}, \\boldsymbol{\\beta}) = \\frac{1}{1 + \\exp(-\\boldsymbol{\\beta}^T \\mathbf{x})}$$\n",
    "\n",
    "### Log-Likelihood\n",
    "$$\\log L(\\boldsymbol{\\beta}) = \\sum_{i=1}^n \\left[y_i \\log p_i + (1-y_i) \\log(1-p_i)\\right]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "072639a3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:12.346960Z",
     "iopub.status.busy": "2026-02-16T16:15:12.346678Z",
     "iopub.status.idle": "2026-02-16T16:15:13.562088Z",
     "shell.execute_reply": "2026-02-16T16:15:13.560700Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Features shape: (200, 3)\n",
      "Class distribution: [100 100]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArkAAAIoCAYAAABtZPNSAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAo0tJREFUeJzt3Xl8U1X6P/BPkqbpmrZAC7RsbdlkFQUE0QqyKl9HxVEZAVEHf8rgwuaCCy6jo+OCKzLOOAOj4ODMuI0zKiACVQEVAYGiQsvegi2F7lva3N8fx6Rpm6RZ7k3uvfm8Xy9epclNenJyk/vcc5/zHIMkSRKIiIiIiHTEGO4GEBERERHJjUEuEREREekOg1wiIiIi0h0GuURERESkOwxyiYiIiEh3GOQSERERke4wyCUiIiIi3WGQS0RERES6wyCXiIiIiHSHQS6Rjt10000wGAwwGAwYO3ZsuJujCY8++qizz3r16hXu5ng1duxYZ1tvuummNvd/+umnuOSSS5CUlOTcLjk5GYA29o1evXo52/joo4+GuzlEpDEMcolUbvPmzc4Dves/k8mEpKQkDB06FHfccQcOHDgQ7qaqkiRJ+PjjjzFr1iz07dsXVqsVZrMZnTt3xvjx4/HHP/4RJ0+eDHczZbd3715ceeWVyM3NRUVFRbib04ZeAtjWn8uoqCgkJCSgW7duuPDCC3HXXXdhx44dsv5NLZ2IEYVTVLgbQESBsdvtqKiowJ49e7Bnzx6sXLkSmzdvxogRI5zbTJ8+HYMGDQIAdO/ePVxNDZvjx4/jhhtuwJdfftnmvuLiYnz++ef4/PPP8cMPP2DVqlWhb2CQ5s6di//7v/8DAOf77PDee++hoaEBAGCxWDB//nx06tQJMTExALSxbzz44IMoLy8HAFx44YVhbo1vmpqaUF1djerqahQWFmLbtm145ZVXMG3aNLzxxhtISUkJdxOJIgaDXCKNuf766zF8+HA0Njbim2++wfvvvw8AqKmpwZNPPokPPvjAue2UKVMwZcqUMLXUNxUVFbBarbI/788//4xLLrkEhw8fdt6WmZmJX/3qV+jcuTPOnj2L7du3uw2AteL666/3eN/Ro0ed/x8xYgSefvrpFvdrYd+49dZbw90EvwwfPhzXX389ampqcPDgQXz00UfOIP29997DkSNH8MUXXyAuLi7MLSWKEBIRqdqmTZskAM5/K1eubHH/oEGDnPf169evxX2zZ8923nfJJZe0uK/1c65fv14aO3asFB8fLyUkJEhTpkyR9u3b16Y9zzzzjHTllVdKffr0kVJSUqSoqCgpKSlJGjFihPTEE09IVVVVbR7T+m998MEH0ujRo6X4+HgpKSlJWrp0qfP+bt26SU1NTS0ev2/fvhbPsX379nb7bfr06S0eM3fuXMlms7XZ7sCBA9Lq1audvz/yyCPOx/Ts2bPFtu+99540c+ZMafDgwVJaWppkNpul+Ph46ZxzzpHmzZsnHT58uM3zl5SUSIsWLZIGDBggxcXFSWazWercubM0YsQIad68edK2bdtabP/hhx9KkydPltLS0qSoqCgpMTFRysrKkq688krpD3/4Q4u+ueSSS5xtnT17tiRJkrRy5coWr7v1P8d23vYNSZKk06dPS48//rh0wQUXSMnJyVJ0dLSUnp4uTZo0SVq7dq1zu9LSUumee+6RLr30Uqlnz55SQkKCZDabpbS0NGnChAnSm2++Kdntduf2rn/X0z+Hnj17Om975JFH2rRxx44d0qxZs6RevXpJFotFio+PlwYOHCgtXLhQOn78eJvtW/fXgQMHpOnTp0sdO3aULBaLNGzYMOmDDz5o8zhv3PWtw9mzZ6UpU6a02Oa+++5rsc1f//pX6dprr5X69+8vdezY0fmeDx06VLr33nulkpIS57atvwvc/XN8Pxw6dEi6++67pYsuukjq1q2bFBcX53wP/+///k/6z3/+49frJNIiBrlEKucpyG1sbJS2bdsmWa1Wj8GKr0HumDFjJIPB0OaA2bFjR6m4uLjF4zp27Oj1IDt48GCpsrLS49+6+OKLW/yelJQkFRUVSWaz2Xnb//73vxaPdw2CBwwY0G6fFRUVtXg95557bpvA2RNvQe4111zj9bVbrVZpz549zu1ra2ulfv36eX2Ma9DTXoAKQKqtrXVur1SQ+80330hdunTx+BxXXnmlc9u9e/e22+abb77Zub1cQe4LL7wgGY1Gj8+RlJQkbdq0qcVjXPtryJAhUmJiYpvHGQwG6bPPPmtnL2nmLciVJEmqrKyUOnfu7NwmISFBqq+vd95//vnne+2LjIwMqbCwUJIk/4Lcjz76qN1tH3vsMZ9fJ5EWMV2BSGNuvvlm3HzzzW1uNxqNuOeeewJ6zq+++gr9+/fHtGnTsHv3bnz88ccAgNLSUvz1r3/F/fff79y2W7duGDduHHr27ImUlBRIkoTDhw/jnXfeQXV1Nfbu3YvXXnsN9957r9u/9cUXX6BTp06YPn06OnbsiLy8PHTt2hXXXHMN1q5dCwB44403cPnllzsf869//avF62/Ppk2bIEmS8/fZs2fDaAx+nm1ycjImTZqEc845BykpKYiOjsbPP/+M999/H8eOHUNFRQXuu+8+Z/9t2rQJP/30EwAgJiYGv/3tb5GRkYFTp04hPz8fW7ZsafH8K1ascP5/xIgR+L//+z80Njbi+PHj+Prrr/HDDz+028YRI0bg2WefxTvvvOOc8JSVlYW5c+cCaJu721plZSV+9atf4dSpU87bLr30UowZMwYVFRVt0juMRiPOOeccjBw5El26dEFycjLq6uqwa9cufPTRR5AkCStXrsTtt9+OkSNHOnOB//CHP+Ds2bMAgIkTJ2LSpEntvjaH3NxcLFy40Pke9+jRA7/5zW9QVVWFlStXoqamBuXl5bjmmmuQn5/vNg92z549SElJwYIFC1BbW4u//OUvaGpqgiRJePbZZzF+/Hif2+NNQkICpk+fjpdeegkAUFVVhR07djhzjNPS0nDFFVcgOzsbHTp0gMlkQmFhId555x2UlpaisLAQTzzxBF577TVkZ2fj2Wefxfr167FhwwYAQEpKCh544AHn33Pk5EdFReHcc8/F8OHDkZqaCqvViurqanz11VfYtGkTAOD3v/+9c58k0iMGuUQ68cQTT2Dq1KkBPbZ79+745ptvkJiYCAA477zzsGvXLgDAt99+22Lb3bt3o7y8HFu3bsWxY8dQXV2Nc845B+effz5yc3MBAOvWrfMY5FqtVnz33Xfo0aNHi9vvvPNOZ5D73//+Fz///DM6d+6Mffv2OYO7qKgozJo1q93XU1hY2OL3/v37t/sYX7zxxhuw2WzYvn07Dh48iIqKCnTr1g3jx4/HypUrAQCff/45bDYbzGYz6urqnI+95JJL8Oqrr7Z4vvr6epw+fdr5u+v2L7/8MkaNGtVi+yNHjiA6OtprGwcOHIiBAwdi3759ziC3e/fuWLx4sU+vcdWqVS0C3CeffLJFEAUAhw4dcv5/wIAB2L9/P44dO4Zvv/0Wp06dgtlsxsUXX4zvvvvO+V6sW7cOI0eOdOYCv/rqq84g98ILL/S5fQCwbNkyZ4CbmJiIb7/9FmlpaQCAqVOnOk+Qzpw5g7///e+YP39+m+cwGAzYuHEjhg0bBkCchLz44osA2u7zwerXr1+L3133z48//hg1NTXYtm0bDh06hKqqKmRmZuKiiy7Chx9+CED0HdD8PlZVVTmDXKvV6rbvHP184MAB7Nq1CyUlJTCbzbj88svx9ddfo6amBo2Njfj88899+kwRaRGDXCKNcUw8a2pqQl5eHv7xj3+gsbERDzzwAGw2G5YuXer3c86aNcsZ4AJA3759nUGuIxABREWH+++/Hy+99JJz5r47J06c8HjfjTfe2CbABUSgc95552Hnzp2w2WxYtWoV7rvvvhajuJdffjk6d+7s12uT05o1azB//vwWgWlrjsC1a9euGDFiBCwWC+rr67Fu3ToMHDgQQ4YMQd++fTFs2DCMHz++xSjaxRdfjD179gAQo5ujR49Gnz59MGDAAOTk5GDw4MGKv0bXkdr
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate synthetic classification data\n",
    "from sklearn.datasets import make_classification\n",
    "\n",
    "X, y = make_classification(n_samples=200, n_features=2, n_redundant=0,\n",
    "                          n_informative=2, random_state=42, n_clusters_per_class=1)\n",
    "\n",
    "# Add intercept\n",
    "X_with_intercept = np.column_stack([np.ones(len(X)), X])\n",
    "\n",
    "print(f\"Features shape: {X_with_intercept.shape}\")\n",
    "print(f\"Class distribution: {np.bincount(y)}\")\n",
    "\n",
    "# Visualize data\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.scatter(X[y == 0, 0], X[y == 0, 1], c='blue', label='Class 0', alpha=0.6, s=50)\n",
    "plt.scatter(X[y == 1, 0], X[y == 1, 1], c='red', label='Class 1', alpha=0.6, s=50)\n",
    "plt.xlabel('Feature 1', fontsize=12)\n",
    "plt.ylabel('Feature 2', fontsize=12)\n",
    "plt.title('Binary Classification Data', fontsize=14, fontweight='bold')\n",
    "plt.legend()\n",
    "plt.grid(alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6163cb52",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:13.565204Z",
     "iopub.status.busy": "2026-02-16T16:15:13.564814Z",
     "iopub.status.idle": "2026-02-16T16:15:14.500460Z",
     "shell.execute_reply": "2026-02-16T16:15:14.499241Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running MCMC for logistic regression...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Posterior estimates:\n",
      "β₀ (intercept): 0.000 ± 0.000\n",
      "β₁: 0.000 ± 0.000\n",
      "β₂: 0.000 ± 0.000\n"
     ]
    }
   ],
   "source": [
    "def log_likelihood_logistic(beta, X, y):\n",
    "    \"\"\"\n",
    "    Log-likelihood for logistic regression.\n",
    "    \"\"\"\n",
    "    z = X @ beta\n",
    "    # Numerically stable sigmoid\n",
    "    p = 1 / (1 + np.exp(-np.clip(z, -500, 500)))\n",
    "    p = np.clip(p, 1e-10, 1 - 1e-10)  # Avoid log(0)\n",
    "    \n",
    "    log_lik = np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))\n",
    "    \n",
    "    # Add weak prior: beta ~ N(0, 10²)\n",
    "    log_prior = -0.5 * np.sum(beta**2) / 100\n",
    "    \n",
    "    return log_lik + log_prior\n",
    "\n",
    "# Prepare data tuple\n",
    "logistic_data = (X_with_intercept, y)\n",
    "\n",
    "# MCMC for logistic regression\n",
    "initial_beta = np.zeros(3)  # [intercept, coef1, coef2]\n",
    "beta_bounds = [(-10, 10)] * 3\n",
    "beta_proposal_std = 0.1\n",
    "\n",
    "print(\"Running MCMC for logistic regression...\")\n",
    "beta_samples = mcmc_sample(\n",
    "    log_likelihood_fn=lambda params: log_likelihood_logistic(params, X_with_intercept, y),\n",
    "    initial_params=initial_beta,\n",
    "    param_bounds=beta_bounds,\n",
    "    proposal_std=beta_proposal_std,\n",
    "    n_samples=15000,\n",
    "    burn_in=1500\n",
    ")\n",
    "\n",
    "print(f\"\\nPosterior estimates:\")\n",
    "print(f\"β₀ (intercept): {beta_samples[:, 0].mean():.3f} ± {beta_samples[:, 0].std():.3f}\")\n",
    "print(f\"β₁: {beta_samples[:, 1].mean():.3f} ± {beta_samples[:, 1].std():.3f}\")\n",
    "print(f\"β₂: {beta_samples[:, 2].mean():.3f} ± {beta_samples[:, 2].std():.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cc9ff7c",
   "metadata": {},
   "source": [
    "### Visualize Decision Boundary with Uncertainty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2eb5d0d5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:14.503549Z",
     "iopub.status.busy": "2026-02-16T16:15:14.503254Z",
     "iopub.status.idle": "2026-02-16T16:15:14.807722Z",
     "shell.execute_reply": "2026-02-16T16:15:14.806187Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1AAAALCCAYAAADQwjLNAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAxopJREFUeJzs3Xl4U1X6B/DvTZom3WlrC7SUrexlEUUQwUrZEXHBGTdAdAYXRB02R1FRcR9F3AcdZ3644Iw67s4oiwwUZRERESwgtJStLC2ldE+bNvf3xyFpkibpTZqbrd/P8/CUJLfJycnN7XnP8h5JlmUZRERERERE1CJNoAtAREREREQUKhhAERERERERKcQAioiIiIiISCEGUERERERERAoxgCIiIiIiIlKIARQREREREZFCDKCIiIiIiIgUYgBFRERERESkEAMoIiIiIiIihRhAUVi55ZZbIEkSJEnCqFGjVH2tQ4cOWV9LkiRs2LBB1dcjdXTt2tX6GT722GOBLk5Y2rBhg9135dChQ149D79zwY2fT5PWnvO++s6Qf/DvSNvDAIo84nhRf/vttwNdJFUEU0PAthyWfxqNBjExMejduzduvfVW/PzzzwErHwWW43dSkiRERkYiISEB3bt3x9ixY7FkyRIcPXo00EVtE5x9HpIkQavVIiEhAYMGDcLdd9+N/fv3B7qoFCCBCo4ee+yxFl/XNhBQuxMyEEKp/RJKZW2LIgJdACJfuuGGG9C/f38AQEZGhqqvlZSUhOeff956OzMzU9XXsyXLMmpqarB//37s378f77//Pj7//HNcfvnlfitDuHjooYdQXl4OALjkkksCXBrfMJlMMJlMqKioQGFhIdatW4cnnngCixcvxuLFi6HR+LfvLDMz0+67kpSU5NXzBPI711pmsxkVFRXYtWsXdu3ahRUrVmDDhg246KKLAl00UoGvznkKDeH4d4TcYwBFYWXixImYOHGiX14rPj4eCxcu9MtrWQwZMgTXX389jEYjtmzZgq+//hqAaDA//PDDQRtAVVdXIyoqyu8NdyVuu+22QBfBp66//noMGTIE5eXl2LFjB1avXo3GxkY0Njbisccew8mTJ7F8+XK/likjI8Mn35VAfOday/J5NDQ0YNu2bfjss88AADU1NXjqqafw+eefB7aAIayiogLx8fGBLoZTvjrnqfX8cZ6E298RUkAm8sD69etlANZ/K1asUPR7x44dkxcuXCj3799fjomJkfV6vdylSxd52rRp8g8//OD0d06fPi3feeedcvv27WWDwSBfeOGF8kcffdSsDIWFhdbfmTlzpvX+yy67zO75du3aJU+bNk3u0qWLHBkZKRsMBjkjI0POycmRH3jgAfnYsWOyLMtyly5d7J7f8Z/leQsLC+3uX79+fbP3sHbtWvm6666TO3fuLOv1ejk+Pl7OysqSZ8+eLZeUlCiqO9vXmDlzpt1jw4YNsz6m1+ud/n5BQYF8zz33yH369JGjo6Nlg8Eg9+3bV77//vtdlmHXrl3yFVdcIcfFxclxcXHyxIkT5Z9//ll+9NFHra/XpUsXu9+xrbdHH31U/u677+QxY8bI8fHxMgC5rKzMeuzOnTvlW2+9Ve7evbtsMBjkmJgY+fzzz5efeuopuaqqqll5Dh06JN9+++1yjx49ZIPBIOv1ejktLU2+5JJL5Hnz5sl79uyxO37FihXyZZddJicnJ8sRERFyu3bt5F69esnXXXed/Prrr7stt6Pt27fLM2bMkLt27Srr9Xo5JiZGzsrKkufPny8fPXq02fGXXXaZ3ee1f/9++YYbbpCTk5NlvV4vDx48WP7888+b/d6KFStcntfutPSd3LNnj9ytWze7Y7755ptmz+PpZyLLslxVVSW/+OKLcnZ2tpyUlCTrdDq5ffv2cnZ2tvzaa6+5LKPte6uqqpKXLFkiDx48WI6NjZUjIiLklJQUedCgQfKsWbPsyqrkO/fxxx/Ll19+udy+fXtZp9PJ7dq1k4cPHy4vXbpUrq6ubna8Y92tWbNGHjVqlBwTEyPHxsbKEydOlH/99dcWPoUmLX0e/fv3tz7Wu3fvZr/f0NAg/+Mf/5BHjx5tPX+TkpLkUaNGyX/7299kk8nk9vUczxtX57fj7xUUFMivv/66PGDAAFmv18spKSnyH//4R/nMmTPNylhdXS3ff//9cqdOnWS9Xi/369dPfu211+SDBw+6/HwOHjwo/+lPf5JHjhwpd+rUSY6OjpYjIyPltLQ0+YorrpC//PLLZq/j+J2orq6WH3zwQblbt25yRESEfM8999id24sWLWr2HAsXLrQ+3rdvXyefWJMzZ87IGo3Genxubq71sTfffNN6/zXXXGP3ecXGxlof++CDD9x+Lu7+tthe4x1//+DBg/Jbb70lDxo0qMXPxxXb67era4zt+eL4N9TxXNq+fbs8efJkOSEhQY6KipJHjhwpf/fdd05f+/Tp0/Ljjz8uDxs2TG7Xrp31sx8/fry1zmxt3LhRvv766+WMjAw5MjJSjouLky+++GL5tddek+vr65sd7/id+/zzz+Xhw4fLMTExckJCgt112dk/279n//jHP+Tf//73cp8+fazfwbi4OHnQoEHyn//8Z6d/N335PVNS1vz8fLtzdfXq1c3KNGTIEOvjd955p9PPhbzHAIo84k0AlZubKycmJrq8GGg0GvmFF16w+52ysjK5T58+To+fMmWKyz8CrgKovLw8OTo62u1FydJQ80UAZTab5VmzZrl9np9//llRnTv742oxdepU62Pp6enNfvfzzz93+77T09ObBR8//vijXYPA8s9gMMjjxo1z+gfHsd6GDx8ua7Vau9+3BFB//etf5YiICJdl6tevn3zixAnr8546dUpOSUlxW5fLly+3Hu/YSHD81759e5fldgygXnzxRbs/Uo7/EhISmjXibf/4DRw4UI6Li2v2e5Ikyd9++63d76kVQMmyLG/bts3umPHjx9s97ulnIssiMO/Zs6fL3xk0aJDLMtq+t1GjRrn9vK6//nrrse6+cw0NDfJ1113n9rn69u0rHz9+3O592D4+YsQIWZKkZr+XnJwsFxcXt+rzaGhokLds2WLtULC9llhUVVXJ2dnZbt/DyJEj5crKSkV1K8vKG3YjR450+nrZ2dl2z1dfXy9feumlTo+dPHmyy8/nq6++cvu+AMhLliyxey3H74Tj6/7pT3+Sn3/+eevttLQ0uaGhweX7f+6551r8/AYPHmw9/qmnnrLeP336dOv9KSkp1vt//PFHuzKdOnXK7efSUh24CqAmTJig6PNxx5cB1NChQ2WdTtesPHq9vtnflG3btskdOnRw+Z6vuuoqu+MffPBBt3V06aWXNuvUcXeeeBpAXXjhhW6PTU9Pl4uKilzWTWu/Z0rLavt9+/3vf29XHsfOjG3btjk5I6g1OIWPVHX27FlMnToVZWVlAICoqCjceuutiI+Px7/+9S8cPnwYZrMZCxcuxIUXXojLLrsMAPDwww9j37591ucZOXIkcnJy8N133+Grr77yuBzvvPMOampqAACdOnXC9OnTERMTg2PHjuHXX3/F1q1brcc+9NBDOHToEJ5++mnrfXfeead1vYWStVVLly7F3//+d+vt5ORkXHfddWjfvj3279+PL774wuP3YMtoNGLz5s1Yu3at9b7rrrvO7pjCwkLceOONqK2tBQBkZWXhmmuugdlsxvvvv4/Dhw+jqKgI1157LXbv3g2tVgtZlvGHP/wBVVVV1ue58cYb0b17d3z00Ud2r+fOli1bEB0djenTpyM9PR0///wztFotNm/ejLvvvhtmsxkAcPHFF2PixImorKzEO++8g9OnT2PPnj24+eabsWbNGgDAJ598gpKSEgBAYmIibr31ViQnJ+P48ePYt28fvvvuO7vXtp2eNnbsWIwaNQrV1dU4evQovv/+e2t9tGTjxo2YP38+ZFkGAHTu3Bk33ngjqqqqsGLFCtTU1KC8vBzXXnst8vPzkZiY2Ow5du3ahcTERMybNw+1tbV466230NjYCFmW8fzzz2PMmDGKytJaF110EQYNGoRffvnF+t4aGxu
      "text/plain": [
       "<Figure size 1000x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot decision boundaries from posterior samples\n",
    "plt.figure(figsize=(10, 8))\n",
    "\n",
    "# Plot data\n",
    "plt.scatter(X[y == 0, 0], X[y == 0, 1], c='blue', label='Class 0', alpha=0.6, s=50, zorder=3)\n",
    "plt.scatter(X[y == 1, 0], X[y == 1, 1], c='red', label='Class 1', alpha=0.6, s=50, zorder=3)\n",
    "\n",
    "# Create grid\n",
    "x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1\n",
    "x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1\n",
    "\n",
    "# Plot decision boundaries from random posterior samples\n",
    "n_boundary_samples = 100\n",
    "indices = np.random.choice(len(beta_samples), n_boundary_samples, replace=False)\n",
    "\n",
    "for idx in indices:\n",
    "    beta = beta_samples[idx]\n",
    "    # Decision boundary: β₀ + β₁x₁ + β₂x₂ = 0\n",
    "    # => x₂ = -(β₀ + β₁x₁) / β₂\n",
    "    if abs(beta[2]) > 0.01:  # Avoid division by zero\n",
    "        x1_line = np.array([x1_min, x1_max])\n",
    "        x2_line = -(beta[0] + beta[1] * x1_line) / beta[2]\n",
    "        plt.plot(x1_line, x2_line, 'gray', alpha=0.02, linewidth=0.5, zorder=1)\n",
    "\n",
    "# Plot mean decision boundary\n",
    "beta_mean = beta_samples.mean(axis=0)\n",
    "if abs(beta_mean[2]) > 0.01:\n",
    "    x1_line = np.array([x1_min, x1_max])\n",
    "    x2_line = -(beta_mean[0] + beta_mean[1] * x1_line) / beta_mean[2]\n",
    "    plt.plot(x1_line, x2_line, 'black', linewidth=3, label='Mean boundary', zorder=2)\n",
    "\n",
    "plt.xlim(x1_min, x1_max)\n",
    "plt.ylim(x2_min, x2_max)\n",
    "plt.xlabel('Feature 1', fontsize=12)\n",
    "plt.ylabel('Feature 2', fontsize=12)\n",
    "plt.title('Logistic Regression: Decision Boundary with Uncertainty', fontsize=14, fontweight='bold')\n",
    "plt.legend()\n",
    "plt.grid(alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c4facb3",
   "metadata": {},
   "source": [
    "## MCMC Diagnostics\n",
    "\n",
    "### Autocorrelation - Check Mixing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0ab8cd82",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-02-16T16:15:14.810825Z",
     "iopub.status.busy": "2026-02-16T16:15:14.810478Z",
     "iopub.status.idle": "2026-02-16T16:15:15.622911Z",
     "shell.execute_reply": "2026-02-16T16:15:15.621002Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABW0AAAHqCAYAAAB/bWzAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZYlJREFUeJzt3Xl8VNXdx/HvTUIWwCwQSAgJJLghBQRJiSgWfEgBtU9BKYrFAqmFRytaxAVoRRSqGKkWURS1Ili1grK4VZSCICiySqOyqBgkARLWJCSEbHOfPyBDJplJJpCZuZN83q/XkMy5Z86ce++E+c1vzj3HME3TFAAAAAAAAADAEgJ83QEAAAAAAAAAwFkkbQEAAAAAAADAQkjaAgAAAAAAAICFkLQFAAAAAAAAAAshaQsAAAAAAAAAFkLSFgAAAAAAAAAshKQtAAAAAAAAAFgISVsAAAAAAAAAsBCStgAAAAAAAABgISRtAaAJ2rt3rwzDsN88acGCBfbn6d+/v0efCwAAAE0D8SyAxi7I1x0A0PSsX79e11xzjUPZN998o5/97GcN0v6aNWu0Zs0aSVKPHj00dOjQBmkXzi1fvlzbt2+XJPXv359AFgAANHrEs40L8SwAKyJpC8DrFixY4LRs1qxZDdL+mjVr9Oijj0qSRo8eTZDrYcuXL9fChQvt96sHuddff73WrVsnSYqIiPBm1wAAADyCeLZxIZ4FYEUkbQF41cmTJ/X222/XKH/99df1xBNPKDAw0Ae9sqaysjKZpqng4GCXdYqLixUSEqKAAOvOdtO2bVu1bdvW190AAABoEMSz7iOeBYBzZ93/FQE0SsuWLVNBQYEk6corr1Tnzp0lSTk5Ofr4449r1O/fv799/qjqIxqqzmG1d+9e+7xWlaMSJGnhwoX2OomJifZy0zS1YMECXXvttWrVqpWaNWum2NhYDRkyRKtWrXLa9+3bt2vMmDHq1KmTQkNDFR4erm7duumBBx5wqFdaWqpnnnlGffr0UUREhIKDg5WQkKDf/va32rp1q0Pd6nNxHTx4UGPGjFHbtm0VEhKiHTt2aM2aNQ778N133+mmm25SVFSUmjdvbj+eFRUVevHFF3XNNdcoKipKwcHB6tixo8aOHavMzEz3TpCkRx55RAMGDFCHDh3UsmVLBQcHKy4uTjfeeKP9Mj1J9n5VHZXw6KOP1pjvq7Y5wOp7Hqoeqx07dugvf/mLOnbsqJCQEHXu3Fmvv/56jcfU9hpyJTEx0f6YqvvszbnTAACANRHPEs9WZdV4VpKWLl2qwYMHKyYmRs2aNXN47qo3RnIDFmUCgBelpqaakkxJ5nPPPWfOmDHDfn/48OE16vfr18++/dVXX3XYVlkuyczMzDQzMzMdyqrfOnbsaJqmaZaXl5tDhw6tte5jjz3m8FwvvfSSGRQU5LRuRESEvV5hYaF55ZVXumw3KCjIXLhwob1+9T5ffPHFDve/+uor89NPP3V4rjZt2jjUOX78uHny5Enz2muvdfm8kZGR5saNG10+b1UxMTEu2zEMw1yyZIlpmqZDv5zd+vXrZ5qmab766qs1ys71PNR2rCpvX3zxhduvIVc6duxof8ynn37q1nEDAABNA/Es8WwlK8ezjzzySK39qnobMmSIW20C8C5G2gLwmqysLK1evVqSFBQUpFtuuUUjR460b3/vvfd0/Pjxc26/Xbt2WrdundLS0uxl1113ndatW6d169bpnXfekSTNnTtXy5cvlyQ1a9ZM06dP17///W/94Q9/sD/uL3/5izZt2iRJ2rFjh+68806Vl5dLOr0YxMKFC7VixQrNnj1bXbp0sT9u6tSp+vLLLyVJLVu21DPPPKMPPvjA/u11eXm5xo0bp6ysLKf7sG/fPk2fPl0ff/yxXnrpJUVHRztsz8/PV1lZmWbPnq1PPvlEzzzzjEJCQvTII4/o008/lSQlJSXp1Vdf1SeffKI77rhDkpSXl6dbb73Vvg+1mTBhghYuXKgPP/xQa9as0ccff6zHH39ckmSapqZOnSpJ6tmzp9atW6frrrvO/ti0tDT78X722WdrfZ76nofq9u/fr6efflrvvvuuunbtai+fM2dOnfsIAABwLohniWersmo8+/3332vGjBn2+9dff73ef/99Pf/88w5z8l555ZVat26dZs6ceV7PB8BDfJ01BtB0PPbYY/Zvc2+44QZ7+VVXXWUvf/755x0eU5+RCZWmTZtmLx89enSNfvTs2dO+/e6773bYlpycbN/2xz/+0TRN07z//vvtZfHx8WZhYaHT/bPZbGbr1q3tdZ966in7tpKSEjMuLs6+7cknnzRNs+YIgTlz5tRot/oIgPfee6/G81YdrfD000+b69ats9/atWtn37ZixQqnz1vVt99+a44aNcpMSkoyQ0JCnH4bX1BQYK8/evRoe/m0adNq9N/VyIT6ngfTdDznlcfQNE3zrbfespdfccUVTs9PfTDSFgAAOEM8SzzrD/Fsenq6va02bdqYxcXF9m3PPfecfVtoaKhZWlp6Xs8FwHMYaQvAa6rOFVV1RELV3+szR9O52rVrl/33vn37Omyrer+y3o4dO+xlgwYNUosWLZy2e/jwYR09etRpW8HBwerdu7fTPlQ1bNiwWvseEhKiX/3qVzWe9/Dhw/b7EydO1DXXXGO/HTx40L7tm2++qbX9r7/+Wr1799Zrr72mzMxMlZSUOK13PiNIKtX3PFQ3YMAA+++tW7e2/37s2LHz7hsAAIAzxLPEs1VZNZ798ccf7b8nJycrNDTUab9OnTqlnJyc83ouAJ5D0haAV3zxxRf67rvv7Pd/+9vf2ie+v+uuu+zlmzZt0s6dO+33qy72VPVSqKpBXWPSrl27WrfHxMSc1wJYhYWFtW5/9tlnVVRUJEm6+OKL9cYbb+izzz5zWLBBkmw22zn3oaG0atXK/ntQUJD9d9M0PfacnmwbAABYG/Gse4hn3eepeLZqW3UJCCAtBFgVf50AvKLqqIT61I2KirL/np2dbf/9/fffd/n4qoGHs2CscoVfSfr8888dtlW9X1mv6hxfn3zyiT0IrFQZVLVp08bhG/KqbZWVlWnz5s1O+1BVXQGss+1t2rRxmCvs448/lmmaNW6FhYWaNm1are3v27fP/vs999yj3/72t7rmmmsUGBjo8jF1HW9X6nsefKXqMan6u0QSFwCApoR4lni2OqvGsz/72c/sv2/dulWnTp1y2q9WrVrVmWQH4Dvuf/0CAOfo1KlTWrRokf3+/fffrwsvvNChztdff63nn39ekvTPf/5Tjz32mAIDA3XJJZfY6/z973/XBRdcoIKCAv3tb39z+XxVA81169bpww8/VEREhGJjY3XRRRdpzJgx+uqrryRJ8+bNU9u2bdWrVy8tW7bMIRAdPXq0JGnMmDH6+9//roqKCmVlZalfv37605/+pJiYGH333Xd66623tH79ehmGoVGjRunvf/+7JGnatGlq1qyZOnXqpFdeeUX79++XdPqSsBEjRpzTsXTGMAylpaVp1qxZkqRRo0Zp8uTJ6tq1qwoLC7Vv3z59+eWX+uCDD1RQUFBrW506dbL//o9//EOJiYk6duyYHnroIZePqXq8//3vf6tv375q3ry5OnbsqISEBJePq+95OFf9+/fX2rVrJUmvvvqqxowZU6/HP/roowoNDVVoaGiNDwlvvPGGBgwYQLALAEAjRzxLPOuMVePZm2++WZMmTdKJEyd06NAhDR8+XHfccYeysrL0l7/8xaFfjLQFLMz70+gCaGrefPNN+2T34eHhZklJSY06x48fN5s1a2av99FHH5mmaZo//vijGRwcXGPhgK5du7pcuGHHjh1mQEBAjcfcfvvtpmmaZnl5uTl06FCnCxJU3v7617869O+FF14wAwMDndaNiIiw1yssLDSvvPJKl+0GBQWZCxcutNd3Z2Grqgs3dOzY0WmdkydPmv379691n6q27+p5MzIyHM5D5a1621WP98cff+z0uWbMmGGapuuFG87lPLjqQ23HqLbFP1ypuhBZ9VtUVJTD/aoLlQEAgMaJeJZ41t/i2aVLlzp93VXerrz
      "text/plain": [
       "<Figure size 1400x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Low autocorrelation at large lags indicates good mixing!\n"
     ]
    }
   ],
   "source": [
    "from statsmodels.graphics.tsaplots import plot_acf\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "plot_acf(samples[:, 0], lags=100, ax=axes[0], alpha=0.05)\n",
    "axes[0].set_title('Autocorrelation: μ', fontsize=13, fontweight='bold')\n",
    "axes[0].set_xlabel('Lag', fontsize=11)\n",
    "\n",
    "plot_acf(samples[:, 1], lags=100, ax=axes[1], alpha=0.05)\n",
    "axes[1].set_title('Autocorrelation: σ', fontsize=13, fontweight='bold')\n",
    "axes[1].set_xlabel('Lag', fontsize=11)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(\"Low autocorrelation at large lags indicates good mixing!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2ed2ccd",
   "metadata": {},
   "source": [
    "## Key Takeaways\n",
    "\n",
    "1. **MCMC samples from complex distributions** using only likelihood evaluations\n",
    "2. **Metropolis-Hastings** uses proposal distribution and accept/reject steps\n",
    "3. **Burn-in period** allows chain to converge to target distribution\n",
    "4. **Diagnostics** (trace plots, autocorrelation) verify convergence and mixing\n",
    "5. **OptimizR provides 50-100x speedup** for likelihood evaluations\n",
    "6. **Bayesian inference** naturally quantifies uncertainty in parameters\n",
    "\n",
    "## Further Reading\n",
    "\n",
    "- Gelman, A., et al. (2013). \"Bayesian Data Analysis\" - Chapter 11\n",
    "- Brooks, S., et al. (2011). \"Handbook of Markov Chain Monte Carlo\"\n",
    "- Betancourt, M. (2017). \"A Conceptual Introduction to Hamiltonian Monte Carlo\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rhftlab",
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}