Skip to content

Commit

Permalink
rename alpha to mu
Browse files Browse the repository at this point in the history
  • Loading branch information
ikrommyd committed Mar 23, 2024
1 parent ca8a1e1 commit ad1be34
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 32 deletions.
41 changes: 20 additions & 21 deletions tests/test_pdf_erfexp.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,48 +9,47 @@

import zfit_physics as zphys

alpha_true = 90.0
mu_true = 90.0
beta_true = 0.08
gamma_true = -1.0
n_true = 0.2

alpha_range = (65.0, 90.0)
mu_range = (65.0, 90.0)
gamma_range = (-10, 10)
beta_range = (0.01, 10)
n_range = (0.1, 0.5)


def _erfexp_numpy(x, alpha, beta, gamma, n):
return special.erfc((x - alpha) * beta) * np.exp(-gamma * (np.power(x, n) - np.power(alpha, n)))
def _erfexp_numpy(x, mu, beta, gamma, n):
return special.erfc((x - mu) * beta) * np.exp(-gamma * (np.power(x, n) - np.power(mu, n)))


erfexp_numpy = np.vectorize(_erfexp_numpy, excluded=["alpha", "beta", "gamma", "n"])
erfexp_numpy = np.vectorize(_erfexp_numpy, excluded=["mu", "beta", "gamma", "n"])


def create_erfexp(alpha, beta, gamma, n, limits):
def create_erfexp(mu, beta, gamma, n, limits):
obs = zfit.Space("obs1", limits=limits)
erfexp = zphys.pdf.ErfExp(alpha=alpha, beta=beta, gamma=gamma, n=n, obs=obs)
erfexp = zphys.pdf.ErfExp(mu=mu, beta=beta, gamma=gamma, n=n, obs=obs)
return erfexp, obs


def test_erfexp_pdf():
# Test PDF here
erfexp, _ = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
erfexp, _ = create_erfexp(mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
assert erfexp.pdf(90.0, norm=False).numpy().item() == pytest.approx(
erfexp_numpy(90.0, alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true), rel=1e-8
erfexp_numpy(90.0, mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true), rel=1e-8
)
assert erfexp.pdf(90.0).numpy().item() == pytest.approx(
erfexp_numpy(90.0, alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true) / 71.18838, rel=1e-8
erfexp_numpy(90.0, mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true) / 71.18838, rel=1e-8
)
np.testing.assert_allclose(
erfexp.pdf(tf.range(50.0, 130, 10_000), norm=False),
erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),
erfexp_numpy(tf.range(50.0, 130, 10_000), mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true),
rtol=1e-8,
)
np.testing.assert_allclose(
erfexp.pdf(tf.range(50.0, 130, 10_000)),
erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true)
/ 71.18838,
erfexp_numpy(tf.range(50.0, 130, 10_000), mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true) / 71.18838,
rtol=1e-8,
atol=1e-8,
)
Expand All @@ -59,41 +58,41 @@ def test_erfexp_pdf():
def test_erfexp_pdf_random_params():
# Test PDF here in a loop with random parameters
for _ in range(1000):
alpha_true = np.random.uniform(*alpha_range)
mu_true = np.random.uniform(*mu_range)
beta_true = np.random.uniform(*beta_range)
gamma_true = np.random.uniform(*gamma_range)
n_true = np.random.uniform(*n_range)

erfexp, __ = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
erfexp, __ = create_erfexp(mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
np.testing.assert_allclose(
erfexp.pdf(tf.range(50.0, 130, 10_000), norm=False),
erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),
erfexp_numpy(tf.range(50.0, 130, 10_000), mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true),
rtol=1e-5,
)


def test_erfexp_integral():
# Test CDF and integral here
erfexp, obs = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
erfexp, obs = create_erfexp(mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
full_interval_numeric = erfexp.numeric_integrate(obs, norm=False).numpy()
true_integral = 71.18838
numpy_full_integral = integrate.quad(erfexp_numpy, 50, 130, args=(alpha_true, beta_true, gamma_true, n_true))[0]
numpy_full_integral = integrate.quad(erfexp_numpy, 50, 130, args=(mu_true, beta_true, gamma_true, n_true))[0]
assert full_interval_numeric == pytest.approx(true_integral, 1e-7)
assert full_interval_numeric == pytest.approx(numpy_full_integral, 1e-7)

numeric_integral = erfexp.numeric_integrate(limits=(80, 100), norm=False).numpy()
numpy_integral = integrate.quad(erfexp_numpy, 80, 100, args=(alpha_true, beta_true, gamma_true, n_true))[0]
numpy_integral = integrate.quad(erfexp_numpy, 80, 100, args=(mu_true, beta_true, gamma_true, n_true))[0]
assert numeric_integral == pytest.approx(numpy_integral, 1e-7)


# register the pdf here and provide sets of working parameter configurations
def erfexp_params_factory():
alpha = zfit.Parameter("alpha", alpha_true)
mu = zfit.Parameter("mu", mu_true)
beta = zfit.Parameter("beta", beta_true)
gamma = zfit.Parameter("gamma", gamma_true)
n = zfit.Parameter("n", n_true)

return {"alpha": alpha, "beta": beta, "gamma": gamma, "n": n}
return {"mu": mu, "beta": beta, "gamma": gamma, "n": n}


tester.register_pdf(pdf_class=zphys.pdf.ErfExp, params_factories=erfexp_params_factory)
22 changes: 11 additions & 11 deletions zfit_physics/models/pdf_erfexp.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,12 +8,12 @@


@z.function(wraps="tensor")
def erfexp_pdf_func(x, alpha, beta, gamma, n):
def erfexp_pdf_func(x, mu, beta, gamma, n):
"""Calculate the ErfExp PDF.
Args:
x: value(s) for which the PDF will be calculated.
alpha: Location parameter.
mu: Location parameter.
beta: Scale parameter.
gamma: Shape parameter.
n: Shape parameter.
Expand All @@ -25,25 +25,25 @@ def erfexp_pdf_func(x, alpha, beta, gamma, n):
Implementation from https://gitlab.cern.ch/cms-muonPOG/spark_tnp/-/blob/Spark3/RooErfExp.cc
The parameters beta and gamma are given in reverse order in this c++ implementation.
"""
return tf.math.erfc((x - alpha) * beta) * znp.exp(-gamma * (znp.power(x, n) - znp.power(alpha, n)))
return tf.math.erfc((x - mu) * beta) * znp.exp(-gamma * (znp.power(x, n) - znp.power(mu, n)))


# Note: There is no analytic integral for the ErfExp PDF
# We tried with sympy, Mathematica, Wolfram Alpha and https://www.integral-calculator.com/
# import sympy as sp
#
# # Define symbols
# x, alpha, beta, gamma, n = sp.symbols('x alpha beta gamma n', real=True)
# x, mu, beta, gamma, n = sp.symbols('x mu beta gamma n', real=True)
#
# # Define the function
# func = sp.erfc((x - alpha) * beta) * sp.exp(-gamma * (x**n - alpha**n))
# func = sp.erfc((x - mu) * beta) * sp.exp(-gamma * (x**n - mu**n))
# sp.integrate(func, x)
class ErfExp(zfit.pdf.BasePDF):
_N_OBS = 1

def __init__(
self,
alpha: ztyping.ParamTypeInput,
mu: ztyping.ParamTypeInput,
beta: ztyping.ParamTypeInput,
gamma: ztyping.ParamTypeInput,
n: ztyping.ParamTypeInput,
Expand All @@ -59,10 +59,10 @@ def __init__(
.. math:
f(x; \\alpha, \\beta, \\gamma, n) = \\text{erfc}(\\beta (x - \\alpha)) \\exp{(-\\gamma (x^n - \\alpha^n))}
f(x; \\mu, \\beta, \\gamma, n) = \\text{erfc}(\\beta (x - \\mu)) \\exp{(-\\gamma (x^n - \\mu^n))}
Args:
alpha: Location parameter.
mu: Location parameter.
beta: Scale parameter.
gamma: Shape parameter.
n: Shape parameter.
Expand All @@ -87,13 +87,13 @@ def __init__(
the PDF for better identification.
Has no programmatical functional purpose as identification. |@docend:pdf.init.name|
"""
params = {"alpha": alpha, "beta": beta, "gamma": gamma, "n": n}
params = {"mu": mu, "beta": beta, "gamma": gamma, "n": n}
super().__init__(obs=obs, params=params, extended=extended, norm=norm)

def _unnormalized_pdf(self, x):
alpha = self.params["alpha"]
mu = self.params["mu"]
beta = self.params["beta"]
gamma = self.params["gamma"]
n = self.params["n"]
x = z.unstack_x(x)
return erfexp_pdf_func(x=x, alpha=alpha, beta=beta, gamma=gamma, n=n)
return erfexp_pdf_func(x=x, mu=mu, beta=beta, gamma=gamma, n=n)

0 comments on commit ad1be34

Please sign in to comment.