%pylab
from __future__ import print_function, division, unicode_literals, absolute_import
import sympy as sy
sy.init_printing()
CCD_calibration_npz = load("CCD_calibration.npz")
CCD_calibration_piles = CCD_calibration_npz['arr_0']
CCD_calibration_npz.close()
%matplotlib inline
Calibration d'un capteur CCD¶
D'après ce que nous avons vu dans le cours précédent, nos observations $Y$ suivent la loi : $$Y \sim G \, \lambda + \sqrt{G^2 \, (\lambda+\sigma_{L}^2)} \, \epsilon \; ,$$ où $G$ est le gain du capteur, $\sigma_{L}^2$ est la « variance de lecture » et $\epsilon$ suit une loi normale centrée réduite. La valeur de $G$ et $\sigma_{L}^2$ est une spécification du capteur, mais l'expérience prouve que les fabriquants ont tendance à être trop optimistes sur les performances de leurs produits. C'est donc une bonne idée de mesurer ces paramètres par des expériences de calibration. Cela va aussi nous permettre de vérifier que notre modèle très simple ci-dessus s'applique.
- notre problème devient donc : comment tester $Y \sim G \, \lambda + \sqrt{G^2 \, (\lambda+\sigma_{L}^2)} \, \epsilon$ ? ou comment faire varier $\lambda$ ;
- considérons un pixel de notre capteur « regardant » un volume fixe d'une lamelle fluorescente, nous avons deux façons de manipuler le paramètre λ :
- faire changer l'intensité $i_{e}$ de la lumière excitant notre colorant ;
- faire varier la durée $\tau$ d'exposition.
Nous pouvons en effet écrire le λ de notre loi comme suit : $\lambda = \phi v c i_{e} \tau \, ,$ où $v$ est le volume de solution « vu » par le pixel, $c$ est la concentration de colorant, $\phi$ est le rendement quantique -- un paramètre qui décrit l'efficacité d'excitation du colorant, les pertes de photons de la lumière excitatrice et de la lumière émise, etc.
En pratique, il est plus facile de faire varier le temps d'exposition τ et c'est ce que nous avons fait dans les expériences qui suivent. J'ai ainsi demandé à mes collaborateurs du laboratoire de Peter Kloppenburg à l'Université de Cologne de : choisir 10 temps d'expositions ; pour chacun des 10 temps effectuer 100 expositions successives ; pour chacune des 10 x 100 expositions, enregistrer la valeur $y_{ij}$ prise par la v.a. $Y_{ij}$ au pixel $i,j$ de la CCD.
La raison pour laquelle on introduit une v.a. $Y_{i,j}$ par pixel et qu'il est très difficile d'avoir une illumination ($i_e$), un volume $v$ et un rendement quantique ($\phi$) uniforment. On a donc pour chaque pixel : $$Y_{i,j} \sim G \, p_{i,j} \tau + \sqrt{G^2 \, (p_{i,j} \tau+\sigma_{L}^2)} \, \epsilon_{i,j}\; ,$$ où $p_{i,j} = c \phi_{i,j} v_{i,j} i_{e,i,j}$.
- si notre modèle est correct alors pour chaque pixel $i,j$ nous devons avoir, pour un temps d'exposition donné, une moyenne empirique : $$\bar{y}_{i,j} = \frac{1}{100} \sum_{k=1}^{100} \quad y_{i,j,k} \approx G \, p_{i,j} \tau \; ;$$
- une variance empirique : $$S_{i,j}^2 = \frac{1}{99} \sum_{k=1}^{100} \quad (y_{i,j,k}-\bar{y}_{i,j})^2 \approx G^2 \, (p_{i,j} \tau+\sigma_{L}^2) \; ;$$
- et le graphe de $S_{i,j}^2$ en fonction de $\bar{y}_{i,j}$ devrait être une droite de pente $G$ et d'ordonnée en zéro $G^2 \sigma_{L}^2$.
Tests sur les données¶
Le graphe de la variance en fonction de la moyenne pour chaque pixel et pour chaque temps d'exposition est obtenu avec :
ADU_m = [pile.mean(2) for pile in CCD_calibration_piles]
ADU_v = [pile.var(2) for pile in CCD_calibration_piles]
for i in range(10):
scatter(ADU_m[i].flatten(),ADU_v[i].flatten(),0.05)
xlabel(r'$\overline{ADU}$')
ylabel("Var(ADU)")
Première estimation des paramètres $G$ et $\sigma_L^2$ basée sur un modèle linéaire¶
L'hétéroscédasdicité (variance non homogène) visible sur le dernier graphe n'est pas une surprise puisque le variance d'une variance d'un échantillon IID, de taille $n$, suivant une loi normale d'espérance μ et de variance $\sigma^2$ est : $$\mathrm{Var}[S^2] = \frac{2\sigma^4}{(n-1)} \; .$$
Dans un premier temps, nous allons ignorer cette hétéroscédasdicité et faire comme si le graphe ci-dessus était une collection d'observations de couples $(X_k,Y_k)$ où $X_k$ est connue et où $Y_k$ est sensé être relié à $X_k$ par une relation du type : $$ Y_k = a + b X_k + \sigma \epsilon_k \, ,$$ où $a$, $b$ et $\sigma$ sont des paramètres à déterminer et où les $\epsilon_k \stackrel{\mathrm{IID}}{\sim} \mathcal{N}(0,1)$. Nous faisons ici deux approximations par-rapport à ce que nous croyons être la réalité du graphe précédent :
- les $X_k$ jouent le rôle des $\bar{y}_{i,j}$, c'est-à-dire qu'ils ne sont pas connus mais résultent de mesures (ce sont des observations de variables aléatoires) ;
- comme nous en avons discuté en début de paragraphe, le « vrai » $\sigma$ dépend de $X_k$. Nous nous permettons ici ces approximations parce-qu'elles ne vont nous servir qu'à obtenir des devinettes initiales pour notre modèle complet et aussi parce-qu'elles nous permettent d'introduire une grand classique de l'analyse statistique : la régression linéaire. Un précision afin d'éviter toute confusion, on parle de régression linaire lorsque la variable dépendante ou la prédiction, ici $Y_k$, est une fonction linéaire des paramètres, ici $(a,b,\sigma^2)$. Nous pourrions ainsi travailler avec le modèle : $Y_k = a_0 + a_1 X_k + a_2 X_k^2 + b \cos X_k + c \exp X_k + \sigma \epsilon_k$ et nous parlerions encore de régression linéaire.
Si la modèle linéaire ci-dessus est correcte (il ne l'est pas mais pour l'instant nous faisons comme si) et que nous disposons d'une collection d'observations IID $Y_k$ nous pouvons écrire leur loi conjointe de la façon suivante : $$\mathrm{Pr}\{(Y_k)\} = \prod_k \frac{1}{\sqrt{2 \pi \sigma^2}} \, \exp\left(- \frac{1}{2} \frac{(Y_k - a - b X_k)^2}{\sigma^2}\right) \, .$$
D'où l'expression de la vraisemblance : $$L(a,b,\sigma^2) = \prod_k \frac{1}{\sqrt{\sigma^2}} \, \exp\left(- \frac{1}{2} \frac{(Y_k - a - b X_k)^2}{\sigma^2}\right) \, ,$$ et celle de la log-vraisemblance : $$l(a,b,\sigma^2) = \sum_{k=1}^K -\frac{1}{2} \log \sigma^2 - \frac{1}{2} \frac{(Y_k - a - b X_k)^2}{\sigma^2} = -\frac{K}{2} \log \sigma^2 - \frac{1}{2 \sigma^2} \sum_{k=1}^K (Y_k - a - b X_k)^2\, .$$ La première composante de la fonction score est : $$\frac{\partial l}{\partial a} = \frac{1}{\sigma^2} \left(Y_{\bullet} - K a - b X_{\bullet}\right) \, ,$$ où $Y{\bullet} = \sum_k Y_k$ et $X{\bullet} = \sum_k X_k$.
La deuxième composante est : $$\frac{\partial l}{\partial b} = \frac{1}{\sigma^2} \left( \sum_k X_k Y_k - a X_{\bullet} - b \sum_k X_k^2\right) \, .$$
La troisième composante est : $$\frac{\partial l}{\partial \sigma^2} = -\frac{K}{2 \sigma^2} + \frac{1}{2 \sigma^4} \sum_{k=1}^K (Y_k - a - b X_k)^2 \, .$$
On voit qu'il est possible de trouver l'EMV, $(\hat{a},\hat{b},\widehat{\sigma^2})$ analytiquement. L'examen des deux premières composantes de la fonction score montre que $\hat{a}$ et $\hat{b}$ sont indépendants de $\widehat{\sigma^2}$. La racine de la première composante de la fonction score nous permet d'exprimer $\hat{a}$ en fonction de $b$ : $$\hat{a} = \frac{Y_{\bullet} - b X_{\bullet}}{K} \, .$$
Par substitution dans l'expression de la racine de la deuxième composante nous obtenons : $$ \hat{b} = \frac{K \sum_k X_k Y_k - X_{\bullet} Y_{\bullet}}{K \sum_k X_k^2 - (X_{\bullet})^2} \, .$$
La troisième composante nous donne immédiatement : $$ \widehat{\sigma^2} = \frac{1}{K} \sum_{k=1}^K (Y_k - \hat{a} - \hat{b} X_k)^2 \, .$$
Une propriété remarquable de cette solution est que $\hat{a}$ et $\hat{b}$ sont des fonctions linéaires des observations $Y_k$. Si ces dernières sont normales on peut conclure, quelque soit la taille de l'échantillon ($K$) que $\hat{a}$ et $\hat{b}$ le sont aussi -- on n'a pas besoin de théorème centrale limite.
Ici pour $\hat{b}$ nous obtenons :
ADU_M = concatenate([m.flatten() for m in ADU_m])
ADU_V = concatenate([v.flatten() for v in ADU_v])
b_chapeau = (len(ADU_M)*sum(ADU_M*ADU_V)-sum(ADU_M)*sum(ADU_V))/(len(ADU_M)*sum(ADU_M**2)-sum(ADU_M)**2)
b_chapeau
Et pour $\hat{a}$ nous obtenons :
a_chapeau = (sum(ADU_v) - b_chapeau*sum(ADU_M))/len(ADU_M)
a_chapeau
Nos premières estimations pour les deux paramètres du capteur CCD sont ainsi :
G_0 = b_chapeau
sigma2_0 = a_chapeau / G_0**2
(G_0,sigma2_0)
Une second estimation prenant en compte l'hétéroscédasdicité¶
Nous pouvons assez simplement raffiner notre modèle en effectuant une nouvelle régression linéaire avec variance inhomogène cette fois ; c'est-à-dire, travailler avec un modèle du type : $$Y_k = a + b X_k + \sigma_k \epsilon_k \, .$$ Nos $Y_k$ ne sont en effet rien de plus que des estimateurs de variance : $Y_k = \hat{\sigma}_k^2$ et leur variance est donnée par la formule du début de section : $\mathrm{Var}[Y_k] = 2 \sigma^4 / (n-1)$ que nous pouvons approcher par l'estimateur /plug-in/ (ou estimateur par substitution) : $$\sigma^2_{Y_k} \approx \frac{2 \hat{\sigma}^4}{(n-1)} = \frac{2 Y_k^2}{(n-1)} \, .$$
Si nous répétons nos calculs précédents en supposant les $\sigma_k^2$ connus nous arrivons rapidement à : $$\hat{a} = \frac{1}{Z} \sum_k \frac{Y_k-b X_k}{\sigma_k^2} \, ,$$ où $$Z = \sum_k \frac{1}{\sigma_k^2}$$ et $$\hat{b} = \left(\sum_k \frac{X_k}{\sigma_k^2} \left(y_k - \frac{1}{Z}\sum_j \frac{y_j}{\sigma_j^2}\right)\right) / \left(\sum_k \frac{X_k}{\sigma_k^2}\left(X_k - \frac{1}{Z}\sum_j \frac{X_j}{\sigma_j^2}\right)\right) \, .$$
Si nous codons cela nous arrivons à :
X_k = concatenate([x.flatten() for x in ADU_m])
y_k = concatenate([x.flatten() for x in ADU_v])
sigma2_k = 2*y_k**2/99
Z = sum(1/sigma2_k)
num1 = sum(X_k/sigma2_k*(y_k-sum(y_k/sigma2_k)/Z))
denom1 = sum(X_k/sigma2_k*(X_k-sum(X_k/sigma2_k)/Z))
b_chapeau = num1/denom1
a_chapeau = sum((y_k-b_chapeau*X_k)/sigma2_k)/Z
(a_chapeau,b_chapeau)
Soit :
G_0 = b_chapeau
sigma2_0 = a_chapeau / G_0**2
(G_0,sigma2_0)
Aventures avec le « modèle exact »¶
Nous avons en fait déjà écrit (implicitement) une version de notre « modèle exact » décrivant la loi de l'observation (compte) au pixel $(i,j) \in \{1,\ldots,I = 60\} \times \{1,\ldots,J = 80\}$ suite à la $k$-ième exposition ($k \in \{1,\ldots,K = 100\}$) de durée $\tau_d$ ($d \in \{1,\ldots,D = 10\}$) : $$Y_{i,j,k,d} \sim G \, p_{i,j} \tau_d + \sqrt{G^2 \, (p_{i,j} \tau_d+\sigma_{L}^2)} \, \epsilon_{i,j,k,d}\; ,$$ où $\epsilon_{i,j,k,d} \stackrel{\mathrm{IID}}{\sim} \mathcal{N}(0,1)$.
La densité de probabilité conjointe correspondante est alors : $$p\left((Y_{i,j,k,d}) = (y_{i,j,k,d})\right) = \prod_{i,j,k,d} \frac{1}{\sqrt{2 \pi}} \frac{1}{G} \sqrt{\frac{1}{p_{i,j} \tau_d+\sigma_{L}^2}} \exp\left(- \frac{1}{2 G^2} \frac{\left(y_{i,j,k,d} - G p_{i,j} \tau_d\right)^2}{\left(p_{i,j} \tau_d+\sigma_{L}^2\right)}\right) \, .$$
Si nous notons : $$\theta = \{G, \sigma_L^2, (p_{i,j})_{i \in \{1,\ldots,I\},j \in \{1,\ldots,J\}}\}$$ les paramètres de notre modèle, alors la log-vraisemblance s'écrit (après les simplifications d'usage) : $$l(\theta) = - 2 I J K D \log G - K \sum_{(i,j,d)} \log \left(p_{i,j} \tau_d+\sigma_{L}^2\right) - \frac{1}{G^2} \sum_{(i,j,k,d)} \frac{ \left(y_{i,j,k,d} - G p_{i,j} \tau_d\right)^2}{\left(p_{i,j} \tau_d+\sigma_{L}^2\right)} \, .$$
Nous pouvons définir : $$l_{i,j,k,d}(\theta) = -2 \log G - \log \left(p_{i,j} \tau_d+\sigma_{L}^2\right) - \frac{1}{G^2} \frac{ \left(y_{i,j,k,d} - G p_{i,j} \tau_d\right)^2}{\left(p_{i,j} \tau_d+\sigma_{L}^2\right)}$$ pour avoir : $$l(\theta) = \sum_{(i,j,k,d)} l_{i,j,k,d}(\theta) \, .$$ Cela va nous faciliter le calcul de la fonction score et de la fonction d'information pour l'application de la méthode de Newton-Raphson. Néanmoins, avant de nous lancer « tête baissée » dans les calculs, nous ne devrions pas oublier que nous travaillons ici avec un modèle à $60 \times 80 + 2 = 4802$ dimensions ! Avec un peu d'expérience, une dimension aussi grande devrait nous inciter à procéder par morceaux. Je m'explique : il est clair lorsque nous regardons les définitions de $l(\theta)$ et $l_{i,j,k,d}(\theta)$ ci-dessus que le bloc de la matrice d'information correspondant aux dérivées partielles par-rapport aux $p_{i,j}$ est diagonal. Il serait ainsi raisonnable de minimiser par-rapport à $(G,\sigma_{L}^2)$ d'une part et par-rapport aux $(p_{i,j})$ d'autre part et d'alterner entre les deux. Cela nous évitera d'avoir à inverser une matrice $4802 \times 4802$.
Mais comme ces calculs sont « un peu lourds », nous allons les faire avec SymPy. Nous commençons par déclarer les variables requises :
G, p_ij, tau_d, s2, y_ijkd = sy.symbols('G, p_ij, tau_d, s2, y_ijkd', real=True, positive=True)
Nous définissons ensuite une fonction de log-vraisemblance « élémentaire » (correspondant aux $l_{i,j,k,d}(\theta)$ ci-dessus) :
l_ijkd = -2*sy.log(G) -sy.log(p_ij*tau_d+s2) - (y_ijkd-G*p_ij*tau_d)**2/G**2/(p_ij*tau_d+s2)
l_ijkd
La dérivée partielle de $l_{i,j,k,d}(\theta)$ par-rapport à $G$ est :
S_ijkd_G = sy.diff(l_ijkd,G)
S_ijkd_G
Et nous pouvons obtenir une expression facile à incorporer dans un code Python numérique avec :
print(S_ijkd_G)
La dérivée partielle de $l_{i,j,k,d}(\theta)$ par-rapport à $s2$ est :
S_ijkd_s2 = sy.diff(l_ijkd,s2)
S_ijkd_s2
Et nous pouvons obtenir une expression facile à incorporer dans un code Python numérique avec :
print(S_ijkd_s2)
La dérivée partielle de $l_{i,j,k,d}(\theta)$ par-rapport à $p_{ij}$ est :
S_ijkd_p_ij = sy.diff(l_ijkd,p_ij)
S_ijkd_p_ij
Et nous pouvons obtenir une expression facile à incorporer dans un code Python numérique avec :
print(S_ijkd_p_ij)
Attention ici parmi les IxJ dérivées partielles de $l_{i,j,k,d}(\theta)$ du type : $\partial l_{i,j,k,d}(\theta) / \partial p_{u,v}$, une seule est non nulle celle pour laquelle $u=i$ et $v=j$.
p_uv = sy.symbols('p_uv', real=True, positive=True)
sy.diff(l_ijkd,p_uv)
Nous pouvons poursuivre en générant les « 6 éléments » de la matrice d'information :
I_ijkd_G_G = - sy.diff(S_ijkd_G,G)
I_ijkd_G_G
print(I_ijkd_G_G)
I_ijkd_G_s2 = - sy.diff(S_ijkd_G,s2)
I_ijkd_G_s2
print(I_ijkd_G_s2)
I_ijkd_G_p_ij = - sy.diff(S_ijkd_G,p_ij)
I_ijkd_G_p_ij
print(I_ijkd_G_p_ij)
La remarque précédente s'applique ici : parmi les IxJ dérivées partielles du type $\partial^2 l_{i,j,k,d}(\theta) / \partial G \partial p_{u,v}$, une seule est non nulle celle pour laquelle $u=i$ et $v=j$.
I_ijkd_s2_s2 = - sy.diff(S_ijkd_s2,s2)
I_ijkd_s2_s2
print(I_ijkd_s2_s2)
I_ijkd_s2_p_ij = - sy.diff(S_ijkd_s2,p_ij)
I_ijkd_s2_p_ij
print(I_ijkd_s2_p_ij)
La remarque précédente s'applique ici : parmi les IxJ dérivées partielles du type $\partial^2 l_{i,j,k,d}(\theta) / \partial s2 \partial p_{u,v}$, une seule est non nulle celle pour laquelle $u=i$ et $v=j$.
I_ijkd_p_ij_p_ij = - sy.diff(S_ijkd_p_ij,p_ij)
I_ijkd_p_ij_p_ij
print(I_ijkd_p_ij_p_ij)
La remarque précédente s'applique ici : parmi les $(I \times J)^2$ dérivées partielles du type $\partial^2 l_{i,j,k,d}(\theta) / \partial p_{i,j} \partial p_{u,v}$, seules sont non nulles celles pour laquelle $u=i$ et $v=j$. Autrement dit, le bloc constitué par les $I \times J$ dernières lignes et par les $I \times J$ dernières colonnes est diagonal.
Définitions des fonctions qui « font le travail »¶
Opposée de la log-vraisemblance¶
Nous commençons par définir une fonction constructrice nous renvoyant l'opposée de la log-vraisemblance dont l'expression est : $$-l(\theta) = 2 I J K D \log G + K \sum_{(i,j,d)} \log \left(p_{i,j} \tau_d+\sigma_{L}^2\right) + \frac{1}{G^2} \sum_{(i,j,k,d)} \frac{ \left(y_{i,j,k,d} - G p_{i,j} \tau_d\right)^2}{\left(p_{i,j} \tau_d+\sigma_{L}^2\right)} \, .$$
def fait_mlv_CCD(liste_de_piles=CCD_calibration_piles,
liste_de_temps=arange(1,11)*10):
D = len(liste_de_temps)
I,J,K = liste_de_piles[0].shape
def mlv_CCD(G,s2,p_ij):
terme1 = 2*I*J*K*D*log(G)
terme2 = 0.
terme3 = 0.
for d in range(D):
tau_d = liste_de_temps[d]
terme2 += sum(log(p_ij*tau_d+s2))
terme3 += sum((liste_de_piles[d]-\
G*tau_d*p_ij.reshape((I,J,1)))**2/(tau_d*p_ij.reshape((I,J,1))+s2))
return terme1+K*terme2+terme3/G**2
return mlv_CCD
Un détail important, dans le code ci-dessus, nous soustrayons des matrices G*tau_d*p_ij de piles liste_de_piles[d], techniquement nous faisons du broadcasting et nous devons spécifier suivant quelle dimension de la pile, « l'extension » de la matrice doit être effectuée. C'est pourquoi nous avons des instructions du type : G*tau_d*p_ij.reshape((I,J,1)).
Obtention de valeurs initiales pour les $p_{ij}$¶
Si nous moyennons nos moyennes à temps d'exposition fixé, pondérées par ce temps d'exposition et que nous divisons par notre valeur initiale pour $G$, nous obtenons une valeur initiale pour $p_{ij}$ : $$p_{i,j} \approx \frac{1}{G} \frac{1}{D} \sum_d \frac{\overline{y}_{i,j,d}}{\tau_d} \, .$$
Soit en Python :
tt = arange(1,11)*10
p_ij_0 = array([m/tau/G_0 for m,tau in zip(ADU_m,tt)]).mean(0)
Nous pouvons immédiatement visualiser les p\_ij\_0 avec :
imshow(transpose(p_ij_0),origin='lower')
set_cmap('gray')
colorbar()
contour(transpose(p_ij_0),colors='red')
Nous pouvons obtenir maintenant l'opposée de notre log-vraisemblance calculée avec le valeur initiale de chacun de nos paramètres :
mlvA = fait_mlv_CCD()
mlvA(G_0,sigma2_0,p_ij_0)
Fonction constructrice du Score et de l'information à $G$ et $\sigma_L^2$ fixés¶
Nous pouvons maintenant définir une fonction « constructrice » de la fonction score et de la fonction d'information pour les paramètres p_ij les deux autres étant considérés fixés :
def fait_S_I_p_ij(liste_de_piles=CCD_calibration_piles,
liste_de_temps=arange(1,11)*10):
I,J,K = liste_de_piles[0].shape
D = len(liste_de_temps)
resultat_S = zeros((I,J,1))
resultat_I = zeros((I,J,1))
fluo = zeros((I,J,1))
def Score(p_ij,G,s2):
resultat_S[:,:,0] = 0.
fluo[:,:,0] = p_ij
for d in range(D):
tau_d = liste_de_temps[d]
y_ijkd = liste_de_piles[d]
resultat_S[:,:,:] += -K*tau_d/(fluo*tau_d + s2) + \
2*tau_d*sum(-G*fluo*tau_d + y_ijkd,2).reshape((I,J,1))/(G*(fluo*tau_d + s2)) + \
tau_d*sum((-G*fluo*tau_d + y_ijkd)**2,2).reshape((I,J,1))/(G**2*(fluo*tau_d + s2)**2)
return resultat_S.reshape((I,J))
def Information(p_ij,G,s2):
resultat_I[:,:,0] = 0.
fluo[:,:,0] = p_ij
for d in range(D):
tau_d = liste_de_temps[d]
y_ijkd = liste_de_piles[d]
resultat_I[:,:,:] += 2*K*tau_d**2/(fluo*tau_d + s2) - \
K*tau_d**2/(fluo*tau_d + s2)**2 + \
4*tau_d**2*sum(-G*fluo*tau_d + y_ijkd,2).reshape((I,J,1))/(G*(fluo*tau_d + s2)**2) + \
2*tau_d**2*sum((-G*fluo*tau_d + y_ijkd)**2,2).reshape((I,J,1))/(G**2*(fluo*tau_d + s2)**3)
return resultat_I.reshape((I,J))
return (Score,Information)
Nous pouvons immédiatement effectuer un petit test :
ScoreA, InfoA = fait_S_I_p_ij()
p_ij = p_ij_0.copy()
p_ij += ScoreA(p_ij,G_0,sigma2_0)/InfoA(p_ij,G_0,sigma2_0)
mlvA(G_0,sigma2_0,p_ij_0) - mlvA(G_0,sigma2_0,p_ij)
Fonction constructrice du Score et de l'information à $p_{ij}$ fixés¶
Nous pouvons maintenant définir une fonction « constructrice » de la fonction score et de la fonction d'information pour les paramètres $G$ et $\sigma_L^2$ les $p_{ij}$ autres étant considérés fixés :
def fait_S_I_G_s2(liste_de_piles=CCD_calibration_piles,
liste_de_temps=arange(1,11)*10):
I,J,K = liste_de_piles[0].shape
D = len(liste_de_temps)
resultat_S = zeros((2,))
resultat_I = zeros((2,2))
fluo = zeros((I,J,1))
def Score(G,s2,p_ij):
resultat_S[0] = -2*I*J*K*D/G
resultat_S[1] = 0.
fluo[:,:,0] = p_ij
p_tau = zeros((I,J,1))
for d in range(D):
tau_d = liste_de_temps[d]
p_tau[:,:,:] = fluo*tau_d
y_ijkd = liste_de_piles[d]
resultat_S[0] += sum(2*p_tau*(-G*p_tau + y_ijkd)/(G**2*(p_tau + s2)) \
+ 2*(-G*p_tau + y_ijkd)**2/(G**3*(p_tau + s2)))
resultat_S[1] += sum(-1/(p_tau + s2) \
+ (-G*p_tau + y_ijkd)**2/(G**2*(p_tau + s2)**2))
return resultat_S
def Information(G,s2,p_ij):
resultat_I[:,:] = 0.
fluo[:,:,0] = p_ij
p_tau = zeros((I,J,1))
non_diagonal = 0.
for d in range(D):
tau_d = liste_de_temps[d]
p_tau[:,:,:] = fluo*tau_d
y_ijkd = liste_de_piles[d]
resultat_I[0,0] += sum(2*p_tau**2/(G**2*(p_tau + s2)) \
- 2/G**2 \
+ 8*p_tau*(-G*p_tau + y_ijkd)/(G**3*(p_tau + s2)) \
+ 6*(-G*p_tau + y_ijkd)**2/(G**4*(p_tau + s2)))
non_diagonal += sum(2*p_tau*(-G*p_tau + y_ijkd)/(G**2*(p_tau + s2)**2) \
+ 2*(-G*p_tau + y_ijkd)**2/(G**3*(p_tau + s2)**2))
resultat_I[1,1] += sum(-1/(p_tau + s2)**2 \
+ 2*(-G*p_tau + y_ijkd)**2/(G**2*(p_tau + s2)**3))
resultat_I[0,1] = non_diagonal
resultat_I[1,0] = non_diagonal
return resultat_I
return (Score,Information)
Et nous effectuons un test :
ScoreB, InfoB = fait_S_I_G_s2()
pGS = array([G_0,sigma2_0])
pGS += solve(InfoB(pGS[0],pGS[1],p_ij),ScoreB(pGS[0],pGS[1],p_ij))
mlvA(G_0,sigma2_0,p_ij) - mlvA(pGS[0],pGS[1],p_ij)
Première optimisation¶
Nous allons commencer par effectuer 100 itérations dont chacune voit se succéder un pas de Newton-Raphson « suivant » les $p_{ij}$ puis un pas de Newton-Raphson suivant le couple $(G,\sigma_L^2)$. Si nous trouvons le temps long, pour nous assurer que tout fonctionne « bien », nous pouvons afficher la valeur de l'opposée de la log-vraisemblance -- qui devrait decroitre -- à chaque itération :
pIJ = p_ij_0
pGS = array([G_0,sigma2_0])
mlv = mlvA(pGS[0],pGS[1],pIJ)
evol = zeros((101,4803))
evol[0,0] = mlv
evol[0,1:3] = pGS
evol[0,3:] = pIJ.flatten()
for i in range(100):
pIJ += ScoreA(pIJ,pGS[0],pGS[1])/InfoA(pIJ,pGS[0],pGS[1])
pGS += solve(InfoB(pGS[0],pGS[1],pIJ),ScoreB(pGS[0],pGS[1],pIJ))
mlv = mlvA(pGS[0],pGS[1],pIJ)
print(mlv)
evol[i+1,0] = mlv
evol[i+1,1:3] = pGS
evol[i+1,3:] = pIJ.flatten()
Nous préparons alors des graphes montrant l'évolution de l'opposée de la log-vraisemblance, de $G$, de $\sigma_L^2$ et d'un $p_{ij}$ central :
subplot(221)
plot(evol[:,0])
grid()
title(u'Opposée de l')
xlabel(u'Itération')
subplot(222)
plot(evol[:,1])
grid()
title(r'$G$')
xlabel(u'Itération')
subplot(223)
plot(evol[:,2])
grid()
title(r'$\sigma_L^2$')
xlabel(u'Itération')
subplot(224)
plot(evol[:,2402])
grid()
title(r'un $p_{ij}$')
xlabel(u'Itération')
On voit que l'opposée de la log-vraisemblance atteint essentiellement un plateau très vite de même que $\hat{\sigma}_L^2$, mais ce paramètre semble atteindre des valeurs incompatibles avec ce que nous montre notre premier graphe (haut de la page) puisqu'après notre optimisation du modèle complet, nous prédisons une ordonnée à l'origine de :
evol[100,2]*evol[100,1]**2
ce qui est très grand comparées aux valeurs de l'ordre de 5 d'où nous étions parties…
Nous poursuivons par la construction du graphe qui doit toujours être construit dans ce type d'analyse : les résidus normalisés (c'est-à-dire divisés par la déviation standard) en fonction de la valeur prédite :
for i in range(10):
scatter(pGS[0]*pIJ.flatten()*tt[i],
(CCD_calibration_piles[i][:,:,0].flatten()-pGS[0]*pIJ.flatten()*tt[i])/sqrt(pGS[0]**2*(pIJ.flatten()*tt[i]+pGS[1])),
0.05,color='black')
xlabel(u'Valeur prédite')
ylabel(u'Résidu normalisé')
grid()
Les ordonnées devraient être des tirages aléatoires suivant une loi normale centrée réduite et ce n'est clairement pas le cas pour les « petites » valeurs prédites… Que se passe-t-il ?
Nos temps d'exposition sont-ils bons ? Continuons par où nous aurions dû commencer !¶
Nous constatons que nos observation sont systématiquement au-dessus de nos prédictions pour les plus petits temps d'exposition. La question de la fiabilité de ces temps se pose. Si nos temps étaient corrects, alors la moyenne des comptes totaux du capteurs (calculée sur les 100 expositions à un temps donné) divisée par ce temps devrait être constante. Est-ce bien le cas ?
ADU_Mn = array([m/tau for m,tau in zip(ADU_m,[10,20,30,40,50,60,70,80,90,100])])
scatter(tt,ADU_Mn[:,0,0])
scatter(tt,ADU_Mn[:,0,79])
scatter(tt,ADU_Mn[:,59,0])
scatter(tt,ADU_Mn[:,59,79])
scatter(tt,ADU_Mn[:,29,39])
scatter(tt,ADU_Mn.mean((1,2)),color='red')
grid()
axhline(mean(ADU_Mn.mean((1,2))[-4:]))
xlabel(u"Temps d'exposition nominal (ms)")
Le graphe montre le compte moyen normalisé (c'est-à-dire divisé par le temps d'exposition nominal) pour les 4 pixels des 4 coins, un pixel central et la moyenne calculée sur le capteur (rouge). Le temps nominal ne peut pas être le temps effectif.
Nouvelle analyse avec des temps d'exposition corrigés¶
Nous pourrions redéfinir le modèle en laissant libre les rapport entre le 9 temps d'exposition les plus longs et le premier temps ; mais nous ne gagnerions pas grand chose (vous pouvez le faire en excercice) et nous allons estimer les temps corriger avant de refaire notre ajustement.
tau2 = tt*ADU_Mn.mean((1,2))/mean(ADU_Mn.mean((1,2))[-4:])
ScoreA2, InfoA2 = fait_S_I_p_ij(liste_de_temps=tau2)
ScoreB2, InfoB2 = fait_S_I_G_s2(liste_de_temps=tau2)
mlv2 = fait_mlv_CCD(liste_de_temps=tau2)
pIJ = p_ij_0
pGS = array([G_0,sigma2_0])
mlv = mlv2(pGS[0],pGS[1],pIJ)
evol2 = zeros((101,4803))
evol2[0,0] = mlv
evol2[0,1:3] = pGS
evol2[0,3:] = pIJ.flatten()
for i in range(100):
pIJ += ScoreA2(pIJ,pGS[0],pGS[1])/InfoA2(pIJ,pGS[0],pGS[1])
pGS += solve(InfoB2(pGS[0],pGS[1],pIJ),ScoreB2(pGS[0],pGS[1],pIJ))
mlv = mlv2(pGS[0],pGS[1],pIJ)
print(mlv)
evol2[i+1,0] = mlv
evol2[i+1,1:3] = pGS
evol2[i+1,3:] = pIJ.flatten()
Nous préparons alors des graphes montrant l'évolution de l'opposée de la log-vraisemblance, de $G$, de $\sigma_L^2$ et d'un $p_{ij}$ central :
subplot(221)
plot(evol2[:,0])
grid()
title(u'Opposée de l')
xlabel(u'Itération')
subplot(222)
plot(evol2[:,1])
grid()
title(r'$G$')
xlabel(u'Itération')
subplot(223)
plot(evol2[:,2])
grid()
title(r'$\sigma_L^2$')
xlabel(u'Itération')
subplot(224)
plot(evol2[:,2402])
grid()
title(r'un $p_{ij}$')
xlabel(u'Itération')
Les valeurs de $\hat{\sigma}_L^2$ sont beaucoup plus raisonnables ! Nous poursuivons par la construction du graphe qui doit toujours être construit dans ce type d'analyse : les résidus normalisés (c'est-à-dire divisés par la déviation standard) en fonction de la valeur prédite :
for i in range(10):
scatter(pGS[0]*pIJ.flatten()*tau2[i],
(CCD_calibration_piles[i][:,:,0].flatten()-\
pGS[0]*pIJ.flatten()*tau2[i])/sqrt(pGS[0]**2*(pIJ.flatten()*tau2[i]+pGS[1])),
0.05,color='black')
xlabel(u'Valeur prédite')
ylabel(u'Résidu normalisé')
grid()
C'est beaucoup mieux !
Les paramètres $\hat{G}$ et $\hat{\sigma}_L^2$ que nous allons (provisoirement) conserver sont donc :
G_chapeau = pGS[0]
S2_chapeau = pGS[1]
(G_chapeau,S2_chapeau)
Rappels sur la stabilisation de la variance¶
Soit $Y$ une v.a. (approximativement) normale d'espérance $\mu_Y$ et de variance $\sigma_Y^2$ : $Y \approx \mu_Y + \sigma_Y \, \epsilon$ avec $\epsilon \sim \mathcal{N}(0,1)$.
Soit $f$ une fonction dérivable et posons $Z = f(Y)$. Un développement limité donne alors : $$Z = f(\mu_Y + \sigma_Y \, \epsilon) \approx f(\mu_Y) + f'(\mu_Y)\, \sigma_Y \, \epsilon + o(\sigma_Y \, \epsilon) \, .$$ D'où $$\mathrm{E}Z \approx f(\mu_Y) = f(\mathrm{E}Y)$$ et $$\mathrm{Var}Z \equiv \mathrm{E}[(Z-\mathrm{E}Z)^2] \approx f'(\mu_Y)^2\, \sigma_Y^2 \, .$$ Soit $$Z \approx f(\mu_Y) + \mid f'(\mu_Y) \mid \, \sigma_Y \, \epsilon \, .$$
Dans le cas qui nous occupe :
- notre modèle de capteur CCD appliqué à un pixel s'écrit : $$Y \sim G \, \lambda + \sqrt{G^2 \, (\lambda+\sigma_{L}^2)} \, \epsilon = \mu_Y + \sqrt{G \, \mu_Y + G^2 \, \sigma_{L}^2} \epsilon \, ;$$
- si $Z = f(Y)$ alors : $$Z \approx f(\mu_Y) + \mid f'(\mu_Y) \mid G \sqrt{\mu_Y / G+\sigma_{L}^2} \, \epsilon \, ;$$
- qu'arrive-t-il si nous choisissons : $f(x) = 2 \, \sqrt{x/G + \sigma_{L}^2}$ ?
- nous avons : $$f'(x) = \frac{1}{G \sqrt{ x / G + \sigma_{L}^2}} \, ;$$
- ce qui donne : $$Z \approx 2 \, \sqrt{\mu_Y / G + \sigma_{L}^2} + \epsilon \, .$$
Appliqué à nos données de calibration cela donne :
ADU_mB = [mean(2*sqrt(pile/G_chapeau+S2_chapeau),2) for pile in CCD_calibration_piles]
ADU_vB = [var(2*sqrt(pile/G_chapeau+S2_chapeau),2) for pile in CCD_calibration_piles]
for i in range(10):
scatter(ADU_mB[i].flatten(),ADU_vB[i].flatten(),0.05)
scatter([mean(m) for m in ADU_mB],[mean(m) for m in ADU_vB],color='red')
grid()
xlabel(u'Moyenne des données transformées')
ylabel(u'Variance des données transformées')
Il y a encore quelque chose qui n'est pas complètement satisfaisant. C'est probablement la façon grossière dont nous avons « corrigé » les temps d'exposition. Pour faire les choses correctement, nous devrions considérer la matrice $p_{ij} \tau_{max}$ comme un paramètre, appelons le $\lambda_{ij}$ et exprimer les autres $p_{ij} \tau_{d}$ sous la forme $\alpha_d \lambda_{ij}$, les $\alpha_d$ devenant des paramètres du modèle. Nous allons néanmoins ici procéder plus rapidement en définissant une fonction de coût à minimiser et utiliser une méthode « automatique » de minimisation :
def cout(par):
G = par[0]
s2 = par[1]
V = sum([sum((var(2*sqrt(pile/G+s2),2)-1)**2) \
for pile in CCD_calibration_piles])
return sum((V-1)**2)
Nous allons effectuer la minimisation avec la fonction minimize du module scipy.optimize et la méthode de Nelder-Mead :
from scipy.optimize import minimize
ajuste = minimize(cout,[G_chapeau,S2_chapeau],method='nelder-mead',options={'disp': True})
Cela nous fournit de nouvelles valeurs pour nos paramètres :
G_chapeau = ajuste.x[0]
S2_chapeau = ajuste.x[1]
(G_chapeau,S2_chapeau)
ainsi qu'une nouvelle version du graphe précédent :
ADU_mB = [mean(2*sqrt(pile/G_chapeau+S2_chapeau),2) for pile in CCD_calibration_piles]
ADU_vB = [var(2*sqrt(pile/G_chapeau+S2_chapeau),2) for pile in CCD_calibration_piles]
for i in range(10):
scatter(ADU_mB[i].flatten(),ADU_vB[i].flatten(),0.05)
scatter([mean(m) for m in ADU_mB],[mean(m) for m in ADU_vB],color='red')
grid()
xlabel(u'Moyenne des données transformées')
ylabel(u'Variance des données transformées')
Détection automatique de région d'intérêt¶
Nous reprenons maintenant notre étude au niveau de la diapo 46 du cours d'« étude de cas » avec lequel nous avions commencé le 5 mai dernier. Les données qui nous intéressent sont téléchargées et chargées avec les commandes :
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()
Nous pouvons les visualiser avec :
imshow(transpose(pile_POMC.mean(2)),origin='lower')
colorbar()
Mise en pratique¶
Sur la pile $ADU_{i,j,k}$ de données POMC :
- après la stabilisation de variance, $Z_{i,j,k} = 2 \, \sqrt{ADU_{i,j,k} / G + \sigma_{L}^2}$, la variance des observations à chaque pixel $(i,j)$, à chaque temps $k$, doit être égale à 1 ;
- si le pixel « ne répond pas à la stimulation », sont intensité de fluorescence doit être constante et la statistique : $$RSS_{i,j} \equiv \sum_{k=1}^{K} (Z_{i,j,k} - \overline{Z}_{i,j})^2 \quad \mathrm{où} \quad \overline{Z}_{i,j} \equiv \frac{1}{K} \sum_{k=1}^{K} Z_{i,j,k}$$ doit suivre une loi du $\chi^2$ à $K-1$ degrés de liberté ;
- nous pouvons donc calculer la valeur de la fonction de survie d'une loi du $\chi_{K-1}^2$ : $$1 - F_{\chi_{K-1}^2}(RSS_{i,j})$$ et chercher les pixels donnant des valeurs très petites -- des probabilités très faibles -- ; l'utilisation d'un échelle log est recommandée !
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')
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é.