Bayesian analysis of biexponential time-decaying signals

Spectrochimica Acta Part A 54 (1998) 559 – 566

Bayesian analysis of biexponential time-decaying signals S.L. Whittenburg Department of Chemistry, Uni6ersity of New Orleans, New Orleans, LA 70148, USA Received 26 April 1996; received in revised form 1 October 1997; accepted 1 October 1997

Abstract In this paper we demonstrate the use of Bayesian analysis methods for the extraction of information from time-domain biexponentially decaying functions. In several forms of time-domain spectroscopy the physically interesting spectrum consists of multiple exponentially decaying functions. Here we demonstrate the use of Bayesian analysis to extract information on both exponentially decaying signals from a biexponentially decaying signal in the presence of noise. Other spectral factors will be considered. © 1998 Elsevier Science B.V. All rights reserved. Keywords: Bayesian analysis; Biexponential decay

1. Introduction In many forms of molecular spectroscopy the physically interesting information under study is extracted from the analysis of an exponentially decaying experimental spectrum [1]. Typically the desired information is extracted from either the amplitude or decay rate of the experimental data which is obtained by conventional curve fitting of the experimental data to a model function consisting of a decaying exponential. Often, however, the experimental data consists of multiple-exponential decays. Removal of DC components or baseline contributions is usually routine. However, extraction of information from biexponential decays proves to be more burdensome, particularly in the presence of noise. In this paper we present the use of Bayesian analysis methods for the extraction of the physical information from a biexponentially decaying time function.

2. Theory Bayes’ theorem is P(A B, I)=

P(A I)× P(B A, I) P(B I)

(1)

where the proposition that A is true given B and any prior information, I, is related to the proposition that A is true given only the prior information, the proposition B is true given A and the prior information divided by the probability B is true given the prior information. To relate this more directly to spectroscopy, consider a general time-domain function, f(t), define as m

f(t)= % Bj Gj (t, {t})

(2)

j=1

where Gj is some function of time which depends parametrically on a set of decay times, t. In this example it is a set of decaying exponential func-

1386-1425/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PII S 1 3 8 6 - 1 4 2 5 ( 9 7 ) 0 0 2 5 6 - 4

S.L. Whittenburg / Spectrochimica Acta Part A 54 (1998) 559–566

560

tions each of which is composed of several relaxation times. Each of these functions has an amplitude, Bj. f(t) is thus a set of model functions which may describe the experimental data, D. Bayes theorem can then be used to determine the probability that the model functions, which may contain adjustable parameters, correctly describe the experimental data P(B, {t} D, I)=

P(B, {t} I) ×P(D B, {t}, I) P(D I) (3)

One advantage of the Bayesian formulation is that linear parameters in the probability may be removed by marginalization, that is, the parameter-dependence may be integrated out of the probability. For exponentially decaying model functions the amplitude may be marginalized so that only the decay rate of each contribution remains. In doing marginalization we have removed the correlated parameters, the amplitude and decay rate, so that non-linear least-squares adjustment of the parameters to optimize the probability now becomes easier. It has previously been shown that marginalization of the amplitude leads to the following expression for the probability [2]

!

P(t D) 1−

"

mh( 2 Nd( 2

function is adjusted. The marginalized parameters, in this case the amplitude of each resonance, may be obtained from the final finals of the estimated parameters. We have previously shown that the Bayesian method provides reliable estimates of frequencies, decay rates and amplitudes of peaks within complicated NMR spectra [3]. In the above formulation of Bayesian probability analysis we can test model functions against a data set. In this example we can test a single exponential decay against a biexponential decay and compute the global likelihood that the model ‘fits’ the data. The model function with the higher global likelihood (related to the sufficient statistic) is more probable in describing the data. Once a biexponential decay model is justified, the physical information for each component of decay may be extracted from the signal. Note that the contribution of each function to the data is not an adjustable parameter in this scheme. Only the decay rates (in this example) are adjusted. The amplitudes of each model function have been marginalized out of the probability and are, therefore, evaluated at the end of the calculation.

m − N/m

(4)

where m is the number of model functions, N is the number of data points, and t is the decay rate, h( 2 is the sufficient statistic, defined by h( 2 =

1 m % h m j = 1 ij

(5)

where hij is the projection of the orthogonalized model function onto the data and d( 2is given by d( 2 =

1 N 2 % di N i=1

(6)

where the orthogonalized model functions are obtained by forming the matrix, gij, of all possible products Gi  Gj of the model functions given in Eq. (2) and obtaining the eigenvectors of this matrix. The Bayesian method is a non-linear leastsquares optimization of the probability, given by Eq. (4), where the decay rate of each model

Fig. 1. Difference between the sufficient statistics for one vs two exponential decays as a function of relative amplitudes of the two components (F) and added noise. Given contours are positive except the uppermost contour which is the zero contour.

S.L. Whittenburg / Spectrochimica Acta Part A 54 (1998) 559–566

561

Fig. 2. (a) Estimated fraction of the larger component for the biexponential model function as a function of the separation between the two decay rates. The true value is F =0.7. The solid line denotes the mean estimate for the decay time and the dashed line denotes one standard deviation. The added peak-to-peak noise is 0.05. (b) The same as (a) except the predicted value is the decay rate of the larger component. The true value is 0.1. (c) The same as (a) except the predicted value is the decay rate of the smaller component. The true value is given by the dashed line.

562

S.L. Whittenburg / Spectrochimica Acta Part A 54 (1998) 559–566

Fig. 2. (Continued)

3. Results and discussion

3.1. Single 6s biexponential decay To test and to display the potential of the method we have generated a synthetic data set similar to one that would be acquired in conventional molecular spectroscopy. The data set is generated in the time-domain. The data is generated with either one or two time-domain contributions, an exponential decay. Data sets were generated as a function of the relative amplitude of the two exponential terms and as a function of added noise. If the fraction, F = 1.0, then the data set is a single exponential decay. F =0.5 means two equal amplitude exponential functions. In each data set 256 points were calculated with a time step of 0.2 and decay rates of 0.1 and 1.0. The total amplitude of the time-domain data was 1.0. To each data set random noise was added with an amplitude (peak-to-peak) in the range 0.05 –0.5. We are interested in the ability of the method to determine if two exponential decays are present in the data set. To this end 100 data

sets were randomly generated with differing random ‘seeds’ at ten discreet values of the signal-tonoise ratio between 0.05 and 0.5. These sets were then Bayesian analyzed as described above to provide 100 estimates of the sufficient statistic with a single-exponential (h 21) and a biexponential (h 22) model function. In the Bayesian method as applied here only the decay rate is an adjustable parameter. Initial estimates of this parameter were randomly chosen for each run in each set in the range 0.05–0.5 as most experimentalists would be able to estimate the decay rate within this range. Shown in Fig. 1 is the difference, h 22 − h 21 vs F and noise-level. Shown in Fig. 1 are contours for which this difference is positive, that is for which the method predicts there are two exponentials in the data set. The last contour at the top of the figure is the zero contour for which there is equal probability of one or two exponentials. Note that much of this figure is as you would expect. As the fraction of the second exponential approaches 0.5 and as the noise level approaches zero the method more easily distinguishes the second exponential contribution producing a large positive peak in

S.L. Whittenburg / Spectrochimica Acta Part A 54 (1998) 559–566

the lower left quadrant of the plot. As the noise level increases and the amplitude of the second (smaller) component decreases (F approaching 1.0) the method has a more difficult time distinguishing the second exponentials. Still, even if the noise level is 30% of the total signal and the second component is at least 10% of the total amplitude, the method correctly predicts a biexponential model is the correct model. Even out to a noise level of 50% if the second exponential is at least 20% of the total amplitude the method can distinguish the second exponential. At this point the noise is approximately 2.5 times the amplitude of the second exponential. Thus, the method is capable of reliably distinguishing single vs biexponential decay.

3.2. Analysis of biexponential decays Since we have demonstrated the utility of the method for predicting biexponential decays we next wish to examine the power of the method in extracting information on the two exponential contributions to a biexponential decay. The information we wish to extract is the relative amplitudes of each contribution (the method provides an estimate of each amplitude separately) and the decay rate of each component. For each of these pieces of information we wish to examine the utility of the method as a function of the added noise and as a function of the relative decay rates of the two components. Clearly, as the noise is increased and as the two decay rates approach the same value, the method should begin to fail. To this end we have generated 100 data sets at F = 0.7 with one decay rate set to 0.1 and the second (smaller component) chosen at ten values between 0.1 and 1.0. The noise level is fixed at 0.05. Fig. 2a shows the estimated F value along with one standard deviation. Down to a separation in decay rates of 0.2 the method reliably estimates the fraction of the two components. Below this separation the method begins to assign all of the probability to a single exponential function. Fig. 2b shows the estimate of the decay rate of the larger (70%) component. Again, down to about 0.2 separation in decay rates the estimate is quite reliable with a very small standard deviation. Fig. 3a shows the estimated decay rate of the smaller component

563

along with the true value (dashed line). Again, above a separation of 0.2 the estimate of the second decay rate is quite reliable. To examine the effect of increasing noise a second data set was generated as described above except now the noise level is fixed at 0.15 or 15% of the total signal. This is near the limit of noisy data set where one would attempt to extract information on biexponential decay. Still, in Fig. 3a and 3b the method reliably estimates the fraction F and the decay rate of the larger component. The difference is now that the standard deviations are significantly larger. To the experimenter this means that estimation of these parameters from a single measurement or a small sample of measurements could lead to a sizeable error is estimating these parameters Fig. 3c shows the estimates of the decay rate of the smaller component along with the true value (dashed line). Here the estimates are always high (predicting a larger decay rate than the true value) and the deviation is sizeable. It is unlikely that on a single measurement reliable estimates of the decay rate will be obtained. We look at this fact with a single data set analysis below.

3.3. Analysis of a single data set To demonstrate the method on a single, noisy data set we have generated a biexponential set with F= 0.7 at the two decay rates of 0.1 and 0.5, with the larger component having the smaller decay rate. The data set was computed as described above with added peak-to-peak noise of 0.15. The results of the Bayesian analysis on this single data set was F=0.71 and k1 = O.10 and k2 = 0.55. As has been noted above the results for this single set are in good agreement with the known values except that the decay rate of the smaller component is slightly overestimated. The results of the procedure are shown in Fig. 4, where the curves show the relative contributions to the data and the sum of the two components.

3.4. Analysis of experimental data As a final test of the method we have analyzed the transient fluorescence spectrum of pyrene excited at 355 nm. 1024 Data points were obtained

564

S.L. Whittenburg / Spectrochimica Acta Part A 54 (1998) 559–566

Fig. 3. (a) The same as Fig. 2a except the added noise is now 0.15. (b) The same as Fig. 2b except added noise is now 0.15. (c) The same as Fig. 2c except added noise is now 0.15.

S.L. Whittenburg / Spectrochimica Acta Part A 54 (1998) 559–566

565

Fig. 3. (Continued)

Fig. 4. Bayesian predicted decay for a single data set with F =0.7, k1 =O.1 and k2 =0.5. The added noise is 0.15. The lower dotted curve shows the contribution of the smaller component while the lower dashed curve shows the contribution of the larger component. The upper dashed curve shows the resulting ‘fit’ to the data.

S.L. Whittenburg / Spectrochimica Acta Part A 54 (1998) 559–566

566

Fig. 5. Bayesian analysis of the transient fluorescence spectrum of pyrene. The dark solid line is the Bayesian result.

with a 5 us window. Conventional non-linear least-squares fitting of the data confirmed two exponential decays with decay rates of 3.67× 106 and 2.01 × 105. The Bayesian method gave identical results. The original data are plotted in Fig. 5. The dark solid line is the Bayesian result.

Due to the generality of the method it may be applied to other cases in which either the experimental time-domain data is either non-exponential. It is also rather straightforward to extend the method to spectroscopic techniques in which the data is not acquired in the time-domain. Such applications await future research. The programs used in this research may be obtained by writing to the author.

4. Conclusions We have demonstrated that the Bayesian probability analysis method may be successfully used to estimate decay rates and amplitudes from biexponentially decaying data in the presence of significant noise and to predict whether a given data set contains single or biexponential decay. We have demonstrated the method using means and standard deviations of a relatively large number of random data sets in order that an estimate of the reliability of a result on a single data set might be predicted. It has been shown that under these circumstances the method successfully extracts reliable estimates in noisy data sets.

Acknowledgements This work is supported by NSF and the Louisiana Educational Quality Support Fund through grant EHR 9108765.

References [1] C.H. Wang, Spectroscopy of Condensed Media, Academic Press, Orlando, FL, 1985. [2] G.L. Bretthorst, Bayesian Spectrum Analysis and Parameter Estimation, Vol. 48, Springer-Verlag, New York, 1988. [3] R.F. Evilia, R. Effiong, S.L. Whittenburg, Spectrosc. Lett. 26 (1993) 1559.