Mise en route de Python¶
Comme d'« habitude », les exemples qui suivent supposent que vous faites tourner IPython, soit avec une interface « toute bête » :
{bash}
ipythonsoit avec la version plus sophistiquée via votre navigateur oueb favori :
{bash}
ipython notebookQuelque soit l'interface choisie, vous employez la commande magique %pylab pour avoir sous la main les graphiques interactifs, les fonctions mathématiques « avancées », etc :
%pylab
Pour la cas, sur lequel je ne compte plus trop, où vous auriez Sympy à disposition, vous chargez aussi le module :
import sympy as sy
sy.init_printing()
Si vous ne l'avez pas, installez une version de Python sur votre machine favorite, par exemple la version anaconda, et refaites les parties dépendant de Sympy en révisant le cours…
Je suppose que Python 2 est employé et pour avoir des appelles à quelques fonctions clés comme la division, compatibles avec Python 3, nous évaluons :
from __future__ import print_function, division, unicode_literals, absolute_import
Exemple 1.2¶
Nous considérons ici des données de fiabilité de roulements à billes avec des durées de vie constatées (en millions de révolutions) :
durees_de_vie = [17.88, 28.92, 33.00, 41.52, 42.12, 45.60, \
48.48, 51.84, 51.96, 54.12, 55.56, 67.80, \
68.64, 68.64, 68.88, 84.12, 93.12, 98.64, \
105.12, 105.84, 127.92, 128.04, 173.40]
Une loi de Weibull est considérée comme un modèle raisonnable pour ces données ; avec la paramétrisation adoptée par Kalbfleisch (p 57), la densité de cette loi s'écrit : $$f(x) = λ β x^{β-1} \exp \{- λ x^β\}, \quad 0 < x < \infty, \quad λ, β > 0\, .$$
Comme les observations sont supposées indépendantes, la loi conjointe est le produit des lois individuelles et la log-vraisemblance s'écrit : $$l(λ,β) = n \log λ + n \log β + (β-1) \sum \log x_i - λ \sum x_i^β \, .$$ Les composants de la fonction score (à valeurs vectorielles) sont : $$S_1(λ,β) = \frac{\partial l}{\partial λ} = \frac{n}{λ} - \sum x_i^β$$ et $$S_2(λ,β) = \frac{\partial l}{\partial β} = \frac{n}{β} + \sum \log x_i - λ \sum x_i^β \log x_i \, .$$ L'équation $S_1(λ,β) = 0$ peut être résolue explicitement pour λ et donne : $$\hat{λ}(β) = \frac{n}{\sum \log x_i}\, .$$ Ceci est l'EMV de λ lorsque β est connu. Pour obtenir $\hat{β}$ nous substituons $\hat{λ}(β)$ dans la définition de $S_2$ et nous cherchons la racine de cette équation : $$g(β) = S_2(\hat{λ}(β),β) = \frac{n}{β} + \sum \log x_i - \frac{ n \sum x_i^β \log x_i}{\sum \log x_i} \, .$$ Nous allons employer la méthode de Newton pour résoudre cette équation : $$β_{nouveau} = β_{ancien} - \frac{g(β_{ancien})}{g'(β_{ancien})} \, .$$ Ici nous avons : $$g'(β) = - \frac{n}{β} - \frac{ n \sum x_i^β (\log x_i)^2}{\sum x_i^β} + \frac{ n (\sum x_i^β \log x_i)^2}{(\sum x_i^β)^2}\, .$$
Nous codons des fonctions Python qui renvoient les fonctions g et g' ci-dessus. En fait, nous allons définir des fermetures) ou clôtures :
def fait_g_dg(donnees):
donnees = array(donnees)
n = len(donnees)
somme = sum(log(donnees))
def g(beta):
terme1 = n/beta
terme2 = n*sum(donnees**beta*log(donnees))/sum(donnees**beta)
return terme1 + somme - terme2
def dg(beta):
terme1 = n/beta**2
terme2 = n*sum(donnees**beta*log(donnees)**2)/sum(donnees**beta)
terme3 = n*(sum(donnees**beta*log(donnees))/sum(donnees**beta))**2
return -terme1 - terme2 + terme3
return (g, dg)
Nous définissons ainsi les fonctions correspondant à notre problème :
g, dg = fait_g_dg(durees_de_vie)
Nous définissons ensuite une fonction mettant en œuvre la méthode de Newton :
def newton(devinette_initiale,
cible,
cible_derivee,
tolerance=1e-6,
iter_max=100):
pos = devinette_initiale
valeur = cible(pos)
idx = 0
while idx <= iter_max and abs(valeur) > tolerance :
pos -= valeur/cible_derivee(pos)
idx += 1
valeur = cible(pos)
return (pos,valeur,idx)
Nous pouvons à présent obtenir $\hat{β}$ numériquement en prenant une devinette initiale de 1 :
(beta_chapeau,cible_a_beta_chapeau,nombre_iter) = newton(1,g,dg)
(beta_chapeau,cible_a_beta_chapeau,nombre_iter)
$\hat{λ}$ est alors directement obtenu avec :
lambda_chapeau = len(durees_de_vie)/sum(array(durees_de_vie)**beta_chapeau)
lambda_chapeau
Exemple 2.1¶
Nous considérons ici une fonction de log-vraisemblance dépendant de deux paramètres : μ~1~ et μ~2~ définie par : $$l(μ_1,μ_2) = -\frac{1}{2}(15.6-μ_1)^2 -\frac{1}{2}(29.3-μ_2)^2 -\frac{1}{2}(45.8-μ_1-μ_2)^2 \, .$$ L'EMV est $(\hat{μ}_1,\hat{μ}_2) = (15.9,29.6)$ et : $$l(\hat{μ}_1,\hat{μ}_2) = -\frac{1}{2} \left(0.3^2+0.3^2+0.3^2\right) = -0.135\, .$$ La fonction de log-vraisemblance relative est donc : $$r((μ_1,μ_2) = l(μ_1,μ_2) + 0.135 \, .$$
Nous souhaitons faire une graphe à lignes de niveau représentant cette fonction avec les niveaux correspondant aux vraisemblances relatives de 0,5, 0,1 et 0,01. Nous commençons par définir une fonction de log-vraisemblance relative :
def r_gogo(mu1,mu2): return -0.5*(15.6-mu1)**2-0.5*(29.3-mu2)**2-0.5*(45.8-mu1-mu2)**2+0.135
Nous souhaitons ensuite évaluer notre fonction en tout point d'une grille allant de 12 à 20 pour μ~1~ et de 26 à 34 pour μ~2~ (nous nous basons sur la figure 10.2.2 de la page 62). Pour cela nous allons faire appel à la fonction meshgrid de numpy (module déjà chargé avec notre commande magique %pylab) :
mu1, mu2 = meshgrid(linspace(12,20,201),linspace(26,34,201))
Comme notre fonction r_gogo est vectorisée, nous pouvons simplement l'évaluer en tout point de notre grille avec :
ma_r_gogo = r_gogo(mu1,mu2)
J'utilise à présent la commande magique permettant de faire apparaître les graphes dans ce document :
%matplotlib inline
Notre graphe est alors construit par un appel à la fonction contour :
contour(mu1,mu2,ma_r_gogo,[log(0.5),log(0.1),log(0.01)])
grid()
xlabel(r'$\mu_1$')
ylabel(r'$\mu_2$')
On peut aussi effectuer cette tâche de façon plus « minimaliste » avec :
x1 = linspace(12,20,201)
x2 = linspace(26,34,201)
ma_r_gogo2 = [r_gogo(x,y) for x in x1 for y in x2]
z = array(ma_r_gogo2)
z.shape = (201,201)
contour(x1,x2,z,log(array([0.5,0.1,0.01])))
grid()
xlabel(r'$\mu_1$')
ylabel(r'$\mu_2$')
Clairement, deux boucles for imbriquées permettraient aussi d'arriver au résultat souhaité (évaluation en tout point de la grille).
Exemple 2.2¶
Nous revenons à l'exemple de fiabilité précédent où nous avons modélisé les durées de bon fonctionnement avec une loi de Weibull. Nous changeons ici la paramétrisation en remplaçant λ par θ avec : λ = θ^-β^, la log-vraisemblance s'écrit alors : $$l(θ,β) = -n β \log θ + n \log β + (β - 1) \sum \log x_i - θ^{-β} \sum x_i^β \, .$$ Nous avons alors $\hat{β} = 2.1021$ et $\hat{θ} = 81.88$ et le maximum de la fonction de vraisemblance est : $$l(\hat{θ},\hat{β}) = -113.691 \, ,$$ d'où la log-vraisemblance relative : $$r(θ,β) = l(θ,β) + 113.691 \, .$$
Pour construire le graphe en lignes de niveau nous procédons comme précédemment :
def fait_r_weibull(donnees):
donnees = array(donnees)
n = len(donnees)
somme = sum(log(donnees))
def r(theta,beta):
return -n*beta*log(theta)+n*log(beta)+(beta-1)*somme-theta**(-beta)*sum(donnees**beta)+113.691
return r
r_weibull = fait_r_weibull(durees_de_vie)
theta_v = linspace(50,120,101)
beta_v = linspace(1,3.5,101)
ma_r_weibull = [exp(r_weibull(theta,beta)) for beta in beta_v for theta in theta_v]
ma_r_weibull = array(ma_r_weibull)
ma_r_weibull.shape = (101,101)
contour(theta_v,beta_v,ma_r_weibull,[0.5,0.1,0.01])
grid()
xlabel(r'$\theta$')
ylabel(r'$\beta$')
On retrouve ainsi la figure 10.2.3 page 64 !
Exemple 5.1¶
Nous considérons ici un test d'insecticide où 5 concentrations différentes (de l'insecticide) sont utilisées sur 5 groupes d'insectes différents. Les concentrations en mg/l sont :
concentrations = [2.6, 3.8, 5.1, 7.7, 10.2]
Les effectifs (nombre d'insectes dans chaque groupe) sont :
effectifs = [50,48,46,49,50]
Le nombre d'insectes tués dans chacun des cinq groupes est :
tues = [6, 16, 24, 42, 44]
Nous allons postuler un modèle logistique à deux inconnues, α et β où la probabilité de décès pour un insecte choisi aléatoirement et soumis à une dose (le log de la concentration d'insecticide) d est donnée par : $$p = 1 - \frac{1}{1+\exp(α+β d)} \, .$$ Ce qui peut se réécrire comme : $$ \log \frac{p}{1-p} = α+β d\, ,$$ où le terme de gauche est nommé transformée logistique de p. Une première façon de juger de la pertinence du modèle logistique est de remplacer p par une estimation $\hat{p}$ -- le nombre de décès divisé par le nombre d'individus soumis à une dose donnée --, d'effectuer les transformations logistiques (des estimations aux différentes doses) et de les tracer en fonction de la dose :
p_chapeau = array(tues)/array(effectifs)
logit_p_chapeau = log(p_chapeau/(1-p_chapeau))
plot(log(array(concentrations)),logit_p_chapeau,'o')
xlabel("Dose")
ylabel("log(y/(n-y))")
ylim([-3,3])
def fait_l_et_plus_logit(n,y,d):
def l(par):
alpha = par[0]; beta = par[1]
p = 1-1/(1+exp(alpha+beta*d))
return sum(y*(alpha+beta*d)+n*log(1-p))
def S(par):
alpha = par[0]; beta = par[1]
mu = n*(1-1/(1+exp(alpha+beta*d)))
return array([sum(y-mu),sum((y-mu)*d)])
def I(par):
alpha = par[0]; beta = par[1]
p = 1-1/(1+exp(alpha+beta*d))
v = n*p*(1-p)
return array([[sum(v),sum(v*d)],[sum(v*d),sum(v*d**2)]])
return (l,S,I)
l_logit, S_logit, I_logit = fait_l_et_plus_logit(n=effectifs,y=tues,d=log(array(concentrations)))
def newton_raphson(devinette_initiale,
Score,
Information,
tolerance=1e-6,
iter_max=100):
pos = devinette_initiale
valeur = Score(pos)
idx = 0
while idx <= iter_max and max(abs(valeur)) > tolerance :
pos += linalg.solve(Information(pos),valeur)
idx += 1
valeur = Score(pos)
return (pos,valeur,idx)
resultat_logit = newton_raphson(array([-5.,3.]),S_logit,I_logit)
resultat_logit