TD7: Test du \( \chi^2 \) (RICM4, 2013)

Application directe

Lancer de dé

Testons si notre dé est pipé ou pas.

freq = c(22, 21, 22, 27, 22, 36)
probs = rep(1/6, 6)
print(chisq.test(freq, p = probs))

## 
##  Chi-squared test for given probabilities
## 
## data:  freq
## X-squared = 6.72, df = 5, p-value = 0.2423

On obtient une p-value de 0.2433, ce qui signifie que dans une loi du \( \chi^2 \) à 5 degrés de liberté, on a 24.33% de chance d'obtenir une valeur supérieure au 6.72 qui correspond à cet échantillon. Rien d'exceptionnel donc et on ne peut donc rejeter pas rejeter l'hypothèse d'uniformité.

En anglais dans le texte

Testons s'il est raisonnable de supposer que la fréquence d'apparition des lettres E, T, N, R, O dans notre texte corresponde à celle de la langue anglaise:

x = c(100, 110, 80, 55, 14)
probs = c(29, 21, 17, 17, 16)/100
print(chisq.test(x, p = probs))

## 
##  Chi-squared test for given probabilities
## 
## data:  x
## X-squared = 55.4, df = 4, p-value = 2.685e-11

On obtient une p-value de 2.685e-11, ce qui signifie que dans une loi du \( \chi^2 \) à 4 degrés de liberté, on n'a quasiment aucune chance d'obtenir une valeur aussi grande que le 55.3955 qui correspond à cet échantillon. Il y a donc peu de chance pour notre texte soit en anglais…

Goodness of fit

Commençons par récupérer notre fonction calculant la dynamique de population.

set.seed(1:6)
city_evolution <- function(n = 500, p1 = 1, p2 = 1) {
    n1 <- c(p1)
    n2 <- c(p2)
    for (t in 1:(n - 1)) {
        if (runif(1) < (n1[t]/(n1[t] + n2[t]))) {
            n1[t + 1] = n1[t] + 1
            n2[t + 1] = n2[t]
        } else {
            n2[t + 1] = n2[t] + 1
            n1[t + 1] = n1[t]
        }
    }
    data.frame(t = (1:n), n1 = n1, n2 = n2)
}

limit_sample <- function(trials = 500, p1 = 1, p2 = 1, n = 500) {
    limit = c()
    for (trial in 1:trials) {
        d = city_evolution(n, p1, p2)
        limit[trial] = d[n, ]$n1/(d[n, ]$n1 + d[n, ]$n2)
    }
    limit
}

Maitenant, calculons la loi limite pour 500 trajectoires

lim = limit_sample(trials = 500, p1 = 1, p2 = 1)

Et faisons notre test du \( chi^2 \):

limc = hist(lim)$counts

r = length(limc)
print(chisq.test(limc, p = rep(1/r, r)))

## 
##  Chi-squared test for given probabilities
## 
## data:  limc
## X-squared = 8.36, df = 9, p-value = 0.4983

Une p-value de 0.4983, c'est énorme. Impossible de rejeter l'hypothèse que ça corresponde à une loi uniforme.

Maitenant, faisons pareil en partant de (10,10)

lim = limit_sample(trials = 500, p1 = 10, p2 = 10)

Et faisons notre test du \( chi^2 \):

limc = hist(lim, xlim = c(0, 1))$counts

r = length(limc)
print(chisq.test(limc, p = rep(1/r, r)))

## 
##  Chi-squared test for given probabilities
## 
## data:  limc
## X-squared = 261.5, df = 11, p-value < 2.2e-16

Une p-value de 1.1018 × 10^-49, là, c'est minuscule et donc on peu rejeter l'hypothèse d'uniformité sans risque. En même temps, vu la tête de l'histogramme, c'est logique.

Un certain nombre d'entre vous ont émis l'hypothèse qu'il s'agisse d'une loi normale centrée en \( 1/2 \) et de variance 0.0115.

br = hist(lim)$breaks

br[1] = -1e+05
br[length(br)] = 1e+05
pb = pnorm(br, mean = 0.5, sd = sd(lim))
print(chisq.test(limc, p = pb[2:length(pb)] - pb[1:length(pb) - 1]))

## Warning: Chi-squared approximation may be incorrect

## 
##  Chi-squared test for given probabilities
## 
## data:  limc
## X-squared = 13.16, df = 11, p-value = 0.2827

Ah et il se plaint car si on regarde limc, on voit que certaines cases sont trop peu remplies.

limc

##  [1]  2 10 30 41 75 84 89 66 54 29 18  2

Je vais donc recommencer en forçant le nombre de breaks à une valeur plus faible.

h = hist(lim, xlim = c(0, 1), breaks = 6)

limc = h$count
br = h$breaks
br[1] = -1e+05
br[length(br)] = 1e+05
pb = pnorm(br, mean = 0.5, sd = sd(lim))
print(chisq.test(limc, p = pb[2:length(pb)] - pb[1:length(pb) - 1]))

## 
##  Chi-squared test for given probabilities
## 
## data:  limc
## X-squared = 4.071, df = 5, p-value = 0.5393

On ne peut toujours pas conclure… En même temps, vu la forme générale et le faible nombre de bins, c'est un peu normal.

Une dernière tentative, cette fois-ci en fusionnant les bins faiblement peuplés:

lim = limit_sample(trials = 1000, p1 = 10, p2 = 10)
h = hist(lim, xlim = c(0, 1))

limc = h$count
br = h$breaks
br = br[2:(length(br) - 1)]
br[1] = -1e+05
br[length(br)] = 1e+05
limc = hist(lim, breaks = br, plot = F)$count
pb = pnorm(br, mean = 0.5, sd = sd(lim))
print(chisq.test(limc, p = pb[2:length(pb)] - pb[1:length(pb) - 1]))

## 
##  Chi-squared test for given probabilities
## 
## data:  limc
## X-squared = 6.868, df = 10, p-value = 0.7378

Bon, et bien toujours pas possible de rejeter l'hypothèse. Vous aviez raison de dire que ça ressemblait à une normale, même si ça n'en est pas une… :) Il faudrait pousser le nombre d'échantillons super loin pour pouvoir commencer à voir la différence.