Compute relative error for the relativistic Breit-Wigner CDF calculation.

check_rel_breitwigner_cdf.py

from mpmath import mp
from mpsci.distributions import rel_breitwigner

import numpy as np
from scipy import stats

import matplotlib.pyplot as plt


def rel_error(value, reference):
    if value == reference:
        return 0.0
    if reference == 0:
        return 1.0
    return abs((value - reference)/reference)


def check_cdf(rho=1, n=10000, start=1e-12, stop=1e50, dps=250):
    # Check the relative error of rel_breitwigner.cdf.
    x = np.geomspace(start, stop, n)
    sp_cdf = stats.rel_breitwigner.cdf(x, rho)
    with mp.workdps(dps):
        mp_cdf = [rel_breitwigner.cdf(t, rho, 1) for t in x]
        sp_relerr = np.array([float(rel_error(p, mp))
                              for p, mp in zip(sp_cdf, mp_cdf)])
    return x, sp_relerr


if __name__ == "__main__":
    rho = 100
    print(f'{rho = }')

    x, relerr = check_cdf(rho=rho, n=50000,
                          start=1e-12, stop=1e20, dps=250)

    alpha = 0.35
    plt.plot(x, relerr, '.', alpha=alpha, label=f'{rho = } (scipy PR)',
             markersize=4)
    plt.title(f'CDF rel. error, {rho = }')
    plt.semilogx()
    plt.grid()

    plt.xlabel('x')
    plt.ylabel('CDF relative error')
    plt.show()