Applications en python
IMAD – Institution genevoise de maintien à domicile, Genève Suisse
October 23, 2023
La statistique est la discipline qui étudie des phénomènes à travers la collecte de données, leur traitement, leur analyse, l’interprétation des résultats et leur présentation afin de rendre ces données compréhensibles par tous. (wikipedia)
\(\Longrightarrow\) Trois professions
La statistique est la discipline qui étudie des phénomènes à travers la collecte de données, leur traitement, leur analyse, l’interprétation des résultats et leur présentation afin de rendre ces données compréhensibles par tous. (wikipedia)
\(\Longrightarrow\) Trois professions
\(\Longrightarrow\) Deux approches
Des chercheurs pointent leur télescope vers le ciel et observent la lumière provenant d’une seule étoile.
Generate data
# Generating some simple photon count data
import numpy as np
from scipy import stats
np.random.seed(42) # for repeatability
F_true = 1000 # true flux, say number of photons measured in 1 second
N = 50 # number of measurements
F = stats.poisson(F_true).rvs(N) # N measurements of the flux
e = np.sqrt(F) # errors on Poisson counts estimated via square root
C’est une première forme d’analyse qui ne renseigne que sur l’échantillon
On considère que les données sont aléatoires mais pas les paramètres.
On cherche à estimer la valeur de flux \(F\) qui maximise la probabilité \(P(D~|~F)\).
\[ P(D_i~|~F_{\rm true}) = \frac{1}{\sqrt{2\pi e_i^2}} \exp{\left[\frac{-(F_i - F_{\rm true})^2}{2 e_i^2}\right]} \]
La fonction de vraisemblance
\[P(D~|~F_{\rm true}) = \prod_{i=1}^N P(D_i~|~F_{\rm true})\]
est ici une fonction de \(F_{\rm true}\) donc définie sur l’espace du paramètre.
On construit sa version log (numérique)
\[\log\mathcal{L} = -\frac{1}{2} \sum_{i=1}^N \left[ \log(2\pi e_i^2) + \frac{(F_i - F_{\rm true})^2}{e_i^2} \right]\]
Pour trouver la valeur de flux qui maximise cette fonction, on cherche la valeur \(F\) du flux telle que
\[d\log\mathcal{L}/dF_{\rm true} = 0\]
Le résultat est le suivant \[ F_{\rm est} = \frac{\sum w_i F_i}{\sum w_i};~~w_i = 1/e_i^2 \]
Remarquez que dans le cas particulier où toutes les erreurs \(e_i\) sont égaux, cela se réduit à \[ F_{\rm est} = \frac{1}{N}\sum_{i=1}^N F_i \]
Conformément à l’intuition, \(F_{\rm est}\) est simplement la moyenne des données observées lorsque les erreurs sont égales.
Quelle est l’erreur de cette estimation ?
En ajustant une approximation gaussienne à la courbe de vraisemblance à son maximum, on obtient comme estimation de son écart-type :
\[ \sigma_{\rm est} = \left(\sum_{i=1}^N w_i \right)^{-1/2} \]
Résultat numérique
On considère que les données ET les paramètres sont probabilistes.
On peut donc écrire la distribution conjointe (impossible en fréquentiste) \[ P(F_{\rm true},~D)\]
On utilise la règle de Bayes pour exprimer la distribution recherchée \[ P(F_{\rm true}~|~D) = \frac{P(D~|~F_{\rm true})~P(F_{\rm true})}{P(D)} \]
Dans la grande majorité des cas on ne sait pas calculer cette probabilité !
Deux solutions :
#| echo: true
# import packages
from jax import random
import numpyro
import numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS, Predictive
# define model
def model(e_obs, F_obs=None):
F = numpyro.sample("F", dist.Uniform(0, 10000)) # prior uniform
numpyro.sample("obs", dist.Normal(F, e_obs), obs=F_obs) # échantillon distribué ~ loi normal
# run sampling
mcmc = MCMC(NUTS(model), num_warmup=500, num_samples=2000, progress_bar=False)
mcmc.run(random.PRNGKey(42), e_obs=e, F_obs=F)
# display result
mcmc.print_summary(0.89)
# import packages
from jax import random
import numpyro
import numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS, Predictive
# define model
def model(e_obs, F_obs=None):
F = numpyro.sample("F", dist.Uniform(0, 10000)) # prior uniform
numpyro.sample("obs", dist.Normal(F, e_obs), obs=F_obs) # échantillon distribué ~ loi normal
# run sampling
mcmc = MCMC(NUTS(model), num_warmup=500, num_samples=2000, progress_bar=False)
mcmc.run(random.PRNGKey(42), e_obs=e, F_obs=F)
# display result
mcmc.print_summary(0.89)
mean std median 5.5% 94.5% n_eff r_hat
F 998.98 4.56 998.85 992.00 1006.73 526.64 1.00
Number of divergences: 0
Question : Une fois que l’animateur a ouvert la porte C, le joueur a-t-il intérêt à changer son choix ?
\[P(H|D \land A) ?\]
\[P(H~|~D \land A) = \frac{P(D~|~H \land A)~P(H~|~A)}{P(D~|~A)}\]
H | \(P(H~|~A)\) | \(P(D~|~H~\land~A)\) | \(P(H~|~D~\land~A)\) |
---|---|---|---|
porte A | 1/3 | ||
porte B | 1/3 | ||
porte C | 1/3 |
\[P(H~|~D \land A) = \frac{P(D~|~H \land A)~P(H~|~A)}{P(D~|~A)}\]
H | \(P(H~|~A)\) | \(P(D~|~H~\land~A)\) | \(P(H~|~D~\land~A)\) |
---|---|---|---|
porte A | 1/3 |
|
|
porte B | 1/3 |
|
|
porte C | 1/3 |
|
\[P(H~|~D \land A) = \frac{P(D~|~H \land A)~P(H~|~A)}{P(D~|~A)}\]
H | \(P(H|A)\) | \(P(D|H \land A)\) | \(P(H|D \land A)\) |
---|---|---|---|
porte A | 1/3 |
|
|
porte B | 1/3 |
|
|
porte C | 1/3 |
|
|
Visualiser les résultats
Boucler entre 3 et 5 jusqu’à obtenir une réponse satisfaisante
Question : Peut-on prédire la taille à partir du poids après 18 ans ?
Préparation des données :
Pour Numpyro
Définition du modèle
\[ \begin{align} W_i \sim Normal(\mu_i, \sigma) \\ mu_i = \alpha + \beta H_i \\ \alpha \sim Normal(0, 10) \\ \beta \sim HalfNormal(1) \\ \sigma \sim HalfNormal(1) \\ \end{align} \]
Implémentation du modèle
# Define the model
def mdl1(height, weight=None):
# Define prior
a = numpyro.sample('a', dist.Normal(0, 1))
b = numpyro.sample('b', dist.Normal(0, 1))
s = numpyro.sample('s', dist.HalfNormal(1))
# Define likelihood
numpyro.sample("obs", dist.Normal(a+b*height, s), obs=weight)
# plot model
numpyro.render_model(mdl1, model_args=dat1.values(), render_distributions=True)
mean std median 5.5% 94.5% n_eff r_hat
a 41.91 0.42 41.93 41.25 42.59 1014.03 1.00
b 0.63 0.03 0.63 0.57 0.68 1297.06 1.00
s 5.09 0.29 5.07 4.63 5.54 935.34 1.00
Number of divergences: 0
Question : Quelle est la différence de poids entre les hommes et les femmes après 18 ans ?
Définition du modèle
\[ \begin{align} W_i \sim Normal(\mu_i, \sigma) \\ \mu_i = \alpha[male] \\ \alpha \sim HalfNormal(0, 10) \\ \sigma \sim HalfNormal(1) \\ \end{align} \]
Implémentation du modèle
# create dict for inference with numpyro
dat2 = {
"male" : data_sup18.male.values,
"weight" : data_sup18.weight.values
}
# Define the model
def mdl2(male, weight=None):
# Define prior
with numpyro.plate("male", 2):
a = numpyro.sample('a', dist.HalfNormal(10))
s = numpyro.sample('s', dist.HalfNormal(1))
diff_a = numpyro.deterministic("diff_a", a[1]- a[0])
# Define likelihood
numpyro.sample("obs", dist.Normal(a[male], s), obs=weight)
# plot model
numpyro.render_model(mdl2, model_args=dat2.values(), render_distributions=True)
Graph du modèle
mean | sd | hdi_5.5% | hdi_94.5% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
---|---|---|---|---|---|---|---|---|---|
a[0] | 41.755 | 0.387 | 41.159 | 42.353 | 0.013 | 0.009 | 856.0 | 552.0 | NaN |
a[1] | 48.511 | 0.439 | 47.904 | 49.291 | 0.013 | 0.009 | 1080.0 | 555.0 | NaN |
diff_a | 6.756 | 0.607 | 5.707 | 7.607 | 0.020 | 0.015 | 874.0 | 658.0 | NaN |
s | 5.315 | 0.196 | 5.028 | 5.645 | 0.007 | 0.005 | 772.0 | 638.0 | NaN |
Question : Quelle est le meilleur modèle pour prédire le poids en fonction de la taille ?
# This time with all the data
df_howell1["height_c"] = center(df_howell1.height)
dat3 = {
"height" : df_howell1.height_c.values,
"weight" : df_howell1.weight.values
}
# Define the model
def mdl3a(height, weight=None):
# Define prior
a = numpyro.sample('a', dist.Normal(0, 10))
b = numpyro.sample('b', dist.Normal(0, 1))
s = numpyro.sample('s', dist.HalfNormal(1))
# Define likelihood
numpyro.sample("obs", dist.Normal(a+b*height, s), obs=weight)
# Define the model
def mdl3b(height, weight=None):
# Define prior
a = numpyro.sample('a', dist.Normal(0, 10))
b1 = numpyro.sample('b1', dist.Normal(0, 1))
b2 = numpyro.sample('b2', dist.Normal(0, 1))
s = numpyro.sample('s', dist.HalfNormal(1))
# Define likelihood
numpyro.sample("obs",
dist.Normal(a+b1*height+b2*height**2, s),
obs=weight)
# Define the model
def mdl3c(height, weight=None):
# Define prior
a = numpyro.sample('a', dist.Normal(0, 10))
b1 = numpyro.sample('b1', dist.Normal(0, 1))
b2 = numpyro.sample('b2', dist.Normal(0, 1))
b3 = numpyro.sample('b3', dist.Normal(0, 1))
s = numpyro.sample('s', dist.HalfNormal(1))
# Define likelihood
numpyro.sample("obs",
dist.Normal(a+b1*height+b2*height**2+b3*height**3, s),
obs=weight)
mean std median 5.5% 94.5% n_eff r_hat
a 35.59 0.20 35.59 35.28 35.91 577.64 1.00
b 0.50 0.01 0.50 0.49 0.51 1082.96 1.00
s 4.91 0.15 4.91 4.68 5.14 1011.70 1.00
Number of divergences: 0
mean std median 5.5% 94.5% n_eff r_hat
a 32.51 0.25 32.52 32.08 32.90 542.73 1.00
b1 0.64 0.01 0.64 0.63 0.66 314.48 1.00
b2 0.00 0.00 0.00 0.00 0.00 417.83 1.00
s 3.91 0.12 3.91 3.73 4.10 389.55 1.00
Number of divergences: 0
mean std median 5.5% 94.5% n_eff r_hat
a 1.76 0.02 1.75 1.73 1.79 3.75 1.53
b1 1.02 0.02 1.02 0.98 1.05 9.46 1.03
b2 0.04 0.00 0.04 0.04 0.04 741.49 1.00
b3 0.00 0.00 0.00 0.00 0.00 75.08 1.00
s 16.30 0.38 16.30 15.74 16.83 117.14 1.00
Number of divergences: 0
rank | elpd_waic | p_waic | elpd_diff | weight | se | dse | warning | scale | |
---|---|---|---|---|---|---|---|---|---|
mdl quardatic | 0 | -1520.630761 | 3.796908 | 0.000000 | 0.945619 | 20.899299 | 0.000000 | False | log |
mdl linear | 1 | -1648.960918 | 3.240282 | 128.330157 | 0.054381 | 14.881861 | 16.194413 | False | log |
mdl cubic | 2 | -2422.694490 | 3.614403 | 902.063729 | 0.000000 | 19.844980 | 22.873266 | True | log |
L’intérêt de ces deux packages est qu’ils s’ouvrent sur les réseaux de neurones…