Probability Theory: The Logic of Science, chapter 9.12 example

janyes_9_12.py

import numpy as np


H1 = [0.499, 0.499, 0.002]
H2 = [0.3333333333, 0.3333333333, 0.3333333333]
obs = [14, 14, 1]


def psi(hypothesis, observed):
    assert np.isclose(sum(hypothesis), 1)

    statistic = 0
    n_obs = sum(observed)
    for o, p in zip(observed, hypothesis):
        statistic += o * np.log10((o / (n_obs * p)))
    return 10 * statistic


def chi2(hypothesis, observed):
    assert np.isclose(sum(hypothesis), 1)

    statistic = 0
    n_obs = sum(observed)
    for o, p in zip(observed, hypothesis):
        statistic += (o - n_obs * p)**2 / (n_obs * p)
    return statistic


print(psi(H1, obs), 'db')
print(psi(H2, obs), 'db')

print(chi2(H1, obs))
print(chi2(H2, obs))