{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "646d6207",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:02.718483Z",
     "iopub.status.busy": "2026-05-17T00:42:02.718275Z",
     "iopub.status.idle": "2026-05-17T00:42:03.180160Z",
     "shell.execute_reply": "2026-05-17T00:42:03.179713Z"
    }
   },
   "outputs": [],
   "source": [
    "import fdfi\n",
    "print('FDFI version:', fdfi.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "495b1ae3",
   "metadata": {},
   "source": [
    "# Confidence Intervals and Statistical Inference\n",
    "\n",
    "This tutorial covers statistical inference with FDFI, including confidence intervals, hypothesis testing, and feature selection.\n",
    "\n",
    "## What You'll Learn\n",
    "\n",
    "1. Computing confidence intervals with `conf_int()`\n",
    "2. One-sided vs two-sided tests\n",
    "3. Variance floor for stable inference\n",
    "4. Practical significance margins\n",
    "5. Statistically-driven feature selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b311614",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:03.181740Z",
     "iopub.status.busy": "2026-05-17T00:42:03.181604Z",
     "iopub.status.idle": "2026-05-17T00:42:03.419859Z",
     "shell.execute_reply": "2026-05-17T00:42:03.419408Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from fdfi.explainers import OTExplainer\n",
    "from fdfi.plots import confidence_interval_plot\n",
    "\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7588803f",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Create a model where we know the true feature importance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b49a0457",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:03.421661Z",
     "iopub.status.busy": "2026-05-17T00:42:03.421519Z",
     "iopub.status.idle": "2026-05-17T00:42:03.451169Z",
     "shell.execute_reply": "2026-05-17T00:42:03.450710Z"
    }
   },
   "outputs": [],
   "source": [
    "n_features = 10\n",
    "n_train = 500\n",
    "n_test = 100\n",
    "\n",
    "# True importance: features 0, 1, 2 are important; rest are noise\n",
    "true_coefs = np.array([2.0, 1.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])\n",
    "\n",
    "def model(X):\n",
    "    return X @ true_coefs\n",
    "\n",
    "X_train = np.random.randn(n_train, n_features)\n",
    "X_test = np.random.randn(n_test, n_features)\n",
    "\n",
    "# Create explainer\n",
    "explainer = OTExplainer(model, data=X_train, nsamples=100)\n",
    "results = explainer(X_test)\n",
    "\n",
    "print(\"True coefficients:\", true_coefs)\n",
    "print(\"Estimated importance:\", results[\"phi_X\"].round(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b5e5198",
   "metadata": {},
   "source": [
    "## Basic Confidence Intervals\n",
    "\n",
    "The `conf_int()` method computes pointwise confidence intervals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8ee9553",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:03.452686Z",
     "iopub.status.busy": "2026-05-17T00:42:03.452574Z",
     "iopub.status.idle": "2026-05-17T00:42:04.114912Z",
     "shell.execute_reply": "2026-05-17T00:42:04.114364Z"
    }
   },
   "outputs": [],
   "source": [
    "# Two-sided 95% confidence intervals\n",
    "ci = explainer.conf_int(alpha=0.05, alternative=\"two-sided\")\n",
    "\n",
    "print(\"Two-sided 95% Confidence Intervals:\")\n",
    "print(\"-\" * 70)\n",
    "print(f\"{'Feature':>8} {'Estimate':>10} {'SE':>10} {'CI Lower':>10} {'CI Upper':>10} {'P-value':>10}\")\n",
    "print(\"-\" * 70)\n",
    "for i in range(n_features):\n",
    "    sig = \"*\" if ci[\"pvalue\"][i] < 0.05 else \"\"\n",
    "    print(f\"{i:>8} {ci['score'][i]:>10.4f} {ci['se'][i]:>10.4f} \"\n",
    "          f\"{ci['ci_lower'][i]:>10.4f} {ci['ci_upper'][i]:>10.4f} {ci['pvalue'][i]:>10.4f} {sig}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c08a2b30",
   "metadata": {},
   "source": [
    "## Visualize Confidence Intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a31af857",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:04.117096Z",
     "iopub.status.busy": "2026-05-17T00:42:04.116800Z",
     "iopub.status.idle": "2026-05-17T00:42:04.220172Z",
     "shell.execute_reply": "2026-05-17T00:42:04.219717Z"
    }
   },
   "outputs": [],
   "source": [
    "feature_names = [f\"X{i}\" for i in range(n_features)]\n",
    "\n",
    "confidence_interval_plot(\n",
    "    ci,\n",
    "    feature_names=feature_names,\n",
    "    show=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b1529ef",
   "metadata": {},
   "source": [
    "## One-Sided Tests\n",
    "\n",
    "For feature importance, we often care if a feature has **positive** importance. Use `alternative=\"greater\"`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7377545b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:04.221797Z",
     "iopub.status.busy": "2026-05-17T00:42:04.221677Z",
     "iopub.status.idle": "2026-05-17T00:42:04.226453Z",
     "shell.execute_reply": "2026-05-17T00:42:04.226070Z"
    }
   },
   "outputs": [],
   "source": [
    "# One-sided test: H0: phi <= 0 vs H1: phi > 0\n",
    "ci_greater = explainer.conf_int(alpha=0.05, alternative=\"greater\")\n",
    "\n",
    "print(\"One-sided test (phi > 0):\")\n",
    "print(\"-\" * 60)\n",
    "print(f\"{'Feature':>8} {'Estimate':>10} {'CI Lower':>10} {'P-value':>10} {'Significant':>12}\")\n",
    "print(\"-\" * 60)\n",
    "for i in range(n_features):\n",
    "    sig = \"Yes\" if ci_greater[\"reject_null\"][i] else \"No\"\n",
    "    print(f\"{i:>8} {ci_greater['score'][i]:>10.4f} \"\n",
    "          f\"{ci_greater['ci_lower'][i]:>10.4f} {ci_greater['pvalue'][i]:>10.4f} {sig:>12}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "685f2536",
   "metadata": {},
   "source": [
    "## Variance Floor\n",
    "\n",
    "When some features have very small variance in their importance estimates, confidence intervals can become too narrow. The **variance floor** adds a minimum standard error.\n",
    "\n",
    "Two methods are available:\n",
    "- `fixed`: Use a constant floor value\n",
    "- `mixture`: Fit a two-component mixture to estimate the floor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50f51917",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:04.228536Z",
     "iopub.status.busy": "2026-05-17T00:42:04.228399Z",
     "iopub.status.idle": "2026-05-17T00:42:04.234505Z",
     "shell.execute_reply": "2026-05-17T00:42:04.234144Z"
    }
   },
   "outputs": [],
   "source": [
    "# Without variance floor\n",
    "ci_no_floor = explainer.conf_int(alpha=0.05, var_floor_c=0)\n",
    "\n",
    "# With fixed variance floor\n",
    "ci_fixed = explainer.conf_int(alpha=0.05, var_floor_method=\"fixed\", var_floor_c=0.1)\n",
    "\n",
    "# With mixture-based floor\n",
    "ci_mixture = explainer.conf_int(alpha=0.05, var_floor_method=\"mixture\", var_floor_quantile=0.95)\n",
    "\n",
    "print(\"Standard errors comparison:\")\n",
    "print(\"-\" * 55)\n",
    "print(f\"{'Feature':>8} {'No Floor':>12} {'Fixed':>12} {'Mixture':>12}\")\n",
    "print(\"-\" * 55)\n",
    "for i in range(n_features):\n",
    "    print(f\"{i:>8} {ci_no_floor['se'][i]:>12.4f} {ci_fixed['se'][i]:>12.4f} {ci_mixture['se'][i]:>12.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d77b0ce",
   "metadata": {},
   "source": [
    "## Practical Significance Margin\n",
    "\n",
    "Instead of testing $H_0: \\phi = 0$, you can test against a **practical threshold** $\\delta$:\n",
    "\n",
    "$$H_0: \\phi \\leq \\delta \\quad \\text{vs} \\quad H_1: \\phi > \\delta$$\n",
    "\n",
    "This identifies features that are not just statistically different from zero, but also practically meaningful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6839f90",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:04.235970Z",
     "iopub.status.busy": "2026-05-17T00:42:04.235875Z",
     "iopub.status.idle": "2026-05-17T00:42:04.240018Z",
     "shell.execute_reply": "2026-05-17T00:42:04.239627Z"
    }
   },
   "outputs": [],
   "source": [
    "# Test with practical margin of 0.5\n",
    "margin = 0.5\n",
    "ci_margin = explainer.conf_int(\n",
    "    alpha=0.05, \n",
    "    alternative=\"greater\",\n",
    "    margin=margin\n",
    ")\n",
    "\n",
    "print(f\"Testing H0: phi <= {margin}\")\n",
    "print(\"-\" * 50)\n",
    "print(f\"{'Feature':>8} {'Estimate':>10} {'P-value':>10} {'Significant':>12}\")\n",
    "print(\"-\" * 50)\n",
    "for i in range(n_features):\n",
    "    sig = \"Yes\" if ci_margin[\"reject_null\"][i] else \"No\"\n",
    "    print(f\"{i:>8} {ci_margin['score'][i]:>10.4f} {ci_margin['pvalue'][i]:>10.4f} {sig:>12}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e159ea1",
   "metadata": {},
   "source": [
    "## Automatic Margin via Mixture Model\n",
    "\n",
    "Use `margin_method=\"mixture\"` to automatically estimate a practical margin:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c949b1f4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:04.241434Z",
     "iopub.status.busy": "2026-05-17T00:42:04.241338Z",
     "iopub.status.idle": "2026-05-17T00:42:04.246882Z",
     "shell.execute_reply": "2026-05-17T00:42:04.246429Z"
    }
   },
   "outputs": [],
   "source": [
    "ci_auto_margin = explainer.conf_int(\n",
    "    alpha=0.05,\n",
    "    alternative=\"greater\",\n",
    "    margin_method=\"mixture\",\n",
    "    margin_quantile=0.95,\n",
    ")\n",
    "\n",
    "print(f\"Automatically selected margin: {ci_auto_margin['margin']:.4f}\")\n",
    "print(f\"Significant features: {np.where(ci_auto_margin['reject_null'])[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f318952",
   "metadata": {},
   "source": [
    "## Feature Selection with Statistical Guarantees\n",
    "\n",
    "Use the confidence intervals to select features with controlled false discovery:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "139638b5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:04.248395Z",
     "iopub.status.busy": "2026-05-17T00:42:04.248265Z",
     "iopub.status.idle": "2026-05-17T00:42:04.259416Z",
     "shell.execute_reply": "2026-05-17T00:42:04.259015Z"
    }
   },
   "outputs": [],
   "source": [
    "def statistical_feature_selection(explainer, X_test, alpha=0.05, margin=0.0):\n",
    "    \"\"\"Select features with statistical significance.\"\"\"\n",
    "    # Compute importance\n",
    "    results = explainer(X_test)\n",
    "    \n",
    "    # Get confidence intervals\n",
    "    ci = explainer.conf_int(\n",
    "        alpha=alpha,\n",
    "        alternative=\"greater\",\n",
    "        margin=margin,\n",
    "        var_floor_method=\"mixture\",\n",
    "    )\n",
    "    \n",
    "    # Select significant features\n",
    "    selected = np.where(ci[\"reject_null\"])[0]\n",
    "    \n",
    "    # Sort by importance\n",
    "    sorted_idx = np.argsort(ci[\"score\"][selected])[::-1]\n",
    "    \n",
    "    return selected[sorted_idx], ci\n",
    "\n",
    "# Run feature selection\n",
    "selected_features, ci_result = statistical_feature_selection(\n",
    "    explainer, X_test, alpha=0.05, margin=0.0\n",
    ")\n",
    "\n",
    "print(\"Selected Features (sorted by importance):\")\n",
    "print(\"-\" * 40)\n",
    "for i, feat in enumerate(selected_features):\n",
    "    print(f\"  {i+1}. Feature {feat} (importance = {ci_result['score'][feat]:.4f})\")\n",
    "\n",
    "print(f\"\\nTrue important features: 0, 1, 2\")\n",
    "print(f\"Correctly identified: {set(selected_features) & {0, 1, 2}}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c94e254",
   "metadata": {},
   "source": [
    "## The `summary()` Method\n",
    "\n",
    "For a quick formatted view, use the built-in `summary()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42d7fc37",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:04.261141Z",
     "iopub.status.busy": "2026-05-17T00:42:04.261011Z",
     "iopub.status.idle": "2026-05-17T00:42:04.265680Z",
     "shell.execute_reply": "2026-05-17T00:42:04.265278Z"
    }
   },
   "outputs": [],
   "source": [
    "# Print formatted summary\n",
    "explainer.summary(\n",
    "    alpha=0.05,\n",
    "    alternative=\"greater\",\n",
    "    var_floor_method=\"mixture\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e67b2d1c",
   "metadata": {},
   "source": [
    "## Multiple Testing Correction\n",
    "\n",
    "When explaining models with many features, it is important to control the false discovery rate (FDR). You can specify a `multitest_method` in both `conf_int()` and `summary()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8d3c5e4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-05-17T00:42:04.267100Z",
     "iopub.status.busy": "2026-05-17T00:42:04.267003Z",
     "iopub.status.idle": "2026-05-17T00:42:04.279563Z",
     "shell.execute_reply": "2026-05-17T00:42:04.279032Z"
    }
   },
   "outputs": [],
   "source": [
    "# Summary with FDR control (Benjamini-Hochberg)\n",
    "explainer.summary(\n",
    "    alpha=0.05,\n",
    "    alternative=\"greater\",\n",
    "    multitest_method=\"fdr_bh\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b40341e9",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ee1b7ff-ef39-446f-9c56-52f8499f3e27",
   "metadata": {},
   "source": [
    "\n",
    "Key takeaways:\n",
    "\n",
    "1. `conf_int()` provides confidence intervals and p-values for feature importance\n",
    "2. Use `alternative=\"greater\"` for one-sided tests of positive importance\n",
    "3. Variance floor (`var_floor_method`) stabilizes inference for small effects\n",
    "4. Practical margin (`margin`) tests against meaningful thresholds\n",
    "5. Use `summary()` for quick formatted output"
   ]
  }
 ],
 "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.10.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}