Calculer un chi 2 sur un tableau de contingence

\chi_2 et tableau de contingence, avec numpy, avec scipy ou sans.

formule

Le test du \chi_2 (wikipedia) sert à comparer deux distributions. Il peut être appliqué sur un tableau de contingence pour comparer la distributions observée avec la distribution qu’on observerait si les deux facteurs du tableau étaient indépendants. On note M=(m_{ij}) une matrice de dimension I \times J. Le test du \chi_2 se calcule comme suit :

  • M = \sum_{ij} m_{ij}

  • \forall i, \; m_{i \bullet} = \sum_j m_{ij}

  • \forall j, \; m_{\bullet j} = \sum_i m_{ij}

  • \forall i,j \; n_{ij} = \frac{m_{i \bullet} m_{\bullet j}}{N}

Avec ces notations :

T = \sum_{ij} \frac{ (m_{ij} - n_{ij})^2}{n_{ij}}

La variable aléatoire T suit asymptotiquement une loi du \chi_2 à (I-1)(J-1) degrés de liberté (table). Comment le calculer avec numpy ?

tableau au hasard

On prend un petit tableau qu’on choisit au hasard, de préférence non carré pour détecter des erreurs de calculs.

[1]:
import numpy

M = numpy.array([[4, 5, 2, 1], [6, 3, 1, 7], [10, 14, 6, 9]])
M
[1]:
array([[ 4,  5,  2,  1],
       [ 6,  3,  1,  7],
       [10, 14,  6,  9]])

calcul avec scipy

Evidemment, il existe une fonction en python qui permet de calculer la statistique T : chi2_contingency.

[2]:
from scipy.stats import chi2_contingency

chi2, pvalue, degrees, expected = chi2_contingency(M)
chi2, degrees, pvalue
[2]:
(6.168598503892621, 6, 0.4045712090580829)

calcul avec numpy

[3]:
N = M.sum()
ni = numpy.array([M[i, :].sum() for i in range(M.shape[0])])
nj = numpy.array([M[:, j].sum() for j in range(M.shape[1])])
ni, nj, N
[3]:
(array([12, 17, 39]), array([20, 22,  9, 17]), 68)

Et comme c’est un usage courant, numpy propose une façon de faire sans écrire une boucle avec la fonction sum :

[4]:
ni = M.sum(axis=1)
nj = M.sum(axis=0)
ni, nj, N
[4]:
(array([12, 17, 39]), array([20, 22,  9, 17]), 68)
[5]:
nij = ni.reshape(M.shape[0], 1) * nj / N
nij
[5]:
array([[ 3.52941176,  3.88235294,  1.58823529,  3.        ],
       [ 5.        ,  5.5       ,  2.25      ,  4.25      ],
       [11.47058824, 12.61764706,  5.16176471,  9.75      ]])
[6]:
d = (M - nij) ** 2 / nij
d.sum()
[6]:
6.168598503892621
[ ]:


Notebook on github