G-test - multinomial goodness-of-fit

2022-03-08 · 3 min read

    \[ \def\Chi{\mathcal{X}} \]

    Overview #

    • $G$ is a test-statistic. Asymptotically the same as the $\Chi^2$ test-statistic.
    • Used for goodness-of-fit test b/w observed and hypothesized multinomial distributions.
    • We use the $\Chi^2$ distribution to compute our desired $G$ threshold from a given $p$-value.
    • It's a log-likelihood ratio. Closely related to KL-divergence and $\Chi^2$ statistic.

    Test Statistic #

    In terms of observed counts, $O_i \ge 0$, and expected counts, $E_i > 0$

    \[ G = 2 \sum_i O_i \cdot \ln\Bigl( \frac{O_i}{E_i} \Bigr) \]

    In terms of probabilities $\hat{p}_i$ and $p_i$ (estimated and expected resp.) and number of samples $n$:

    \[ \begin{align*} G &= 2 \sum_i \hat{p}_i \cdot n \cdot \ln\Bigl( \frac{\hat{p}_i \cdot n}{p_i \cdot n} \Bigr) \\ &= 2 \sum_i \hat{p}_i \cdot n \cdot \ln\Bigl( \frac{\hat{p}_i}{p_i} \Bigr) \\ &= 2 n \sum_i \hat{p}_i \cdot \ln\Bigl( \frac{\hat{p}_i}{p_i} \Bigr) \\ &= 2 \, n \; D_{\text{KL}}(\hat{p} \;||\; p) \end{align*} \]

    where $D_{\text{KL}}(\hat{p} \;||\; p)$ is the Kullback-Leibler divergence.

    Example Implementation #

    Source: https://github.com/phlip9/kcddice/blob/master/kcddice/src/stats.rs

    use claim::{debug_assert_ge, debug_assert_le};
    use ndarray::{ArrayView1, Zip};
    use statrs::distribution::{ChiSquared, ContinuousCDF};
    
    const EPS: f64 = 1e-10;
    
    /// Return true iff `supp(p) ⊆ supp(q)` for dense PMFs `p` and `q`.
    fn is_pmf_subset(p: ArrayView1<f64>, q: ArrayView1<f64>) -> bool {
        Zip::from(p).and(q).all(|&p_i, &q_i| {
            // A = (q_i == 0.0)
            // B = (p_i == 0.0)
            // (A ==> B) <==> (¬A ∨ B)
            (q_i > 0.0) || (p_i <= 0.0)
        })
    }
    
    /// Compute the [KL-divergence](https://www.wikiwand.com/en/Kullback%E2%80%93Leibler_divergence).
    /// between dense PMFs `p` and `q`.
    ///
    /// `D_{KL}(p || q) = \sum_i p_i * \ln(p_i / q_i)`
    ///
    /// Note: p's support must be a strict subset of q's support!
    ///       this is also referred to as the "absolute continuity" assumption,
    ///       where `q_i = 0` implies `p_i = 0`.
    fn kl_divergence(p: ArrayView1<f64>, q: ArrayView1<f64>) -> f64 {
        // caller should check this before
        debug_assert!(is_pmf_subset(p, q));
    
        Zip::from(p)
            .and(q)
            .fold(0.0, |sum, &p_i, &q_i| sum + kl_div_term(p_i, q_i))
    }
    
    #[inline]
    fn kl_div_term(p_i: f64, q_i: f64) -> f64 {
        if q_i > EPS {
            if p_i > EPS {
                p_i * (p_i / q_i).ln()
            } else {
                0.0
            }
        } else {
            if p_i > EPS {
                debug_assert!(false);
                f64::INFINITY
            } else {
                0.0
            }
        }
    }
    
    /// The G-test statistic.
    ///
    /// * Used for comparing observed multinomial distribution with expected
    ///   hypothesis multinomial distribution.
    /// * Asmyptotically approximates chi^2-test statistic.
    ///
    /// `n`: the number of samples
    /// `p`: the expected PMF
    /// `p_hat`: the observed PMF
    ///
    /// G-test: https://www.wikiwand.com/en/G-test
    fn g_test(n: usize, p: ArrayView1<f64>, p_hat: ArrayView1<f64>) -> f64 {
        (n as f64) * (2.0 * kl_divergence(p_hat, p))
    }
    

    See also: https://github.com/msberends/AMR/blob/main/R/g.test.R