Power Curve

The following is a simple code to generate a power curve.


In [2]:
from __future__ import division

from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import prettyplotlib as ppl
import brewer2mpl

mpl.rcParams['lines.linewidth'] = 2
mpl.rcParams['lines.color'] = 'r'
mpl.rcParams['axes.titlesize'] = 32
mpl.rcParams['axes.labelsize'] = 24 
mpl.rcParams['axes.labelsize'] = 24 
mpl.rcParams['xtick.labelsize'] = 24 
mpl.rcParams['ytick.labelsize'] = 24 

colors = brewer2mpl.get_map('Set1', 'qualitative', 9).mpl_colors

def plot_power_curve(alpha, stddev, sample_sizes, alternative=">"):
    if alternative == ">":
        z_crit = norm.interval(1 - alpha * 2)[1]
        x = np.arange(0, 14, 0.01)
    elif alternative == "<":
        x = - np.arange(0, 14, 0.01)
        z_crit = norm.interval(1 - alpha * 2)[0]
    plt.figure(figsize=(12, 8))
    for i, n in enumerate(sample_sizes):
        z = (-x / (stddev / np.sqrt(n))) + z_crit
        z = -z if alternative == "<" else z
        ppl.plot(x, 1 - norm.cdf(z), color=colors[i], linewidth=3, label="n = %d" % n)
    plt.legend(loc=4, fontsize=24)
    plt.xlabel("Difference ($\mu_{hyp} - \mu_{true}$)")
    plt.ylabel("Power")
    plt.title("Power curve for $\sigma = %d$ and $\\alpha = %.2f$" % (stddev, alpha))
    plt.tight_layout()
    plt.ylim([0, 1.1])
    return plt

In [4]:
plt = plot_power_curve(0.05, 15, [20, 30, 50, 100, 1000])


Beta analysis


In [19]:
x = np.arange(80, 200, 1)
plt.plot(x, norm.pdf(x, 160, 34/np.sqrt(10)), color=colors[0], label="hypothesized")
plt.plot(x, norm.pdf(x, 110, 34/np.sqrt(10)), color=colors[1], label="true")
plt.legend()
x = np.arange(140.9, 180)
plt.fill_between(x, norm.pdf(x, 110, 34/np.sqrt(10)), color='grey')


Out[19]:
<matplotlib.collections.PolyCollection at 0x10a02cb50>