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é.

Préliminaires

Charegement des données, etc…

Nous chargeons les données et les modules dont nous allons avoir besoin :

In [1]:
%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()    
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib
In [2]:
%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 :

In [3]:
G_chapeau = 0.147397980373
S2_chapeau = 268.741290284

La version stabilisée de pile_POMC est simplement obtenue avec :

In [4]:
pile_POMC_stab = 2*sqrt(pile_POMC/G_chapeau+S2_chapeau)

La matrice des $RSS_{i,j}$ est obtenues avec :

In [5]:
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 :

In [6]:
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')
/home/xtof/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:2: RuntimeWarning: divide by zero encountered in log
  from IPython.kernel.zmq import kernelapp as app

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 :

In [7]:
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 :

In [8]:
bons_pxl = [x[:(len(x)-1)] for x in bons_pxl]; bons_pxl
Out[8]:
[array([25, 25, 25, 25, 25, 26, 26, 26, 26, 26, 26, 27, 27, 27, 27, 27, 27,
        28, 28, 28, 28, 28, 28, 29, 29, 29, 29, 29, 29, 30, 30, 30, 30, 31,
        31]),
 array([35, 36, 37, 38, 39, 35, 36, 37, 38, 39, 40, 35, 36, 37, 38, 39, 40,
        35, 36, 37, 38, 39, 40, 35, 36, 37, 38, 39, 40, 37, 38, 39, 40, 38,
        39])]

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) :

In [9]:
bons_signaux = pile_POMC[bons_pxl[0],bons_pxl[1],:]
In [10]:
ligne_de_base = bons_signaux[:,:5].mean(1)
In [11]:
f_0 = mean(bons_signaux/ligne_de_base.reshape(35,1),0)
In [12]:
b_0 = 100
In [13]:
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()]
In [14]:
def var_stab(x,G=G_chapeau,S2=S2_chapeau): return 2*sqrt(x/G+S2)
In [15]:
par_0 = zeros((1+35+167,))
par_0[0] = b_0
par_0[1:36] = phi_0
par_0[36:] = f_0[1:]
In [16]:
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
In [17]:
bons_signaux_stab = var_stab(bons_signaux)    
rss0 = fait_rss(bons_signaux_stab)
rss0(log(par_0))
Out[17]:
$$1113689.08989$$
In [18]:
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})
Warning: Maximum number of function evaluations has been exceeded.
Warning: Maximum number of function evaluations has been exceeded.
Warning: Desired error not necessarily achieved due to precision loss.
         Current function value: 6015.136848
         Iterations: 300
         Function evaluations: 155596
         Gradient evaluations: 759
Optimization terminated successfully.
         Current function value: 6015.136848
         Iterations: 1660
         Function evaluations: 3212
In [19]:
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:])
In [20]:
plot(f_1)
Out[20]:
[<matplotlib.lines.Line2D at 0x7fe2ac8210d0>]
In [ ]: