WO2010089363A1 - Method for estimating parameters of the distribution of particle response times of a system, in particular applied to fluorescence measurements - Google Patents

Abstract

The present invention relates to a method for estimating the distribution of particle response times of a system, in particular for the purpose of characterising said particles. The method is applicable, for example, to time-resolved fluorescence measurements. The method at least comprises: a step (31) of determining a function approximating the log-likelihood of the target distribution ?o corresponding to a distribution with no pile-up effect based on measurements of response times of particles detected with a pile-up effect, that is, measurements of the distribution g?o; a step (32) of maximising the approximated log-likelihood of all the parameters (?) in order to determine the estimator as the vector ? that maximises the approximated log-likelihood.

Description

 METHOD OF ESTIMATING PARAMETERS OF DISTRIBUTION OF PARTICLE RESPONSE TIMES

OF A SYSTEM APPLIES IN PARTICULAR TO FLUORESCENCE MEASUREMENTS

The present invention relates to a method for estimating the distribution of particle response times of a system, in particular for the purpose of characterizing it. It applies, for example, to time-resolved fluorescence measurements.

Fluorescence measurements make it possible to characterize molecules in particular. More particularly, molecules are illuminated by lasers, and the molecules thus excited release photons. The times between the moment of illumination by the laser and the instants of arrival of the photons on a detector make it possible to characterize these molecules in a known manner. The number of photons emitted during each illumination is indeterminate, as an indication, the stronger the pulse emitted by the laser, the more photons there are. There is also an additional hazard on response times. A general technical problem is the measurement of the distribution of the response times of a given system in order to characterize it in the image of the preceding example.

temps de réponse, θo désignant un vecteur de paramètres appartenant à un ensemble Θ de vecteurs θ. A tout élément θ dans Θ on associe une densité de probabilité et ainsi on définit une famille de distribution

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

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.The distribution of these response times is generally parameterized by a set of numerical parameters and it is these parameters that one wants to measure. This dependence on a parameter set will be noted later by calling the probability density of the distribution of

response time, θo designating a vector of parameters belonging to a set Θ of vectors θ. To every element θ in Θ one associates a density of probability and thus one defines a family of distribution

From a statistical point of view, this approach is a parametric estimation problem. For example, in the context of time-resolved fluorescence, where the studied system is then a set of fluorescent molecules, it is assumed that the distribution function, parameterized by a set of parameters, is a mixture of exponentials and

is the set of exponential density mixtures. Modeling is more complex when it comes to the distribution of the response times of a optical system or a transparent medium like the atmosphere. However, the extension of this parametric approach to domains other than time-resolved fluorescence is possible.

The measurements of the response times are obtained by means of pulsed excitation of the system to be characterized. One emits a very brief excitation, one measures the time put by the system to answer this excitation and one records this time. These three operations are repeated a very large number of times until thousands or millions of response times are available. Finally, it is from this collection of independently collected response times that the parameters of their common distribution are estimated.

The field of the invention relates to contexts where the system can issue a random number of responses to excitation and where only the first response is detected and recorded. In the context of time-resolved fluorescence, the excitation is a laser pulse sent on the fluorescent sample and the response of this sample is generally composed of several photons, because there is a very large number of fluorescent molecules present. Only the response time corresponding to the first photon is recorded because this first photon blinds the electronic detection and recording apparatus. This phenomenon of forgetting the second, third and subsequent photons is called stacking or "stack-up" in the Anglo-Saxon literature. The invention thus relates to the contexts which comprise in particular the four following characteristics: - Measurement of the response time distribution of a system,

- Parametric distribution,

- Provision of a generally large collection of independently measured response times,

- The studied system emits a random number of responses to each excitation but only the first response is seen by the measuring equipment.

These four conditions are met in the Time Correlated Single Photon Counting (TCSPC) measurement method in time-resolved fluorescence. Data of the "stack up" type are encountered in time-resolved fluorescence. This involves estimating the distribution of fluorescence lifetimes. The lifetime constant of a molecule is the mean time between the moment of excitation of the molecule and the instant of emission of a photon. It contains valuable information about molecules and their environment, such as pH, polarization, or viscosity. It is therefore of great interest in biology, medicine or physics in particular. Typically, the distribution of the lifetime of a molecule is modeled by an exponential law. Now, the distribution of the arrival times of a mixture of molecules corresponds to a mixture of exponential laws. In addition, the number of photons emitted by excitation laser pulse is modeled by a Poisson law.

The observations available for the estimation of the exponential parameters can be obtained by the TCSPC method, also called the SPT (Single Photon Timing) method, as shown in particular by the documents of B.Value "Molecular Fluorescence", WILEY-VCH, Weinheim, 2002, JRLakowicz "Principles of Fluorescence Spectroscopy", Academic / Plenium, New York, 1999, or D. O'Connors et al "Time-Correlated Single Photon Counting", Academy Press, London, 1984. In the method TCSPC a laser shot excites a random number of photons, but for technical reasons only the response time of the photon that first arrives at the detector can be recorded. Other fluorescent photons are not observable. In TCSPC the setting of the laser intensity, or the attenuation of the optical chain, determines the average number of photons detected by excitation.

En effet, les estimateurs classiques, comme l'estimateur du maximum de vraisemblance ou les méthodes des moments, ne sont pas faisable.A technical problem is the opposite problem. More particularly, we want to estimate the parameters of a probability density M), but we do not have direct observations of the density M>. However, we observe the minimum of a random number of random variables independent of the density hh. It is obvious that the density of the observations, denoted by ->'^, is not the same as the density of probability M'. Indeed, $ & $ results from h'h by a nonlinear distortion, effect of stacking up. The same estimators can not be applied to the "stack up" data from the density M "as to the stack up data from the density". The estimation problem comes especially the complicated form of density

Indeed, classical estimators, such as the maximum likelihood estimator or the moment methods, are not feasible.

There are at least two ways to deal with the opposite problem where an observation represents the minimum of a random number of random variables. In TCSPC it is common to work with a laser intensity so weak that there is at most one photon emitted by excitation with a very high probability. So, the problem of a loss of photons, that is to say not to observe all the photons arriving on the detector, is avoided. In this case, the observed arrival times can be considered as independent achievements of an exponential mixture. Exponential mixing parameter estimators are obtained by conventional statistical methods such as the least squares method described in B.Valeur et al., Or the entropy maximum method described in JCBrochon's document "Maximum entropy method of data analysis in time-resolved spectroscopy Methods Enzymology, 240: 262-311, 1994.

However, a very low laser intensity leads to a long acquisition time, because on 100 excitations there is only one to three photons detected in general. The other excitations are not followed by the detection of a photon, so they contain no information on the parameters to estimate. It then takes a large number of excitations to obtain a sufficient amount of arrival time for a reliable estimate of the parameters. A study of the Cramér-Rao kiosk shows that a tiny laser intensity leads to a sharp increase in the variability of the estimators and that it is preferable to work with a laser stronger than those usually used to increase the information on the parameters that is contained in the data. However, with a higher intensity, the data undergo the "stack up" effect and statistical methods adapted to the "stack up" model are needed.

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 « Itérative 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é 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 ^ιeme 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.The literature on the statistical model of "stack up", where only the minimum of a random number of random variables is observed, is very limited and, in TCSPC, these methods are rarely used in practice. A first method, described in the PBCoates paper "The correction of photon 'stack-up' in the measurement of relative lifetimes" Journal of Physics E: Scientific Instruments, 2 (1): 878-879, 1968, is a correction of the histogram for TCSPC data. It is a question of reweighting the groups or classes of the histogram, also called "bins", in order to obtain a histogram of the target distribution

A generalization of the PBCoates approach to a large class of distributions is presented in A.Souloumiac's document "Apparatus and method for characterizing a system for counting elementary particles", patent application WO2007088173 (A1). In addition, the paper by A.Souloumiac expands the method of reweighting the distribution function to correct the "stack up" effect. In the JGWalker document "Iterative correction for 'stack-up' in single-photon lifetime measurement," Optics Communications, 201: 271-27, 2002, the method of PBCoates is adapted to the case of an unstable laser. A disadvantage of these reweighting approaches is that they do not provide parameter estimators. After reweighting the histogram or distribution function, another statistical method must be used to arrive at estimated parameter values. Reweighting is only a first step in the estimation. Direct parameter estimation methods for the "stack up" model have been developed. The document by Okano et al, "Abstracting reliable parameters from time-correlated single photon counting experiments", Journal of Photochemistry and Photobiology A: Chemistry, 175, 221-225, 2005 proposes a least-squares method for directly adjusting the function of density to data with stacking. The methods of

least squares are "large sample methods" which means that we need a very large number of response times in order to obtain reliable estimators. As we seek to reduce the acquisition time, this method is not suitable. Another direct method is presented in T. Rebafka's "An MCMC approach for estimating a fluorescence lifetime with stack-up distortion", Proceedings of the 21st ^GRETSI Symposium on Signal and Image Processing (GRETSI 2007), Troyes, France. , where a Gibbs sampler is developed which is well adapted to the exponential mixing model and to a Poisson number of photons by excitation. However, this method has several disadvantages, in particular it requires a long calculation time and the difficulty in including the device feature in the model. It can not therefore be used in TCSPC where the measurements usually contain additive noise.

Finally the two previous methods are not very flexible in the sense that they are restricted to the multi-exponential case. Even if the exponential law is the most used in fluorescence, there are applications with Kohlrausch distributions, as described in the paper by MNBerberan-Santos et al "Mathematical functions for the analysis of luminescence decays with underlying distribution: 1. Kohlrausch decay function (stretched exponential), Chemical Physics, 315: 171-182, 2005, or Becquerel distributions as described in MNBerberan-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, 2005.

An object of the invention is in particular to overcome the aforementioned drawbacks by making it possible to obtain a direct estimate of the parameters of the time distribution, adapted to the "stack-up" effect, thus not requiring the attenuation of the power of the laser used. For this purpose, the subject of the invention is a method for estimating the parameters θ 0 of the distribution of the response times of particles reemitted by a system, in response to a series of particle emissions to said system, a response time a corresponding particle being measured between an emission instant and the instant of detection of said particle, the method comprising at least:

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

a step of determining a function approaching the log-likelihood of the target distribution corresponding to a

non-stacking distribution from measured particle response time measurements with stacking effect, i.e., measurements of the distribution

afin de déterminer l'estimateur comme le vecteur θ qui maximise la log-vraisemblance approchée.a step of maximizing the approximate log-likelihood on all the parameters

to determine the estimator as the vector θ that maximizes the approximate log-likelihood.

distribuées selon la fonction de distribution ^ = avec effet d'empilement. Supposons que Zi,...,Z_n représentent les n mesures de temps de réponse. La fonction approchée est obtenue en pondérant chaque terme log(f_θ(Zi)) de la log-vraisemblance par un poids qui compense le fait que Z, suit une loi ^- > et non à cause du phénomène du « pile-up ». On peut prendre un poids

de la forme 1

où m note la fonction génératrice de N,The approximate function is an estimate of the log-likelihood of the target distribution function based on observed measurements

distributed according to the distribution function ^ = with stacking effect. Suppose that Zi, ..., Z _n represent the n response time measurements. The approximate function is obtained by weighting each term log (f _θ (Zi)) of the log-likelihood by a weight which compensates for the fact that Z follows a law ^ -> and not because of the "stack-up" phenomenon. . We can take a weight

of form 1

where m denotes the generating function of N,

et m' la dérivée de m et &u 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

alors le poids

correspondant au terme log(f_θ(Z_(l))) s'écrit

and m 'the derivative of m and & u is an estimator of the distribution function. If we choose the empirical estimator of the distribution function for the latter and if we reorder the observations in increasing order in the form

then the weight

corresponding to the term log (f _θ (Z _(l) )) is written

est par exemple considérée comme étant un mélange de densités exponentielles.The target distribution function

is for example considered to be a mixture of exponential densities.

Advantageously, in a particular mode of implementation, the target distribution function takes into account the measuring device function, the target function being the product of convolution of the target function without device function (h \ ->) with the probability density of the device function η.

The particles are, for example, photons emitted by a laser, the system re-emitting photons.

The system can be a set of molecules to be characterized by the target distribution function of the response times.

Other characteristics and advantages of the invention will become apparent with the aid of the following description made with reference to appended drawings which represent: FIG. 1, an illustration of the effect of "stack-up" on the distributions of response time measurements; FIGS. 2a, 2b and 2c illustrate histograms of measurements with "stack-up" effect for different intensities of photon emission;

- Figure 3, an overview of the possible steps for the implementation of a method according to the invention.

et évoquées

précédemment.FIG. 1 illustrates by two curves 1, 2 the "stack-up" effect mentioned previously in a system of axes, where the abscissas represent the measured response times and the ordinate axis the number of detected photons corresponding to these response times. The values indicated are standardized. The first curve 1 and the second curve 2 respectively illustrate the probability densities

and evoked

previously.

In the TCSPC method, a sample of molecules is illuminated with a laser beam, for example. The molecules excited by photons from the laser in turn emit photons which are detected by appropriate means, including using suitable filters. The arrival time of the photons is measured and recorded. This gives a distribution of the response times characterizing the molecules. The times being of the order of several nanoseconds for a required accuracy of the order of one picosecond.

In reality, there is not a single photon. There are actually two hazards. The first hazard concerns the number of photons re-emitted by the molecules and the second random relates to the response times, because each photon has its associated response time which is specific and therefore indeterminate. The stronger the laser emits, the more photons re-emitted by the molecules. The times being very short, the chronometer used can measure only one duration, the duration corresponding to the time of the first photon received. The time the timer resets to engage a second is then too long compared to the observed physical phenomenon. The fluorescence is indeed over and all the photons already received. At each laser pulse, we have a number of photons arriving but only the response time of the first is measured. The distribution obtained is therefore constructed only for the first photon received, and therefore for the shortest response time, thus leading to a compressed image and finally to a measurement distortion. This is the stacking effect, where only the first photon is detected and measured. The density function associated with the measurements is illustrated by the first curve 1. This curve 1 has a high concentration of measurements corresponding to the very short response times.

A known solution is to lower the power of the emitted laser beam so that there is only one photon reemitted at each laser pulse. The distribution obtained is satisfactory because it is devoid of distortion. The second curve 2, of exponential form, illustrates the density associated with this solution. The difference of the first curve 1 with respect to this second curve 2 illustrates the distortion caused by the "stack-up" effect. The advantage of the measurements corresponding to curve 2 is that we know how to extract parameters θo sought with statistical methods, since it is an exponential density or in more complex cases of an exponential mixture . However it follows that in the vast majority of cases, in almost 95% pulses, there are no re-emitted photons. Only 2% to 5% of the emitted laser pulses give rise to the re-emission of a photon. Obtaining a distribution according to the second curve 2 is therefore with a significant loss of time. There are applications where the distribution must be obtained quickly and where the loss of time of the method illustrated by this second curve 2 is not admissible. In particular, in cases where the observed system is unstable and the longer the measurement, the more possibilities there is for the state of the system to evolve making the set of measurements unsuitable for characterizing the system. There are also applications that require a high flow of measurements, especially in the medical field.

Figures 2a, 2b and 2c illustrate the measurement histograms 21, 22, 23 with the effect of "stack-up" for different intensities. As indicated previously, in TCSPC the adjustment of the laser intensity, or the attenuation of the optical chain, determines the average number of photons detected by excitation. The effect of different intensity choices on the distribution of observations is illustrated by the three histograms of Figures 2a, 2b and 2c corresponding to 1000 observations in each case. The target density J ^f A> is an exponential of average 1 and each histogram is composed of 1000 observations, that is to say 1000 laser pulses, obtained by Monte-Carlo simulation. The last bin, noted + °°, represented on the abscissa scale, contains the number of observations where no photon is detected. The intensity of the laser is expressed by the average number λ of photons detected by excitation. More precisely, λ is the parameter of the Poisson's law which is the distribution of the number of photons detected by excitation. When the intensity is low, λ = 0.02, there are only 13 photons detected in 1000 excitations, as illustrated in Figure 2a. When the intensity is very high, λ = 20, the distribution is almost degenerate because almost all observations are in the first bin, as illustrated in Figure 2b. In these two extreme cases, it is clear that the histograms contain very little information on the exponential parameter, the distribution of the response times, to be estimated. It seems favorable to use an average intensity. Figure 2c shows that λ = 2 leads to a "usual" histogram because the observations where a photon has been detected are very numerous and the arrival times are quite dispersed.

An object of the invention is to obtain a distribution of the response times according to the first curve 1, without loss of time, that is to say obtained much faster than the distribution according to the second curve 2. In particular , the power of the emitted laser beams is not reduced as much as for the second curve 2. This is in particular the resolution of the inverse problem mentioned above, leading to obtaining the parameters of the distribution according to the second curve 2. In the case of a simple exponential, for example, the parameter obtained is the parameter τ which characterizes the exponential density

le vecteur des paramètres qui caractérise la famille des distributions à laquelle appartient la densité cible

. Puis dans une deuxième étape 32, on maximise cette log-vraisemblance approchée en θ. Rechercher les paramètres θo de la fonction î

revient en fait à rechercher le vecteur θ des paramètres où la log-vraisemblance atteint son maximum. La maximisation est faite sur tout l'ensemble Θ des paramètres θ qui détermine la famille des distributions à laquelle appartient la densité cible

For this purpose, the invention comprises at least two steps as illustrated by FIG. 3. In a first step 31, a simple function approximating the log-likelihood of the target distribution, which is the desired function, illustrated by FIG. second curve 2 of Figure 1. The log-likelihood, depending on the response time measurements, is a function in

the parameter vector that characterizes the family of distributions to which the target density belongs

. Then, in a second step 32, this approximate log-likelihood is maximized at θ. Find the parameters θo of the function î

come back in search for the vector θ of the parameters where the log-likelihood reaches its maximum. The maximization is done on the whole set Θ of the parameters θ which determines the family of distributions to which the target density belongs

de la densité fβ afin de corriger la distorsion, W^_n représentant le poids associé à la donnée Z_(i). Dans cette phase de repondération, on donne un poids aux données Z_(i) 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

Thus, in the first step 31, we construct an estimator of the log-likelihood function E [log (fβ (Y))] based on data according to the density,, where E is the expected expectation function and Y is a variable. random density of density h, that is E [log (fθ (Y))] = J log (fθ (y)) fβ (y) dy. Let Zi, ..., Z _n photons measures response times. The same measures are ordered in ascending order by Z ₍₁₎ , ..., Z _(n ). That is, Zf ₁ ) denotes the smallest measured response time, Z ₍₂ ) the second smallest, and so on. The estimation of the log-likelihood is obtained by a reweighting of the log-likelihood

the density fβ in order to correct the distortion, W ^ _n representing the weight associated with the datum Z _(i) . In this reweighting phase, we give a weight to the data Z _(i) which depends on their temporal positions. In other words, the data Z ^ are ordered according to their value and their weight depends on their order, that is to say that the weight of the detected photons Z, depends on their order of arrival. In particular, the highest arrival time Z _(n) is associated with the greatest weight in the function

In the second step 32, the log-likelihood thus corrected is maximized by using a suitable algorithm, such as the EM algorithm, abbreviation for "Expectation-Maximization", presented in the document by APDempster et al, "Maximum likelihood of incomplete data". via the EM algorithm ", Journal of the Royal Statistical Society, Series B, 39 (1), 1-38, 1977.

An example of an estimation method is then presented for a "stack-up" application in the simple case of a mixture of exponential densities. Then we present an estimation method for an application with a device function.

In a statistical approach with "stack up" for exponential mixtures, it is assumed that the distribution of life durations in Fluorescence represents a mixture of exponential laws. We consider the density family {f _θ , θeΘ} given by the following relation (1):

et

\_e vecteur de paramètres et Θ l'ensemble de paramètres θ. Les constantes de durées de vie en fluorescence sont données par les paramètres exponentiels

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

la fonction de répartition associée à la densité donnée par

where>κ> 0 are the exponential parameters and the weights satisfy

and

\ _E parameter vector Θ and all parameters θ. The constants of fluorescence lifetimes are given by the exponential parameters

In addition, K denotes the order of the exponential mixture which corresponds to the number of different molecules. We assume that K is known and fixed. We notice

the distribution function associated with the density given by

est le paramètre inconnu à estimer. Les

représentent les temps d'arrivée des photons de fluorescence qui frappent le détecteur.We denote by Y1 Y2, ... a series of random variables independent of density where

is the unknown parameter to estimate. The

represent the arrival times of the fluorescence photons striking the detector.

. Les probabilités sont données par la relation suivante :The number N of photons striking the detector by laser pulse is modeled by a Poisson distribution of parameter λ (restricted to

. The probabilities are given by the following relation:

On suppose que N est indépendant de la suite Y1 Y2, .... On définit alors l'observation avec « pile up » par

la mesure retenue Z correspond, comme indiqué précédemment, au temps de réponse le plus court. La fonction de répartition avec « pile-up » '"^:! de la variable aléatoire Z est donnée par la relation suivante :

We assume that N is independent of the sequence Y1 Y2, .... We then define the observation with "stack up" by

the retained measurement Z corresponds, as indicated above, to the shortest response time. The distribution function with "stack-up"'" ^:! Of the random variable Z is given by the following relation:

et la densité avec « pile up » -^ par la relation suivante

and the density with "stack up" - ^ by the following relation

par un estimateur basé sur des observations avec « pile up ». Ensuite, on obtient un estimateur du paramètre θo 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.To carry out the estimation, the method according to the invention approaches the log-likelihood of

by an estimator based on observations with "stack up". Then, we obtain an estimator of the parameter θo by maximizing this approximate log-likelihood. The simple form of this log-likelihood allows the application of the EM algorithm, defined later, to determine its maximum.

est donnée par la relation suivante :The log-likelihood of fe is then approximated. The log-likelihood of the function f _θ associated with a sample Yi, ..., Y _m of independent and identically distributed random variables (iid) of density

is given by the following relation:

et l'estimateur du maximum de vraisemblance classique Θ_ML est défini par la relation suivante :

and the classical maximum likelihood estimator Θ _ML is defined by the following relation:

II découle de la loi forte des grands nombres que L₀(θ) est un estimateur consistent de E[log(fe(Y))]. Si Z est une variable aléatoire de densité -^ * , on établit la propriété suivante en utilisant relation (3),

It follows from the strong law of large numbers that L ₀ (θ) is a consistent estimator of E [log (fe (Y))]. If Z is a random variable of density - ^ *, we establish the following property using relation (3),

On en déduit que la quantité

est un estimateur consistant de E[log(fo(Y))] basé sur un échantillon Z₁, ... ,Z_n i.i.d. de densité

Pourtant, la quantité L₁(O; θ₀) ne peut pas être utilisée pour estimer le paramètre inconnu θo, car elle-même dépend de θo. Or, une autre modification est indispensable. On déduit de la relation (2) que :

We deduce that the quantity

is a consistent estimator of E [log (fo (Y))] based on a sample Z ₁ , ..., Z _n iid of density

However, the quantity L ₁ (O; θ ₀ ) can not be used to estimate the unknown parameter θo, since it itself depends on θo. However, another modification is essential. We deduce from relation (2) that:

Λ; i, par la fonction de répartition empirique

on établit un nouvel estimateur deReplacing

Λ; i, by the empirical distribution function

a new estimator of

E [log (f _θ (Y))] similar to L i (θ; θo) by

où Z(ι) désigne la /-ème statistique d'ordre de Z₁, ... ,Z_n. On introduit les poids

On définit un nouvel estimateur de θo par la relation suivante

where Z (ι) denotes the / -st statistical order of Z ₁ , ..., Z _n . We introduce the weights

We define a new estimator of θo by the following relation

Donc, on cherche le paramètre où la log-vraisemblance approchée atteint son maximum. On peut considérer comme un estimateur de maximum de vraisemblance approché.

So, we look for the parameter where the approximate log-likelihood reaches its maximum. We can consider as an approximate maximum likelihood estimator.

Cela changera les poids qui sont associés

aux données mais le principe restera le même.

Note that other estimators of the distribution function can also be used

This will change the weights that are associated

to the data but the principle will remain the same.

Subsequently, a maximization is performed by the EM algorithm.

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).We find that the problem of maximizing the relation (5) differs from the problem of the relation (4) only by positive constants w

However, the maximization of the relation (4) can be done by the EM algorithm. Due to the similar structure of the problems, the EM algorithm also applies to relation (5).

We introduce the missing variable S from an observation Y of the exponential mixture fe being the label of the component of the mixture that generated the observation Y. More precisely, S is a random variable with values in the attached law of (Y ₁ S) is given by

As in the classical EM algorithm, we define the function

Ce problème de maximisation est facile car il se décompose en plusieurs problèmes de maximisation indépendants:and iteratively calculates the sequence (θ ^(t) ) _t by maximizing Q (θ, (θ ^(t) ), leading to the following relation:

This problem of maximization is easy because it breaks down into several independent maximization problems:

où la dernière maximisation est faite sous la contrainte que tous les ctk sont à valeurs dans ]0,1[ et ∑kCik = 1. On peut montrer qu'avec chaque itération la log-vraisemblance approchée L(θ) augmente. Plus précisément on a pour toute valeur initiale θ^<1) e θ et pour tout t ≥1,

where the last maximization is made under the constraint that all ctk are with values in] 0,1 [and ΣkCik = 1. It can be shown that with each iteration the approximate log-likelihood L (θ) increases. More precisely, for any initial value θ ^<1) e θ and for all t ≥1,

critique de L(θ).Moreover, we can show that the sequence converges towards a point

critical of L (θ).

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(θ^(t)) à toute itération et on s'arrête si l'accroissement par rapport à l'étape précédente est inférieur à un certain seuil. Un autre critère étudie la variation des valeurs du paramètre θ^m. A toute itération on évalue la différence

désigne la norme infinie. Si la différence est inférieure à un certain seuil, on considère la dernière valeur de θ^m comme limite de la suite

There are several stop criteria to determine when the next

has converged and where we can stop the iterative calculation. One criterion is based on the evolution of the approximate log-likelihood. The value of L (θ ^(t) ) is calculated at any iteration and stopped if the increase from the previous step is below a certain threshold. Another criterion studies the variation of the values of the parameter θ ^m . At any iteration we evaluate the difference

designates the infinite norm. If the difference is less than a certain threshold, consider the last value of θ ^m as the limit of the following

The EM algorithm to maximize the L (O) function. in the case of an exponential mixture, can be defined by the following phases:

We set a threshold σ> 0 for the stopping criterion.

i.i.d. de densité

On note Z

les observations ordonnées, n étant le nombre total de mesures avec photons détectés, et / étant l'ordre d'arrivée du photon

1. On pose t = 1.Let λ be the known value of the Poisson parameter during the experiment. Let _{n be} the observations with "pile-up", that is to say a sample

density iid

We note Z

the ordered observations, n being the total number of measurements with photons detected, and / being the order of arrival of the photon

1. We put t = 1.

On calcule, pour / = 1, ...,n

We choose an initial value

We calculate, for / = 1, ..., n

2. The following steps are repeated until convergence:

(a) Calculate the conditional probabilities, for i = 1, ..., n and s = 1, ..., K,

(b) Determine the new parameters, for

(c) D'après le critère d'arrêt choisi, on évalueWe pose

(c) According to the chosen stopping criterion,

(d) If σ _{t +} i <σ, we go to phase 3 below which gives the estimator,

(e) if not, set t = t + 1 and go back to step 2 (a). 3. The estimator of the parameter θo is given by

The previous model is only close to reality for measuring devices with a laser that emits light instantly, that is to say if the laser pulse can be modeled by a Dirac law with all the mass in zero. However, the light source of most measuring devices emits light for a short time. We call function of device the light remission distribution which is different for any device, but which can be easily registered before performing the experiment. So a TCSPC observation is the time that the molecule remains in excited state, the lifetime, to which is added a random duration that comes from the device. The EM method takes into account the hazards caused by the device.

donnée par la relation (1 ).Η is the probability density of the device function. It will be assumed in the following that η is known. Let (E ₁ , i≥1} and (Y ₁ , i≥1) be two independent sequences of random variables independent of law η respectively of the mmlaellgegege eexxppoonneennttiieell

given by relation (1).

We pose

The density of is the convoi ution of

par ' Or, une observation en

TCSPC a la même distribution que : .

We note the distribution function of

by 'Or, an observation in

TCSPC has the same distribution as:.

et la densité avec « pile-up » 5Kest donnée par

De la même manière que précédemment, on déduit l'estimateur comme le vecteur θ qui maximise la log-vraisemblance approchée,

où

désigne des observations avec « pile-up » de densité et

la log-vraisemblance approchée est donnée parThe distribution function with "stack-up" is given by:

and the density with "stack-up" 5K is given by

In the same way as above, we deduce the estimator as the vector θ that maximizes the approximate log-likelihood,

or

designates observations with "stack-up" density and

the approximate log-likelihood is given by

Pour l'observation

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

est donnée par :

For observation

we consider the two missing variables (S, E) where S is the label of the component of the mixture that generated Y, and E is the duration added to the life of the device. The attached law of

is given by:

dans le cas d'une fonction d'appareil est défini par les phases suivantes :The EM algorithm to maximize

in the case of a device function is defined by the following phases:

We set a threshold σ> 0 for the stopping criterion.

Let λ be the known value of the Poisson parameter during the experiment. Let η be a piecewise constant known density given by

densité .9^ . On note ^ les observations ordonnées

So observations with stack up, a sample iid of

density .9 ^. We note the ordered observations

On calcule, pour i=1, ...,n ,1. We set t = 1 We choose an initial value ^v

We calculate, for i = 1, ..., n,

2. The following steps are repeated until convergence: (a) We calculate, for i = 1, ..., n and s = 1, ..., K,

où aΛb dénote le minimum de a et b, et on calcule

et

where aΛb denotes the minimum of a and b, and we calculate

and

et(b) The new parameters are determined, for s = 1, ..., K,

and

(c) D'après le critère d'arrêt choisi, on évalueWe put ^t

(c) According to the chosen stopping criterion,

passer à la phase (3) ci-dessous.(d) If

proceed to phase (3) below.

(e) In the opposite case we set t = t + 1 and return to step 2 (a).

3. The estimator of the parameter θ ₀ is given by

une famille de densités quelconque telle que l'algorithme EM s'applique à un échantillon i.i.d. de densité

Soient

une suite de variables aléatoires indépendantes de densité On note

la fonction de répartition associée à la densité

The invention uses the preceding algorithm for the general inverse problem mentioned above, that is to say for density families other than the family of exponential mixtures and for other discrete laws than the Poisson's law. This new EM algorithm, used to solve the general inverse problem, is thus not restricted to the exponential mixture and Poisson's law. We consider

any density family such that the EM algorithm applies to a density iid sample

Let

a sequence of random variables independent of density We note

the distribution function associated with the density

On note m la fonction génératrice de A/^*,

et m' la dérivée de m. On définit l'observation avec « pile-up » par :

We consider N ^* a random variable with values in {1,2, ...} independent of the sequence i

We denote by m the generating function of A / ^* ,

and m 'the derivative of m. We define the observation with "pile-up" by:

_{avec (<} pjie-up » :Approximate log-likelihood based on observations is established

_{with (<} pjie-up ":

OÙ

OR

où S dénote les variables manquante associées à la variable

et désigne la densité jointe de ^ - ^ k

We define

where S denotes the missing variables associated with the variable

and denotes the attached density of ^ - ^ k

suivantes :The general EM algorithm for maximizing involves phases

following:

2. On réitère jusqu'à la convergence les étapes suivantes : (a) On poser t = t+1. (b) On calcule

We set a threshold σ> 0 for the stopping criterion. 1. We choose an initialization value We set t = 1.

2. The following steps are repeated until convergence: (a) We put t = t + 1. (b) We calculate

(c) According to the chosen stopping criterion,

(d) If σt <σ we go to phase (3) below, otherwise we go back to step 2 (a).

3. The estimator of the parameter θo is given by

An estimation method according to the invention can provide very accurate measurements when applied to real data. It requires little computing time, takes into account the device function, is easy to implement and can handle the case of a mixture of several exponentials.

In comparison with existing methods, it is particularly adapted to the nonlinear distortion of the general problem. In particular, the results reported in Rebafka et al., 2008, can be applied to the optimal setting of laser intensity. In other words, the experiment can be carried out under optimal conditions and thus a reduction of the acquisition time, by a factor of 10 for example, can be achieved compared to the techniques currently used.

The system to be characterized may be composed of molecules as described above. The invention can nevertheless be applied to other systems. A system to be characterized may be, for example, an optical fiber in which photons are injected and whose lifetimes are measured. Other types of particles than photons can also be used.

Claims

A method for estimating the parameters θ 0 of the distribution of the response times of particles reissued by a system, in response to a series of particle emissions to said system, a response time of a corresponding particle being measured between a instant of emission and the instant of detection of said particle, characterized in that the method comprises at least: a step (31) of determining a function approaching the log-likelihood of the target distribution <M> corresponding to a distribution without stacking effect from measurements of response time of particles detected with stacking effect, that is, measures of the distribution ^ o; a step (32) for maximizing the approximate log-likelihood on the set of parameters Θ in order to determine the estimator as the vector θ that maximizes the approximate log-likelihood. 2. Method according to claim 1, characterized in that the approximate function is an estimate of the log-likelihood of the target distribution function <M> based on observed measurements distributed according to the distribution function ^ with stacking effect. 3. Method according to claim 2, characterized in that the approximate function is obtained by a weighting of the log-likelihood of the target distribution function J **. 4. Method according to claim 3, characterized in that the weighting of the log-likelihood is done so as to order the observed response time measurements Zi, ..., Z _n according to the order of arrival of the corresponding particles, the weight assigned to each measurement depending on its order. 5. Method according to claim 4, characterized in that the weight assigned to a measurement increases with its associated order of arrival. 6. Method according to any one of the preceding claims, characterized in that n being the total number of particles detected and / being the order of arrival of a particle, the weight w ^* _hn associated with the measurement of the time of Z _(l) response of a particle of order / is given by the following relation:

where m denotes the generating function of N,

and m 'the derivative of m and is an estimator of the distribution function.

7. Method according to any one of the preceding claims, characterized in that the target distribution function is considered

as a mixture of exponential densities. 8. Method according to any one of the preceding claims, characterized in that the target distribution function takes into account the measurement apparatus function, the target function being the product of convolution of the target function without device function.

with the probability density of the device function η. 9. Method according to any one of the preceding claims, characterized in that the particles are photons emitted by a laser, the system reemitting photons. 10. Method according to any one of the preceding claims, characterized in that the system is a set of molecules to be characterized by the target distribution function of the response times.