diff --git a/.zenodo.json b/.zenodo.json index c635095..450823e 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -1,12 +1,12 @@ { "description": "pyhf plug-in for spey package", "license": "MIT", - "title": "SpeysideHEP/spey-pyhf: v0.1.6", - "version": "v0.1.6", + "title": "SpeysideHEP/spey-pyhf: v0.1.7", + "version": "v0.1.7", "upload_type": "software", "creators": [ { - "affiliation": "Thomas Jefferson National Accelerator Facility", + "affiliation": "Stony Brook University", "name": "Araz, Jack Y.", "orcid": "0000-0001-8721-8042" } @@ -29,7 +29,7 @@ }, { "scheme": "url", - "identifier": "https://github.com/SpeysideHEP/spey-pyhf/tree/v0.1.6", + "identifier": "https://github.com/SpeysideHEP/spey-pyhf/tree/v0.1.7", "relation": "isSupplementTo" }, { diff --git a/docs/releases/changelog-v0.1.md b/docs/releases/changelog-v0.1.md index cdc9eb0..7121bc2 100644 --- a/docs/releases/changelog-v0.1.md +++ b/docs/releases/changelog-v0.1.md @@ -7,6 +7,9 @@ and usage can be found in the online documentation. ([#2](https://github.com/SpeysideHEP/spey-pyhf/pull/2)) +* Full statistical model mapping to effective sigma model have been implemented. + ([#14](https://github.com/SpeysideHEP/spey-pyhf/pull/14)) + ## Improvements * Sampler functionality has been extended to isolate auxiliary data. diff --git a/src/spey_pyhf/_version.py b/src/spey_pyhf/_version.py index 5260485..0e31141 100644 --- a/src/spey_pyhf/_version.py +++ b/src/spey_pyhf/_version.py @@ -1,3 +1,3 @@ """Version of the spey - pyhf plugin""" -__version__ = "0.1.6" +__version__ = "0.1.7" diff --git a/src/spey_pyhf/simplify.py b/src/spey_pyhf/simplify.py index a15589b..0b01d10 100644 --- a/src/spey_pyhf/simplify.py +++ b/src/spey_pyhf/simplify.py @@ -3,13 +3,19 @@ import logging import warnings from contextlib import contextmanager +from pathlib import Path from typing import Callable, List, Literal, Optional, Text, Union import numpy as np import spey import tqdm -from scipy.stats import moment, multivariate_normal -from spey.backends.default_pdf import CorrelatedBackground, ThirdMomentExpansion +from scipy.stats import moment, multivariate_normal, norm +from spey.backends.default_pdf import ( + CorrelatedBackground, + EffectiveSigma, + ThirdMomentExpansion, +) +from spey.helper_functions import covariance_to_correlation from spey.optimizer.core import fit from . import WorkspaceInterpreter @@ -78,7 +84,8 @@ class Simplify(spey.ConverterBase): fittype (``Text``, default ``"postfit"``): what type of fitting should be performed ``"postfit"`` or ``"prefit"``. convert_to (``Text``, default ``"default_pdf.correlated_background"``): conversion type. Should - be either ``"default_pdf.correlated_background"`` or ``"default_pdf.third_moment_expansion"``. + be either ``"default_pdf.correlated_background"``, ``"default_pdf.third_moment_expansion"`` + or ``"default_pdf.effective_sigma"``. number_of_samples (``int``, default ``1000``): number of samples to be generated in order to estimate contract the uncertainties into a single value. control_region_indices (``List[int]`` or ``List[Text]``, default ``None``): indices or names of the control and @@ -171,13 +178,15 @@ def __call__( statistical_model: spey.StatisticalModel, fittype: Literal["postfit", "prefit"] = "postfit", convert_to: Literal[ - "default_pdf.correlated_background", "default_pdf.third_moment_expansion" + "default_pdf.correlated_background", + "default_pdf.third_moment_expansion", + "default_pdf.effective_sigma", ] = "default_pdf.correlated_background", number_of_samples: int = 1000, control_region_indices: Optional[Union[List[int], List[Text]]] = None, include_modifiers_in_control_model: bool = False, save_model: Optional[Text] = None, - ) -> Union[CorrelatedBackground, ThirdMomentExpansion]: + ) -> Union[CorrelatedBackground, ThirdMomentExpansion, EffectiveSigma]: assert statistical_model.backend_type == "pyhf", ( "This method is currently only available for `pyhf` full statistical models." @@ -390,6 +399,13 @@ def __call__( # in the full statistical model! background_yields = np.mean(samples, axis=0) + save_kwargs = { + "covariance_matrix": covariance_matrix, + "background_yields": background_yields, + "data": data, + "channel_order": stat_model_pyhf.config.channels, + } + third_moments = [] if convert_to == "default_pdf.correlated_background": backend = CorrelatedBackground( @@ -400,6 +416,7 @@ def __call__( ) elif convert_to == "default_pdf.third_moment_expansion": third_moments = moment(samples, moment=3, axis=0) + save_kwargs.update({"third_moments": third_moments}) backend = ThirdMomentExpansion( signal_yields=signal_yields, @@ -408,21 +425,37 @@ def __call__( covariance_matrix=covariance_matrix, third_moment=third_moments, ) + elif convert_to == "default_pdf.effective_sigma": + # Get 68% quantiles + q = (1.0 - (norm.cdf(1.0) - norm.cdf(-1.0))) / 2.0 + absolute_uncertainty_envelops = np.stack( + [np.quantile(samples, q, axis=0), np.quantile(samples, 1 - q, axis=0)], + axis=1, + ) + save_kwargs.update( + {"absolute_uncertainty_envelops": absolute_uncertainty_envelops} + ) + + backend = EffectiveSigma( + signal_yields=signal_yields, + background_yields=background_yields, + data=data, + correlation_matrix=covariance_to_correlation( + covariance_matrix=covariance_matrix + ), + absolute_uncertainty_envelops=absolute_uncertainty_envelops, + ) else: raise ConversionError( "Currently available conversion methods are " - + "'default_pdf.correlated_background', 'default_pdf.third_moment_expansion'" + + "'default_pdf.correlated_background', 'default_pdf.third_moment_expansion'," + + " 'default_pdf.effective_sigma'" ) if save_model is not None: - if save_model.endswith(".npz"): - np.savez_compressed( - save_model, - covariance_matrix=covariance_matrix, - background_yields=background_yields, - third_moments=third_moments, - data=data, - channel_order=stat_model_pyhf.config.channels, - ) + save_path = Path(save_model) + if save_path.suffix != ".npz": + save_path = save_path.with_suffix(".npz") + np.savez_compressed(str(save_path), **save_kwargs) return backend