Plotting Confidance Intervals

example.py

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy import stats

# thanks to
# https://codereview.stackexchange.com/questions/84414/obtaining-prediction-bands-for-regression-model


# an example of using matplotlib to draw confidence bands

def confidence_band(x, y, covMat, dfdp, confprob):
    """ given a value for x calculate the error df in y = model(p,x)
    returns lower and upper bounds for each value in x

    dfdp - list of derivatives
    confprob - the confidence probability
    """

    # if the confidence probability is given out of 100 then we convert to a percentage
    if confprob > 1.0:
        confprob = confprob / 100.00

    # we derive alpha = 1 - confprob
    alpha = 1.0 - confprob
    prb = 1.0 - (alpha / 2.0)
    tval = stats.t.ppf(prb, len(dfdp))

    df2 = np.zeros(len(x))
    for j in range(len(dfdp)):
        for k in range(len(dfdp)):
            df2 += dfdp[j] * dfdp[k] * covMat[j, k]

    df = np.sqrt(df2)
    delta = tval * df

    upper_band = y + delta
    lower_band = y - delta

    return lower_band, upper_band


def exp_func(x, a, beta, c):
    '''Exponential 3-param function.'''
    return a * np.exp(-beta * x) + c


def exp_fit(ax, x, y, yerr):

    # it is important that we use the absolute sigma
    # the manner in which we construct the confidence band assumes that our
    # covariance matrix is the 'absolute' version
    popt, pcov = curve_fit(exp_func, x, y, sigma=yerr, absolute_sigma=True)

    # generate a range of x values for a smooth fit
    x_min, x_max = ax.get_xlim()
    new_x = np.linspace(x_min, x_max, num=200)

    # do the fit
    fit = exp_func(new_x, *popt)
    ax.plot(new_x, fit, '--', color='red', label="Least Squares exponential fit")

    # for the exponential function we need to figure out the derivatives by hand
    dfdp = [
        # df/da = e^{-b * x}
        np.exp(-popt[1] * new_x),
        # df/db = -x * a * e^{-b * x}
        -new_x * popt[0] * np.exp(-popt[1] * new_x),
        # df/dc = 1
        1,
    ]

    # draw the 1 sigma interval
    confprob = 68.27
    fill_label = fr'''$1\sigma$ or ${confprob}\%$ confidence interval'''
    lb, ub = confidence_band(new_x, fit, pcov, dfdp, confprob)
    ax.fill_between(new_x, lb, ub, color='k', alpha=.30, label=fill_label)

    # draw the 2 sigma interval
    confprob = 95.45
    fill_label = fr'''$2\sigma$ or ${confprob}\%$ confidence interval'''
    lb, ub = confidence_band(new_x, fit, pcov, dfdp, confprob)
    ax.fill_between(new_x, lb, ub, color='tab:cyan', alpha=.10, label=fill_label)


def poly_func(x, a, b, c):
    '''Polynomial 3-param function.'''
    return a + (b * x) + (c * x ** 2)


def poly_fit(ax, x, y, yerr, degree=2):

    # it is important that we use the absolute sigma
    # the manner in which we construct the confidence band assumes that our
    # covariance matrix is the 'absolute' version
    popt, pcov = curve_fit(poly_func, x, y, sigma=yerr, absolute_sigma=True)

    # generate a range of x values for a smooth fit
    x_min, x_max = ax.get_xlim()
    new_x = np.linspace(x_min, x_max, num=200)

    fit = poly_func(new_x, *popt)

    poly_label = f"Least Squares polynomial fit of degree {degree:d}"
    ax.plot(new_x, fit, '--', color='red', label=poly_label)

    # here I explicitly show the derivatives for this example
    dfdp = [
        # df/da = 1
        np.exp(-popt[1] * new_x),
        # df/db = x
        -new_x * popt[0] * np.exp(-popt[1] * new_x),
        # df/dc = x^2
        1,
    ]

    # alternatively we could use list comprehension to generate the derivatives
    # for any polynomial function
    dfdp = [new_x**n for n in range(0, degree + 1)]

    # draw the 1 sigma interval
    confprob = 68.27
    fill_label = fr'''$1\sigma$ or ${confprob}\%$ confidence interval'''
    lb, ub = confidence_band(new_x, fit, pcov, dfdp, confprob)
    ax.fill_between(new_x, lb, ub, color='k', alpha=.30, label=fill_label)

    # draw the 2 sigma interval
    confprob = 95.45
    fill_label = fr'''$2\sigma$ or ${confprob}\%$ confidence interval'''
    lb, ub = confidence_band(new_x, fit, pcov, dfdp, confprob)
    ax.fill_between(new_x, lb, ub, color='tab:cyan', alpha=.20, label=fill_label)


def generate_fake_data(n_data_points):

    # a few seeds that seem to give reasonable random data on my machine
    reasonable_seeds = [
        2860442663,
        2968839790,
        3251284835,
        1253373634,
    ]

    np.random.seed(1253373634)

    # generate some fake data
    x = np.sort(np.random.random_sample(n_data_points))
    y = np.random.random_sample(n_data_points)
    yerr = np.random.random_sample(n_data_points)

    return x, y, yerr


def plot_example():

    N = 35
    x, y, yerr = generate_fake_data(N)

    # create a new figure and axis upon which we will plot the polynomial fit
    fig, ax = plt.subplots(1, 1)

    # plot the fake data
    ax.errorbar(x, y, xerr=None, yerr=yerr,
                marker='o', mfc="None", color='orange', ecolor='k',
                label="Data we are fitting")

    # fit the fake data
    poly_fit(ax, x, y, yerr)
    ax.legend()
    mpl.interactive(True)
    plt.show()

    # generate some data from the exponential distribution
    y = np.sort(np.random.exponential(size=N))

    # create a new figure and axis upon which we will plot the exponential fit
    fig, ax = plt.subplots(1, 1)

    # plot the fake data
    ax.errorbar(x, y, xerr=None, yerr=yerr,
                marker='o', mfc="None", color='orange', ecolor='k',
                label="Data we are fitting")

    # fit the fake data
    exp_fit(ax, x, y, yerr)
    ax.legend()
    mpl.interactive(False)
    plt.show()


if (__name__ == "__main__"):
    plot_example()