Merge pull request #86 from XENONnT/asymp_n_dof

Add optional argument `degree_of_freedom` for `asymptotic_critical_value`

alea/examples/configs/unbinned_wimp_running.yaml

@@ -49,6 +49,7 @@ computation_options:
         n_batch: 40,
         compute_confidence_interval: True,
         limit_threshold: "thresholds.json",
+        asymptotic_dof: 1,
         toydata_mode: "generate_and_store",
         toydata_filename: "toyfile_wimp_mass_{wimp_mass:d}_poi_expectation_{poi_expectation:.2f}.h5",
       }

alea/model.py

@@ -72,6 +72,7 @@ def __init__(
         confidence_level: float = 0.9,
         confidence_interval_kind: str = "central",  # one of central, upper, lower
         confidence_interval_threshold: Optional[Callable[[float], float]] = None,
+        asymptotic_dof: Optional[int] = 1,
         data: Optional[Union[dict, list]] = None,
         **kwargs,
     ):
@@ -93,6 +94,7 @@ def __init__(
             raise ValueError("confidence_interval_kind must be one of central, upper, lower")
         self._confidence_interval_kind = confidence_interval_kind
         self.confidence_interval_threshold = confidence_interval_threshold
+        self.asymptotic_dof = asymptotic_dof
         self._define_parameters(parameter_definition)
 
         self._check_ll_and_generate_data_signature()
@@ -355,6 +357,7 @@ def _confidence_interval_checks(
         confidence_level: Optional[float] = None,
         confidence_interval_kind: Optional[str] = None,
         confidence_interval_threshold: Optional[Callable[[float], float]] = None,
+        asymptotic_dof: Optional[int] = None,
         **kwargs,
     ) -> Tuple[str, Callable[[float], float], Tuple[float, float]]:
         """Helper function for confidence_interval that does the input checks and return bounds.
@@ -400,8 +403,22 @@ def _confidence_interval_checks(
                 confidence_interval_threshold = self.confidence_interval_threshold
             else:
                 # use asymptotic thresholds assuming the test statistic is Chi2 distributed
+                if asymptotic_dof is None:
+                    if self.asymptotic_dof is not None:
+                        degree_of_freedom = self.asymptotic_dof
+                    else:
+                        degree_of_freedom = 1
+                else:
+                    if self.asymptotic_dof is not None:
+                        if asymptotic_dof != self.asymptotic_dof:
+                            warnings.warn(
+                                f"You set asymptotic_dof as {asymptotic_dof}, "
+                                f"which is different from the default value "
+                                f"{self.asymptotic_dof}. Be careful!"
+                            )
+                    degree_of_freedom = asymptotic_dof
                 critical_value = asymptotic_critical_value(
-                    confidence_interval_kind, confidence_level
+                    confidence_interval_kind, confidence_level, degree_of_freedom
                 )
 
                 def confidence_interval_threshold(_):
@@ -427,6 +444,7 @@ def confidence_interval(
         confidence_interval_threshold: Optional[Callable[[float], float]] = None,
         confidence_interval_args: Optional[dict] = None,
         best_fit_args: Optional[dict] = None,
+        asymptotic_dof: Optional[int] = None,
     ) -> Tuple[float, float]:
         """Uses self.fit to compute confidence intervals for a certain named parameter. If the
         parameter is a rate parameter, and the model has expectation values implemented, the bounds
@@ -454,6 +472,7 @@ def confidence_interval(
                 profile likelihood-- mainly used for 1-D slices of higher-dimensional confidence
                 volumes, where the global best-fit may not be along the profile.
                 If None, will be set to confidence_interval_args.
+            asymptotic_dof (int, optional (default=None)): Degrees of freedom for asymptotic
 
         """
         if confidence_interval_args is None:
@@ -466,6 +485,7 @@ def confidence_interval(
             confidence_level,
             confidence_interval_kind,
             confidence_interval_threshold,
+            asymptotic_dof,
             **confidence_interval_args,
         )
         (

alea/models/blueice_extended_model.py

@@ -91,7 +91,7 @@ def data(self, data: Union[dict, list]):
                 If data is a list, it must be a list of length len(self.likelihood_names) + 1.
 
         Raises:
-            ValueError: If data is not a list of length len(self.likelihood_names) + 1.
+            Warning: If data is not a list of length len(self.likelihood_names) + 1.
 
         Caution:
             The self._data is read-only, so you can not change the data after it is set.

alea/runner.py

@@ -156,6 +156,7 @@ def __init__(
             confidence_interval_kind,
             confidence_level,
             statistical_model_args.get("limit_threshold_interpolation", False),
+            statistical_model_args.get("asymptotic_dof", 1),
         )
 
     @property

alea/submitter.py

@@ -453,6 +453,10 @@ def update_statistical_model_args(
             runner_args["statistical_model_args"][
                 "limit_threshold_interpolation"
             ] = runner_args.pop("limit_threshold_interpolation")
+        if "asymptotic_dof" in runner_args:
+            runner_args["statistical_model_args"]["asymptotic_dof"] = runner_args.pop(
+                "asymptotic_dof"
+            )
 
     @staticmethod
     def check_redunant_arguments(runner_args, allowed_special_args: List[str] = []):

alea/submitters/local.py

@@ -7,7 +7,7 @@
 import operator
 from functools import reduce
 from copy import deepcopy
-from typing import List, Dict, Any, cast
+from typing import List, Dict, Any, Optional, Callable, cast
 
 import numpy as np
 from scipy.interpolate import interp1d, RegularGridInterpolator
@@ -349,6 +349,7 @@ def get_confidence_interval_thresholds(
         confidence_interval_kind,
         confidence_level,
         limit_threshold_interpolation,
+        asymptotic_dof: Optional[int] = 1,
     ):
         """Get confidence interval threshold function from limit_threshold file. If the
         limit_threshold file does not contain the threshold, it will interpolate the threshold from
@@ -364,6 +365,8 @@ def get_confidence_interval_thresholds(
             confidence_level (float): confidence level
             limit_threshold_interpolation (bool): whether to interpolate the threshold from the
                 existing threshold, if the limit_threshold file does not contain the threshold
+            asymptotic_dof (int, optional (default=1)):
+                degrees of freedom for asymptotic critical value
 
         """
 
@@ -372,7 +375,7 @@ def get_confidence_interval_thresholds(
 
         threshold = load_json(limit_threshold)
 
-        func_list = []
+        func_list: List[Optional[Callable]] = []
         for i_hypo in range(len(hypotheses_values)):
             hypothesis = hypotheses_values[i_hypo]
 
@@ -454,7 +457,9 @@ def get_confidence_interval_thresholds(
                 poi_values,
                 threshold_values,
                 bounds_error=False,
-                fill_value=asymptotic_critical_value(confidence_interval_kind, confidence_level),
+                fill_value=asymptotic_critical_value(
+                    confidence_interval_kind, confidence_level, asymptotic_dof
+                ),
             )
             func_list.append(func)
 

alea/utils.py

@@ -283,25 +283,37 @@ def get_template_folder_list(likelihood_config, extra_template_path: Optional[st
     return template_folder_list
 
 
-def asymptotic_critical_value(confidence_interval_kind: str, confidence_level: float):
+def asymptotic_critical_value(
+    confidence_interval_kind: str, confidence_level: float, degree_of_freedom: Optional[int] = None
+):
     """Return the critical value for the confidence interval.
 
     Args:
         confidence_interval_kind (str): confidence interval kind, either 'lower', 'upper' or
             'central'
         confidence_level (float): confidence level
+        degree_of_freedom (int, optional (default=None)): degree of freedom
 
     Returns:
         float: critical value
 
-    Caution:
-        The critical value is calculated from chi2 distribution with 1 degree of freedom.
+    Raises:
+        ValueError: if confidence_interval_kind is not 'lower', 'upper' or 'central'
+        ValueError: if degree_of_freedom is not None and not 1, when confidence_interval_kind is
+            'lower' or 'upper'
 
     """
     if confidence_interval_kind in {"lower", "upper"}:
+        if (degree_of_freedom is not None) and (degree_of_freedom != 1):
+            raise ValueError(
+                f"degree_of_freedom must be 1 for {confidence_interval_kind} confidence interval"
+            )
         critical_value = chi2(1).isf(2 * (1.0 - confidence_level))
     elif confidence_interval_kind == "central":
-        critical_value = chi2(1).isf(1.0 - confidence_level)
+        if degree_of_freedom is None:
+            critical_value = chi2(1).isf(1.0 - confidence_level)
+        else:
+            critical_value = chi2(degree_of_freedom).isf(1.0 - confidence_level)
     else:
         raise ValueError(
             f"confidence_interval_kind must be either 'lower', 'upper' or 'central', "

tests/test_utils.py

@@ -40,6 +40,15 @@ def test_asymptotic_critical_value(self):
         critical_value = chi2(1).isf(2 * (1.0 - confidence_level))
         self.assertEqual(asymptotic_critical_value("lower", confidence_level), critical_value)
         self.assertEqual(asymptotic_critical_value("upper", confidence_level), critical_value)
+        with self.assertRaises(ValueError):
+            asymptotic_critical_value("invalid", confidence_level)
+        with self.assertRaises(ValueError):
+            asymptotic_critical_value("lower", confidence_level, 2)
+
+        critical_value_central = chi2(2).isf((1.0 - confidence_level))
+        self.assertEqual(
+            asymptotic_critical_value("central", confidence_level, 2), critical_value_central
+        )
 
     def test_within_limits(self):
         """Test of the within_limits function."""