From b18ec29afe7f54798a86f330bfd344d3e567d35b Mon Sep 17 00:00:00 2001 From: Michael-Howes Date: Fri, 20 Mar 2026 14:45:03 -0700 Subject: [PATCH 1/3] fix test_mean.py - Fix broadcasting error in test_ppi_mean_multid. - Change lambd_optim_mode to lam_optim_mode in test_ppi_mean_elem. - Change confromal baseline to use np.inf instead of np.infty. --- ppi_py/baselines.py | 2 +- tests/test_mean.py | 17 ++++++++++------- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/ppi_py/baselines.py b/ppi_py/baselines.py index 5e0c17c..799b7df 100644 --- a/ppi_py/baselines.py +++ b/ppi_py/baselines.py @@ -144,7 +144,7 @@ def conformal_mean_ci(Y, Yhat, Yhat_unlabeled, alpha=0.1, bonferroni=True): else (1 - alpha) * (1 + 1 / n) ) if level >= 1: - return -np.infty, np.infty + return -np.inf, np.inf conformal_quantile = np.quantile(scores, level, method="higher") imputed_estimate = Yhat_unlabeled.mean() return ( diff --git a/tests/test_mean.py b/tests/test_mean.py index 784135f..2f51998 100644 --- a/tests/test_mean.py +++ b/tests/test_mean.py @@ -46,7 +46,10 @@ def test_ppi_mean_multid(): included = (ci[0] <= 0) & (ci[1] >= 0) includeds[j] += included.astype(int) - failed = np.any(includeds / trials < 1 - alphas - epsilon) + print(includeds/trials) + failures = (includeds / trials).T < 1 - alphas - epsilon + print(failures) + failed = np.any(failures) assert not failed @@ -56,21 +59,21 @@ def test_ppi_mean_elem(): Yhat = np.random.normal(-2, 1, 10000) Yhat_unlabeled = np.random.normal(-2, 1, 10000) - ppi_mean_pointestimate(Y, Yhat, Yhat_unlabeled, lambd_optim_mode="element") + ppi_mean_pointestimate(Y, Yhat, Yhat_unlabeled, lam_optim_mode="element") ppi_mean_ci( - Y, Yhat, Yhat_unlabeled, alpha=alpha, lambd_optim_mode="element" + Y, Yhat, Yhat_unlabeled, alpha=alpha, lam_optim_mode="element" ) - ppi_mean_pval(Y, Yhat, Yhat_unlabeled, lambd_optim_mode="element") + ppi_mean_pval(Y, Yhat, Yhat_unlabeled, lam_optim_mode="element") Y = np.random.normal(0, 1, (10000, 5)) Yhat = np.random.normal(-2, 1, (10000, 5)) Yhat_unlabeled = np.random.normal(-2, 1, (10000, 5)) - ppi_mean_pointestimate(Y, Yhat, Yhat_unlabeled, lambd_optim_mode="element") + ppi_mean_pointestimate(Y, Yhat, Yhat_unlabeled, lam_optim_mode="element") ppi_mean_ci( - Y, Yhat, Yhat_unlabeled, alpha=alpha, lambd_optim_mode="element" + Y, Yhat, Yhat_unlabeled, alpha=alpha, lam_optim_mode="element" ) - ppi_mean_pval(Y, Yhat, Yhat_unlabeled, lambd_optim_mode="element") + ppi_mean_pval(Y, Yhat, Yhat_unlabeled, lam_optim_mode="element") def test_ppi_mean_pval(): From 1f2f63f28d7adcfb60a16910fa7c26ca78b6a8b6 Mon Sep 17 00:00:00 2001 From: Michael-Howes Date: Fri, 20 Mar 2026 16:20:06 -0700 Subject: [PATCH 2/3] fix test_ppi_logistic_pointestimate_debias --- tests/test_cross.py | 4 ++-- tests/test_logistic.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_cross.py b/tests/test_cross.py index 43ab3fd..5c0ddce 100644 --- a/tests/test_cross.py +++ b/tests/test_cross.py @@ -117,7 +117,7 @@ def test_crossppi_logistic_pointestimate_debias(): X = np.random.randn(n, d) beta = np.random.randn(d) beta_prediction = beta + np.random.randn(d) + 2 - Y = expit(X.dot(beta) + np.random.randn(n)) + Y = np.random.binomial(1, expit(X.dot(beta))) Yhat = expit(X.dot(beta_prediction) + np.random.randn(n)) # Make a synthetic unlabeled data set with predictions Yhat X_unlabeled = np.random.randn(N, d) @@ -127,7 +127,7 @@ def test_crossppi_logistic_pointestimate_debias(): # Compute the point estimate beta_ppi_pointestimate = crossppi_logistic_pointestimate( X, - (Y > 0.5).astype(int), + Y, Yhat, X_unlabeled, Yhat_unlabeled, diff --git a/tests/test_logistic.py b/tests/test_logistic.py index 00c977c..ad1b65e 100644 --- a/tests/test_logistic.py +++ b/tests/test_logistic.py @@ -17,7 +17,7 @@ def test_ppi_logistic_pointestimate_debias(): X = np.random.randn(n, d) beta = np.random.randn(d) beta_prediction = beta + np.random.randn(d) + 2 - Y = expit(X.dot(beta) + np.random.randn(n)) + Y = np.random.binomial(1, expit(X.dot(beta))) Yhat = expit(X.dot(beta_prediction) + np.random.randn(n)) # Make a synthetic unlabeled data set with predictions Yhat X_unlabeled = np.random.randn(N, d) @@ -27,7 +27,7 @@ def test_ppi_logistic_pointestimate_debias(): # Compute the point estimate beta_ppi_pointestimate = ppi_logistic_pointestimate( X, - (Y > 0.5).astype(int), + Y, Yhat, X_unlabeled, Yhat_unlabeled, From c991e38bb39eaddc7abb4019f0bb377dc827958c Mon Sep 17 00:00:00 2001 From: Michael-Howes Date: Fri, 20 Mar 2026 16:23:52 -0700 Subject: [PATCH 3/3] format with black --- examples/tree_cover_ptd.ipynb | 602 ++++++++++++++++++++++------------ tests/test_mean.py | 10 +- 2 files changed, 403 insertions(+), 209 deletions(-) diff --git a/examples/tree_cover_ptd.ipynb b/examples/tree_cover_ptd.ipynb index afadc12..fd3cbd5 100644 --- a/examples/tree_cover_ptd.ipynb +++ b/examples/tree_cover_ptd.ipynb @@ -28,6 +28,7 @@ "outputs": [], "source": [ "import os, sys\n", + "\n", "sys.path.append(os.path.abspath(os.path.join(os.getcwd(), os.pardir)))\n", "\n", "import numpy as np\n", @@ -61,6 +62,7 @@ "from statsmodels.genmod.families import Binomial\n", "from statsmodels.genmod.families.links import Logit\n", "\n", + "\n", "def classical_logistic_regression_ci(X, Y, w=None, alpha=0.05):\n", " \"\"\"\n", " Computes confidence intervals for logistic regression coefficients using the classical method.\n", @@ -74,10 +76,13 @@ " Returns:\n", " tuple: lower and upper bounds of classical confidence intervals for the coefficients\n", " \"\"\"\n", - " regression = GLM(endog=Y, exog=X, freq_weights=w, family=Binomial(link=Logit())).fit()\n", + " regression = GLM(\n", + " endog=Y, exog=X, freq_weights=w, family=Binomial(link=Logit())\n", + " ).fit()\n", " ci = regression.conf_int(alpha=alpha).T\n", " return ci\n", "\n", + "\n", "def classical_linear_regression_ci(X, Y, w=None, alpha=0.05):\n", " \"\"\"\n", " Computes confidence intervals for linear regression coefficients using the classical method.\n", @@ -361,11 +366,11 @@ "metadata": {}, "outputs": [], "source": [ - "truth_Y = np.array(data['truth_tree']).reshape(-1, 1)\n", - "preds_Y = np.array(data['preds_tree']).reshape(-1, 1)\n", + "truth_Y = np.array(data[\"truth_tree\"]).reshape(-1, 1)\n", + "preds_Y = np.array(data[\"preds_tree\"]).reshape(-1, 1)\n", "\n", - "truth_X = np.array(data[['truth_elevation', 'truth_population']])\n", - "preds_X = np.array(data[['truth_elevation', 'preds_population']])" + "truth_X = np.array(data[[\"truth_elevation\", \"truth_population\"]])\n", + "preds_X = np.array(data[[\"truth_elevation\", \"preds_population\"]])" ] }, { @@ -388,13 +393,19 @@ "outputs": [], "source": [ "np.random.seed(seed=100)\n", - "calibration_indices = np.random.choice(np.arange(0, len(data)), size=500, replace=False)\n", + "calibration_indices = np.random.choice(\n", + " np.arange(0, len(data)), size=500, replace=False\n", + ")\n", "X = statsmodels.tools.add_constant(truth_X[calibration_indices])\n", "Xhat = statsmodels.tools.add_constant(preds_X[calibration_indices])\n", - "Xhat_unlabeled = statsmodels.tools.add_constant(np.delete(preds_X, calibration_indices, axis=0)) # all predicted datapoints except calibration indices\n", + "Xhat_unlabeled = statsmodels.tools.add_constant(\n", + " np.delete(preds_X, calibration_indices, axis=0)\n", + ") # all predicted datapoints except calibration indices\n", "Y = truth_Y[calibration_indices]\n", "Yhat = preds_Y[calibration_indices]\n", - "Yhat_unlabeled = np.delete(preds_Y, calibration_indices, axis=0) # all predicted datapoints except calibration indices" + "Yhat_unlabeled = np.delete(\n", + " preds_Y, calibration_indices, axis=0\n", + ") # all predicted datapoints except calibration indices" ] }, { @@ -420,16 +431,29 @@ "metadata": {}, "outputs": [], "source": [ - "true_coeff = ptd.algorithm_linear_regression(data=[statsmodels.tools.add_constant(truth_X), truth_Y], w=None)\n", - "\n", - "tuning_matrix, ptd_pointestimate, ptd_ci = ptd.ptd_linear_regression(X, Xhat, Xhat_unlabeled, Y, Yhat, Yhat_unlabeled, \n", - " B=2000, alpha=0.05, tuning_method='optimal')\n", + "true_coeff = ptd.algorithm_linear_regression(\n", + " data=[statsmodels.tools.add_constant(truth_X), truth_Y], w=None\n", + ")\n", + "\n", + "tuning_matrix, ptd_pointestimate, ptd_ci = ptd.ptd_linear_regression(\n", + " X,\n", + " Xhat,\n", + " Xhat_unlabeled,\n", + " Y,\n", + " Yhat,\n", + " Yhat_unlabeled,\n", + " B=2000,\n", + " alpha=0.05,\n", + " tuning_method=\"optimal\",\n", + ")\n", "\n", "classical_ci = classical_linear_regression_ci(X, Y, alpha=0.05)\n", - "classical_pointestimate = (classical_ci[0]+classical_ci[1])/2\n", + "classical_pointestimate = (classical_ci[0] + classical_ci[1]) / 2\n", "\n", - "naive_ci = classical_linear_regression_ci(statsmodels.tools.add_constant(preds_X), preds_Y, alpha=0.05)\n", - "naive_pointestimate = (naive_ci[0]+naive_ci[1])/2" + "naive_ci = classical_linear_regression_ci(\n", + " statsmodels.tools.add_constant(preds_X), preds_Y, alpha=0.05\n", + ")\n", + "naive_pointestimate = (naive_ci[0] + naive_ci[1]) / 2" ] }, { @@ -464,57 +488,83 @@ "source": [ "plt.rcParams[\"font.sans-serif\"] = \"Arial\"\n", "\n", - "plt.rc('font', size=43) \n", - "plt.rc('axes', titlesize=43) \n", - "plt.rc('axes', labelsize=43) \n", - "plt.rc('xtick', labelsize=43) \n", - "plt.rc('ytick', labelsize=43) \n", - "plt.rc('figure', titlesize=43)\n", + "plt.rc(\"font\", size=43)\n", + "plt.rc(\"axes\", titlesize=43)\n", + "plt.rc(\"axes\", labelsize=43)\n", + "plt.rc(\"xtick\", labelsize=43)\n", + "plt.rc(\"ytick\", labelsize=43)\n", + "plt.rc(\"figure\", titlesize=43)\n", "\n", - "covariates = ['Intercept', 'Elevation', 'Population']\n", + "covariates = [\"Intercept\", \"Elevation\", \"Population\"]\n", "pd.set_option(\"display.precision\", 8)\n", "\n", "fig = plt.figure(figsize=(32, 6.5))\n", - "colors = ['blue', 'green', 'red']\n", + "colors = [\"blue\", \"green\", \"red\"]\n", "for i, covariate in enumerate(covariates):\n", - " \n", + "\n", " data_dict = {}\n", " classical_ci_width = classical_ci[1][i] - classical_ci[0][i]\n", " ptd_ci_width = ptd_ci[1][i] - ptd_ci[0][i]\n", - " data_dict['category'] = ['PTD','Classical','Naive']\n", - " data_dict['lower'] = [ptd_ci[0][i], classical_ci[0][i], naive_ci[0][i]]\n", - " data_dict['upper'] = [ptd_ci[1][i], classical_ci[1][i], naive_ci[1][i]]\n", - " data_dict['pointestimate'] = [ptd_pointestimate[i], classical_pointestimate[i], naive_pointestimate[i]]\n", + " data_dict[\"category\"] = [\"PTD\", \"Classical\", \"Naive\"]\n", + " data_dict[\"lower\"] = [ptd_ci[0][i], classical_ci[0][i], naive_ci[0][i]]\n", + " data_dict[\"upper\"] = [ptd_ci[1][i], classical_ci[1][i], naive_ci[1][i]]\n", + " data_dict[\"pointestimate\"] = [\n", + " ptd_pointestimate[i],\n", + " classical_pointestimate[i],\n", + " naive_pointestimate[i],\n", + " ]\n", " dataset = pd.DataFrame(data_dict)\n", - " \n", - " subplot = fig.add_subplot(141+i)\n", - " subplot.spines['top'].set_visible(False)\n", - " subplot.spines['right'].set_visible(False)\n", - " subplot.spines['left'].set_visible(False)\n", - " subplot.spines['bottom'].set_linewidth(2)\n", - " \n", + "\n", + " subplot = fig.add_subplot(141 + i)\n", + " subplot.spines[\"top\"].set_visible(False)\n", + " subplot.spines[\"right\"].set_visible(False)\n", + " subplot.spines[\"left\"].set_visible(False)\n", + " subplot.spines[\"bottom\"].set_linewidth(2)\n", + "\n", " # show true coefficient (yellow dotted line)\n", - " plt.axvline(true_coeff[i], color='tab:olive', linestyle='--', lw=4, label='True \\ncoefficient')\n", - " \n", + " plt.axvline(\n", + " true_coeff[i],\n", + " color=\"tab:olive\",\n", + " linestyle=\"--\",\n", + " lw=4,\n", + " label=\"True \\ncoefficient\",\n", + " )\n", + "\n", " # plot confidence intervals\n", - " for lower,upper,pointestimate,y in zip(dataset['lower'],dataset['upper'],dataset['pointestimate'],range(len(dataset))):\n", + " for lower, upper, pointestimate, y in zip(\n", + " dataset[\"lower\"],\n", + " dataset[\"upper\"],\n", + " dataset[\"pointestimate\"],\n", + " range(len(dataset)),\n", + " ):\n", " subplot.scatter([pointestimate], [y], color=colors[y], s=400)\n", - " subplot.scatter([lower, upper], [y, y], color=colors[y], marker='|', s=400, lw=5)\n", - " subplot.plot((lower,upper),(y,y),color=colors[y], lw=6)\n", + " subplot.scatter(\n", + " [lower, upper], [y, y], color=colors[y], marker=\"|\", s=400, lw=5\n", + " )\n", + " subplot.plot((lower, upper), (y, y), color=colors[y], lw=6)\n", " if y == 0:\n", " # compute PTD effective sample size improvement over classical method\n", - " ptd_effective_n_improvement = np.round((classical_ci_width/ptd_ci_width)**2, 1)\n", - " subplot.text((lower+upper)/2, 0.2, f'{ptd_effective_n_improvement}x', fontsize = 40, c='blue', weight='bold')\n", - " \n", - " subplot.tick_params(axis=u'both', which=u'both',length=0)\n", + " ptd_effective_n_improvement = np.round(\n", + " (classical_ci_width / ptd_ci_width) ** 2, 1\n", + " )\n", + " subplot.text(\n", + " (lower + upper) / 2,\n", + " 0.2,\n", + " f\"{ptd_effective_n_improvement}x\",\n", + " fontsize=40,\n", + " c=\"blue\",\n", + " weight=\"bold\",\n", + " )\n", + "\n", + " subplot.tick_params(axis=\"both\", which=\"both\", length=0)\n", " if i == 0:\n", - " subplot.set_yticks(range(len(dataset)),list(dataset['category']))\n", + " subplot.set_yticks(range(len(dataset)), list(dataset[\"category\"]))\n", " else:\n", " subplot.set_yticks([])\n", - " subplot.set_xlabel(f'{covariate}')\n", + " subplot.set_xlabel(f\"{covariate}\")\n", " subplot.set_ylim(-0.2)\n", - " \n", - "plt.tight_layout() \n", + "\n", + "plt.tight_layout()\n", "plt.show()" ] }, @@ -543,11 +593,11 @@ "metadata": {}, "outputs": [], "source": [ - "truth_Y = np.array(data['truth_tree'] > 10, dtype=float).reshape(-1, 1)\n", - "preds_Y = np.array(data['preds_tree'] > 10, dtype=float).reshape(-1, 1)\n", + "truth_Y = np.array(data[\"truth_tree\"] > 10, dtype=float).reshape(-1, 1)\n", + "preds_Y = np.array(data[\"preds_tree\"] > 10, dtype=float).reshape(-1, 1)\n", "\n", - "truth_X = np.array(data[['truth_elevation', 'truth_population']])\n", - "preds_X = np.array(data[['truth_elevation', 'preds_population']])" + "truth_X = np.array(data[[\"truth_elevation\", \"truth_population\"]])\n", + "preds_X = np.array(data[[\"truth_elevation\", \"preds_population\"]])" ] }, { @@ -566,13 +616,19 @@ "outputs": [], "source": [ "np.random.seed(seed=100)\n", - "calibration_indices = np.random.choice(np.arange(0, len(data)), size=1000, replace=False)\n", + "calibration_indices = np.random.choice(\n", + " np.arange(0, len(data)), size=1000, replace=False\n", + ")\n", "X = statsmodels.tools.add_constant(truth_X[calibration_indices])\n", "Xhat = statsmodels.tools.add_constant(preds_X[calibration_indices])\n", - "Xhat_unlabeled = statsmodels.tools.add_constant(np.delete(preds_X, calibration_indices, axis=0)) # all predicted datapoints except calibration indices\n", + "Xhat_unlabeled = statsmodels.tools.add_constant(\n", + " np.delete(preds_X, calibration_indices, axis=0)\n", + ") # all predicted datapoints except calibration indices\n", "Y = truth_Y[calibration_indices]\n", "Yhat = preds_Y[calibration_indices]\n", - "Yhat_unlabeled = np.delete(preds_Y, calibration_indices, axis=0) # all predicted datapoints except calibration indices" + "Yhat_unlabeled = np.delete(\n", + " preds_Y, calibration_indices, axis=0\n", + ") # all predicted datapoints except calibration indices" ] }, { @@ -598,16 +654,29 @@ "metadata": {}, "outputs": [], "source": [ - "true_coeff = ptd.algorithm_logistic_regression(data=[statsmodels.tools.add_constant(truth_X), truth_Y], w=None)\n", - "\n", - "tuning_matrix, ptd_pointestimate, ptd_ci = ptd.ptd_logistic_regression(X, Xhat, Xhat_unlabeled, Y, Yhat, Yhat_unlabeled, \n", - " B=2000, alpha=0.05, tuning_method='optimal')\n", + "true_coeff = ptd.algorithm_logistic_regression(\n", + " data=[statsmodels.tools.add_constant(truth_X), truth_Y], w=None\n", + ")\n", + "\n", + "tuning_matrix, ptd_pointestimate, ptd_ci = ptd.ptd_logistic_regression(\n", + " X,\n", + " Xhat,\n", + " Xhat_unlabeled,\n", + " Y,\n", + " Yhat,\n", + " Yhat_unlabeled,\n", + " B=2000,\n", + " alpha=0.05,\n", + " tuning_method=\"optimal\",\n", + ")\n", "\n", "classical_ci = classical_logistic_regression_ci(X, Y, alpha=0.05)\n", - "classical_pointestimate = (classical_ci[0]+classical_ci[1])/2\n", + "classical_pointestimate = (classical_ci[0] + classical_ci[1]) / 2\n", "\n", - "naive_ci = classical_logistic_regression_ci(statsmodels.tools.add_constant(preds_X), preds_Y, alpha=0.05)\n", - "naive_pointestimate = (naive_ci[0]+naive_ci[1])/2" + "naive_ci = classical_logistic_regression_ci(\n", + " statsmodels.tools.add_constant(preds_X), preds_Y, alpha=0.05\n", + ")\n", + "naive_pointestimate = (naive_ci[0] + naive_ci[1]) / 2" ] }, { @@ -642,57 +711,83 @@ "source": [ "plt.rcParams[\"font.sans-serif\"] = \"Arial\"\n", "\n", - "plt.rc('font', size=43) \n", - "plt.rc('axes', titlesize=43) \n", - "plt.rc('axes', labelsize=43) \n", - "plt.rc('xtick', labelsize=43) \n", - "plt.rc('ytick', labelsize=43) \n", - "plt.rc('figure', titlesize=43)\n", + "plt.rc(\"font\", size=43)\n", + "plt.rc(\"axes\", titlesize=43)\n", + "plt.rc(\"axes\", labelsize=43)\n", + "plt.rc(\"xtick\", labelsize=43)\n", + "plt.rc(\"ytick\", labelsize=43)\n", + "plt.rc(\"figure\", titlesize=43)\n", "\n", - "covariates = ['Intercept', 'Elevation', 'Population']\n", + "covariates = [\"Intercept\", \"Elevation\", \"Population\"]\n", "pd.set_option(\"display.precision\", 8)\n", "\n", "fig = plt.figure(figsize=(32, 6.5))\n", - "colors = ['blue', 'green', 'red']\n", + "colors = [\"blue\", \"green\", \"red\"]\n", "for i, covariate in enumerate(covariates):\n", - " \n", + "\n", " data_dict = {}\n", " classical_ci_width = classical_ci[1][i] - classical_ci[0][i]\n", " ptd_ci_width = ptd_ci[1][i] - ptd_ci[0][i]\n", - " data_dict['category'] = ['PTD','Classical','Naive']\n", - " data_dict['lower'] = [ptd_ci[0][i], classical_ci[0][i], naive_ci[0][i]]\n", - " data_dict['upper'] = [ptd_ci[1][i], classical_ci[1][i], naive_ci[1][i]]\n", - " data_dict['pointestimate'] = [ptd_pointestimate[i], classical_pointestimate[i], naive_pointestimate[i]]\n", + " data_dict[\"category\"] = [\"PTD\", \"Classical\", \"Naive\"]\n", + " data_dict[\"lower\"] = [ptd_ci[0][i], classical_ci[0][i], naive_ci[0][i]]\n", + " data_dict[\"upper\"] = [ptd_ci[1][i], classical_ci[1][i], naive_ci[1][i]]\n", + " data_dict[\"pointestimate\"] = [\n", + " ptd_pointestimate[i],\n", + " classical_pointestimate[i],\n", + " naive_pointestimate[i],\n", + " ]\n", " dataset = pd.DataFrame(data_dict)\n", - " \n", - " subplot = fig.add_subplot(141+i)\n", - " subplot.spines['top'].set_visible(False)\n", - " subplot.spines['right'].set_visible(False)\n", - " subplot.spines['left'].set_visible(False)\n", - " subplot.spines['bottom'].set_linewidth(2)\n", - " \n", + "\n", + " subplot = fig.add_subplot(141 + i)\n", + " subplot.spines[\"top\"].set_visible(False)\n", + " subplot.spines[\"right\"].set_visible(False)\n", + " subplot.spines[\"left\"].set_visible(False)\n", + " subplot.spines[\"bottom\"].set_linewidth(2)\n", + "\n", " # show true coefficient (yellow dotted line)\n", - " plt.axvline(true_coeff[i], color='tab:olive', linestyle='--', lw=4, label='True \\ncoefficient')\n", - " \n", + " plt.axvline(\n", + " true_coeff[i],\n", + " color=\"tab:olive\",\n", + " linestyle=\"--\",\n", + " lw=4,\n", + " label=\"True \\ncoefficient\",\n", + " )\n", + "\n", " # plot confidence intervals\n", - " for lower,upper,pointestimate,y in zip(dataset['lower'],dataset['upper'],dataset['pointestimate'],range(len(dataset))):\n", + " for lower, upper, pointestimate, y in zip(\n", + " dataset[\"lower\"],\n", + " dataset[\"upper\"],\n", + " dataset[\"pointestimate\"],\n", + " range(len(dataset)),\n", + " ):\n", " subplot.scatter([pointestimate], [y], color=colors[y], s=400)\n", - " subplot.scatter([lower, upper], [y, y], color=colors[y], marker='|', s=400, lw=5)\n", - " subplot.plot((lower,upper),(y,y),color=colors[y], lw=6)\n", + " subplot.scatter(\n", + " [lower, upper], [y, y], color=colors[y], marker=\"|\", s=400, lw=5\n", + " )\n", + " subplot.plot((lower, upper), (y, y), color=colors[y], lw=6)\n", " if y == 0:\n", " # compute PTD effective sample size improvement over classical method\n", - " ptd_effective_n_improvement = np.round((classical_ci_width/ptd_ci_width)**2, 1)\n", - " subplot.text((lower+upper)/2, 0.2, f'{ptd_effective_n_improvement}x', fontsize = 40, c='blue', weight='bold')\n", - " \n", - " subplot.tick_params(axis=u'both', which=u'both',length=0)\n", + " ptd_effective_n_improvement = np.round(\n", + " (classical_ci_width / ptd_ci_width) ** 2, 1\n", + " )\n", + " subplot.text(\n", + " (lower + upper) / 2,\n", + " 0.2,\n", + " f\"{ptd_effective_n_improvement}x\",\n", + " fontsize=40,\n", + " c=\"blue\",\n", + " weight=\"bold\",\n", + " )\n", + "\n", + " subplot.tick_params(axis=\"both\", which=\"both\", length=0)\n", " if i == 0:\n", - " subplot.set_yticks(range(len(dataset)),list(dataset['category']))\n", + " subplot.set_yticks(range(len(dataset)), list(dataset[\"category\"]))\n", " else:\n", " subplot.set_yticks([])\n", - " subplot.set_xlabel(f'{covariate}')\n", + " subplot.set_xlabel(f\"{covariate}\")\n", " subplot.set_ylim(-0.2)\n", - " \n", - "plt.tight_layout() \n", + "\n", + "plt.tight_layout()\n", "plt.show()" ] }, @@ -743,34 +838,51 @@ "metadata": {}, "outputs": [], "source": [ - "def ptd_quantile(X, Xhat, Xhat_unlabeled, quantile, B=2000, alpha=0.05, tuning_method='optimal_diagonal'):\n", + "def ptd_quantile(\n", + " X,\n", + " Xhat,\n", + " Xhat_unlabeled,\n", + " quantile,\n", + " B=2000,\n", + " alpha=0.05,\n", + " tuning_method=\"optimal_diagonal\",\n", + "):\n", " \"\"\"\n", - " Computes tuning matrix, point estimates, and confidence intervals for quantile estimation using the Predict-then-Debias bootstrap algorithm. \n", - " \n", + " Computes tuning matrix, point estimates, and confidence intervals for quantile estimation using the Predict-then-Debias bootstrap algorithm.\n", + "\n", " Args:\n", " X (ndarray): ground truth values in labeled data (dimensions n x 1)\n", " Xhat (ndarray): predicted values in labeled data (dimensions n x 1)\n", " Xhat_unlabeled (ndarray): predicted values in unlabeled data (dimensions N x 1)\n", " quantile (scalar): the desired quantile probability. Must be in the range [0, 1].\n", " B (int, optional): number of bootstrap steps\n", - " alpha (float, optional): error level (must be in the range (0, 1)). The PTD confidence interval will target a coverage of 1 - alpha. \n", - " tuning_method (str, optional): method used to create the tuning matrix: \"optimal_diagonal\", \"optimal\", or None. (If tuning_method is None, the identity matrix is used.) \n", - " \n", + " alpha (float, optional): error level (must be in the range (0, 1)). The PTD confidence interval will target a coverage of 1 - alpha.\n", + " tuning_method (str, optional): method used to create the tuning matrix: \"optimal_diagonal\", \"optimal\", or None. (If tuning_method is None, the identity matrix is used.)\n", + "\n", " Returns:\n", " ndarray: the tuning matrix computed from the selected tuning method (1 x 1)\n", - " ndarray: PTD point estimate of the quantile \n", + " ndarray: PTD point estimate of the quantile\n", " tuple: lower and upper bounds of PTD confidence intervals with (1-alpha) coverage\n", " \"\"\"\n", + "\n", " def algorithm_quantile(data, weights):\n", " pointestimate = np.quantile(data, quantile)\n", " return np.array([pointestimate])\n", - " \n", + "\n", " # ptd_bootstrap requires input data arrays to be in list form (List[ndarray])\n", " data_truth = [X]\n", " data_pred = [Xhat]\n", " data_pred_unlabeled = [Xhat_unlabeled]\n", - " \n", - " return ptd.ptd_bootstrap(algorithm_quantile, data_truth, data_pred, data_pred_unlabeled, B=B, alpha=alpha, tuning_method=tuning_method)" + "\n", + " return ptd.ptd_bootstrap(\n", + " algorithm_quantile,\n", + " data_truth,\n", + " data_pred,\n", + " data_pred_unlabeled,\n", + " B=B,\n", + " alpha=alpha,\n", + " tuning_method=tuning_method,\n", + " )" ] }, { @@ -792,6 +904,7 @@ "source": [ "from scipy.stats.mstats import mquantiles_cimj\n", "\n", + "\n", "def classical_quantile_ci(X, quantile, alpha=0.05):\n", " ci = mquantiles_cimj(X, prob=quantile, alpha=alpha)\n", " return ci" @@ -812,8 +925,8 @@ "metadata": {}, "outputs": [], "source": [ - "truth_X = np.array(data['truth_tree']).reshape(-1, 1)\n", - "preds_X = np.array(data['preds_tree']).reshape(-1, 1)" + "truth_X = np.array(data[\"truth_tree\"]).reshape(-1, 1)\n", + "preds_X = np.array(data[\"preds_tree\"]).reshape(-1, 1)" ] }, { @@ -834,10 +947,14 @@ "outputs": [], "source": [ "np.random.seed(seed=100)\n", - "calibration_indices = np.random.choice(np.arange(0, len(data)), size=500, replace=False)\n", + "calibration_indices = np.random.choice(\n", + " np.arange(0, len(data)), size=500, replace=False\n", + ")\n", "X = truth_X[calibration_indices]\n", "Xhat = preds_X[calibration_indices]\n", - "Xhat_unlabeled = np.delete(preds_X, calibration_indices, axis=0) # all predicted datapoints except calibration indices" + "Xhat_unlabeled = np.delete(\n", + " preds_X, calibration_indices, axis=0\n", + ") # all predicted datapoints except calibration indices" ] }, { @@ -866,14 +983,21 @@ "quantile = 0.8\n", "true_coeff = [np.quantile(truth_X, quantile)]\n", "\n", - "tuning_matrix, ptd_pointestimate, ptd_ci = ptd_quantile(X, Xhat, Xhat_unlabeled, quantile, \n", - " B=2000, alpha=0.05, tuning_method='optimal')\n", + "tuning_matrix, ptd_pointestimate, ptd_ci = ptd_quantile(\n", + " X,\n", + " Xhat,\n", + " Xhat_unlabeled,\n", + " quantile,\n", + " B=2000,\n", + " alpha=0.05,\n", + " tuning_method=\"optimal\",\n", + ")\n", "\n", "classical_ci = classical_quantile_ci(X, quantile, alpha=0.05)\n", - "classical_pointestimate = (classical_ci[0]+classical_ci[1])/2\n", + "classical_pointestimate = (classical_ci[0] + classical_ci[1]) / 2\n", "\n", "naive_ci = classical_quantile_ci(preds_X, quantile, alpha=0.05)\n", - "naive_pointestimate = (naive_ci[0]+naive_ci[1])/2" + "naive_pointestimate = (naive_ci[0] + naive_ci[1]) / 2" ] }, { @@ -908,57 +1032,83 @@ "source": [ "plt.rcParams[\"font.sans-serif\"] = \"Arial\"\n", "\n", - "plt.rc('font', size=43) \n", - "plt.rc('axes', titlesize=43) \n", - "plt.rc('axes', labelsize=43) \n", - "plt.rc('xtick', labelsize=43) \n", - "plt.rc('ytick', labelsize=43) \n", - "plt.rc('figure', titlesize=43)\n", + "plt.rc(\"font\", size=43)\n", + "plt.rc(\"axes\", titlesize=43)\n", + "plt.rc(\"axes\", labelsize=43)\n", + "plt.rc(\"xtick\", labelsize=43)\n", + "plt.rc(\"ytick\", labelsize=43)\n", + "plt.rc(\"figure\", titlesize=43)\n", "\n", - "variables = [f'Tree Cover {100*quantile}% quantile']\n", + "variables = [f\"Tree Cover {100*quantile}% quantile\"]\n", "pd.set_option(\"display.precision\", 8)\n", "\n", "fig = plt.figure(figsize=(35, 6.5))\n", - "colors = ['blue', 'green', 'red']\n", + "colors = [\"blue\", \"green\", \"red\"]\n", "for i, variable in enumerate(variables):\n", - " \n", + "\n", " data_dict = {}\n", " classical_ci_width = classical_ci[1][i] - classical_ci[0][i]\n", " ptd_ci_width = ptd_ci[1][i] - ptd_ci[0][i]\n", - " data_dict['category'] = ['PTD','Classical','Naive']\n", - " data_dict['lower'] = [ptd_ci[0][i], classical_ci[0][i], naive_ci[0][i]]\n", - " data_dict['upper'] = [ptd_ci[1][i], classical_ci[1][i], naive_ci[1][i]]\n", - " data_dict['pointestimate'] = [ptd_pointestimate[i], classical_pointestimate[i], naive_pointestimate[i]]\n", + " data_dict[\"category\"] = [\"PTD\", \"Classical\", \"Naive\"]\n", + " data_dict[\"lower\"] = [ptd_ci[0][i], classical_ci[0][i], naive_ci[0][i]]\n", + " data_dict[\"upper\"] = [ptd_ci[1][i], classical_ci[1][i], naive_ci[1][i]]\n", + " data_dict[\"pointestimate\"] = [\n", + " ptd_pointestimate[i],\n", + " classical_pointestimate[i],\n", + " naive_pointestimate[i],\n", + " ]\n", " dataset = pd.DataFrame(data_dict)\n", - " \n", - " subplot = fig.add_subplot(141+i)\n", - " subplot.spines['top'].set_visible(False)\n", - " subplot.spines['right'].set_visible(False)\n", - " subplot.spines['left'].set_visible(False)\n", - " subplot.spines['bottom'].set_linewidth(2)\n", - " \n", + "\n", + " subplot = fig.add_subplot(141 + i)\n", + " subplot.spines[\"top\"].set_visible(False)\n", + " subplot.spines[\"right\"].set_visible(False)\n", + " subplot.spines[\"left\"].set_visible(False)\n", + " subplot.spines[\"bottom\"].set_linewidth(2)\n", + "\n", " # show true coefficient (yellow dotted line)\n", - " plt.axvline(true_coeff[i], color='tab:olive', linestyle='--', lw=4, label='True \\ncoefficient')\n", - " \n", + " plt.axvline(\n", + " true_coeff[i],\n", + " color=\"tab:olive\",\n", + " linestyle=\"--\",\n", + " lw=4,\n", + " label=\"True \\ncoefficient\",\n", + " )\n", + "\n", " # plot confidence intervals\n", - " for lower,upper,pointestimate,y in zip(dataset['lower'],dataset['upper'],dataset['pointestimate'],range(len(dataset))):\n", + " for lower, upper, pointestimate, y in zip(\n", + " dataset[\"lower\"],\n", + " dataset[\"upper\"],\n", + " dataset[\"pointestimate\"],\n", + " range(len(dataset)),\n", + " ):\n", " subplot.scatter([pointestimate], [y], color=colors[y], s=400)\n", - " subplot.scatter([lower, upper], [y, y], color=colors[y], marker='|', s=400, lw=5)\n", - " subplot.plot((lower,upper),(y,y),color=colors[y], lw=6)\n", + " subplot.scatter(\n", + " [lower, upper], [y, y], color=colors[y], marker=\"|\", s=400, lw=5\n", + " )\n", + " subplot.plot((lower, upper), (y, y), color=colors[y], lw=6)\n", " if y == 0:\n", " # compute PTD effective sample size improvement over classical method\n", - " ptd_effective_n_improvement = np.round((classical_ci_width/ptd_ci_width)**2, 1)\n", - " subplot.text((lower+upper)/2, 0.2, f'{ptd_effective_n_improvement}x', fontsize = 40, c='blue', weight='bold')\n", - " \n", - " subplot.tick_params(axis=u'both', which=u'both',length=0)\n", + " ptd_effective_n_improvement = np.round(\n", + " (classical_ci_width / ptd_ci_width) ** 2, 1\n", + " )\n", + " subplot.text(\n", + " (lower + upper) / 2,\n", + " 0.2,\n", + " f\"{ptd_effective_n_improvement}x\",\n", + " fontsize=40,\n", + " c=\"blue\",\n", + " weight=\"bold\",\n", + " )\n", + "\n", + " subplot.tick_params(axis=\"both\", which=\"both\", length=0)\n", " if i == 0:\n", - " subplot.set_yticks(range(len(dataset)),list(dataset['category']))\n", + " subplot.set_yticks(range(len(dataset)), list(dataset[\"category\"]))\n", " else:\n", " subplot.set_yticks([])\n", - " subplot.set_xlabel(f'{variable}')\n", + " subplot.set_xlabel(f\"{variable}\")\n", " subplot.set_ylim(-0.2)\n", - " \n", - "plt.tight_layout() \n", + "\n", + "plt.tight_layout()\n", "plt.show()" ] }, @@ -987,11 +1137,11 @@ "metadata": {}, "outputs": [], "source": [ - "truth_Y = np.array(data['truth_tree']).reshape(-1, 1)\n", - "preds_Y = np.array(data['preds_tree']).reshape(-1, 1)\n", + "truth_Y = np.array(data[\"truth_tree\"]).reshape(-1, 1)\n", + "preds_Y = np.array(data[\"preds_tree\"]).reshape(-1, 1)\n", "\n", - "truth_X = np.array(data[['truth_elevation', 'truth_population']])\n", - "preds_X = np.array(data[['truth_elevation', 'preds_population']])" + "truth_X = np.array(data[[\"truth_elevation\", \"truth_population\"]])\n", + "preds_X = np.array(data[[\"truth_elevation\", \"preds_population\"]])" ] }, { @@ -1256,14 +1406,14 @@ "source": [ "np.random.seed(seed=100)\n", "# map each point to its longitude quartile (1, 2, 3, or 4)\n", - "perc = np.percentile(data['lon'], [25,50,75])\n", - "data['quartile_number'] = 1 + np.digitize(data['lon'], perc)\n", + "perc = np.percentile(data[\"lon\"], [25, 50, 75])\n", + "data[\"quartile_number\"] = 1 + np.digitize(data[\"lon\"], perc)\n", "\n", "# assign quartile scores q(i) = 5^i\n", - "data['quartile_score'] = np.array([5**i for i in data['quartile_number']])\n", + "data[\"quartile_score\"] = np.array([5**i for i in data[\"quartile_number\"]])\n", "\n", - "# compute probability of sampling each point by normalizing the quartile scores to sum to 1 \n", - "data['probability'] = data['quartile_score']/data['quartile_score'].sum()\n", + "# compute probability of sampling each point by normalizing the quartile scores to sum to 1\n", + "data[\"probability\"] = data[\"quartile_score\"] / data[\"quartile_score\"].sum()\n", "\n", "display(data)" ] @@ -1285,14 +1435,20 @@ "source": [ "# sample calibration set using the given probabilities\n", "n = 1000\n", - "calibration_indices = np.random.choice(np.arange(0, len(data)), size=n, replace=False, p=data['probability'])\n", + "calibration_indices = np.random.choice(\n", + " np.arange(0, len(data)), size=n, replace=False, p=data[\"probability\"]\n", + ")\n", "\n", "X = statsmodels.tools.add_constant(truth_X[calibration_indices])\n", "Xhat = statsmodels.tools.add_constant(preds_X[calibration_indices])\n", - "Xhat_unlabeled = statsmodels.tools.add_constant(np.delete(preds_X, calibration_indices, axis=0)) # all predicted datapoints except calibration indices\n", + "Xhat_unlabeled = statsmodels.tools.add_constant(\n", + " np.delete(preds_X, calibration_indices, axis=0)\n", + ") # all predicted datapoints except calibration indices\n", "Y = truth_Y[calibration_indices]\n", "Yhat = preds_Y[calibration_indices]\n", - "Yhat_unlabeled = np.delete(preds_Y, calibration_indices, axis=0) # all predicted datapoints except calibration indices" + "Yhat_unlabeled = np.delete(\n", + " preds_Y, calibration_indices, axis=0\n", + ") # all predicted datapoints except calibration indices" ] }, { @@ -1313,11 +1469,13 @@ "outputs": [], "source": [ "# to compute sample weights, we need to rescale the probabilities so they sum to n (the number of calibration samples)\n", - "data['probability'] = n*data['probability'] \n", + "data[\"probability\"] = n * data[\"probability\"]\n", "\n", "# sample weights (Inverse Probability Weighting)\n", - "w = (1/data['probability'])[calibration_indices].to_numpy()\n", - "w_unlabeled = np.delete(1/(1-data['probability']), calibration_indices, axis=0) # all predicted datapoints except calibration indices" + "w = (1 / data[\"probability\"])[calibration_indices].to_numpy()\n", + "w_unlabeled = np.delete(\n", + " 1 / (1 - data[\"probability\"]), calibration_indices, axis=0\n", + ") # all predicted datapoints except calibration indices" ] }, { @@ -1337,17 +1495,31 @@ "metadata": {}, "outputs": [], "source": [ - "true_coeff = ptd.algorithm_linear_regression(data=[statsmodels.tools.add_constant(truth_X), truth_Y], w=None)\n", - "\n", - "tuning_matrix, ptd_pointestimate, ptd_ci = ptd.ptd_linear_regression(X, Xhat, Xhat_unlabeled, Y, Yhat, Yhat_unlabeled, \n", - " w=w, w_unlabeled=w_unlabeled,\n", - " B=2000, alpha=0.05, tuning_method='optimal')\n", + "true_coeff = ptd.algorithm_linear_regression(\n", + " data=[statsmodels.tools.add_constant(truth_X), truth_Y], w=None\n", + ")\n", + "\n", + "tuning_matrix, ptd_pointestimate, ptd_ci = ptd.ptd_linear_regression(\n", + " X,\n", + " Xhat,\n", + " Xhat_unlabeled,\n", + " Y,\n", + " Yhat,\n", + " Yhat_unlabeled,\n", + " w=w,\n", + " w_unlabeled=w_unlabeled,\n", + " B=2000,\n", + " alpha=0.05,\n", + " tuning_method=\"optimal\",\n", + ")\n", "\n", "classical_ci = classical_linear_regression_ci(X, Y, alpha=0.05, w=w)\n", - "classical_pointestimate = (classical_ci[0]+classical_ci[1])/2\n", + "classical_pointestimate = (classical_ci[0] + classical_ci[1]) / 2\n", "\n", - "naive_ci = classical_linear_regression_ci(statsmodels.tools.add_constant(preds_X), preds_Y, alpha=0.05)\n", - "naive_pointestimate = (naive_ci[0]+naive_ci[1])/2" + "naive_ci = classical_linear_regression_ci(\n", + " statsmodels.tools.add_constant(preds_X), preds_Y, alpha=0.05\n", + ")\n", + "naive_pointestimate = (naive_ci[0] + naive_ci[1]) / 2" ] }, { @@ -1571,58 +1743,84 @@ "source": [ "plt.rcParams[\"font.sans-serif\"] = \"Arial\"\n", "\n", - "plt.rc('font', size=43) \n", - "plt.rc('axes', titlesize=43) \n", - "plt.rc('axes', labelsize=43) \n", - "plt.rc('xtick', labelsize=43) \n", - "plt.rc('ytick', labelsize=43) \n", - "plt.rc('figure', titlesize=43)\n", + "plt.rc(\"font\", size=43)\n", + "plt.rc(\"axes\", titlesize=43)\n", + "plt.rc(\"axes\", labelsize=43)\n", + "plt.rc(\"xtick\", labelsize=43)\n", + "plt.rc(\"ytick\", labelsize=43)\n", + "plt.rc(\"figure\", titlesize=43)\n", "\n", - "covariates = ['Intercept', 'Elevation', 'Population']\n", + "covariates = [\"Intercept\", \"Elevation\", \"Population\"]\n", "pd.set_option(\"display.precision\", 8)\n", "\n", "fig = plt.figure(figsize=(32, 6.5))\n", - "colors = ['blue', 'green', 'red']\n", + "colors = [\"blue\", \"green\", \"red\"]\n", "for i, covariate in enumerate(covariates):\n", - " \n", + "\n", " data_dict = {}\n", " classical_ci_width = classical_ci[1][i] - classical_ci[0][i]\n", " ptd_ci_width = ptd_ci[1][i] - ptd_ci[0][i]\n", - " data_dict['category'] = ['PTD (weighted)', 'Classical (weighted)','Naive']\n", - " data_dict['lower'] = [ptd_ci[0][i], classical_ci[0][i], naive_ci[0][i]]\n", - " data_dict['upper'] = [ptd_ci[1][i], classical_ci[1][i], naive_ci[1][i]]\n", - " data_dict['pointestimate'] = [ptd_pointestimate[i], classical_pointestimate[i], naive_pointestimate[i]]\n", + " data_dict[\"category\"] = [\"PTD (weighted)\", \"Classical (weighted)\", \"Naive\"]\n", + " data_dict[\"lower\"] = [ptd_ci[0][i], classical_ci[0][i], naive_ci[0][i]]\n", + " data_dict[\"upper\"] = [ptd_ci[1][i], classical_ci[1][i], naive_ci[1][i]]\n", + " data_dict[\"pointestimate\"] = [\n", + " ptd_pointestimate[i],\n", + " classical_pointestimate[i],\n", + " naive_pointestimate[i],\n", + " ]\n", " dataset = pd.DataFrame(data_dict)\n", " display(dataset)\n", - " \n", - " subplot = fig.add_subplot(141+i)\n", - " subplot.spines['top'].set_visible(False)\n", - " subplot.spines['right'].set_visible(False)\n", - " subplot.spines['left'].set_visible(False)\n", - " subplot.spines['bottom'].set_linewidth(2)\n", - " \n", + "\n", + " subplot = fig.add_subplot(141 + i)\n", + " subplot.spines[\"top\"].set_visible(False)\n", + " subplot.spines[\"right\"].set_visible(False)\n", + " subplot.spines[\"left\"].set_visible(False)\n", + " subplot.spines[\"bottom\"].set_linewidth(2)\n", + "\n", " # show true coefficient (yellow dotted line)\n", - " plt.axvline(true_coeff[i], color='tab:olive', linestyle='--', lw=4, label='True \\ncoefficient')\n", - " \n", + " plt.axvline(\n", + " true_coeff[i],\n", + " color=\"tab:olive\",\n", + " linestyle=\"--\",\n", + " lw=4,\n", + " label=\"True \\ncoefficient\",\n", + " )\n", + "\n", " # plot confidence intervals\n", - " for lower,upper,pointestimate,y in zip(dataset['lower'],dataset['upper'],dataset['pointestimate'],range(len(dataset))):\n", + " for lower, upper, pointestimate, y in zip(\n", + " dataset[\"lower\"],\n", + " dataset[\"upper\"],\n", + " dataset[\"pointestimate\"],\n", + " range(len(dataset)),\n", + " ):\n", " subplot.scatter([pointestimate], [y], color=colors[y], s=400)\n", - " subplot.scatter([lower, upper], [y, y], color=colors[y], marker='|', s=400, lw=5)\n", - " subplot.plot((lower,upper),(y,y),color=colors[y], lw=6)\n", + " subplot.scatter(\n", + " [lower, upper], [y, y], color=colors[y], marker=\"|\", s=400, lw=5\n", + " )\n", + " subplot.plot((lower, upper), (y, y), color=colors[y], lw=6)\n", " if y == 0:\n", " # compute PTD effective sample size improvement over classical method\n", - " ptd_effective_n_improvement = np.round((classical_ci_width/ptd_ci_width)**2, 1)\n", - " subplot.text((lower+upper)/2, 0.2, f'{ptd_effective_n_improvement}x', fontsize = 40, c='blue', weight='bold')\n", - " \n", - " subplot.tick_params(axis=u'both', which=u'both',length=0)\n", + " ptd_effective_n_improvement = np.round(\n", + " (classical_ci_width / ptd_ci_width) ** 2, 1\n", + " )\n", + " subplot.text(\n", + " (lower + upper) / 2,\n", + " 0.2,\n", + " f\"{ptd_effective_n_improvement}x\",\n", + " fontsize=40,\n", + " c=\"blue\",\n", + " weight=\"bold\",\n", + " )\n", + "\n", + " subplot.tick_params(axis=\"both\", which=\"both\", length=0)\n", " if i == 0:\n", - " subplot.set_yticks(range(len(dataset)),list(dataset['category']))\n", + " subplot.set_yticks(range(len(dataset)), list(dataset[\"category\"]))\n", " else:\n", " subplot.set_yticks([])\n", - " subplot.set_xlabel(f'{covariate}')\n", + " subplot.set_xlabel(f\"{covariate}\")\n", " subplot.set_ylim(-0.2)\n", - " \n", - "plt.tight_layout() \n", + "\n", + "plt.tight_layout()\n", "plt.show()" ] } diff --git a/tests/test_mean.py b/tests/test_mean.py index 2f51998..85cd98e 100644 --- a/tests/test_mean.py +++ b/tests/test_mean.py @@ -46,7 +46,7 @@ def test_ppi_mean_multid(): included = (ci[0] <= 0) & (ci[1] >= 0) includeds[j] += included.astype(int) - print(includeds/trials) + print(includeds / trials) failures = (includeds / trials).T < 1 - alphas - epsilon print(failures) failed = np.any(failures) @@ -60,9 +60,7 @@ def test_ppi_mean_elem(): Yhat_unlabeled = np.random.normal(-2, 1, 10000) ppi_mean_pointestimate(Y, Yhat, Yhat_unlabeled, lam_optim_mode="element") - ppi_mean_ci( - Y, Yhat, Yhat_unlabeled, alpha=alpha, lam_optim_mode="element" - ) + ppi_mean_ci(Y, Yhat, Yhat_unlabeled, alpha=alpha, lam_optim_mode="element") ppi_mean_pval(Y, Yhat, Yhat_unlabeled, lam_optim_mode="element") Y = np.random.normal(0, 1, (10000, 5)) @@ -70,9 +68,7 @@ def test_ppi_mean_elem(): Yhat_unlabeled = np.random.normal(-2, 1, (10000, 5)) ppi_mean_pointestimate(Y, Yhat, Yhat_unlabeled, lam_optim_mode="element") - ppi_mean_ci( - Y, Yhat, Yhat_unlabeled, alpha=alpha, lam_optim_mode="element" - ) + ppi_mean_ci(Y, Yhat, Yhat_unlabeled, alpha=alpha, lam_optim_mode="element") ppi_mean_pval(Y, Yhat, Yhat_unlabeled, lam_optim_mode="element")