{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# P-value Function Diagnostics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## INTRO & SETTINGS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `p-values` calibration method fits a monotonic probabilistic classifier that estimates the\n", "rejection probability $P(T \\le \\tau \\mid \\theta)$ -- i.e. a p-value/CDF function of the test statistic,\n", "conditional on $\\theta$. Notebooks 2 and 4 use this fitted function to build confidence sets, but never\n", "ask *how good the fit itself is*.\n", "\n", "This notebook is about that second question: diagnosing the fit of the p-value function, independent\n", "of the confidence sets built from it. We reuse the same model as notebook 2 (`Posterior` +\n", "`calibration_method='p-values'` on `GaussianMean`) so the diagnostics here are directly comparable to\n", "that notebook's construction." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SETTINGS\n", "\n", "LIKELIHOOD_COV = 0.01\n", "PRIOR_LOC = 0\n", "PRIOR_COV = 0.1\n", "\n", "PARAM_DIM = 2\n", "DATA_DIM = 2\n", "BATCH_SIZE = 1\n", "PARAM_SPACE_BOUNDS = {'low': -1.5, 'high': 1.5}\n", "PARAM_GRID_SIZE = 1_000\n", "\n", "CONFIDENCE_LEVEL = 0.90\n", "\n", "B = 20_000\n", "B_PRIME = 10_000\n", "MONTE_CARLO_SIZE = 2_000 # MC draws per grid point for the diagnostics themselves" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## SIMULATE" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import torch\n", "\n", "from lf2i.simulator.gaussian import GaussianMean\n", "\n", "gm = GaussianMean(\n", " likelihood_cov=LIKELIHOOD_COV,\n", " prior='gaussian',\n", " prior_kwargs={'loc': PRIOR_LOC, 'cov': PRIOR_COV},\n", " poi_space_bounds=PARAM_SPACE_BOUNDS,\n", " poi_grid_size=PARAM_GRID_SIZE,\n", " poi_dim=PARAM_DIM,\n", " data_dim=DATA_DIM,\n", " batch_size=BATCH_SIZE,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## FIT THE P-VALUE FUNCTION" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from lf2i.inference import LF2I\n", "from lf2i.test_statistics import Posterior\n", "from sbi.inference import SNPE\n", "\n", "posterior_ts = Posterior(poi_dim=PARAM_DIM, estimator=SNPE())\n", "lf2i = LF2I(test_statistic=posterior_ts)\n", "\n", "x_obs = gm(param=torch.tensor([[0.5, -0.3]])).reshape(1, DATA_DIM)\n", "\n", "_ = lf2i.inference(\n", " x=x_obs,\n", " evaluation_grid=gm.poi_grid,\n", " confidence_level=CONFIDENCE_LEVEL,\n", " calibration_method='p-values',\n", " calibration_model='nn',\n", " simulator=gm,\n", " b=B,\n", " b_prime=B_PRIME,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## DIAGNOSE THE FIT\n", "\n", "For each $\\theta$ on the evaluation grid, `monte_carlo_pvalue_diagnostics` draws fresh Monte Carlo\n", "samples, evaluates the test statistic, and compares the fitted p-value/CDF model's predictions to the\n", "*empirical* CDF at that $\\theta$ -- returning, per grid point, a normalized CRPS (0 = perfect fit,\n", "1 = no better than ignoring $\\theta$ and using the marginal CDF) and pinball losses at a range of\n", "quantile levels." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from lf2i.diagnostics.monte_carlo_methods import monte_carlo_pvalue_diagnostics\n", "\n", "evaluation_grid_out, estimation_errors = monte_carlo_pvalue_diagnostics(\n", " test_statistic=posterior_ts,\n", " calibration_model=lf2i.calibration_model,\n", " simulator=gm,\n", " evaluation_grid=gm.poi_grid,\n", " monte_carlo_size=MONTE_CARLO_SIZE,\n", ")\n", "list(estimation_errors.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Calibration score heatmap" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from lf2i.plot.calibration_diagnostics import calibration_score_plot\n", "\n", "calibration_score_plot(\n", " parameters=gm.poi_grid.numpy(),\n", " scores=estimation_errors['crps'],\n", " score_label='CRPS (normalized)',\n", " param_dim=PARAM_DIM,\n", " title='Normalized CRPS of the fitted p-value function across the parameter grid',\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### CDF comparison at a single theta" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from lf2i.plot.calibration_diagnostics import plot_cdf_comparison\n", "\n", "plot_cdf_comparison(\n", " test_statistic=posterior_ts,\n", " calib_model=lf2i.calibration_model[f'{CONFIDENCE_LEVEL:.2f}'],\n", " theta_eval=torch.tensor([[0.5, -0.3]]),\n", " simulator=gm,\n", " monte_carlo_size=MONTE_CARLO_SIZE,\n", " title='Fitted vs. empirical CDF at theta = (0.5, -0.3)',\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Combined diagnostic panel" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from lf2i.plot.calibration_diagnostics import calibration_cdf_panel\n", "\n", "calibration_cdf_panel(\n", " evaluation_grid=gm.poi_grid,\n", " estimation_errors=estimation_errors,\n", " test_statistic=posterior_ts,\n", " calibration_model=lf2i.calibration_model,\n", " simulator=gm,\n", " param_dim=PARAM_DIM,\n", " score_key='crps',\n", " monte_carlo_size=MONTE_CARLO_SIZE,\n", ")" ] } ], "metadata": { "kernelspec": { "display_name": "tsi", "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.11.9" } }, "nbformat": 4, "nbformat_minor": 5 }