Correction de l'exercice 1 de la section 9.2, page 16

Partie (a)

• On considère une population dont la fréquence d'un gène pour le daltonisme (color blindness) est θ dans la population mâle (♂) et θ² dans la population femelle (♀).
• Dans un échantillon aléatoire de M individus ♂, on en trouve m daltoniens et dans un échantillon aléatoire de N individus ♀, on en trouve n daltoniens.
• On cherche la fonction de vraisemblance basée sur les deux échantillons.
• Nous devons montrer que l'EMV $\hat{θ}$ est solution d'une équation quadratique.

Solution

• Pour chacun des deux échantillons, un modèle probabiliste binomiale s'applique avec une probabilité de succés θ pour l'échantillon de ♂ et de θ² pour l'échantillon de ♀.
• Soit $\mathbb{P}$(m daltoniens parmi M|♂) = M!/(m! (M-m)!) θ^m(1-θ)^M-m et $\mathbb{P}$(n daltoniennes parmi N|♀) = N!/(n! (N-n)!) θ²ⁿ(1-θ²)^N-n.
• Et : L(θ) = c x $\mathbb{P}$(m daltoniens parmi M|♂) x $\mathbb{P}$(n daltoniennes parmi N|♀)
• D'où : L(θ) = θ^m+2n(1-θ)^M-m(1-θ²)^N-n.
• Soit : l(θ) = (m+2n) log θ + (M-m) log (1-θ) + (N-n) log (1-θ²).
• D'où : S(θ) = (m+2n) / θ - (M-m) / (1-θ) - 2 θ (N-n) / (1-θ²).

Utilisation d'un système de calcul formel

Nous allons profiter de la résolution de ce problème pour nous familiariser avec un système de calcul formel : le module SymPy du logiciel que nous allons utiliser dans ce cours, Python.

Nous commençons par « lancer » Python de la façon la plus simple, en tapant en ligne de commande :

python

Nous chargeons ensuite le module Sympy :

import sympy as sy

Nous définissons ensuite les variables M, m, N, n et theta comme des symboles avec :

M, m, N, n, theta = sy.symbols('M m N n theta')

Pour avoir l'impression la plus « jolie » possible avec notre terminal nous évaluons la fonction init_printing avec :

sy.init_printing()

Nous pouvons maintenant définir notre fonction de vraisemblance :

L = theta**(m+2*n)*(1-theta)**(M-m)*(1-theta**2)**(N-n)

Un appel à la version symbolique de la fonction log nous donne après un « petit » traitement supplémentaire une version « agréable » de la log-vraisemblance :

l = sy.expand_log(sy.log(L),force=True)

Le petit traitement est effectué par la fonction expand_log et le paramètre formel optionnel force=True dit au système de ne pas se poser de questions (sur le valeur de theta) et de faire le développement comme si il avait le droit de le faire.

Un appel à la fonction diff nous donne alors la fonction score :

S = sy.diff(l,theta)

La réponse à la première partie de la question est alors obtenue en combinant simplify et collect :

S = sy.collect(sy.simplify(S),theta)

Pour trouver l'expression analytique de $\hat{θ}$ nous devons extraire le numérateur avec la fonction fraction :

S_num = sy.fraction(S)[0]

Partie (b)

L'expression analytique de $\hat{θ}$ est finalement obtenue avec un appel à la fonction solve :

sy.solve(S_num,theta)

Dans l'énoncer de la partie (b) on nous dit que M = N = 100, m = 11 et n = 2. La méthode subs nous permet de substituer ces valeurs numériques dans nos expressions :

solution = [exp.subs([(M,100),(m,11),(N,100),(n,2)])
            for exp in sy.solve(S_num,theta)]

Ici, nous avons employé une liste en compréhension une possibilité très pratique de la syntaxe de Python. Pour « forcer » une conversion des expressions rationnelles obtenues vers des réels représenté en double précision, il faut utiliser la méthode evalf :

solution[1].evalf()

Si nous n'utilisons que les données de l'échantillon ♂, nous obtenons :

(m/M).subs([(M,100),(m,11)]).evalf()

Alors que si nous n'utilisons que les données de l'échantillon ♀, nous obtenons :

sy.sqrt(n/N).subs([(N,100),(n,2)]).evalf()