Introduction¶
Aujourd'hui, nous finissons l'étude de cas entamée il y a un mois avec la solution de l'excercice suivant basé sur les « données POMC »…
Rappel : énoncé de l'exercice¶
- nous allons être « prudents » et garder comme pixels de notre région d'intérêt les 35 pixels formant une région connexe;
- nous allons modéliser l'intensité de fluorescence de chaque pixel par : $$S_{i,j}(t) = \phi_{i,j} \, f(t) + b \; ,$$ où $f(t)$ est la dynamique du signal _pour chaque pixel_, $\phi_{i,j}$ est un paramètre spécifique de chaque pixel et $b$ décrit l'auto-fluorescence supposée identique pour chaque pixel ;
- le temps $t$ prend en fait des valeurs discrètes, $t = \delta \, k$ ($\delta$ = 150 ms) et nous cherchons un ensemble d'estimateurs ponctuels : $\{f_1,f_2,\ldots,f_K\}$ ($K$ = 168) où $f_k = f(\delta \, k)$ ;
- nous avons alors 35 ($\phi_{i,j}$) + 168 ($f_k$) + 1 ($b$) = 204 paramètres pour 35 x 168 = 5880 observations.
- nous devons encore ajouter une contrainte puisqu'avec notre spécification du modèle : $$S_{i,j,k} = \phi_{i,j} \, f_k + b \; ,$$ nous pouvons multiplier tous les $\phi{i,j}$ par 2 et diviser les $f_k$ par 2 sans changer la prédiction (_le modèle n'est pas identifiable) ;
- nous allons poser $f_1=1$ et nos estimateurs correspondent alors à ce que les neurophysiologistes font lorsqu'ils analysent ce type de données -- ils travaillent avec $\Delta S(t) / S(1)$ -- car : $$\Delta S(t) / S(1) = \frac{S(t) - S(1)}{S(1)} = f(t) - 1 + \mathrm{bruit} \, ;$$
- remarquez qu'aucune mesure indépendante de l'auto-fluorescence n'est nécessaire.
Vous définirez une fonction de coût correspondant à ce modèle comme la somme des résidus aux carrés ; les résidus étant calculés à partir des données et de la prédiction stabilisés. Vous utiliserez ensuite minimize et la méthode BFGS. Il faudra réfléchir un peu pour trouver les valeurs initiales. En utilisant les paramètres estimés, vous simulerez entre 500 et 1000 jeux de données en utilisant le « modèle exacte », c'est-à-dire :
- vous effectuerez d'abord des tirages suivant des lois de Poisson : $$N_{i,j,k} \sim \mathcal{P}(\hat{\phi_{i,j}} \, \hat{f}_k + \hat{b}) \,;$$
- vous effectuerez ensuite des tirages suivant une loi normale : $L_{i,j,k} \sim \mathcal{N}(0,\sigma_L^2)$ ce que nous avons aussi noté : $L_{i,j,k} \sim \sigma_L \epsilon_{i,j,k}$ ;
- vous formerez les sommes : $M_{i,j,k} = N_{i,j,k} + L_{i,j,k}$ que vous multiplierez par $G$ pour obtenir vos observations simulez : $$ADU_{i,j,k} = G M_{i,j,k} \, .$$ Vous ajusterez le modèle sur chaque jeu de données simulés ce qui vous donnera une version simulée (on parle de bootstrap paramétrique) de la loi d'échantillonage de votre estimateur. Vous comparerez la valeur de la somme des résidus au carré que vous avez obtenus sur le vrai jeu de données à la distribution de cette statistique résultant des simulations. Que concluez vous ? Nous comparerons en cours les erreurs types de chaque paramètre à l'élément correspondant de l'espérance de la matrice d'information.
Quelques indices :
- la spécification ci-dessus revient à dire que les observations $adu_{i,j,k}$ sont modélisées comme des réalisations de variables aléatoires normales et indépendantes de lois : $$ADU_{i,j,k} \sim G \left(\phi_{i,j} \, f_k + b\right) + G \sqrt{\phi_{i,j} \, f_k + b + \sigma_L^2} \epsilon_{i,j,k} \, ;$$
- si nous travaillons avec les versions stabilisées des observations : $z_{i,j,k} = 2 \sqrt{adu_{i,j,k}/G + \sigma_L^2}$, nous pouvons modélisées les $z_{i,j,k}$ comme des réalisations de variables aléatoires indépendantes normales de variance 1 dont les lois sont : $$Z_{i,j,k} \sim 2 \sqrt{\phi_{i,j} \, f_k + b + \sigma_L^2} + \epsilon_{i,j,k} \, ;$$
- nous aboutissons alors à une fonction de cout (l'opposée de la log-vraisemblance) qui ressemble à celle de notre modèle complet utilisé pour les données de calibration du capteur CCD, c'est-à-dire sur une somme des résidus au carré.
%pylab
from __future__ import print_function, division, unicode_literals, absolute_import
import sympy as sy
sy.init_printing()
import urllib
urllib.urlretrieve("http://xtof.disque.math.cnrs.fr/data/Data_POMC.npz","Data_POMC.npz")
POMC = load("Data_POMC.npz")
temps_POMC = POMC['arr_0']
pile_POMC = POMC['arr_1']
POMC.close()
%matplotlib inline
Détection automatique de la région d'intérêt¶
Le gain et la variance de lecture du capteur CCD que nous avons déterminés sont :
G_chapeau = 0.147397980373
S2_chapeau = 268.741290284
La version stabilisée de pile_POMC est simplement obtenue avec :
pile_POMC_stab = 2*sqrt(pile_POMC/G_chapeau+S2_chapeau)
La matrice des $RSS_{i,j}$ est obtenues avec :
pile_POMC_RSS = sum((pile_POMC_stab - pile_POMC_stab.mean(2).reshape((60,80,1)))**2,2)
Et le graphe des pixels tels que $\log\left(1 - F_{\chi_{K-1}^2}(RSS_{i,j})\right)$ est proche de $-\infty$ est construit avec :
from scipy.stats import chi2
pile_POMC_logCFRSS = log(1-chi2.cdf(pile_POMC_RSS,167))
imshow(transpose(pile_POMC_logCFRSS == -Inf),origin='lower')
set_cmap('gray')
Les pixels pour lesquels $\log\left(1 - F_{\chi_{K-1}^2}(RSS_{i,j})\right)$ est numériquement égale à $-\infty$ sont obtenus par un appel à la fonction where :
bons_pxl = where(pile_POMC_logCFRSS == -Inf)
Là, une impression du résultat montre que le pixel isolé est le dernier trouvé et les coordonnées des pixels connexes sont donc obtenues avec :
bons_pxl = [x[:(len(x)-1)] for x in bons_pxl]; bons_pxl
Valeurs initiales des paramètres¶
Nous allons obtenir une valeur intiale du vecteur $f_k$ en moyennant une version normalisée du décours temporel obtenu en chacun des pixels de notre région d'intérêt. Nous estimons la ligne de base de chaque pixel en moyennant les observations obtenues au cours des cinq premiers temps. Ici il faut faire un peu attention car notre région connexe ne constitue pas un rectangle, nous ne pouvons donc pas extraire les pixels intéressants et conserver une structure de pile. En extrayant les signaux des pixels d'intérêts, nous allons donc réduire la dimension des données (consultez la documentation de Advanced indexing pour plus de détails) :
bons_signaux = pile_POMC[bons_pxl[0],bons_pxl[1],:]
ligne_de_base = bons_signaux[:,:5].mean(1)
f_0 = mean(bons_signaux/ligne_de_base.reshape(35,1),0)
b_0 = 100
phi_0 = [(len(f_0)*sum(f_0*y)-sum(f_0)*sum(y))/(len(f_0)*sum(f_0**2)-sum(f_0)**2) for y in bons_signaux.tolist()]
def var_stab(x,G=G_chapeau,S2=S2_chapeau): return 2*sqrt(x/G+S2)
par_0 = zeros((1+35+167,))
par_0[0] = b_0
par_0[1:36] = phi_0
par_0[36:] = f_0[1:]
def fait_rss(donnees):
npix,K = donnees.shape
f = ones((K,))
def rss(par):
par = exp(par)
b = par[0]
phi = par[1:(npix+1)]
f[1:] = par[(npix+1):]
pred = var_stab(outer(phi,f)+b)
return sum((donnees-pred)**2)
return rss
bons_signaux_stab = var_stab(bons_signaux)
rss0 = fait_rss(bons_signaux_stab)
rss0(log(par_0))
from scipy.optimize import minimize
ajuste0 = minimize(rss0,log(par_0),method='nelder-mead',options={'disp': True})
ajuste0 = minimize(rss0,ajuste0.x,method='nelder-mead',options={'disp': True})
ajuste0 = minimize(rss0,ajuste0.x,method='BFGS',options={'disp': True})
ajuste0 = minimize(rss0,ajuste0.x,method='nelder-mead',options={'disp': True})
b_1 = exp(ajuste0.x[0])
phi_1 = exp(ajuste0.x[1:36])
f_1 = ones((bons_signaux.shape[1],))
f_1[1:] = exp(ajuste0.x[36:])
plot(f_1)
