Brevet FR2942040A1 — publié le 20100813

PROCEDE D'ESTIMATION DES PARAMETRES DE LA DISTRIBUTION DES TEMPS DE REPONSE DE PARTICULES D'UN SYSTEME, APPLIQUE NOTAMMENT AUX MESURES DE FLUORESCENCE La présente invention concerne un procédé d'estimation de la distribution des temps de réponse de particules d'un système, dans le but notamment de le caractériser. Elle s'applique par exemple aux mesures de fluorescence résolues en temps. Des mesures de fluorescence permettent de caractériser notamment des molécules. Plus particulièrement, on éclaire des molécules par des lasers, puis les molécules ainsi excitées libèrent des photons. Les temps entre l'instant d'éclairement par le laser et les instants d'arrivée des photons sur un détecteur permettent de caractériser ces molécules de façon connue. Le nombre de photons émis lors de chaque éclairement est indéterminé, à titre indicatif plus l'impulsion émise par le laser est forte, plus il y a de photons. II y a par ailleurs un aléa supplémentaire sur les temps de réponse. Un problème technique général est la mesure de la distribution des temps de réponse d'un système donné dans le but de le caractériser à l'image de l'exemple précédent. La distribution de ces temps de réponse est généralement paramétrée par un ensemble de paramètres numériques et ce sont ces paramètres que l'on veut mesurer. On notera par la suite cette dépendance envers un jeu de paramètre en appelant fo., la densité de probabilité de la distribution des temps de réponse, 90 désignant un vecteur de paramètres appartenant à un ensemble O de vecteurs O. A tout élément 9 dans O on associe une densité de probabilité fe, et ainsi on définit une famille de distribution {fe, 06O}. D'un point de vue statistique, cette approche est un problème d'estimation paramétrique. Par exemple, dans le contexte de la fluorescence résolue en temps, où le système étudié est alors un ensemble de molécules fluorescentes, on suppose que la fonction de distribution, paramétrée par un jeu de paramètres, est un mélange d'exponentielles et {fe, 9e0} est l'ensemble des mélanges de densités exponentielles. La modélisation est plus complexe lorsqu'il s'agit de la distribution des temps de réponse d'un système optique ou d'un milieu transparent comme l'atmosphère. Cependant l'extension de cette approche paramétrique à d'autres domaines que la fluorescence résolue en temps est possible. Les mesures des temps de réponse sont obtenues au moyen d'une excitation pulsée du système à caractériser. On émet une excitation très brève, on mesure le temps mis par le système à répondre à cette excitation et on enregistre ce temps. Ces trois opérations sont répétées un très grand nombre de fois jusqu'à ce qu'on dispose de milliers ou de millions de temps de réponse. Enfin, c'est à partir de cette collection de temps de réponse collectés indépendamment que l'on estime les paramètres de leur distribution commune. Le domaine de l'invention concerne les contextes où le système peut émettre un nombre aléatoire de réponses à l'excitation et où seule la première réponse est détectée et enregistrée. Dans le contexte de la fluorescence résolue en temps, l'excitation est une impulsion laser envoyée sur l'échantillon fluorescent et la réponse de cet échantillon est constituée en général de plusieurs photons, car il y a un très grand nombre de molécules fluorescentes présentes. On enregistre seulement le temps de réponse correspondant au premier photon car ce premier photon aveugle l'appareillage électronique de détection et d'enregistrement. Ce phénomène d'oubli des deuxième, troisième photons et des suivants est appelé empilement ou � pile-up � dans la littérature anglo-saxonne. L'invention concerne donc les contextes qui comportent notamment les quatre caractéristiques suivantes : Mesure de la distribution de temps de réponse d'un système, Distribution paramétrée, Mise à disposition d'une collection généralement grande de temps de réponse mesurés indépendamment, Le système étudié émet un nombre aléatoire de réponses à chaque 30 excitation mais seule la première réponse est vue par l'appareillage de mesure. Ces quatre conditions sont rencontrées dans la méthode de mesure dite � Time Correlated Single Photon Counting � (TCSPC) en fluorescence résolue en temps. Des données du type � pile up � sont rencontrées en fluorescence résolue en temps. II s'agit d'estimer la distribution des durées de vie en fluorescence. La constante de durée de vie d'une molécule est le temps moyen entre l'instant d'excitation de la molécule et l'instant d'émission d'un photon. Elle contient des précieuses informations sur les molécules et leur milieu, par exemple le pH, la polarisation, ou la viscosité. Elle est donc de grand intérêt en biologie, en médecine ou en physique notamment. Typiquement, la distribution de la durée de vie d'une molécule est modélisée par une loi exponentielle. Or, la distribution des temps d'arrivée d'un mélange des molécules correspond à un mélange des lois exponentielles. En outre, le nombre de photons émis par impulsion de laser d'excitation est modélisé par une loi de Poisson. Les observations disponibles pour l'estimation des paramètres exponentiels peuvent être obtenues par la méthode TCSPC, aussi appelée la méthode SPT (Single Photon Timing), comme le montrent notamment les documents de B.Valeur � Molecular Fluorescence �, WILEY-VCH, Weinheim, 2002, de J.R.Lakowicz � Principles of Fluorescence Spectroscopy �, Academic/Plenium, New-York, 1999, ou de D.O'Connors et al � Time-Correlated Single Photon Counting �, Academic Press, London, 1984. Dans la méthode TCSPC un coup de laser excite un nombre aléatoire de photons, mais pour des raisons techniques seulement le temps de réponse du photon qui arrive le premier au détecteur peut être enregistré. Les autres photons fluorescents ne sont pas observables. En TCSPC le réglage de l'intensité de laser, ou l'atténuation de la chaîne optique, détermine le nombre moyen de photons détectés par excitation. Un problème technique est notamment le problème inverse. Plus particulièrement, on veut estimer les paramètres d'une densité de probabilité foi, , mais on ne dispose pas d'observations directes de la densité foi: . Pourtant on observe le minimum d'un nombre aléatoire des variables aléatoires indépendantes de la densité foo . II est évident que la densité des observations, notée 90,, n'est pas la même que la densité de probabilité foot En effet, yo0 résulte de foo par une distorsion non linéaire, effet de l'empilement � pile up �. On ne peut pas appliquer les mêmes estimateurs aux données sans � pile up �, issues de la densité foo , qu'aux données avec pile up, issues de la densité go o. Le problème d'estimation vient notamment de la forme compliquée de la densité 9or;. En effet, les estimateurs classiques, comme l'estimateur du maximum de vraisemblance ou les méthodes des moments, ne sont pas faisable. II existe au moins deux façons d'affronter le problème inverse où une observation représente le minimum d'un nombre aléatoire de variables aléatoires. En TCSPC il est courant de travailler avec une intensité de laser tellement faible qu'il y a au plus un photon émis par excitation avec une très forte probabilité. Donc, le problème d'une perte de photons, c'est-à-dire de ne pas observer tous les photons qui arrivent sur le détecteur, est évité. Dans ce cas, les temps d'arrivée observés peuvent être considérés comme des réalisations indépendantes d'un mélange exponentiel. Des estimateurs de paramètres du mélange exponentiel sont obtenues par des méthodes statistiques classiques comme la méthode des moindres carrées décrite dans le document de B.Valeur et al précité, ou la méthode du maximum d'entropie décrite dans le document de J.C.Brochon � Maximum entropy method of data analysis in time-resolved spectroscopy � Methods Enzymology, 240: 262-311, 1994. Or, une très basse intensité de laser entraîne un long temps d'acquisition, car sur 100 excitations il n'y a qu'un à trois photons détectés en général. Les autres excitations ne sont pas suivies par la détection d'un photon, donc elles ne contiennent aucune information sur les paramètres à estimer. Il faut alors un grand nombre d'excitations pour obtenir une quantité suffisante de temps d'arrivée pour une estimation fiable des paramètres. Une étude de la borne de Cramér-Rao montre qu'une minuscule intensité de laser entraîne une forte augmentation de la variance des estimateurs et qu'il est préférable de travailler avec un laser plus fort que ceux utilisés habituellement afin d'augmenter l'information sur les paramètres qui est contenue dans les données. Cependant, avec une plus forte intensité, les données subissent l'effet de � pile up � et il faut des méthodes statistiques adaptées au modèle avec � pile up �. La littérature sur le modèle statistique de � pile up �, où uniquement le minimum d'un nombre aléatoire de variables aléatoires est observé, est très limitée et, en TCSPC, ces méthodes sont rarement utilisées dans la pratique. Une première méthode, décrite dans le document de P.B.Coates � The correction of photon 'pile-up' in the measurement of relative lifetimes � Journal of Physics E : Scientific Instruments, 2(1) : 878-879, 1968, est une correction de l'histogramme pour des données en TCSPC. Il s'agit de repondérer les groupes ou classes de l'histogramme, encore appelés � bins �, afin d'obtenir un histogramme de la distribution cible fo() . Une généralisation de l'approche de P.B.Coates à une grande classe de distributions est présentée dans le document de A.Souloumiac � Apparatus and method for characterizing a system for counting elementary particles �, demande de brevet WO2007088173(A1). En plus, le document de A.Souloumiac étend la méthode de la repondération à la fonction de répartition afin de corriger l'effet de � pile up �. Dans le document de J.G.Walker � Iterative correction for 'pile-up' in single-photon lifetime measurement �, Optics Communications, 201 : 271-277, 2002, la méthode de P.B.Coates est adaptée au cas d'un laser non stable. Un inconvénient de ces approches de repondération est qu'elles ne fournissent pas d'estimateurs de paramètres. Après la repondération de l'histogramme ou de la fonction de répartition, il faut appliquer une autre méthode statistique pour arriver à des valeurs estimées des paramètres. La repondération ne représente qu'une première étape dans l'estimation. Des méthodes d'estimation de paramètres directe pour le modèle � pile up � ont été développées. Le document de L.Okano et al � Abstracting reliable parameters from time-correlated single photon counting experiments �, Journal of Photochemistry and Photobiology A : Chemistry, 175, 221-225, 2005 propose une méthode des moindres carrées pour ajuster directement la fonction de densité ge aux données avec empilement. Les méthodes des moindres carrées sont des � large sample methods � ce qui veut dire que l'on a besoin d'un très grand nombre de temps de réponse afin d'obtenir des estimateurs fiables. Comme nous cherchons à diminuer le temps d'acquisition, ce procédé n'est pas adapté. Une autre méthode directe est présentée dans le document de T. Rebafka � An MCMC approach for estimating a fluorescence lifetime with pile-up distorsion �, Actes 21 ième Colloque GRETSI sur le Traitement du Signal et des Images (GRETSI 2007), Troyes, France, où un échantillonneur de Gibbs est développé qui est bien adaptée au modèle de mélange exponentiel et à un nombre Poissonien de photons par excitation. Cette méthode présente cependant plusieurs inconvénients, en particulier elle nécessite un long temps de calcul et la difficulté d'inclure la fonction d'appareil dans le modèle. Elle ne peut donc pas être utilisée en TCSPC où les mesures contiennent généralement du bruit additif. Enfin les deux méthodes précédentes ne sont pas très souples dans le sens qu'elles sont restreintes au cas multi-exponentiel. Même si la loi exponentielle est la plus utilisée en fluorescence, il y a des applications avec des distribution de Kohlrausch, comme décrit dans le document de M.N.Berberan-Santos et al � Mathematical functions for the analysis of luminescence decays with underlying distribution : 1. Kohlrausch decay ~o function (stretched exponential), Chemical Physics, 315: 171-182, 2005, ou des distributions de Becquerel comme décrit dans le document de M.N.Berberan-Santos et al � Mathematical functions for the analysis of luminescence decays with underlying distribution : 2. Becquerel (compressed hyperbola) and related decay functions, Chemical Physics, 317: 57-62, 15 2005. Un but de l'invention est notamment de pallier les inconvénients précités en permettant d'obtenir une estimation directe des paramètres de la distribution des temps, adapté à l'effet de � pile-up �, ne nécessitant donc pas 20 l'atténuation de la puissance du laser utilisé. A cet effet, l'invention a pour objet un procédé d'estimation des paramètres 6o de la distribution des temps de réponse de particules réémises par un système, en réponse à une série d'émissions de particules vers ledit système, un temps de réponse d'une particule correspondant étant mesuré entre un instant d'émission et l'instant 25 de détection de ladite particule, le procédé comportant au moins : une étape de détermination d'une fonction approchant la log- vraisemblance de la distribution cible foo correspondant à une distribution sans effet d'empilement à partir des mesures de temps de réponse de particules détectées avec effet d'empilement, c'est-à-dire des 30 mesures de la distribution goo ; - une étape de maximisation de la log-vraisemblance approchée sur l'ensemble des paramètres O afin de déterminer l'estimateur comme le vecteur 6 qui maximise la log-vraisemblance approchée. La fonction approchée est une estimation de la log-vraisemblance de la 35 fonction de distribution cible foo basée sur des mesures observées distribuées selon la fonction de distribution 900 avec effet d'empilement. Supposons que Z1,...,Zn représentent les n mesures de temps de réponse. La fonction approchée est obtenue en pondérant chaque terme log(fe(Z;)) de la log-vraisemblance par un poids qui compense le fait que Z; suit une loi 9o0 5 et non foo à cause du phénomène du � pile-up �. On peut prendre un poids w*;,,, de la forme 1 rn in'(in-1 (1 - Gä(Z7))) où m note la fonction génératrice de N, ,x, in(u) = E P(V = `12.~21n. 2l E [0. 1] n-1 1 o et m' la dérivée de m et en est un estimateur de la fonction de répartition. Si l'on choisit l'estimateur empirique de la fonction de répartition pour ce dernier et si l'on réordonne les observations dans un ordre croissant sous la forme Z(1) La figure 1 illustre par deux courbes 1, 2 l'effet de � pile-up � évoqué précédemment dans un système d'axes, où les abscisses représentent les temps de réponse mesurés et l'axe des ordonnées le nombre de photons détectés correspondant à ces temps de réponse. Les valeurs indiquées sont normalisées. La première courbe 1 et la deuxième courbe 2 illustrent respectivement les densités de probabilité 9ûo et f�. évoquées précédemment. Dans la méthode TCSPC on éclaire donc un échantillon de molécules avec un rayon laser par exemple. Les molécules excitées par des photons issus du laser émettent à leur tour des photons qui sont détectés par des moyens appropriés, utilisant notamment des filtres adaptés. On mesure et on enregistre le temps d'arrivée des photons. On obtient ainsi une distribution des temps de réponse caractérisant les molécules. Les temps étant de l'ordre de plusieurs nanosecondes pour une précision exigée de l'ordre d'une picoseconde. En réalité, il n'y a pas un seul photon. II y a en fait deux aléas. Le premier aléa porte sur le nombre de photons réémis par les molécules et le deuxième aléa porte sur les temps de réponse, car chaque photon a son temps de réponse associé qui lui est propre et donc indéterminé. Plus le laser émet des impulsions fortes, plus il y a de photons réémis par les molécules. Les temps étant très courts, le chronomètre utilisé ne peut mesurer qu'une seule durée, la durée correspondant au temps du premier photon reçu. Le temps que le chronomètre se réinitialise pour engager une deuxième est alors trop long en regard du phénomène physique observé. La fluorescence est en effet terminée et tous les photons déjà reçus. A chaque impulsion laser, on a donc un certain nombre de photons qui arrivent mais seul le temps de réponse du premier est mesuré. La distribution obtenue n'est donc construite que pour le premier photon reçu, et donc pour le temps de réponse le plus court, conduisant ainsi à une image compressée et finalement à une distorsion de mesure. C'est l'effet d'empilement ou � pileup �, où seul le premier photon est détecté et mesuré. La fonction de densité associée aux mesures est illustrée par la première courbe 1. Cette courbe 1 présente une forte concentration de mesures correspondant aux temps de réponse très courts. Une solution connue consiste à baisser la puissance du faisceau laser émis de façon à ce qu'il n'y ait qu'un seul photon réémis à chaque impulsion laser. La distribution obtenue est satisfaisante car dénuée de distorsion. La deuxième courbe 2, de forme exponentielle, illustre la densité associée à cette solution. La différence de la première courbe 1 par rapport à cette deuxième courbe 2 illustre la distorsion provoquée par l'effet de � pile-up �. L'avantage des mesures correspondant à la courbe 2 est que l'on sait comment en extraire des paramètres e0 recherchés avec des méthodes statistiques, puisqu'il s'agit d'une densité exponentielle ou dans des cas plus complexes d'un mélange exponentiel. Cependant il s'ensuit que dans la grande majorité des cas, pratiquement dans environ 95% impulsions, il n'y a pas de photons réémis. Seules 2% à 5% des impulsions laser émises donnent lieu à la réémission d'un photon. L'obtention d'une distribution selon la deuxième courbe 2 se fait donc avec une perte de temps important. Il existe des applications où la distribution doit être obtenue rapidement et où les pertes de temps de la méthode illustrée par cette deuxième courbe 2 n'est pas admissible. En particulier, dans les cas où le système observé est instable et où plus la mesure est longue, plus il y a de possibilités que l'état du système évolue rendant l'ensemble des mesures inaptes à caractériser le système. II y a aussi des applications qui nécessitent un haut débit de mesures, notamment dans le domaine médical. Les figures 2a, 2b et 2c illustrent les histogrammes de mesures 21, 22, 23 avec l'effet de � pile-up � pour des intensités différentes. Comme indiqué précédemment, en TCSPC le réglage de l'intensité de laser, ou l'atténuation de la chaîne optique, détermine le nombre moyen de photons détectés par excitation. L'effet des choix différents de l'intensité sur la distribution des observations est illustré par les trois histogrammes des figures 2a, 2b et 2c correspondant à 1000 observations dans chacun des cas. La densité cible foo est une exponentielle de moyenne 1 et chaque histogramme est composé de 1000 observations, c'est-à-dire 1000 impulsions de laser, obtenues par simulation Monte-Carlo. Le dernier bin, noté +00, représenté sur l'échelle des abscisses, contient le nombre d'observations où aucun photon n'est détecté. L'intensité du laser est exprimée par le nombre moyen A de photons détectés par excitation. Plus précisement, A est le paramètre de la loi de Poisson qui est la distribution du nombre de photons détectés par excitation. Quand l'intensité est faible, A = 0,02, il n'y a que 13 photons détectés sur 1000 excitations, cas illustré par la figure 2a. Quand l'intensité est très élevée, A = 20, la distribution est presque dégénérée car presque toutes les observations se trouvent dans le premier bin, cas illustré par la figure 2b. Dans ces deux cas extrêmes, il est clair que les histogrammes contiennent très peu d'information sur le paramètre exponentiel, de la distribution des temps de réponse, à estimer. Il semble favorable d'utiliser une intensité moyenne. La figure 2c montre que A = 2 mène à un histogramme � habituel � car les observations où un photon a été détecté sont très nombreuses et les temps d'arrivée sont assez dispersés. Un but de l'invention est d'obtenir une distribution des temps de réponse conforme à la première courbe 1, sans perte de temps, c'est-à-dire obtenue beaucoup plus rapidement que la distribution selon la deuxième courbe 2. En particulier, la puissance des faisceaux laser émis n'est pas diminuée autant que pour la deuxième courbe 2. C'est notamment la résolution du problème inverse évoqué précédemment, conduisant à l'obtention des paramètres de la distribution selon la deuxième courbe 2. Dans le cas d'une simple exponentielle, par exemple, le paramètre obtenu est le paramètre r 1 e-,/r qui caractérise la densité exponentielle ~- A cet effet, l'invention comporte au moins deux étapes telles qu'illustrées par la figure 3. Dans une première étape 31 on détermine une fonction simple approchant la log-vraisemblance de la distribution cible, qui est la fonction f�o recherchée, illustrée par la deuxième courbe 2 de la figure 1. La log-vraisemblance, dépendant des mesures de temps de réponse, est une fonction en 0, le vecteur des paramètres qui caractérise la famille des distributions à laquelle appartient la densité cible f�.:, Puis dans une deuxième étape 32, on maximise cette log-vraisemblance 35 approchée en O. Rechercher les paramètres 00 de la fonction fo() revient en fait à rechercher le vecteur e des paramètres où la log-vraisemblance atteint son maximum. La maximisation est faite sur tout l'ensemble O des paramètres 9 qui détermine la famille des distributions à laquelle appartient la densité cible f�o Ainsi dans la première étape 31, on construit un estimateur de la fonction de log-vraisemblance E[log(fe(Y))J basée sur des données selon la densité 900, E étant la fonction d'espérance mathématique et Y étant une variable aléatoire de densité f�o c'est-à-dire E[log(fe(Y))J = $ log(fe(y)) fe(y)dy. Soient Z,,...,Zä des mesures des temps de réponse des photons. On note les mêmes mesures ordonnées par ordre croissant par Z(l),...,Z(0) C'est-à-dire que Z(l) désigne le plus petit temps de réponse mesuré, Z(2) le deuxième plus petit et ainsi de suite. L'estimation de la log-vraisemblance est obtenue par une repondération de la log-vraisemblance Ei=1 U'i,-z log fe(Z(,)) de la densité fe afin de corriger la distorsion, w;,,, représentant le poids associé à la donnée Zco. Dans cette phase de repondération, on donne un poids aux données Z(;) qui dépend de leurs positions temporelles. En d'autres termes, les données Z(;) sont ordonnées en fonction de leur valeur et leur poids dépend de leur ordre, c'est-à-dire que le poids des photons détectés Z; dépend de leur ordre d'arrivée. En particulier, le temps d'arrivée Z(n) le plus élevé est associé au plus grand poids dans la fonction Eiù1 uà.n logfe(Z(i) Dans la deuxième étape 32, on maximise la log-vraisemblance ainsi corrigée en utilisant un algorithme adapté, comme l'algorithme EM, abbréviation pour � Expectation-Maximization �, presenté dans le document de A.P.Dempster et al, � Maximum likelihood from incomplete data via the EM algorithm �, Journal of the Royal Statistical Society, Series B, 39(1), 1-38, 1977. On présente par la suite un exemple de méthode d'estimation pour une application avec effet de � pile-up � dans le cas simple d'un mélange de densités exponentielles. Puis on présente une méthode d'estimation pour une application avec une fonction d'appareil. Dans une approche statistique avec � pile up � pour des mélanges exponentiels, on suppose que la distribution des durées de vie en fluorescence représente un mélange de lois exponentielles. On considère la famille de densités (fe, ece) donnée par la relation (1) suivante : K fo(y) = cake lTk y e IR+ (1) ~ Tk où Tk>0 sont les paramètres exponentiels et les poids satisfont ak>0 et 10 Ek=i = 1.. On note 8 = (01 0k', T1, • Th')T le vecteur de paramètres et 0 l'ensemble de paramètres O. Les constantes de durées de vie en fluorescence sont données par les paramètres exponentiels Tk . En plus, K désigne l'ordre du mélange exponentiel qui correspond au nombre de molécules différentes. On suppose que K est connue et fixe. On note Fe la 15 fonction de répartition associée à la densité fe donnée par y Fo(y) = fo(t)dt o On note Y1,Y2,... une suite de variables aléatoires indépendantes de densité foo où e0cO est le paramètre inconnu à estimer. Les Yk représentent les temps d'arrivée des photons de fluorescence qui frappent le détecteur. 20 Le nombre N de photons qui frappent le détecteur par impulsion laser est modélisé par une loi de Poisson de paramètre A (restreinte à {1,2,...}). Les probabilités sont données par la relation suivante : n !( 1 ù eù.a ) n=1.2 25 On suppose que N est indépendant de la suite Y,, Y2, .... On définit alors l'observation avec � pile up � par Z =min{Y,,...,YN} la mesure retenue Z correspond, comme indiqué précédemment, au temps de réponse le plus court. Àne-A =n)= La fonction de répartition avec � pile-up � G00 de la variable aléatoire Z est donnée par la relation suivante : z R+ G00(z) = 1 ù c-a 1ùexp{ùÀFoo(z)} (2) et la densité avec � pile up � go() par la relation suivante : ÀfBo (z) 9B°(~) 1ùeù~ zER+ (3) Pour effectuer l'estimation, la méthode selon l'invention approche la log-vraisemblance de fe par un estimateur basé sur des observations avec � pile up �. Ensuite, on obtient un estimateur du paramètre eo en maximisant cette log-vraisemblance approchée. La forme simple de cette log-vraisemblance permet l'application de l'algorithme EM, défini par la suite, pour déterminer son maximum. On procède ensuite à l'approximation de la log-vraisemblance de fo. La log-vraisemblance de la fonction fo associée à un échantillon Y1,...,Ym de variables aléatoires indépendantes et identiquement distribuées (i.i.d) de densité foo est donnée par la relation suivante : L(0) = -Eloy (foO`Ç)) et l'estimateur du maximum de vraisemblance classique °ML est défini par la 20 relation suivante : m °ML = arg max Lo(8) = arg ma.x E log fo(Y) . BEN oee i=1 (4) Il découle de la loi forte des grands nombres que Lo(9) est un estimateur consistent de E[log(fe(Y))]. Si Z est une variable aléatoire de densité g0 , on établit la propriété suivante en utilisant relation (3), 1 IE [log(ff(Y))} = log( f�( y))foO(y)dy 0 log(fo(y))eÀ!ootngo0(y)dy IE lo- g(ff,(Z)) a~~(1 . A ' On en déduit que la quantité _ -a n L1(0; Op) = 1 nA E log fo(Z;)eakoo(Z,) est un estimateur consistant de E[log(fe(Y))] basé sur un échantillon Z1,...,Zä i.i.d. de densité goo. Pourtant, la quantité L,(8; 80) ne peut pas être utilisée pour estimer le paramètre inconnu e0, car elle-même dépend de 80. Or, une autre modification est indispensable. On déduit de la relation (2) que : e-ArGo(2) = 1 ù (1 ù c-`)Coo(') En remplaçant Cou( ) par la fonction de répartition empiriqueGn (zù n =i L { Z, 1 ùc-À log fo(Z(,)) 1 ùe-a A(1 ù (1 ù c-')n) où 4 désigne la i-ème statistique d'ordre de Z1,...,Z,,. On introduit les poids 1ùc-A u'i.n = À(1ù(1ùe-a)l) On définit un nouvel estimateur de 80 par la relation suivante : r, 891 = arg max L(y) = arg ruax> u}i,, log f0(Zri)) o{e o{e i=1 (5) Donc, on cherche le paramètre où la log-vraisemblance approchée atteint son maximum. On peut considérer On comme un estimateur de maximum de vraisemblance approché. On remarque que l'on peut également utiliser d'autres estimateurs de la fonction de répartition que Cn. Cela changera les poids w;,,, qui sont associés aux données Z(,), mais le principe restera le même. On effectue par la suite une maximisation par l'algorithme EM. On constate que le problème de maximisation de la relation (5) diffère du problème de la relation (4) seulement par des constantes positives w,.,,. Or, la maximisation de la relation (4) peut se faire par l'algorithme EM. Grâce à la structure semblable des problèmes, l'algorithme EM s'applique également à la relation (5). On introduit la variable manquante S d'une observation Y du mélange exponentiel fe étant l'étiquette de la composante du mélange qui a généré l'observation Y. Plus précisément, S est une variable aléatoire à valeurs dans (1,...,K) et 1F(S = /-) = aÀ. La loi jointe de (Y,S) est donnée par 'rd(y` s) = ase-u T, Tg y>O,se{l K}. Comme dans l'algorithme EM classique, on définit la fonction n et on calcule itérativement la suite (9('`))r en maximisant Q(e,(e(�) , conduisant à la relation suivante : 0(1+1) = argmaxQ(0.0t~') . oee Q(0. A') = E 2Ci nJE0 [log 7r�(Z(.i), S)IY = Z( i=1 n K _ E E(log(ak) ù log Tk ù Z. )P01(S = kilt = Z(o) Ty iù1 k=1 (6) Ce problème de maximisation est facile car il se décompose en plusieurs problèmes de maximisation indépendants: (t+1) Ti = arg max E u'i.n T,>o (_1og() 1Poa>(S = 11Y = Z(il) i=1 n K (Q t+1) ... ah+1) ) = arg ma,x E u'i.n E log( ak)1P�ut,(S = kil = Z(il) ni nh i=1 k=1 où la dernière maximisation est faite sous la contrainte que tous les ak sont à valeurs dans 10,1[ et Ikak = 1. On peut montrer qu'avec chaque itération la log-vraisemblance approchée L(0) augmente. Plus précisément on a pour toute valeur initiale el) e e et pour tout t ?1, L(0(t+1)) > L(8(t)) . Par ailleurs, on peut montrer que la suite (e(t))t converge vers un point 10 critique de L(6). Il y a plusieurs critères d'arrêt pour déterminer le moment quand la suite (6(0)t a convergé et où on peut arrêter le calcul itératif. Un critère repose sur l'évolution de la log-vraisemblance approchée. On calcule la valeur de L(e(t)) à toute itération et on s'arrête si l'accroissement par rapport à l'étape 15 précédente est inférieur à un certain seuil. Un autre critère étudie la variation des valeurs du paramètre e"°. A toute itération on évalue la différence Ilo(t) où Il . Il Y désigne la norme infinie. Si la différence est inférieure à un certain seuil, on considère la dernière valeur de 6(t) comme limite de la suite (e(�)t . 20 L'algorithme EM pour maximiser la fonction L(6).dans le cas d'un mélange exponentiel, peut être défini par les phases suivantes : On fixe un seuil a > 0 pour le critère d'arrêt. Soit A la valeur connue du paramètre Poissonien lors de l'expérience. 25 Soit Z1,...,Z, les observations avec � pile-up �, c'est-à-dire un échantillon i.i.d. de densité 9o,). On note Z(l),...,Z(n) les observations ordonnées, n étant le nombre total de mesures avec photons détectés, et i étant l'ordre d'arrivée du photon Z(;). 5 1. On pose t = 1. On choisit une valeur initiale 8(1) �ail) On calcule, pour i = 1,...,n 1 ù u'i.n= À(1-(1-e-,) 2. On réitère jusqu'à la convergence les étapes suivantes : (a) Calculer les probabilités conditionnelles, pour i = 1,...,n et s = 1,...,K , a(1) T(1) T(.1)) E e p,.i = Pr3w(S = s'Y _ Z(t}) = tt) an -z TF- T EK n' e z,, ,r`Tk if) . k=1Titt k , i Tgtt; (b) Déterminer les nouveaux paramètres, pour s = 1,...,K , (1+1) tVi.nPs.t Qg - En1 10 i=1 u t.n et T(t+1) = f::=1 ZIt)u'i.nPs.i s ii n i=1 u'i.nPs.i (r+1} ù t+1} (t+1) (1+1) it+1)) On pose 0(t+1) = (0'i ..,ah- T1 ,...,Th. (c) D'après le critère d'arrêt choisi, on évalue L(8(t+1 -) ù L(8('}) 0"r+1 ù 110(t+1) ù 0(t)jj : . (d) Si Qt+, _1) deux suites indépendantes de variables aléatoires indépendantes de loi q respectivement du mélange exponentiel foo donnée par la relation (1). On pose Yi = Yi Ei. La densité foo de Y est la convolution de foo et q, foo = foo On note la fonction de répartition de 1' par Foo. Or, une observation en TCSPC a la même distribution que : Z = mii1{11. .. ,1-'.v } . La fonction de répartition avec � pile-up � Goode Z est donnée par : Ger ( ) = 1 ù exp { ù À È00 (z) } ti E R+ 1 ù (ùa et la densité avec � pile-up � 9oäest donnée par : Af�n(z)eûaF't)tz, 9no(z)= zER+. 1 - De la même manière que précédemment, on déduit l'estimateur comme le vecteur e qui maximise la log-vraisemblance approchée, 8n = arg max L(A) . oee où Z1, .. . . Zn désigne des observations avec � pile-up � de densité gon et la log-vraisemblance approchée est donnée par n L(8) = z E "wi.n log fo(Z(il) iù125 Pour l'observation 1 = Y -1- E,, on considère les deux variables manquantes (S, E) où S est l'étiquette de la composante du mélange qui a généré Y, et E désigne la durée ajouté à la durée de vie par l'appareil. La loi jointe de ( S, E) est donnée par : os ù(üùu?~Ts 1ro(y• s,'u) = q(u)ùe TSO 0 pour le critère d'arrêt. Soit A la valeur connue du paramètre Poissonien lors de l'expérience. Soit ri une densité connue constante par morceaux donnée par p q(u) = E ci L {t'{[h3 1.~,)} .1=1 Soient Z1 Zä les observations avec pile up, un échantillon i.i.d. de densité gon. On note Z(1)- ... Ztn) les observations ordonnées 1. On pose t = 1 On choisit une valeur initiale e(1) _ ( 1), • • .0( 2), . • • • TFl� E H On calcule, pour i=1,...,n , 1 ù a(1ù(1ùc-~)) 2. On réitère jusqu'à convergence les étapes suivantes : (a) On calcule, pour i=1,...,n et s=1,...,K, (t) qs.i _ - r](u )e-(t' u) ré, du ù (t) Ts 0 p ((i EL {Z(t) ? hj_1}c j=1 = Q`(t)e_j!t E,(h,nlltT6t ù ~.i où a/1b dénote le minimum de a et b, et on calcule Ps.i = (S _ sii. = 2(i� = _ Qs.i E K k=1 qk.t t_l-~(Zl)eùfÎ~_,iùu}/re (iu P t[1 E rr jù1 1=1 tai x (Tvt) ù hi_1) e"' ni-t + ((h1 A Z(j),l ù TSt) (b) On détermine les nouveaux paramètres, pour s=1,...,K, n ~(t+1) _ Ei=1 Z~'i.nPs.i c ~` n Lei=1 et En w1.n Li_1 Et _1 (t+1) = ~_~ ek,t El=1 ui.nPs.t (t+1) ù (t+1) (t+l) (fil) (t+i) On pose 00+1) = (al ,...,Q.- ,Tl . .TK ). (c) D'après le critère d'arrêt choisi, on évalue Ï (oit+1)) _ L(B(t) Qt+1 = 0(t+1) _ y(t)~I (d) Si Qt+, 1 } une suite de variables aléatoires indépendantes de densité fco. On note Fo la fonction de répartition associée à la densité f�0. On considère N* une variable aléatoire à valeurs dans (1,2,...} indépendante de la suite { 1k, k ? 11. On note m la fonction génératrice de N*, m(2) _ E IP(N* = n)u.n. u E [0, 11 n-1 et m' la dérivée de m. On définit l'observation avec � pile-up � par : Z* = min{ i *, } . On établit la log-vraisemblance approchée basée sur des observations Z, ... , Zn i.i.d. avec � pile-up � : L*(0) = 1 E Z(1 n log(f�(Z(*0)) n iù1 n où 1 172.-1(1 ù n)) On définit n `2*(0, 0� = E uni ilE0' [1og(7r (Z(4), S))IY* = Z~1] i-1 où S dénote les variables manquante associées à la variable Y*, et ' désigne la densité jointe de (1*, S)• L'algorithme EM général pour maximiser Lx(0) comporte les phases suivantes : On fixe un seuil Q>0 pour le critère d'arrêt. 1. On choisit une valeur d'initialisation 0(') e e. On pose t = 1. 2. On réitère jusqu'à la convergence les étapes suivantes : (a) On poser t = t+1. (b) On calcule Ou) = arg 1Ï1âX e Q* (8 B't-1) ) (c) D'après le critère d'arrêt choisi, on évalue L(00)) - L*(O`tû1)) = 110(t) _ (d) Si ai Le système à caractériser peut être composé de molécules comme décrit 20 précédemment. L'invention peut néanmoins s'appliquer à d'autres systèmes. Un système à caractériser peut être par exemple une fibre optique dans laquelle on injecte des photons et dont mesure les durées de vie. D'autres types de particules que les photons peuvent aussi être utilisés. 22 25 REVENDICATIONS1. Procédé d'estimation des paramètres e0 de la distribution des temps de réponse de particules réémises par un système, en réponse à une série d'émissions de particules vers ledit système, un temps de réponse d'une particule correspondant étant mesuré entre un instant d'émission et l'instant de détection de ladite particule, caractérisé en ce que le procédé comporte au moins : une étape (31) de détermination d'une fonction approchant la log-vraisemblance de la distribution cible foo correspondant à une distribution sans effet d'empilement à partir des mesures de temps de réponse de particules détectées avec effet d'empilement, c'est-à-dire des mesures de la distribution gon ; une étape (32) de maximisation de la log-vraisemblance approchée sur 15 l'ensemble des paramètres O afin de déterminer l'estimateur comme le vecteur e qui maximise la log-vraisemblance approchée. 2. Procédé selon la revendication 1, caractérisé en ce que la fonction approchée est une estimation de la log-vraisemblance de la fonction de 20 distribution cible foä basée sur des mesures observées distribuées selon la fonction de distribution gr),, avec effet d'empilement. 3. Procédé selon la revendication 2, caractérisé en ce que la fonction approchée est obtenue par une pondération de la log-vraisemblance de la 25 fonction de distribution cible fo . 4. Procédé selon la revendication 3, caractérisé en ce que la pondération de la log-vraisemblance se fait de façon à ordonner les mesures de temps de réponse observées Z,,...,Zä en fonction de l'ordre d'arrivée des particules 30 correspondante, le poids affecté à chaque mesure dépendant de son ordre. 5. Procédé selon la revendication 4, caractérisé en ce que le poids affecté à une mesure croît avec son ordre d'arrivée associé. 6. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que n étant le nombre total de particules détectées et i étant l'ordre d'arrivée d'une particule, le poids w*;,,, associé à la mesure du temps de réponse Z(0 d'une particule d'ordre i est donné par la relation suivante : 1 t ni'(m-1(1 ù C (Z(=)))) où m note la fonction génératrice de N, rz(u) = E P(N = n)un, u E [0, 1] nù1 et m' la dérivée de m et Gn est un estimateur de la fonction de répartition. 7. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que la fonction de distribution cible fo. est considérée comme étant un mélange de densités exponentielles. 15 8. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que la fonction de distribution cible prend en compte la fonction d'appareil de mesures, la fonction cible étant le produit de convolution de la fonction cible sans fonction d'appareil (fo() ) avec la densité 20 de probabilité de la fonction d'appareil q. 9. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que les particules sont des photons émis par un laser, le système réémettant des photons. 10. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que le système est un ensemble de molécules à caractériser par la fonction de distribution cible des temps de réponse. 25 30