Jeg åpnet en pose twist og fikk resultatet ovenfor. Ved første øyekast er det mistenkelig få typer av enkelte sjokolader.

Vi teller opp antall av hver type sjokolade:

På Wikipedias artikkel om twist står det at “variasjonen av hvor mange av hver bit, er helt tilfeldig.” Det er uklart om det er (1) tilfeldig i hver pose, men med ulik sannsynlighet for hver type, eller (2) om artikkelforfatteren mener at det er samme sannsynlighet for hver type.

Gir observasjonene ovenfor grunnlag til å påstå at det er ulik sannsynlighet for hver type? For å svare på dette, kan vi utføre en G-test. G-verdien er gitt av en log-likelihood ratio test, der den underliggende sannsynlighetsfordelingen er multinomial. Formelen for \(G\) er

der \(O_i\) er observationer og \(E_i\) er forventede observasjoner under nullhypotesen \(H_0\). Observert \(G\) for vårt datasett er \(8.3\).

Vi simulerer ti tusen ganger fra en multinomial fordeling for å få distribusjonen til \(G\) under \(H_0\).

• Laveste mulige \(G\) er her \(0.93\), som oppstår når vi får \(4\) observasjoner av \(6\) typer sjokolade og \(3\) observasjoner av de resterende \(7\) typene.
• Høyeste mulige \(G\) er her \(230.85\), som oppstår når vi får \(45\) observasjoner av \(1\) type sjokolade og \(0\) observasjoner av de resterende \(12\) typene.

En figur som viser distribusjon av \(G\) under \(H_0\) ser slik ut:

Her er \(p\)-verdien arealet til høyre for observert \(G=8.3\), som gir oss en \(p\)-verdi på \(0.7984\). Det er altså ingen grunn til å forkaste \(H_0\), det er faktisk svært sannsynlig å få mer ekstreme fordelinger (og \(G\)-verdier) enn det vi observerte.

Observasjonene gir ikke grunnlag til å påstå at det er ulik sannsynlighet for hver type.

import numpy as np
from scipy.special import rel_entr
from scipy.stats import multinomial


def G_statistic(observed, expected):
    return 2 * rel_entr(observed, expected)


# The observed counts
observed = np.array([1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5, 6])

# Compute observed test statistic
expected = np.ones_like(observed) * np.mean(observed)
G_observed = G_statistic(observed, expected).sum()

# Parameters for multinomial (under H_0)
n = np.sum(observed)
p = np.ones_like(observed) / len(observed)

# Simulate draws from multinomial under H_0
simulations = 10_000
simulation_observed = multinomial(n, p, seed=42).rvs(simulations)
simulation_expected = np.outer(np.ones(simulations), expected)

# Compute test statistic for each draw under H_0
G_simulated = G_statistic(simulation_observed, simulation_expected).sum(axis=1)

print(
    "Probability of observed or more extreme, given H_0:",
    (G_simulated >= G_observed).mean(),
)

Jeg fikk i etterkant tilsendt resultatene fra en ny pose twist:

La oss slå sammen tellingene ved å endre koden til:

bag1 = np.array([1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5, 6])
bag2 = np.array([3, 2, 2, 2, 3, 5, 4, 3, 3, 5, 4, 4, 4])
observed = bag1 + bag2

Nå får vi en \(p\)-verdi på \(0.6638\). Det er fremdeles ingen grunn til å påstå at det er ulik sannsynlighet for hver type.