Closed form A/B/C/D tests

binary_outcomes.py

from numpy import log, exp
from scipy.special import betaln as logbeta


def probability_B_beats_A(α_A, β_A, α_B, β_B):
    total = 0.0
    for i in range(α_B):
        total += exp(
            logbeta(α_A + i, β_B + β_A)
            - log(β_B + i)
            - logbeta(1 + i, β_B)
            - logbeta(α_A, β_A)
        )
    return total


def probability_C_beats_A_and_B(α_A, β_A, α_B, β_B, α_C, β_C):
    total = 0.0
    for i in range(α_A):
        for j in range(α_B):
            total += exp(
                logbeta(α_C + i + j, β_A + β_B + β_C)
                - log(β_A + i)
                - log(β_B + j)
                - logbeta(1 + i, β_A)
                - logbeta(1 + j, β_B)
                - logbeta(α_C, β_C)
            )
    return (
        1
        - probability_B_beats_A(α_C, β_C, α_A, β_A)
        - probability_B_beats_A(α_C, β_C, α_B, β_B)
        + total
    )


def probability_D_beats_ABC(α_A, β_A, α_B, β_B, α_C, β_C, α_D, β_D):
    total = 0.0
    for i in range(α_A):
        for j in range(α_B):
            for k in range(α_C):
                total += exp(
                    logbeta(α_D + i + j + k, β_A + β_B + β_C + β_D)
                    - log(β_A + i)
                    - log(β_B + j)
                    - log(β_C + k)
                    - logbeta(1 + i, β_A)
                    - logbeta(1 + j, β_B)
                    - logbeta(1 + k, β_C)
                    - logbeta(α_D, β_D)
                )

    return (
        1
        - probability_B_beats_A(α_A, β_A, α_D, β_D)
        - probability_B_beats_A(α_B, β_B, α_D, β_D)
        - probability_B_beats_A(α_C, β_C, α_D, β_D)
        + probability_C_beats_A_and_B(α_A, β_A, α_B, β_B, α_D, β_D)
        + probability_C_beats_A_and_B(α_A, β_A, α_C, β_C, α_D, β_D)
        + probability_C_beats_A_and_B(α_B, β_B, α_C, β_C, α_D, β_D)
        - total
    )

via Evan Miller