diff --git a/appletree/parameter.py b/appletree/parameter.py index dd52e91a..19402b60 100644 --- a/appletree/parameter.py +++ b/appletree/parameter.py @@ -4,6 +4,7 @@ import numpy as np from appletree.randgen import TwoHalfNorm +from appletree.utils import errors_to_two_half_norm_sigmas class Parameter: @@ -81,10 +82,13 @@ def sample_prior(self): val = np.random.normal(**kwargs) self._parameter_dict[par_name] = np.clip(val, *setting["allowed_range"]) elif prior_type == "twohalfnorm": + # We need to convert errors to sigmas + # See the docstring of errors_to_two_half_norm_sigmas for details + sigmas = errors_to_two_half_norm_sigmas([args["sigma_pos"], args["sigma_neg"]]) kwargs = { "mu": args["mu"], - "sigma_pos": args["sigma_pos"], - "sigma_neg": args["sigma_neg"], + "sigma_pos": sigmas[0], + "sigma_neg": sigmas[1], } val = TwoHalfNorm.rvs(**kwargs) self._parameter_dict[par_name] = np.clip(val, *setting["allowed_range"]) @@ -150,14 +154,15 @@ def log_prior(self): std = args["std"] log_prior += -((val - mean) ** 2) / 2 / std**2 elif prior_type == "twohalfnorm": + # We need to convert errors to sigmas + # See the docstring of errors_to_two_half_norm_sigmas for details + sigmas = errors_to_two_half_norm_sigmas([args["sigma_pos"], args["sigma_neg"]]) mu = args["mu"] - sigma_pos = args["sigma_pos"] - sigma_neg = args["sigma_neg"] log_prior += TwoHalfNorm.logpdf( x=val, mu=mu, - sigma_pos=sigma_pos, - sigma_neg=sigma_neg, + sigma_pos=sigmas[0], + sigma_neg=sigmas[1], ) elif prior_type == "free": pass diff --git a/appletree/utils.py b/appletree/utils.py index 5d42c7ec..c5340ef6 100644 --- a/appletree/utils.py +++ b/appletree/utils.py @@ -10,6 +10,9 @@ import matplotlib as mpl from matplotlib.patches import Rectangle from matplotlib import pyplot as plt +from scipy.special import erf +from scipy.optimize import root +from scipy.stats import chi2 import GOFevaluation from appletree.share import _cached_configs @@ -602,8 +605,35 @@ def cum_integrate_midpoint(x, y): return x_mid, np.cumsum(dx * y_mid) +@export def check_unused_configs(): """Check if there are unused configs.""" unused_configs = set(_cached_configs.keys()) - _cached_configs.accessed_keys if unused_configs: warn(f"Detected unused configs: {unused_configs}, you might set the configs incorrectly.") + + +@export +def errors_to_two_half_norm_sigmas(errors): + """This function solves the sigmas for a two-half-norm distribution, such that the 16 and 84 + percentile corresponds to the given errors. + + In the two-half-norm distribution, the positive and negative errors are assumed to be + the std of the glued normal distributions. While we interpret the 16 and 84 percentile as + the input errors, thus we need to solve the sigmas for the two-half-norm distribution. + The solution is determined by the following conditions: + - The 16 percentile of the two-half-norm distribution should be the negative error. + - The 84 percentile of the two-half-norm distribution should be the positive error. + - The mode of the two-half-norm distribution should be 0. + + """ + + def _to_solve(x, errors, p): + return [ + x[0] / (x[0] + x[1]) * (1 - erf(errors[0] / x[0] / np.sqrt(2))) - p / 2, + x[1] / (x[0] + x[1]) * (1 - erf(errors[1] / x[1] / np.sqrt(2))) - p / 2, + ] + + res = root(_to_solve, errors, args=(errors, 1 - chi2.cdf(1, 1))) + assert res.success, f"Cannot solve sigmas of TwoHalfNorm for errors {errors}!" + return res.x