Automatic peak identification method

[0001] The present application claims priority from U.S. Provisional Application Ser. No. 60/664,010, filed Mar. 22, 2005.

[0002] This invention was made with government support under Grant No. 5R44 CA101479-03 awarded by NIH. The government has certain rights in the invention.

[0003] The invention is related to the field of analyzing data which includes multiple peak data and, in particular, pertains to an automatic method for analyzing data which includes multiple peak data resulting in high resolution estimates.

[0004] Work in chromatography, spectroscopy, engineering, pharmacology, and other fields frequently requires analysis of data sets exhibiting multiple peaks. Analysis of multi-peak data is particularly difficult when peaks overlap, or when data are noisy. Peak finding algorithms of many types are used to identify and interpret peaks and peak positions.

[0005] Present state-of-the-art peak finders are either based upon 1] simple matched filters, or wavelets, or searches for local maxima in the time series (resulting in less-than-optimal resolution), 2] parameter fitting to peaks using optimization methods, sometimes on the entire time series (such methods are slow), or 3] direct methods which involve significant operator interaction (i.e. by placing cursors over peaks on a computer screen, estimating baselines by eye, etc.). All of these methods have disadvantages. The method described herein runs essentially like a filter, results in much higher resolution estimates of peak positions and amplitudes than the existing state-of-the-art, and is almost completely automatic.

[0006] The invention described herein details a protocol to improve analysis and peak identification in spectroscopic data. Bayesian methods are used to automatically identify peaks in data sets. After identifying peak shapes, the method tests the hypothesis that a given number of peaks is found within any given data window. If a peak is identified within a given window, then the likelihood function is maximized in order to estimate peak position and amplitude.

[0007] This process yields a spectrum with high resolution and minimal artifacts. Furthermore, the process is not overly computational intensive, with low delay times required. A method that simply maximizes the likelihood function is likely to generate artifacts, while a method that only tests the hypothesis that a peak is present, without the subsequent maximization of the likelihood function, will yield a spectrum with decreased resolution.

[0008] The peak identification method described herein is potentially useful in any data set wherein peaks of interest are distinguished. The method described herein is useful for identifying peaks in data sets obtained from spectroscopy, including but not limited to various forms of absorption and emission spectroscopy, including absorption spectroscopy, UV/VIS spectroscopy, IR spectroscopy, X-ray spectroscopy, X-ray crystallography, Nuclear Magnetic Resonance spectroscopy (NMR), and raman spectroscopy. In particular, the methods described herein are useful for applications in mass spectroscopy, including TOF-SIMS and SELDI-MS. For example, the methods described herein are useful for analyzing TOF-SIMS mass spectroscopic data related to proteomics. The methods described herein are not limited to spectroscopic applications; for example, the methods are particularly useful when used for imaging applications.

[0012] The following symbols are accorded the following meanings:

• [0013] N total number of data points in the window
• [0014] s_j=s(t_j) j=1 . . . N observed signal prior to any processing
• [0015] x_j=x(t_j) ideal wavelet shape
• [0016] η_j=η(t_j) noise
• [0017] σ_η²=ηη^T noise co-variance

[0018] Herein we describe the use of Bayesian methods to automatically identify peaks in data sets. We start with a brief overview of the logic. The first step is to compare the hypotheses that there are a given number of peaks (i.e., 0, 1, 2, . . . n) in the observation window. If a peak is present, we then estimate its position and amplitude via parameter fitting.

[0019] After this summary of the logic, we develop algebras for two simple cases where Gaussian or Poisson noise is involved. Both of these cases are simple enough that all results can be arrived at in closed form. We first present how to derive a matched filter for a known peak shape that is embedded in additive white (Gaussian) noise of unknown amplitude. We then show the similar procedure can be applied to Poisson noised data. Some preliminary results on SELDI-MS (Surface Enhanced Laser Desorption Ionization Mass Spectrometer) and TOF-SIMS (Time of Flight Secondary Ionization Mass Spectrometer) data are shown. It is easy to add color (finite time correlations) to the noise. Once these simple situations are understood, extending to other types of noise, such combinations of Poisson, Gaussian, and coherent ‘noise’ processes, is logically straightforward, though algebraically more involved. All of the methods described herein require knowledge of the ‘target’ peak shape.

[0020] FIG. 1 shows a sample mass spectrum, in which the number of ions is was plotted versus the time that ions take to fly to detector. Ideally, ions of a specific m/z would exhibit the same time of flight, resulting in sharp peak line shape like a delta function. However, in reality, they enter the drift region with velocity and time distributions of finite width, which results in a finite peak width for ions of a specific m/z.

[0021] 1.1 Description of the Underlying Assumptions

[0022] 1.1.1 The Relationship Between Peak Signal and Probability Density

[0023] In counting experiments, ideally only one ion of a given mass arrives, or at most a few ions of the same mass arrive, at the detector during a given temporal sampling interval and are counted individually. The measurement is repeated many times under nominally identical conditions and one gradually builds up a portrait of the probability distribution of arrival times for each ion mass. This means that the final time series s(t) that is to be analyzed after all the data has been gathered, is proportional to the conditional probability distribution
s(t_k)Δt∝p(t_k|ƒ₀(ν*(m),t*))Δt,  (1)
wherein s(t)Δt is the signal amplitude at time t_k, Δt is the sampling interval, and p(t_k|ƒ₀(ν*(m),t*)) is the probability that an ion of mass m will strike the detector between t_k and t_k+Δt given that it entered the drift region at time t* with a velocity close to ν*(m). If only one ion species is present, the initial velocity distribution, ƒ₀(v*(m),t*) is just a sharply peaked function of ν of the form

[0024] $\begin{array}{cc}{f}_0\left(v_{*}\left(m\right),t_{*}\right)=g\left(\frac{v-v_{*}\left(m\right)}{\sigma \left(m\right)}\right),\int d\mathrm{sg}\left(s\right)=1.& \left(2\right)\end{array}$
If a variety of different ions is present, ƒ₀ will be a superposition of terms like (2), properly normalized so that ƒ₀ will have unit total area. This assumes that the shape of the velocity spread is universal, up to translation and rescaling. We will assume ν*(m) and σ(m) to be weak functions of the mass, and to scale like 1/√{square root over (m)}. Fixing their functional form is part of the calibration process. In practice, the velocity distributions for each peak will contain information about the dynamics of the ion formation process, which depend upon the ion species. However, if we assume that the ion optics are designed so that it will only pass a narrow range of most possible velocities for any given mass, then maximizing the entropy of the probability distribution ƒ₀ for a single peak, given v* and σ, leads to the choice of a Gaussian for g:

[0025] $\begin{array}{cc}g\left(\frac{v-v_{*}}{\sigma }\right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\left(v-v_{*}\right)}^2}{2{\sigma }^2}\right)& \left(3\right)\end{array}$
The choice of the maximum-entropy distribution is founded upon the principle that it maximizes the number of possible outcomes of repeated observations that are consistent with this assignment of the probability. Hence, it is the least biasing assignment of the probability that is consistent with our limited knowledge of the initial distribution. Transforming from the initial velocity distribution g to the temporal peak shape requires solution of a simple Fokker-Planck equation, and is described later.

[0026] 1.1.2 The Assumption of Independence

[0027] The observed peaks for a portion of a spectrum show in FIG. 1 have finite width in time. One of the causes is the finite width of the velocity distribution upon entering free flight. However, it is important to realize that ions that arrive between time t_k and t_k+Δt are independent of those that arrive at any other time. This is because, as stated in the previous section, each event/measurement only results in one or at most a few ions of a given m/z. The resulting spectrum is an accumulation of many measurements under the same experiment condition. Therefore, the joint probability of arrival of two ions at two different times satisfies
p(t_k,t_j|ƒ₀)=p(t_k|ƒ₀)p(t_j|ƒ₀), j≠k.  (4)
even for counts associated with the same mass peak. Any correlations in the signal are assumed to be due to the electronics. If such correlations are present (as they are in MALDI and SELDI spectra) then this must be included in the analysis as part of the model used, as will be discussed later. The assumption of independence, however, underlies the entire analysis that we pursue below in the first few examples.

[0028] Because the measurement is repeated many times under ‘identical’ conditions, producing the same ƒ₀ at each repetition (with the ‘clock time’ t reset to zero at the start of each measurement), if a total number of n ions are counted over the entire course of the measurement, then the total number of counts between t_k and t_k+Δt is very nearly
n_k≈p(t_k|ƒ₀) n  (5)

[0029] 1.2 Model Comparisons (i.e. Counting Peaks in the Window, Comparing Noise Models)

[0030] We now consider the peak of a particular of mass m. We introduce a window that includes only N points in the time series. The width of the window is chosen to be a slowly varying function of the mass so that it typically will run only from the half-maxima to the left and right of a ‘typical’ peak of mass m. The change in window width is due to the nature of time-of-flight apparatus where heavier ions will arrive later in time and show up with wider peak shape. This is a rough-and-ready compromise between the desire to include as much data as possible in the window to improve the sampling statistics, and the realization that peaks will overlap and that our peak shape model is probably not very good out on the tails of the peak.

[0031] Focus attention now on a particular window position. We label the window location, t₀ in such a way that if a mass peak appeared in the window with zero width then the likelihood function (defined momentarily) will be maximized at t₀. Because of the mass-dependent asymmetry in peak line shape, t₀ is certainly not at the leading or trailing edge of the window. It is also not located at the center of the window, nor at the maximum of the ‘target’ temporal peak shape, but lies somewhere between these two points at a position which is weakly dependent upon m. In this way, we correct for systematic errors in the estimation of the peak shape due to its mass-dependent asymmetry.

[0032] We define the target shape of the temporal peak to be x(t−t₀) which is then sampled at discrete time intervals, giving the N-component vector
x=(x₁, . . . x_N)≡(x(t₁−t₀), . . . x(t_N−t₀))  (6)
This ‘slides’ with the window and for convenience is normalized to have unit area:

[0033] $\begin{array}{cc}\sum _{k=1}^N{x}_k=1& \left(7\right)\end{array}$
This will only introduce a constant correction for the peak amplitude computed later. We now assume that the observed signal s(t) within the window is given by
s_k=αx_k+ξ_k, k=1,2, . . . N  (8)
where a is an unknown amplitude and ξ=(ξ₁,ξ₂, . . . ξ_N) is a random process of some appropriate type. If more than one peak is present in the window, this complicates the analysis because more parameters must be fit, but the logic is similar to what we describe below. For counting experiments, we do not observe the signal s(t) directly, but instead count how many times it crosses some detection threshold in a time interval Δt. This is because a single ion striking a micro-channel plate, or a channeltron, will initiate a cascade of electron emission that amplifies the original micro-pulse into one that is macroscopically detectable. The noise process will be Poisson with some appropriately chosen rate, while for analog experiments it is probably a combination of Poisson and Gaussian processes. For the moment, we leave the choice of noise model aside and simply assert that choosing a noise model is part of the model selection process, and choosing between them is part of the hypothesis testing we now describe.

[0034] Under the assumption that the ion counts at different times are independent, the probability of observing the particular count sequence n=(n₁,n₂,K n_N) is simply

[0035] $\begin{array}{cc}\begin{array}{c}p\left(n\text{❘}a,\lambda ,t_0,M\right)=p\left(\left\{n_1,n_2,K{n}_N\right\}\text{❘}a,\lambda ,t_0,M\right)\\ =\prod _{k=1}^{N}p\left(n_k\text{❘}a,\lambda ,t_0,M\right)\end{array}& \left(9\right)\end{array}$
By the notation λ we indicate a parameter that characterizes the noise process (e.g. the variance σ for a Gaussian process or the ‘dark’ rate r₀ for a Poisson process). By M we mean to emphasize a particular choice of model class, including a peak shape x(t) and a noise model. The probability (9) is the likelihood function of observing the data in the window labeled t₀, given the model class M and the parameters.

[0036] We first want to decide whether there is a peak in the window or not. This is a comparison between two hypotheses:

• • [0037] H₀=There is no peak in the widow t₀. The data are noise of the assumed type. Let's call the associated pure-noise model M₀.
  • [0038] H₁=There is a single peak in the widow t₀ with the shape x. Deviations from this shape in the data are due to noise of the assumed type. Let's call the associated peak-plus-noise model M₁.
    We wish to compute the probability of which each of these hypotheses was true, given the data, and then take the ratio

[0039] $\begin{array}{cc}\frac{p\left(H_1\text{❘}n\right)}{p\left(H_0\text{❘}n\right)}=\frac{p\left(M_1\text{❘}n,t_0\right)}{p\left(M_0\text{❘}n,t_0\right)}& \left(10\right)\end{array}$
If this ratio is large compared to one, we can be confident that there is probably a peak in the window, while if it is approximately one we interpret that as saying there is only weak evidence of a peak in the window (because α=0 is a possible estimate of the peak amplitude, which we interpret as ‘no peak’). Therefore, one may set an appreciable threshold for peak detection. In order to derive (10) from (9), we first invoke Bayes' theorem which states:

[0040] $\begin{array}{cc}p\left(A\text{❘}B\right)p\left(B\right)=p\left(B\text{❘}A\right)p\left(A\right)\Rightarrow p\left(A\text{❘}B\right)=\frac{p\left(B\text{❘}A\right)p\left(A\right)}{p\left(B\right)}& \left(11\right)\end{array}$
Therefore, identifying A as the model class M_k, and B as the observed data in the time window t₀:

[0041] $\begin{array}{cc}p\left(M_k\text{❘}n,t_0\right)=\frac{p\left(n\text{❘}{M}_k,t_0\right)p\left(M_k\text{❘}{t}_0\right)}{p\left(n\text{❘}{t}_0\right)}& \left(12\right)\end{array}$
Since we are given that we observed the data n in the window t₀, the denominator is unity. If we have no reason to prefer one model class over another, we should assign them all equal prior probabilities. For example, if we are comparing two types of models (a single peak vs. no peak), then p(M₀|t₀)=p(M₁|t₀)=½. Therefore, we have the simple result that

[0042] $\begin{array}{cc}p\left(M_k\text{❘}n,t_0\right)=\frac{1}{2}p\left(n\text{❘}{M}_k,t_0\right)& \left(13\right)\end{array}$
Hence, we need to calculate the likelihood that we observed the data given the model class, p(n|M_k,t₀), a quantity called the evidence. We get this by integrating the likelihood function (9) over the model parameters using an appropriate prior for the parameters:
p(n|M_k,t₀)=∫dadλp(n|α,λ,M_k,t₀)p(α,λ|M_k,t₀)  (14)
where p(α,λ|M_k,t₀) is prior distributions of parameters (α,λ). For each model class, there will be a prior probability distribution for the parameters. For example, when no peak is present we choose as the prior:

[0043] $\begin{array}{cc}p\left(a,\lambda \text{❘}{t}_0,M_0\right)=\frac{1}{2}\delta \left(a\right)p\left(\lambda \text{❘}{t}_0,M_0\right)& \left(15\right)\end{array}$
where the ½ factor appears because the δ-function will be integrated only over the positive values of α. If we know nothing about the values of λ, then we choose a uniform prior, or some other prior that is very broad in λ-space on the grounds that when we integrate against (9) only the neighborhood of the maximum likelihood value of λ will contribute.

[0044] The likelihood of the data given the model (14) will be a bit less mysterious when we calculate for specific models momentarily. Before doing so, we summarize that to compare the hypotheses that there is, or is not, a peak in the window, we need to compute the ratio (10), which requires computation of (14) for each model. Setting the threshold for the level of evidence we require for ‘peak detection’ is a separate issue. Given that we have detected that a peak lies within a certain region of the time axis, we then fix the position and amplitude of the peak by maximizing the likelihood over the parameters (α,λ,t₀), as discussed in next section. Notice, however, that maximizing over t₀ has a different logical character than the other parameters, because we are comparing different data sets as we slide the window across the peak. The justification is based upon the physical reasonableness of the approach: the width of the window is large compared to the uncertainty in the position of the peak, hence near the maximum of the likelihood, most of the data being compared comes from overlapping windows.

[0045] 1.3 Parameter Fitting

[0046] As the window sliding across the spectrum one point at a time, the ratio (10) is calculated for each window position. We can then justify in which region in the spectrum we are confident that there is a peak. We then look into particular regions of interest, fit parameters (α, λ) by maximizing equation (9), i.e. solve the following equations, for each window t₀:

[0047] $\begin{array}{cc}\frac{\partial L}{\partial a}=0\frac{\partial L}{\partial \lambda }=0& \left(16\right)\end{array}$
where L(t₀), the natural logarithm of likelihood function (9):
L(t₀)=1n(p(n|α,λ,M_k,t₀)  (17)

[0048] Solving equations (16) gives the maximum likelihood estimations of parameters (α*,λ*). If the data is very informative, then the likelihood would sharply peak around the point (α*,λ*) in the parameter space. It is natural to Taylor expand (17) around (α*,λ*):

[0049] $\begin{array}{cc}\begin{array}{c}L\approx L\left(a^{*},{\lambda }^{*}\right)+\frac{1}{2}\left[\begin{array}{c}\frac{{\partial }^{2}L}{\partial a^2}{|}_{a^{*},{\lambda }^{*}}{\left(a-a^{*}\right)}^2+\frac{{\partial }^{2}L}{\partial a\partial \lambda }{|}_{a^{*},{\lambda }^{*}}\\ \left(a-a^{*}\right)\left(\lambda -{\lambda }^{*}\right)+\frac{{\partial }^{2}L}{\partial {\lambda }^2}{|}_{a^{*},{\lambda }^{*}}{\left(\lambda -{\lambda }^{*}\right)}^2\end{array}\right]\\ =L\left(a^{*},{\lambda }^{*}\right)+\frac{1}{2}{\left(X-X^{*}\right)}^{\prime }\nabla \nabla L\left(a^{*},{\lambda }^{*}\right)\left(X-X^{*}\right)\end{array}\mathrm{where}X=\left[\begin{array}{c}a\\ \lambda \end{array}\right],\mathrm{and}& \left(18\right)\\ \nabla \nabla L\left[\begin{array}{cc}\frac{{\partial }^{2}L}{\partial a^2}& \frac{{\partial }^{2}L}{\partial a\partial \lambda }\\ \frac{{\partial }^{2}L}{\partial a\partial \lambda }& \frac{{\partial }^{2}L}{\partial {\lambda }^2}\end{array}\right]& \left(19\right)\end{array}$
is the Hessian matrix. It follows from (18) that the likelihood function (9) is approximately:

[0050] $\begin{array}{cc}\begin{array}{c}p\left(n|a,\lambda ,M,t_0\right)=\exp \left[L\left(t_0\right)\right]\\ \approx \exp \left(L\left(a^{*},{\lambda }^{*}\right)+\frac{1}{2}{\left(X-X^{*}\right)}^{\prime }\nabla \nabla L\left(a^{*},{\lambda }^{*}\right)\left(X-X^{*}\right)\right)\\ =e^{L\left(a^{*},{\lambda }^{*}\right)}{e}^{\frac{1}{2}{\left(X-X^{*}\right)}^{\prime }\nabla \nabla L\left(a^{*},{\lambda }^{*}\right)\left(X-X^{*}\right)}\end{array}& \left(20\right)\end{array}$
This implies that the likelihood function (9) looks like a multivariate normal distribution in parameter space, with uncertainties:

[0051] $\begin{array}{cc}{\sigma }_a={\left(-\frac{{\partial }^{2}L}{\partial a^2}{|}_{a^{*},{\lambda }^{*}}\right)}^{-1/2}{\sigma }_{\lambda }={\left(-\frac{{\partial }^{2}L}{\partial {\lambda }^2}{|}_{a^{*},{\lambda }^{*}}\right)}^{-1/2}& \left(21\right)\end{array}$

[0052] Moreover, when computing the evidence, i.e. marginalizing the likelihood function (9) over the prior p(α,λ|M_k,t₀), as seen in equation (14), if the prior is independent of (α, λ), for example, α and λ are uniformly distributed in some region (α_min,α_max) and (λ_min,λ_max), substituting equation (20) into equation (14) results in the following:

[0053] $\begin{array}{cc}\begin{array}{c}p\left(n|M_k,t_0\right)=\int dad\lambda p\left(n|a,\lambda ,M_k,t_0\right)p\left(a,\lambda |M_k,t_0\right)\\ ={\int }_{a_{\min }}^{a_{\max }}{\int }_{{\lambda }_{\min }}^{{\lambda }_{\max }}dad{\mathrm{\lambda e}}^{L(a^{*},{\lambda }^{*})}\\ e^{\frac{1}{2}{\left(X-X^{*}\right)}^{\prime }\nabla \nabla L(a^{*},{\lambda }^{*})\left(X-X^{*}\right)}\\ \frac{1}{a_{\max }-a_{\min }}\frac{1}{{\lambda }_{\max }-{\lambda }_{\min }}\\ =\frac{1}{a_{\max }-a_{\min }}\frac{1}{{\lambda }_{\max }-{\lambda }_{\min }}{e}^{L\left(a^{*},{\lambda }^{*}\right)}\frac{{\left(2\pi \right)}^{m/2}}{\sqrt{\det \left[\nabla \nabla L\left(a^{*},{\lambda }^{*}\right)\right]}}\end{array}& \left(22\right)\end{array}$
where m in equation (22) is the dimension of parameter space. In the last step of the integration, the lower and upper boundaries of integration are extended to infinity. This is doable on the basis that (α_min,α_max) and (λ_min,λ_max) are large enough such that contributions from outside these regions are negligible. Note det[∇∇L(α*,λ*)] is the determinant of the Hessian matrix evaluated at (α*,λ*), and 1/det[∇∇L(α*,λ*)] is proportional to the volume within σ_α and σ_λ around (α*,λ*) in parameter space.

[0054] 2. Finding a Peak in a Time-of-flight Mass Spectrum.

[0055] We now consider some specific examples, following the discussion above. The first concerns the detection of a peak in white Gaussian noise. We show that the approach described above leads to what is called the ‘matched filter’ in the engineering literature. After that, we move on the applications in ‘counting’ problems, where the noise is due to the discreteness of the underlying process and a Poisson noise model is more appropriate. The later case is of particular interest as it is related to TOF-SIMS (Time of Flight Secondary Ionization Mass Spectrometer). The Gaussian noise case is also important as it may relate to other mass spectrometer, such as SELDI (Surface Enhanced Laser Desorption Ionization). In this section, we do not specify the peak shape for both Gaussian noise case and Poisson noise case, we simply use x_i to denote peak shape in general. In section 3, before we present some preliminary results, we will choose a particular peak shape, motivated by the kinematics of TOF-MS instruments.

[0056] Before we go to detailed calculations for Gaussian and Poisson noise cases, let us summarize that, in order to compare model M₀ v.s M₁, we need to evaluate the ratio (10), in which we need to compute (14) by marginalizing (9) over prior distribution p(α,λ|M_k,t₀)). Parameters are fitted by maximizing the likelihood function (9), via solving equation (16).

[0057] For convenience, from now on, except for t₀, when we refer to parameters associated with model M₀, a subscript 0 will be used, and a subscript 1 will be used for parameters associated with model M₁. For both models, a superscript of star * will be used for best estimations of parameters.

[0058] 2.1 Some Preliminary Results for Gaussian Noise Models

[0059] 2.1.1 Model Comparison

[0060] We assume the noise process is Gaussian and white:

[0061] $\begin{array}{cc}\langle \eta \left(t_j\right)\rangle =\mu ,\langle \eta \left(t_j\right)\eta \left(t_k\right)\rangle ={\delta }_{\mathrm{jk}},p_1\left(\eta \left(t_k\right)\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{\eta }}\exp \left(-\frac{{\left(\eta \left(t_k\right)-\mu \right)}^2}{2{\sigma }_{\eta }^2}\right)& \left(23\right)\end{array}$
wherein the subscript on the probability density indicates that this is the PDF for a single time step. We use this single-step PDF to compute the joint probability that we would have observed the N-step noise sequence η(t_k), k=1, 2 . . . N.

[0062] We want to compare the models having a peak in the window (M₁) to those lacking a peak in the window (M₀). Therefore, we need to evaluate (14) for both M₁ and M₀.

[0063] We first assume there is no peak present, only noise, and that each time step is independent of the others:

[0064] $\begin{array}{cc}{p}_N\left(\eta |{\mu }_0,{\sigma }_{\eta 0},M_0\right)=\prod _{j=1}^N\frac{1}{\sqrt{2\pi }{\sigma }_{\eta 0}}\exp \left(-\frac{1}{2{\sigma }_{\eta 0}^2}{\left(\eta \left(t_j\right)-{\mu }_0\right)}^2\right).& \left(24\right)\end{array}$
The notation η in the argument of the likelihood function represents an N-dimensional vector of noise values. Collecting terms, we find

[0065] $\begin{array}{cc}{p}_N\left(\eta |{\mu }_0,{\sigma }_{\eta 0},M_0\right)=\frac{1}{{\left(2\pi \right)}^{N/2}{\sigma }_{\eta 0}^N}\exp \left(-\frac{1}{2{\sigma }_{\eta 0}^2}\sum _{j=1}^N{\left(\eta \left(t_j\right)-{\mu }_0\right)}^2\right)& \left(25\right)\end{array}$
The exponent is a quadratic function of μ₀, therefore we complete the square and write it as

[0066] $\begin{array}{cc}\sum _{j=1}^N{\left(\eta \left(t_j\right)-{\mu }_0\right)}^2={N\left({\mu }_0-{\mu }_{*}\right)}^2+N{\sigma }_{*}^2\mathrm{where}& \left(26\right)\\ {\mu }_{*}=\frac{1}{N}\sum _{j=1}^N{\eta }_j,\mathrm{and}{\sigma }_{*}^2=\frac{1}{N}\sum _{j=1}^N{\left({\eta }_j-{\mu }_{*}\right)}^2& \left(27\right)\end{array}$
The likelihood function now becomes:

[0067] $\begin{array}{cc}{p}_N\left(\eta |{\mu }_0,{\sigma }_{\eta 0},M_0\right)=\frac{1}{{\left(2\pi \right)}^{N/2}{\sigma }_{\eta 0}^N}\exp \left(-\frac{N{\sigma }_{*}^2}{2{\sigma }_{\eta 0}^2}\right)\exp \left(-\frac{{N\left({\mu }_0-{\mu }_{*}\right)}^2}{2{\sigma }_{\eta 0}^2}\right)& \left(28\right)\end{array}$
From this form it is clear that the likelihood is maximized if we choose μ₀*=μ* and σ_η0*=σ*, with maximum:

[0068] $\begin{array}{cc}{p}_N^{*}\left(\eta \right)\equiv p_N^{*}\left(\eta |{\mu }_0^{*},{\sigma }_{\eta 0}^{*},M_0\right)=\frac{1}{{\left(2\pi e{\sigma }_{\eta 0}^{*2}\right)}^{N/2}}.& \left(29\right)\end{array}$
Comparing this value for different window positions can be used to measure the likelihood that the signal in each window is a Gaussian white noise process. Comparing different models using those parameter values in each model that maximizse the likelihood may not the most robust way to proceed. Instead, it may be preferable to marginalize over all possible values of the model parameters (in this case μ₀ and σ_η0) to get the probability for the model that there is no peak in the window (M₀) given the signal observed in the window, which is parameter-free. Doing the same marginalization for the case where a peak is assumed in the window, gives a parameter free probability for model M₁. These two parameter-free probabilities are directly comparable.

[0069] In evaluating equation (14), we now assume that we know nothing about the prior probabilities for the mean and variance of the noise process, aside from the requirement that they be probabilities in the parameters μ₀ and σ_η0. The simplest choice is to assume that they are constant over some range. Marginalizing, we have:

[0070] $\begin{array}{cc}{p}_N\left(\eta |M_0\right)={\int }_0^{\infty }\frac{d{\sigma }_{\eta 0}}{{\sigma }_{\eta 0_{\max }}}{\int }_{-\infty }^{\infty }\frac{d{\mu }_0}{{\mu }_{0_{\max }}}{p}_N\left(\eta |{\mu }_0,{\sigma }_{\eta 0},M_0\right)& \left(30\right)\end{array}$
Performing the integral over the parameter μ₀ first we are left with

[0071] $\begin{array}{cc}{p}_N\left(\eta |M_0\right)=\frac{1}{{\sigma }_{\eta 0_{\max }}{\mu }_{0_{\max }}}\frac{1}{{\left(2\pi \right)}^{\left(N-1\right)/2}\sqrt{N}}{\int }_0^{\infty }\frac{d{\sigma }_{\eta 0}}{{\sigma }_{\eta 0}^{N-1}}\exp \left(-\frac{N{\sigma }_{\eta 0}^{*2}}{2{\sigma }_{\mathrm{\eta 0}}^2}\right)& \left(31\right)\end{array}$
Changing variables to

[0072] $\begin{array}{cc}t=\frac{N{\sigma }_{\eta 0}^{*2}}{2{\sigma }_{\eta 0}^2}& \left(32\right)\end{array}$
we find

[0073] $\begin{array}{cc}{p}_N\left(\eta |M_0\right)=\frac{1}{2\sqrt{2}{\pi }^{\left(N-1\right)/2}}\frac{\Gamma \left(\frac{N}{2}-1\right)}{N^{\frac{N-1}{2}}{\sigma }_{\eta 0}^{*N-2}}\frac{1}{{\sigma }_{\eta 0_{\max }}{\mu }_{0_{\max }}}& \left(33\right)\end{array}$
Using Stirling's asymptotic series for the Γ-function (Γ(x+1)≈√{square root over (2πx)}x^xe^−x(1+ . . . ))

[0074] $\begin{array}{cc}\begin{array}{c}\Gamma \left(\frac{N}{2}-1\right)\approx {e^{-\frac{N}{2}+1}\left(\frac{N}{2}-1\right)}^{\frac{N}{2}-1-\frac{1}{2}}\sqrt{2\pi }\\ \approx {e^{-\frac{N}{2}+1}\left(\frac{N}{2}\right)}^{\frac{N-3}{2}}\sqrt{2\pi }\end{array}& \left(34\right)\end{array}$
we find:

[0075] $\begin{array}{cc}{p}_N\left(\eta |M_0\right)=\frac{1}{{\left(2\mathrm{\pi e}{\sigma }_{\eta 0}^{*2}\right)}^{N/2}}\frac{\pi \sqrt{2}{\sigma }_{\eta 0}^{*2}}{N}=p_N^{*}\frac{\pi \sqrt{2}{\sigma }_{\eta 0}^{*2}}{N}\frac{e}{{\sigma }_{\max }{\mu }_{\max }}& \left(35\right)\end{array}$
Notice that this still depends upon the window position. This result could have been anticipated by the following argument: the integration of the likelihood over μ₀ and σ_η0 will be dominated by the region of μ₀σ_η0-space in the immediate neighborhood of the maximum likelihood values μ₀* and σ_η0*. Over most of parameter space the likelihood is vanishingly small (if the data is ‘informative’). The volume of space over which the likelihood is not vanishingly small is given by the product of the half-widths of the likelihood function in each direction (which are essentially the uncertainties in the parameter estimates). The uncertainties can be estimated from the second derivative of the ln of the likelihood function, and are of the order of Δμ₀□σ_η0*/√{square root over (N)} and Δσ_η0□σ_η0*/√{square root over (2N)}. Therefore, we have the appealing result that (35) is simply

[0076] $\begin{array}{cc}{p}_N\left(\eta |M_0\right)=2\pi p_N^{*}\left(\frac{{\mathrm{\Delta \mu }}_0}{{\mu }_{0_{\max }}}\right)\left(\frac{{\mathrm{\Delta \sigma }}_{\eta 0}}{{\sigma }_{\eta 0_{\max }}}\right)e& \left(36\right)\end{array}$
The results (35) and (29) depend only upon the noise model and are independent of the rest of the model (e.g. peak shapes and positions, if there are any). Equation (33), (35) and (36) are different expressions for evidence.

[0077] The computation for p_N(η|M₀) is complete. However, we want to point out that there is a way that arrives at (35) and (36) while avoiding above tedious integrals. That is to utilize (22) and (29). The Hessian matrix, as defined in (19), can be easily calculated in case:

[0078] $\begin{array}{cc}\nabla \nabla L\left({\mu }_0^{*},{\sigma }_{\eta 0}^{*}\right)=\left[\begin{array}{cc}\frac{N}{{\sigma }_{\eta 0}^{*2}}& 0\\ 0& \frac{2N}{{\sigma }_{\eta 0}^{*2}}\end{array}\right]& \left(37\right)\end{array}$

[0079] Substitute (29) and (37) into (22):

[0080] $\begin{array}{cc}\begin{array}{c}{p}_N^{\prime }\left(\eta |M_0\right)\approx \frac{1}{{\mu }_{0_{\max }}}\frac{1}{{\sigma }_{\eta 0_{\max }}}\frac{2{\mathrm{\pi \sigma }}_{\eta 0}^{*2}}{\sqrt{2}N}\frac{1}{{\left(2\pi \right)}^{N/2}{\sigma }_{\eta 0}^{*N}}{e}^{-N/2}\\ =\frac{1}{{\mu }_{0_{\max }}}\frac{1}{{\sigma }_{\eta 0_{\max }}}\frac{{\sigma }_{\eta 0}^{*}}{\sqrt{N}}\frac{{\sigma }_{\eta 0}^{*}}{\sqrt{2N}}2\pi \frac{1}{{\left(2\pi \right)}^{N/2}{\sigma }_{\eta 0}^{*N}}{e}^{-N/2}\\ =\frac{1}{{\mu }_{0_{\max }}}\frac{1}{{\sigma }_{\eta 0_{\max }}}\frac{{\sigma }_{\eta 0}^{*}}{\sqrt{N}}\frac{{\sigma }_{\eta 0}^{*}}{\sqrt{2N}}2\pi p_N^{*}\left(\eta |{\mu }_0^{*},{\sigma }_{\eta 0}^{*},M_0\right)\\ =\frac{{\mathrm{\Delta \mu }}_0}{{\mu }_{0_{\max }}}\frac{{\mathrm{\Delta \sigma }}_{\eta 0}}{{\sigma }_{\eta 0_{\max }}}2\pi p_N^{*}\left(\eta |{\mu }_0^{*},{\sigma }_{\eta 0}^{*},M_0\right)\end{array}& \left(38\right)\end{array}$

[0081] Comparing (38) with (36), they differ by a factor e which comes from Stirling series expansion of

[0082] $\Gamma \left(\frac{N}{2}-1\right)$
in (34).

[0083] We now consider there is a peak in the window. Suppose the peak has a line shape described by x_i, which peaks near the center of the window. The observed data within the window may be written as:
s(t_j)=α₁x(t_j−t₀)+η(t_j) j=1 . . . N  (39)
where the amplitude (α₁) and arrival time (t₀) of the wavelet, as well as the noise amplitude σ_η1 and mean μ₁, are assumed to be unknown. The marginalization in this case looks like:

[0084] $\begin{array}{cc}\begin{array}{c}p\left(\eta |M_1\right)=\int \int \int d{a}_{1}d{\mu }_{1}d{\sigma }_{\eta 1}{p}_N\left(\eta |a_1,{\mu }_1,{\sigma }_{\mathrm{\eta 1}},t_0,M_1\right)\\ \frac{1}{a_{1_{\max }}}\frac{1}{{\mu }_{1_{\max }}}\frac{1}{{\sigma }_{{\mathrm{\eta 1}}_{\max }}}\end{array}& \left(40\right)\end{array}$
where:

[0085] $\begin{array}{cc}\begin{array}{c}{p}_N\left(\eta |a_1,{\mu }_1,{\sigma }_{\mathrm{\eta 1}},t_0,M_1\right)=\frac{1}{{\left(\sqrt{2\pi }{\sigma }_{\mathrm{\eta 1}}\right)}^N}\exp \\ \left(-\frac{1}{2{\sigma }_{\mathrm{\eta 1}}^2}\sum _{i=1}^N{\left(s_i-a_1{x}_i-u_1\right)}^2\right)\end{array}& \left(41\right)\end{array}$
is the realization of likelihood function (9). As in the case of no peak, complete the square, let:

[0086] $u_{*}=\frac{1}{N}\sum \left(s_i-a_1{x}_i\right)$ ${\sigma }_{{\eta }^{*}}^2=\frac{1}{N}\sum {\left(s_i-a_1{x}_i-u_{*}\right)}^2$
one may get:

[0087] $\begin{array}{cc}\begin{array}{c}p\left(\eta |M_1\right)=\int \int d{a}_{1}d{\mu }_{1}d{\sigma }_{\mathrm{\eta 1}}\frac{1}{{\left(\sqrt{2\pi }{\sigma }_{\mathrm{\eta 1}}\right)}^N}\\ e^{-\frac{1}{2{\sigma }_{\mathrm{\eta 1}}^2}{N\left({\mu }_1-{\mu }_{*}\right)}^2}{e}^{-\frac{1}{2{\sigma }_{\mathrm{\eta 1}}^2}N{\sigma }_{{\eta }^{*}}^2}\\ \frac{1}{a_{1_{\max }}}\frac{1}{{\mu }_{1_{\max }}}\frac{1}{{\sigma }_{{\mathrm{\eta 1}}_{\max }}}\\ =\frac{1}{{\left(\sqrt{2\pi }\right)}^{N-1}\sqrt{N}}\int d{a}_{1}d{\sigma }_{\mathrm{\eta 1}}\frac{1}{{\sigma }_{\mathrm{\eta 1}}^{N-1}}{e}^{-\frac{1}{2{\sigma }_{\mathrm{\eta 1}}^2}N{\sigma }_{{\eta }^{*}}^2}\\ \frac{1}{a_{1_{\max }}}\frac{1}{{\mu }_{1_{\max }}}\frac{1}{{\sigma }_{{\mathrm{\eta 1}}_{\max }}}\\ =\int d{a}_1\frac{1}{2\sqrt{2}{\left(N\pi \right)}^{\frac{N-1}{2}}{\sigma }_{{\eta }^{*}}^{N-2}}\Gamma \left(\frac{N}{2}-1\right)\\ \frac{1}{a_{1_{\max }}}\frac{1}{{\sigma }_{{\mathrm{\eta 1}}_{\max }}}\frac{1}{{\mu }_{1_{\max }}}\end{array}& \left(42\right)\end{array}$
Note, if we let (α₁*,μ₁*,σ_η1*) be the point where the likelihood is maximized, the derivation of (α₁*,μ₁*,σ_η1*) will be given momentarily, then

[0088] $\begin{array}{cc}\begin{array}{c}{\sigma }_{{\eta }^{*}}^2=\frac{1}{N}\sum {\left(s_i-a_1{x}_i-u_{*}\right)}^2\\ =\frac{1}{N}\sum {\left(s_i-a_1^{*}{x}_i-u_1^{*}-\left(a_1-a_1^{*}\right)x_i-\left(u_{*}-u_1^{*}\right)\right)}^2\\ ={\sigma }_{\mathrm{\eta 1}}^{*}+{{\sigma }_x^2\left[a_1-a_1^{*}-\frac{\sum \left(s_i-a_1^{*}{x}_i-u_1^{*}\right)\left(x_i-\stackrel{\_}{x}\right)}{{\sigma }_x^{2}N}\right]}^2-\\ \frac{{\left[\sum \left(s_i-a_1^{*}{x}_i-u_1^{*}\right)\left(x_i-\stackrel{\_}{x}\right)\right]}^2}{{\sigma }_x^2{N}^2}\\ ={A\left(a-a_1^{\prime }\right)}^2+B\end{array}& \left(43\right)\end{array}$
where

[0089] $\begin{array}{cc}\stackrel{\_}{x}=\frac{\sum x_i}{N}{\sigma }_x^2=\frac{1}{N}\sum {\left(x_i-\stackrel{\_}{x}\right)}^{2}A={\sigma }_x^{2}B={\sigma }_{\mathrm{\eta 1}}^{*2}-\frac{{\left[\sum \left(s_i-a_1^{*}{x}_i-u_1^{*}\right)\left(x_i-\stackrel{\_}{x}\right)\right]}^2}{{\sigma }_x^2{N}^2}{a}_1^{\prime }=a_1^{*}+\frac{\sum \left(s_i-a_1^{*}{x}_i-u_1^{*}\right)\left(x_i-\stackrel{\_}{x}\right)}{{\sigma }_x^{2}N}& \left(44\right)\end{array}$
then marginalization may be written as:

[0090] $\begin{array}{cc}\begin{array}{c}p\left(\eta |M_1\right)=\int da\frac{1}{2\sqrt{2}{{\left(N\pi \right)}^{\frac{N-1}{2}}\left[{A\left(a-a_1^{\prime }\right)}^2+B\right]}^{\frac{N-2}{2}}}\\ \Gamma \left(\frac{N}{2}-1\right)\frac{1}{a_{\max }}\frac{1}{{\sigma }_{{\eta }^{\max }}}\frac{1}{{\mu }_{\max }}\end{array}& \left(45\right)\end{array}$
Changing variable t=(α₁−α₁′)² gives:

[0091] $\begin{array}{cc}p\left(\eta |M_1\right)=\frac{1}{4\sqrt{2}{\pi }^{\frac{N}{2}-1}{N}^{\frac{N-1}{2}}}\frac{B^{-\frac{N}{2}+\frac{3}{2}}}{\sqrt{A}}\Gamma \left(\frac{N}{2}-\frac{3}{2}\right)& \left(46\right)\end{array}$
Now, having p(η|M₀) and p(η|M₁) computed, one may substitute them into (10) to see which model is more preferable.

[0092] Up to this point, the model comparison is complete, as we derived in equations (33) and (46), and equation (10) can be used to compare two hypotheses. Yet, as in parallel with what we did in the M₀ case, we can also arrive at an equation similar to equation (46) by computing the Hessian matrix, but this is going to need (α₁*,μ₁*,σ_η1*), the best estimation of (α₁,μ₁,σ_η1).

[0093] 2.1.2 Parameter Fitting for Gaussian Noise

[0094] Once it is determined which model is more appealing, parameters may be fitted via maximizing appropriate likelihood.

[0095] For the case where there is no peak in the window, the likelihood (9) is simply:

[0096] $\begin{array}{cc}\begin{array}{c}{p}_N\left(\eta |{\mu }_0,{\sigma }_{\mathrm{\eta 0}},M_0\right)=\frac{1}{{\left(2\pi \right)}^{N/2}{\sigma }_{\mathrm{\eta 0}}^N}\exp \\ \left(-\frac{1}{2{\sigma }_{\mathrm{\eta 0}}^2}\sum _{j=1}^N{\left(\eta \left(t_j\right)-{\mu }_0\right)}^2\right)\end{array}& \left(47\right)\end{array}$
where
η(t_j)=s(t_j)  (48)
as only noise is assumed in presence. As it is pointed out in previous section, (47) is maximized at

[0097] $\begin{array}{cc}{\mu }_0^{*}=\frac{1}{N}\sum s_i{\sigma }_{\mathrm{\eta 0}}^{*}=\frac{1}{N}\sum {\left(s_i-{\mu }_0\right)}^2& \left(49\right)\end{array}$

[0098] The uncertainties are, as one can see from the Hessian matrix (37):

[0099] $\begin{array}{cc}{\mathrm{\Delta \mu }}_0=\frac{{\sigma }_{\mathrm{\eta 0}}^{*}}{\sqrt{N}}{\mathrm{\Delta \sigma }}_{\mathrm{\eta 0}}=\frac{{\sigma }_{\mathrm{\eta 0}}^{*}}{\sqrt{2N}}& \left(50\right)\end{array}$

[0100] If it is preferred that there is a peak in the window, then the likelihood function (9) looks like:

[0101] $\begin{array}{cc}\begin{array}{c}{p}_N\left(\eta |a_1,{\mu }_1,{\sigma }_{\mathrm{\eta 1}},M_1\right)=\frac{1}{{\left(\sqrt{2\pi }{\sigma }_{\mathrm{\eta 1}}\right)}^N}\exp \left(-\frac{1}{2{\sigma }_{\mathrm{\eta 1}}^2}\sum _{i=1}^N{\eta }_i^2\right)\\ =\frac{1}{{\left(\sqrt{2\pi }{\sigma }_{\mathrm{\eta 1}}\right)}^N}\exp \\ \left(-\frac{1}{2{\sigma }_{\mathrm{\eta 1}}^2}\sum _{i=1}^N{\left(s_i-a_1{x}_i-u_1\right)}^2\right)\end{array}& \left(51\right)\end{array}$
the natural logarithm of the likelihood is:

[0102] $\begin{array}{cc}\begin{array}{c}L=\log \left[p_N\left(\eta |a_1,{\mu }_1,{\sigma }_{\mathrm{\eta 1}},M_1\right)\right]\\ =-\frac{N}{2}\log \left(2\pi \right)-N\log \left({\sigma }_{\mathrm{\eta 1}}\right)-\frac{1}{2{\sigma }_{\mathrm{\eta 1}}^2}\sum _{i=1}^N{\left(s_i-a_1{x}_i-{\mu }_1\right)}^2\end{array}& \left(52\right)\end{array}$
In order to maximize this value, let:

[0103] $\begin{array}{cc}\frac{\partial L}{\partial a_1}=-\frac{1}{{\sigma }_{\mathrm{\eta 1}}^2}\sum _{i=1}^N\left(s_i-a_1{x}_i-{\mu }_1\right)\left(-x_i\right)=0\frac{\partial L}{\partial {\mu }_1}=-\frac{1}{{\sigma }_{\mathrm{\eta 1}}^2}\sum _{i=1}^N\left(s_i-a_1{x}_i-{\mu }_1\right)\left(-1\right)=0\frac{\partial L}{\partial {\sigma }_{\mathrm{\eta 1}}}=-\frac{N}{{\sigma }_{\mathrm{\eta 1}}}+\frac{1}{{\sigma }_{\mathrm{\eta 1}}^3}\sum _{i=1}^N{\left(s_i-a_1{x}_i-{\mu }_1\right)}^2=0& \left(53\right)\end{array}$
Solving these equations gives:

[0104] $\begin{array}{cc}\begin{array}{c}{a}_1^{*}=\frac{\stackrel{\_}{\left(\mathrm{sx}\right)}-\stackrel{\_}{s}\stackrel{\_}{x}}{x^2-{\stackrel{\_}{x}}^2}\\ =\frac{\stackrel{\_}{\left(s-\stackrel{\_}{s}\right)\left(x-\stackrel{\_}{x}\right)}}{\stackrel{\_}{{\left(x-\stackrel{\_}{x}\right)}^2}}\\ =\frac{\stackrel{\_}{\left(s-\stackrel{\_}{s}\right)\left(x-\stackrel{\_}{x}\right)}}{{\sigma }_x^2}\end{array}{\mu }_1^{*}=\stackrel{\_}{s}-a_1^{*}\stackrel{\_}{x}{\sigma }_{\mathrm{\eta 1}}^{*}=\stackrel{\_}{{\left(s_i-\stackrel{\_}{s}-a_1\left(x_i-\stackrel{\_}{x}\right)\right)}^2}={\sigma }_s^2-a_1^2{\sigma }_x^2& \left(54\right)\end{array}$

[0105] The log of the likelihood function at (α₁*,μ₁*,σ_η1*) is simply:

[0106] $\begin{array}{cc}L\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\mathrm{\eta 1}}^{*}\right)=-N\log \left({\sigma }_{\mathrm{\eta 1}}^{*}\right)-\frac{N}{2}\log \left(2\pi \right)-\frac{N}{2}& \left(55\right)\end{array}$

[0107] The elements in Hessian matrix ∇∇L(α₁,μ₁,σ_η1) evaluated at (α₁*,μ₁*,σ_η1*) are:

[0108] $\begin{array}{cc}-\frac{{\partial }^{2}L}{\partial a_1^2}{|}_{\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\mathrm{\eta 1}}^{*}\right)}=\frac{1}{{\sigma }_{\mathrm{\eta 1}}^{*2}}\sum _{i=1}^N{x}_i^2=\frac{N\left({\sigma }_x^2+{\stackrel{\_}{x}}^2\right)}{{\sigma }_{\mathrm{\eta 1}}^{*2}}-\frac{{\partial }^{2}L}{\partial {\mu }_1^2}{|}_{\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\mathrm{\eta 1}}^{*}\right)}=\frac{1}{{\sigma }_{\mathrm{\eta 1}}^{*2}}\sum _{i=1}^N\left(1\right)=\frac{N}{{\sigma }_{\mathrm{\eta 1}}^{*2}}\begin{array}{c}-\frac{{\partial }^{2}L}{\partial {\sigma }_{\mathrm{\eta 1}}^2}{|}_{\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\mathrm{\eta 1}}^{*}\right)}=-\frac{N}{{\sigma }_{\mathrm{\eta 1}}^{*2}}+\frac{3}{{\sigma }_{\mathrm{\eta 1}}^{*4}}\sum _{i=1}^N{\left(s_i-a_1^{*}{x}_i-{\mu }_1^{*}\right)}^2\\ =\frac{2N}{{\sigma }_{\mathrm{\eta 1}}^{*2}}\end{array}-\frac{{\partial }^{2}L}{\partial a_1\partial {\mu }_1}{|}_{\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\mathrm{\eta 1}}^{*}\right)}=\frac{\sum _{i=1}^N{x}_i}{{\sigma }_{\mathrm{\eta 1}}^{*2}}-\frac{{\partial }^{2}L}{\partial {\mu }_1\partial {\sigma }_{\mathrm{\eta 1}}}{|}_{\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\mathrm{\eta 1}}^{*}\right)}=\frac{2}{{\sigma }_{\mathrm{\eta 1}}^{*3}}\sum _{i=1}^N\left(s_i-a_1^{*}{x}_i-{\mu }_1^{*}\right)=0-\frac{{\partial }^{2}L}{\partial a_1\partial {\sigma }_{\mathrm{\eta 1}}}{|}_{\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\mathrm{\eta 1}}^{*}\right)}=\frac{2}{{\sigma }_{\mathrm{\eta 1}}^{*3}}\sum _{i=1}^N\left(s_i-a_1^{*}{x}_i-{\mu }_1^{*}\right)x_i=0& \left(56\right)\end{array}$

[0109] The last 2^nd term is zero obviously because μ₁*= s−α₁* x. The last term is also zero because it is proportional to the correlation between the residue, which is presumably white noise with zero mean, and the target x_i, which is zero. Nevertheless, it can be shown by following:

[0110] $\begin{array}{cc}\begin{array}{c}\sum _{i=1}^N\left(s_i-a_1^{*}{x}_i-{\mu }_1^{*}\right)x_i=\sum _{i=1}^N\left(s_i-a_1^{*}{x}_i-\stackrel{\_}{s}+a_1^{*}\stackrel{\_}{x}\right)x_i\\ =\sum _{i=1}^N{s}_i{x}_i-a_1^{*}\sum _{i=1}^N{x}_i^2-\\ \stackrel{\_}{s}\sum _{i=1}^N{x}_i+a_1^{*}\stackrel{\_}{x}\sum _{i=1}^N{x}_i\\ =N\left(\stackrel{\_}{\mathrm{sx}}-a_1^{*}\stackrel{\_}{x^2}-\stackrel{\_}{s}\stackrel{\_}{x}+a_1^{*}{\stackrel{\_}{x}}^2\right)\\ =N\left(\stackrel{\_}{\mathrm{sx}}-\stackrel{\_}{s}\stackrel{\_}{x}-a_1^{*}\left(\stackrel{\_}{x^2}-{\stackrel{\_}{x}}^2\right)\right)\\ =N\left(\stackrel{\_}{\left(s-\stackrel{\_}{s}\right)\left(x-\stackrel{\_}{x}\right)}-a_1^{*}{\sigma }_x^2\right)\\ =0\end{array}& \left(57\right)\end{array}$

[0111] With all elements of Hessian matrix computed, we can now write out the Hessian matrix ∇∇L(α₁*,μ₁*,σ_η1*):

[0112] $\begin{array}{cc}\nabla \nabla L\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\mathrm{\eta 1}}^{*}\right)=\left[\begin{array}{ccc}\frac{N\left({\sigma }_x^2+{\stackrel{\_}{x}}^2\right)}{{\sigma }_{\mathrm{\eta 1}}^{*2}}& \frac{\sum _{i=1}^N{x}_i}{{\sigma }_{\mathrm{\eta 1}}^{*2}}& 0\\ \frac{\sum _{i=1}^N{x}_i}{{\sigma }_{\mathrm{\eta 1}}^{*2}}& \frac{N}{{\sigma }_{\mathrm{\eta 1}}^{*2}}& 0\\ 0& 0& \frac{2N}{{\sigma }_{\mathrm{\eta 1}}^{*2}}\end{array}\right].& \left(58\right)\end{array}$

[0113] From which, one can get the uncertainties:

[0114] $\begin{array}{cc}\Delta a_1=\frac{{\sigma }_{\mathrm{\eta 1}}^{*}}{\sqrt{N\left({\sigma }_x^2+{\stackrel{\_}{x}}^2\right)}}\Delta {\mu }_1=\frac{{\sigma }_{\mathrm{\eta 1}}^{*}}{\sqrt{N}}\Delta {\sigma }_{\mathrm{\eta 1}}=\frac{{\sigma }_{\mathrm{\eta 1}}^{*}}{\sqrt{2N}},& \left(59\right)\end{array}$

[0115] From (22), we have:

[0116] $\begin{array}{cc}\begin{array}{c}{p}^{\prime }\left(\eta |M_1\right)\approx \frac{1}{a_{1_{\max }}}\frac{1}{{\mu }_{1_{\max }}}\frac{1}{{\sigma }_{{\mathrm{\eta 1}}_{\max }}}\\ \frac{{\left(2\pi \right)}^{3/2}}{\sqrt{\text{det}\left(\nabla \nabla L\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\eta 1}^{*}\right)\right)}}{e}^{L\left(a_1^{*},{\mu }_1^{*},{\sigma }_{\eta 1}^{*}\right)}\end{array}& \left(60\right)\end{array}$