Probability interval boundaries with Python

Probability interval boundaries with Python

In the previous post we calculated proabilities and z values with Python. Here we continue with finding the interval boundaries belonging to a particular tail or two-sided probabiltiy.

Libraries and helper functions

import math as m
import numpy as np

def to_string(number: float, column_width: int = 20) -> str:
    # return str(number).ljust(column_width)
    return f"{str(number): <{column_width}}"

def calc_normal_cdf(x: float, mu: float = 0, sigma: float = 1) -> float:
    return (1 + m.erf((x - mu) / m.sqrt(2) / sigma)) / 2

def calc_inverse_normal_cdf(p: float, mu:float = 0, sigma: float = 1, tolerance: float = 1E-5, show_steps=False) -> float:

    if p == 0: return -np.inf
    if p == 1: return np.inf

    # In case it is not a standard normal distribution, calculate the standard normal first and then rescale
    if mu != 0 or sigma != 1:
        return mu + sigma * calc_inverse_normal_cdf(p, tolerance=tolerance)

    low_z = -10
    hi_z = 10

    if show_steps: print(f"{'':<19}".join(['low_z', 'mid_z', 'hi_z']), "\n")

    while hi_z - low_z > tolerance:
        mid_z = (low_z + hi_z) / 2
        mid_p = calc_normal_cdf(mid_z)

        if mid_p < p:
            low_z = mid_z
        else:
            hi_z = mid_z

        if show_steps: print("\t".join(map(to_string, [low_z, mid_z, hi_z])))

    return mid_z

Upper bound

def normal_upper_bound(probabilty: float, mu: float = 0, sigma: float = 1) -> float:
    """Return z for which P(Z <= z) = probability"""
    return calc_inverse_normal_cdf(probabilty, mu, sigma)

def normal_lower_bound(probabilty: float, mu: float = 0, sigma: float = 1) -> float:
    """Return z for which P(Z >= z) = probability"""
    return calc_inverse_normal_cdf(1 - probabilty, mu, sigma)

for i in np.arange(0, 1.1, .1):
    print(f"{i:.1f}: {normal_upper_bound(i): >10.4f}, {normal_lower_bound(i): >10.4f}")

0.0:       -inf,        inf
0.1:    -1.2816,     1.2816
0.2:    -0.8416,     0.8416
0.3:    -0.5244,     0.5244
0.4:    -0.2533,     0.2533
0.5:    -0.0000,    -0.0000
0.6:     0.2533,    -0.2533
0.7:     0.5244,    -0.5244
0.8:     0.8416,    -0.8416
0.9:     1.2816,    -1.2816
1.0:        inf,       -inf

Two sided bounds

Two variations

def normal_two_sided_bounds(probability: float, mu: float = 0, sigma: float = 1) -> float:
    if probability == 0: return 0, 0

    tail_probability = (1 - probability) / 2
    
    lower_bound = normal_upper_bound(tail_probability, mu, sigma)
    upper_bound = normal_lower_bound(tail_probability, mu, sigma)

    return lower_bound, upper_bound

def normal_two_sided_bounds2(probability: float, mu: float = 0, sigma: float = 1) -> float:
    if probability == 0: return 0, 0

    if mu != 0:
        raise ValueError("Requires standard normal distribution")

    tail_probability = (1 - probability) / 2
    z = normal_upper_bound(tail_probability, mu, sigma)
    
    return -z, z

By default they do not produce the same results at the extremes because of how the inverse cdf function defines the theoretical uppern and lower boundaries (as in these cases it should produce +/- infinity).

Accordingly, it is useful to define extreme values (i.e. probablities of 0 and 1)

for i in np.arange(0, 1.1, .1):
    print(f"{i:.1f}:  {normal_two_sided_bounds(i) == normal_two_sided_bounds2(i)}, {normal_two_sided_bounds(i)}, {normal_two_sided_bounds2(i)}")

0.0:  True, (0, 0), (0, 0)
0.1:  False, (-0.12566566467285156, 0.12566566467285156), (0.12566566467285156, -0.12566566467285156)
0.2:  False, (-0.2533435821533203, 0.2533435821533203), (0.2533435821533203, -0.2533435821533203)
0.3:  False, (-0.3853130340576172, 0.3853130340576172), (0.3853130340576172, -0.3853130340576172)
0.4:  False, (-0.5243968963623047, 0.5243968963623047), (0.5243968963623047, -0.5243968963623047)
0.5:  False, (-0.6744861602783203, 0.6744861602783203), (0.6744861602783203, -0.6744861602783203)
0.6:  False, (-0.8416271209716797, 0.8416271209716797), (0.8416271209716797, -0.8416271209716797)
0.7:  False, (-1.0364246368408203, 1.0364246368408203), (1.0364246368408203, -1.0364246368408203)
0.8:  False, (-1.2815570831298828, 1.2815570831298828), (1.2815570831298828, -1.2815570831298828)
0.9:  False, (-1.6448497772216797, 1.6448497772216797), (1.6448497772216797, -1.6448497772216797)
1.0:  False, (-inf, inf), (inf, -inf)