import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
Dans un plan d'experience, on a souvent les resultat qui tombent les uns apres les autres. Supposons qu'on ait prévu de faire une expérience sur 100 personnes pour demontrer que l'hypothese nulle n'est pas vérifiée. En pratique, on effectue les expériences par groupes de 5 personnes, et on commence a calculer une p-valeur a partir de 20 personnes minimum, sans quoi on a pas convergence dans le TCL ou le Chi2.
Comment évolue la p-valeur au fil des expérience de 20 a 100 personnes?
# afficher l'evolution de la p valeur de 1 a 100
pi = np.ones(10)/10 # loi uniforme a 10 valeurs
samples_1 = np.random.multinomial(20, pi) # les 20 premiers echantilons
p_val_1 = scipy.stats.chisquare(samples_1)[1]
samples_2 = samples_1 + np.random.multinomial(10, pi) # on ajoute les 10 suivants
p_val_2 = scipy.stats.chisquare(samples_2)[1]
samples_3 = samples_2 + np.random.multinomial(10, pi) # on ajoute les 10 suivants
p_val_3 = scipy.stats.chisquare(samples_3)[1] # on ajoute les 10 suivants
# etc...
p_val_1, p_val_2, p_val_3
(0.7399182920946538, 0.4685950118243071, 0.5341462169096916)
On peut etre tenté de s'arreter dès que la p-valeur est inférieure au niveau alpha fixé, apres avoir fait un minimum de test (30 par exemple) mais avant d'avoir fait les experiences sur les 100 personnes.
Quelle est l'influence sur le niveau du test lorsqu'on choisit quand on s'arrete ? A-t-on toujours correspondance entre alpha et la frequence de rejet de l'hypothese nulle alors qu'elle est vérifiée ?
# mode 'normal'
# on répète le test N fois, et on regarde le pourcentage de rejet de H_0, qui doit etre proche de alpha.
N = 1000 # nombre de tests d'independance
n = 100 # nombre d'echantillon par test
alpha = 0.05
pi = np.ones(10)/10 # loi uniforme a 10 valeurs
n_reject = 0
for i in range(N):
samples = np.random.multinomial(n, pi)
p_val = scipy.stats.chisquare(samples)[1]
n_reject += p_val < alpha
print(f'N = {N} experiments with {n} samples, alpha={100*alpha}%, %reject {100*n_reject/N:.2f}')
N = 1000 experiments with 100 samples, alpha=5.0%, %reject 5.20
# mode 'arret quand ca nous arrange'
# on répète le test N fois
N = 10000
# on prend au minimum 30 patients et on ajoute par group de 5 jusqu'a 100.
n_min = 30
n_per_group = 5
n = 100
n_reject = 0
for i in range(N):
# A VOUS DE COMPLETER
# on commence par faire tirer une multinomiale avec n=n_min, et on calcule la p-valeur
samples = np.random.multinomial(n_min, pi)
p_val = scipy.stats.chisquare(samples)[1]
n_samples = n_min
while p_val >= alpha and n_samples < n:
samples += np.random.multinomial(n_per_group, pi)
n_samples += n_per_group
p_val = scipy.stats.chisquare(samples)[1]
if p_val < alpha:
n_reject += 1
print(f'N = {N} experiments with {n} samples, alpha={100*alpha}%, %reject {100*n_reject/N:.2f}')
N = 10000 experiments with 100 samples, alpha=5.0%, %reject 16.81
# probabilites observees par le medecin
p_cancer = 0.0005
p_covid = 0.5
p_gastro = 0.2
p_toux_s_cancer = 0.95
p_toux_s_gastro = 0.003
p_toux_s_covid = 0.89
k = 1
# ratio entre les p(H| donnees) ?
p_cancer_s_toux = k * p_cancer * p_toux_s_cancer
p_gastro_s_toux = k * p_gastro * p_toux_s_gastro
p_covid_s_toux = k * p_covid * p_toux_s_covid
p_covid_s_toux / p_cancer_s_toux, p_covid_s_toux / p_gastro_s_toux
(936.8421052631579, 741.6666666666666)
# generer des v.a. uniforme sur [0,1]
X = np.random.random(1000)
# valeurs theoriques
eps = 0.1
sigma_X = np.sqrt(1/12)
sigma_Y = np.sqrt((1+eps**2)/12)
a_th = 1
b_th = 0.5 * eps
corr_XY = sigma_X / sigma_Y
MSE_th = eps**2 / 12
print(f'theoretical values a={a_th}, b={b_th}, MSE={MSE_th:.6f}')
theoretical values a=1, b=0.05, MSE=0.000833
# valeurs estimmées par méthode de Monte-Carlo
N = 1000
X = np.random.random(N)
Y = X + eps* np.random.random(N)
E_X, sigma_X = np.mean(X), np.std(X)
E_Y, sigma_Y = np.mean(Y), np.std(Y)
cov_XY = np.mean(X*Y) - E_X*E_Y
a = cov_XY / sigma_X**2
b = E_Y - a*E_X
MSE = np.mean((Y - (a*X +b))**2)
print(f'estimated values a={a}, b={b}, MSE={MSE:.6f}')
estimated values a=1.0016593020905673, b=0.04790184878505627, MSE=0.000824