Spectral baseline estimation using penalized least squares with weights derived from the Bayesian method

Spectral baseline estimation using penalized least squares with weights derived from the Bayesian method 1 2 corresp first-author † wangqian2016@impcas.ac.cn † wangqian2016@impcas.ac.cn wangqian2016@impcas.ac.cn https://orcid.org/0000-0002-0193-0301 1 3 https://orcid.org/0000-0003-3513-8870 1 1 2 https://orcid.org/0000-0003-2408-8089 1 2 Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000, China Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000, China University of Chinese Academy of Sciences, Beijing 100049, China University of Chinese Academy of Sciences, Beijing 100049, China Van Swinderen Institute, University of Groningen, Groningen 9747 AA, the Netherlands Van Swinderen Institute, University of Groningen, Groningen 9747 AA, the Netherlands The penalized least squares (PLS) method with appropriate weights has proven to be a successful baseline estimation method for various spectral analyses. It can extract the baseline from the spectrum while retaining the signal peaks in the presence of random noise. The algorithm is implemented by iterating over the weights of the data points. In this study, we propose a new approach for assigning weights based on the Bayesian rule. The proposed method provides a self-consistent weighting formula and performs well, particularly for baselines with different curvature components. This method was applied to analyze Schottky spectra obtained in 86 Kr projectile fragmentation measurements in the experimental Cooler Storage Ring (CSRe) at Lanzhou. It provides an accurate and reliable storage lifetime with a smaller error bar than existing PLS methods. It is also a universal baseline-subtraction algorithm that can be used for spectrum-related experiments, such as precision nuclear mass and lifetime measurements in storage rings. Penalized least squares Baseline correction Bayesian rule Spectrum analysis Introduction 1 A typical one-dimensional spectrum obtained from a spectroscopy experiment is an array that encodes intensity information at a continuous equidistant interval. This can be assumed to be the superposition of a slowly varying baseline, random noise, and a few narrow signal peaks. These signal peaks directly or indirectly reveal the physical and chemical properties of the measured objects. To correctly determine the peak position, intensity, and area, many methods have been developed to eliminate the influence of the baseline and noise, which are commonly observed in spectra. The main goal of these methods is to subtract the baseline from the spectrum while retaining the signal peaks. Traditional methods for estimating baselines are based on simple or modified polynomial functions that fit the data. Even with automatic optimization in linear programming [ 1 1 , 2 2 ], the performance of the fitting procedure is highly dependent on the knowledge of the baseline, i.e., both the semi-empirical function and initial guess of its fitting coefficients. Its accuracy depends on user experience [ 3 3 ]. As the degrees of freedom of nonlinear fitting increase, fitting instability may become a fatal problem with this method [ 4 4 ]. Modern methods avoid any assumptions about the baseline-specific functional form by using digital filtering [ 5 5 , 6 6 ] and clipping [ 7 7 ]. For example, the rolling ball method [ 5 5 ] operates as a frequency-differentiated nonlinear digital filter by gradually increasing the diameter of the ball, rolling it below the spectrum, and marking the baseline as the locus of the ball at each point. Peak distortion occasionally occurs after baseline correction and its effect is sometimes not negligible. Another example is the wavelet transform method, which is effective in suppressing random noise in data [ 6 6 ]. This method assumes that the baseline is properly separated from the signal in the transformed domain. However, real-world signals do not always fulfill this assumption. The third example is the sensitive nonlinear iterative peak clipping method, which is included in CERN-ROOT software packages [ 7 7 , 8 8 ]. This method is primarily used in the baseline correction of the proton-induced X-ray emission (PIXE) spectrum [ 9 9 ]. The algorithm compares the values of the data points with a local baseline approximation based on binomial expansions, iteratively stripping sharp peaks and preserving the baseline. One limitation of this approach is its low correction efficiency in regions with relatively large baseline slopes. In the past two decades, penalized least squares (PLS) baseline estimation has been widely used in the analysis of various types of spectra to overcome the drawbacks mentioned above. The concept of PLS was first proposed in 1922 by Whittaker in his paper on constructing effective smoothers for data [ 10 10 ]. In the 1990s, a similar procedure known as the Hodrick-Prescott filter became popular in the econometric literature [ 11 11 ]. In 2003, Eilers revived the algorithm [ 12 12 , 13 13 ] and applied it to peak alignment and baseline correction [ 14 14 ]. The core of the algorithm is to balance two conflicting goals: fidelity to the data, and smoothness of the recommended baseline estimates. Using carefully constructed weights for each data point, the algorithm can effectively extract the baseline shape from the data [ 14 14 - 16 16 14-16 ]. Owing to the efficient processing of sparse matrices by modern computers, PLS-based algorithms are fast and flexible in terms of implementation. These algorithms do not require knowledge of the specific functions of the baseline and are capable of dealing with complicated baseline shapes. However, the weights of existing methods are empirically determined and may be inconsistent under different conditions, leading to unreasonable results. In this study, we derived a self-consistent function of weights that adheres to the Bayesian rule for a PLS-based baseline estimation algorithm. In the subsequent sections, the theory of the algorithm and its implementation procedure are introduced. The proposed method is reliable and even more effective than existing PLS-based algorithms in the analysis of synthetic and real Schottky experimental data. In particular, this method stands out in estimating rapidly changing baselines in the presence of signal peaks and strong noise. Theory 1 PLS algorithm 2 The concept of PLS originated from the Whittaker smoother [ 10 10 ], which was initially designed to smooth one-dimensional equally spaced data. Consider a data sequence y of length N , sampled at equal time intervals. The algorithm smoothes a noisy sequence y into a smooth sequence z by minimizing the following objective function: <math display="block" id="M1"><mrow><mi>S</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mo stretchy="false">)</mo><mo>=</mo><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>y</mtext></mstyle><mo>−</mo><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><msup><mo stretchy="false">)</mo><mtext>T</mtext></msup><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>y</mtext></mstyle><mo>−</mo><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mo stretchy="false">)</mo><mo>+</mo><mi>λ</mi><msup><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mtext>T</mtext></msup><msup><mstyle mathvariant="bold"><mtext>D</mtext></mstyle><mtext>T</mtext></msup><mstyle mathvariant="bold"><mtext>DZ</mtext></mstyle><mn>,</mn></mrow></math> where D is an m -th-order difference matrix with N - m rows and N columns. In this study, the discussion is limited to the 2nd-order difference ( m = 2). Thus, D is expressed as <math display="block" id="M2"><mrow><mstyle mathvariant="bold"><mtext>D</mtext></mstyle><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mn>1</mn></mtd><mtd><mrow><mo>−</mo><mn>2</mn></mrow></mtd><mtd><mn>1</mn></mtd><mtd><mn>0</mn></mtd><mtd><mo>⋯</mo></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd><mtd><mrow><mo>−</mo><mn>2</mn></mrow></mtd><mtd><mn>1</mn></mtd><mtd><mo>⋯</mo></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋱</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mo>⋯</mo></mtd><mtd><mn>1</mn></mtd><mtd><mrow><mo>−</mo><mn>2</mn></mrow></mtd><mtd><mn>1</mn></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>.</mn></mrow></math> The first term in Eq. (1) represents the fidelity of z to data y, whereas the second term represents the smoothness of z. The balance between fidelity and smoothness is controlled by the parameter λ . A larger λ leads to a smoother z. In a typical spectrum, some data points have only baseline and noise components, whereas the rest of the points contain additional signal components. For baseline estimation, it is optimal to obtain a slowly varying baseline if there is a way to distinguish between both cases by considering only the former and ignoring the latter. Therefore, combining the Whittaker smoother with the concept of asymmetric least squares [ 17 17 ], Eilers and Boelens [ 14 14 ] first introduced a weight for each point to obtain a smooth baseline. Let W be a diagonal matrix with a weight array w of data points on its diagonal. Thus, the objective function in Eq. (1) becomes <math display="block" id="M3"><mrow><mi>S</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mo stretchy="false">)</mo><mo>=</mo><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>y</mtext></mstyle><mo>−</mo><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><msup><mo stretchy="false">)</mo><mtext>T</mtext></msup><mstyle mathvariant="bold"><mtext>W</mtext></mstyle><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>y</mtext></mstyle><mo>−</mo><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mo stretchy="false">)</mo><mo>+</mo><mi>λ</mi><msup><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mtext>T</mtext></msup><msup><mstyle mathvariant="bold"><mtext>D</mtext></mstyle><mtext>T</mtext></msup><mstyle mathvariant="bold"><mtext>DZ</mtext></mstyle><mn>.</mn></mrow></math> Setting the partial derivatives <math display="block" id="M4"><mrow><mfrac><mrow><mo>∂</mo><mi>S</mi></mrow><mrow><mo>∂</mo><msup><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mtext>T</mtext></msup></mrow></mfrac><mo>=</mo><mo>−</mo><mn>2</mn><mstyle mathvariant="bold"><mtext>W</mtext></mstyle><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>y</mtext></mstyle><mo>−</mo><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mo stretchy="false">)</mo><mo>+</mo><mn>2</mn><mi>λ</mi><msup><mstyle mathvariant="bold"><mtext>D</mtext></mstyle><mtext>T</mtext></msup><mstyle mathvariant="bold"><mtext>DZ</mtext></mstyle></mrow></math> to zero, we obtain the solution for minimizing Eq. (3): <math display="block" id="M5"><mrow><mstyle mathvariant="bold"><mtext>z</mtext></mstyle><mo>=</mo><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>W</mtext></mstyle><mo>+</mo><mi>λ</mi><msup><mstyle mathvariant="bold"><mtext>D</mtext></mstyle><mtext>T</mtext></msup><mstyle mathvariant="bold"><mtext>D</mtext></mstyle><msup><mo stretchy="false">)</mo><mrow><mo>−</mo><mn>1</mn></mrow></msup><mstyle mathvariant="bold"><mtext>W</mtext></mstyle><mstyle mathvariant="bold"><mtext>y</mtext></mstyle><mn>.</mn></mrow></math> Intuitively, the weight should be 0 if the data point is definitely in the peak region and 1 if the data point does not contain any peak component of the signal. However, the regions containing no signal are difficult to determine because of the presence of noise and peaks. Thus, several methods have been proposed using iterative algorithms to circumvent peak finding. A pioneering method is the asymmetric least squares (AsLS) method developed by Eilers and Boelens [ 14 14 ]. In this method, an asymmetry parameter p , recommended to be set between 0.001 and 0.01, is introduced, and the weights are assigned as <math display="block" id="M6"><mrow><msub><mi>w</mi><mi>i</mi></msub><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mtext>  </mtext><mi>p</mi></mrow></mtd><mtd columnalign="left"><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>></mo><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mn>1</mn><mo>−</mo><mi>p</mi></mrow></mtd><mtd columnalign="left"><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>⩽</mo><msub><mi>z</mi><mi>i</mi></msub><mn>.</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> The algorithm first sets w to an array of all 1s and solves Eq. (5) to obtain the initial z. Then, the data points are divided into two zones: positive deviations y i > z i and negative residuals <math id="M7"><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>⩽</mo><msub><mi>z</mi><mi>i</mi></msub></mrow></math> . By applying the asymmetric weight function given by Eq. (6), the positive deviations, which are supposed to contain signals, gain fewer weights than the negative residuals in the next iteration. The baseline, i.e., z, is updated by solving Eqs. (5) using the new w. By applying this new z, we can again obtain the updated weight array w. The iteration ends when the weight array w remains unchanged or the number of iterations reaches the maximum pre-set value. The final updated baseline z is the result of the estimation method. Given that the AsLS method only assigns two values to the weights of each point, this assignment is unreasonable for data containing different signal peak heights. To improve the weight assignment, the adaptive iteratively reweighted penalized least squares (airPLS) [ 15 15 ] and asymmetrically reweighted penalized least squares (arPLS) [ 16 16 ] methods were successively developed. The airPLS algorithm considers the effects of signal peak heights. It uses an exponential function to ensure that the weights are small for large peak heights and vice versa. This method also adapts to two-segment weight assignment. The subtraction d = y - z, d - is part of d in the region where <math id="M8"><mrow><mo>{</mo><msub><mi>d</mi><mi>i</mi></msub><mo>|</mo><msub><mi>d</mi><mi>i</mi></msub><mo>=</mo><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub><mo><</mo><mn>0,</mn><msub><mi>d</mi><mi>i</mi></msub><mo>∈</mo><mi>d</mi><mo>}</mo></mrow></math> , and d + is part of the region where <math id="M9"><mrow><mo>{</mo><msub><mi>d</mi><mi>i</mi></msub><mo>|</mo><msub><mi>d</mi><mi>i</mi></msub><mo>=</mo><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub><mo>></mo><mn>0,</mn><msub><mi>d</mi><mi>i</mi></msub><mo>∈</mo><mi>d</mi><mo>}</mo></mrow></math> . The weight array w for each iteration step t in the airPLS method is constructed as follows: <math display="block" id="M10"><mrow><msub><mi>w</mi><mi>i</mi></msub><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mn>0</mn></mtd><mtd columnalign="left"><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>⩾</mo><msub><mi>z</mi><mi>i</mi></msub></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mtext>exp</mtext><mrow><mo>[</mo><mrow><mfrac><mrow><mi>t</mi><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><msub><mi>d</mi><mi>j</mi></msub><mo>∈</mo><msub><mi>d</mi><mo>−</mo></msub></mrow></msub><mrow><mo>|</mo><msub><mi>d</mi><mi>j</mi></msub><mo>|</mo></mrow></mstyle></mrow></mfrac></mrow><mo>]</mo></mrow></mrow></mtd><mtd columnalign="left"><mrow><msub><mi>y</mi><mi>i</mi></msub><mo><</mo><msub><mi>z</mi><mi>i</mi></msub><mn>.</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> The step number t in Eq. (7) is introduced to speed up the convergence of the iterations. The iteration of this method stops when the maximum number of iterations or the following termination criterion is fulfilled: <math display="block" id="M11"><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><msub><mi>d</mi><mi>j</mi></msub><mo>∈</mo><msub><mi>d</mi><mo>−</mo></msub></mrow></munder></mstyle><mo>|</mo><msub><mi>d</mi><mi>j</mi></msub><mo>|</mo><mo><</mo><mn>0.001</mn><mo>×</mo><mo>|</mo><mi>y</mi><mo>|</mo><mn>.</mn></mrow></math> Given that both the airPLS and AsLS methods attempt to set the weights of the data points to zero, or very close to zero, they tend to result in a baseline lower than the accurate baseline. To overcome this problem, the arPLS method is proposed, in which an asymmetric weight array is constructed as follows: <math display="block" id="M12"><mrow><msub><mi>w</mi><mi>i</mi></msub><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mtext>exp</mtext><mrow><mo>[</mo><mrow><mfrac><mn>2</mn><mrow><msub><mi>σ</mi><mo>−</mo></msub></mrow></mfrac><mo stretchy="false">(</mo><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>−</mo><mo stretchy="false">(</mo><mo>−</mo><msub><mi>m</mi><mo>−</mo></msub><mo>+</mo><mn>2</mn><msub><mi>σ</mi><mo>−</mo></msub><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow><mo>]</mo></mrow></mrow></mfrac></mrow></mtd><mtd><mrow><mtext> </mtext><msub><mi>y</mi><mi>i</mi></msub><mo>⩾</mo><msub><mi>z</mi><mi>i</mi></msub></mrow></mtd></mtr><mtr><mtd><mn>1</mn></mtd><mtd><mrow><mtext> </mtext><msub><mi>y</mi><mi>i</mi></msub><mo>⩽</mo><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn></mrow></mtd></mtr></mtable></mrow></mrow><mn>,</mn></mrow></math> where m - and σ - are respectively the mean and standard deviation of d - . The arPLS method performs baseline correction as successfully as the AsLS and airPLS methods without a priori peak detection. Indeed, it obtains a more satisfactory baseline in most cases. Similar to the previous method, the weight assignment is still constructed in different forms for both zones (d - and d + ). In addition, it does not consider a specific form of the weights at the points corresponding to different signal peak heights. When dealing with higher and wider peaks in a rapidly changing baseline, the results are no longer reliable. In the next section, we propose a unified method to construct weights for the PLS method by invoking the Bayesian rule. The new algorithm outperforms arPLS when it comes to handling rapidly changing baselines and varying peak heights under different noise levels. Bayesian approach of weight assignment 2 Generally, experimental data, denoted by y i , consist of three components: baseline z i , noise ϵ i , and possible signal s i : <math display="block" id="M13"><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>=</mo><msub><mi>s</mi><mi>i</mi></msub><mo>+</mo><msub><mi>z</mi><mi>i</mi></msub><mo>+</mo><msub><mi>ϵ</mi><mi>i</mi></msub><mn>.</mn></mrow></math> The baseline and noise are commonly present throughout the spectrum, whereas signals are only present in certain parts. Therefore, some data points in different spectral regions have signal contributions, whereas others do not. As an analogy to the good-and-bad data model [ 18 18 ], the likelihood <math id="M14"><mrow><mtext>prob</mtext><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><mi>β</mi><mn>,</mn><mi mathvariant="script">I</mi><mo stretchy="false">)</mo></mrow></math> of the data point i can be considered as a weighted combination of case B i , which contains no signal but only baseline and noise, and case <math id="M15"><mrow><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub></mrow></math> , which contains also a signal, as follows: <math display="block" id="M16"><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><mi>β</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>=</mo><mi>β</mi><mtext> prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>+</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mtext> prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mtext>B</mtext><mi>i</mi></msub><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> where 0 < β < 1 and the nuisance parameter β represents the probability that a datum contains a signal component. All relevant information <math id="M17"><mi mathvariant="script">I</mi></math> includes the nature of the physical situation and knowledge of the experiment. The prior probability of case B i is <math display="block" id="M18"><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mtext>B</mtext><mi>i</mi></msub><mo>|</mo><mi>β</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>=</mo><mn>1</mn><mo>−</mo><mi>β</mi><mn>.</mn></mrow></math> The task is to obtain the appropriate weight using Eq. (5) for each point to achieve a smooth baseline z i . According to the Bayesian theorem [ 19 19 ], the posterior probability for case B i containing no signal is <math display="block" id="M19"><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mtext>B</mtext><mi>i</mi></msub><mo>|</mo><msub><mi>y</mi><mi>i</mi></msub><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mtext>B</mtext><mi>i</mi></msub><mn>,</mn><mi>β</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mtext> prob</mtext><mrow><mo>(</mo><mrow><msub><mtext>B</mtext><mi>i</mi></msub><mo>|</mo><mi>β</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mrow><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><mi>β</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mrow></mfrac><mn>,</mn></mrow></math> where <math id="M20"><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mtext>B</mtext><mi>i</mi></msub><mn>,</mn><mi>β</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mrow></math> is the likelihood function of the data point i in the absence of a signal. Starting with the common assumption of Gaussian noise with zero mean and standard deviation σ , the likelihood function of case B i containing no signal can be approximated by a Gaussian probability density function: <math display="block" id="M21"><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mtext>B</mtext><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msqrt><mrow><mn>2</mn><mi>π</mi></mrow></msqrt><mi>σ</mi></mrow></mfrac><mtext>exp</mtext><mrow><mo>[</mo><mrow><mo>−</mo><mfrac><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow><mrow><mn>2</mn><msup><mi>σ</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><mn>.</mn></mrow></math> When the datum contains a signal contribution s i , the likelihood of <math id="M22"><mrow><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub></mrow></math> is given by <math display="block" id="M23"><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><msub><mi>s</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mrow><msqrt><mrow><mn>2</mn><mi>π</mi></mrow></msqrt><mi>σ</mi></mrow></mfrac><mtext>exp</mtext><mrow><mo>[</mo><mrow><mo>−</mo><mfrac><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub><mo>−</mo><msub><mi>s</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow><mrow><mn>2</mn><msup><mi>σ</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><mn>.</mn></mrow></math> Because s i is unknown before the analysis, it must be marginalized by <math display="block" id="M24"><mtable><mtr><mtd><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>=</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>∞</mi></msubsup></mrow></mstyle><mi>d</mi><msub><mi>s</mi><mi>i</mi></msub><mtext> prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><msub><mi>s</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mtext> prob</mtext><mrow><mo>(</mo><mrow><msub><mi>s</mi><mi>i</mi></msub><mo>|</mo><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mn>.</mn></mtd></mtr></mtable></math> As a convention, s i is always positive in the task of baseline correction. Here, a common scale factor μ is introduced to describe the signal; let μ be the mean value of all signal heights in the spectrum. The maximum entropy principle leads to a prior probability of s i [ 19 19 ] under a Poisson distribution <math display="block" id="M25"><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>s</mi><mi>i</mi></msub><mo>|</mo><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mn>1</mn><mi>μ</mi></mfrac><mtext>exp</mtext><mrow><mo>[</mo><mrow><mo>−</mo><mfrac><mrow><msub><mi>s</mi><mi>i</mi></msub></mrow><mi>μ</mi></mfrac></mrow><mo>]</mo></mrow><mn>.</mn></mrow></math> Substituting Eqs. (15) and (17) into Eqs. (16) for s i , we obtain the likelihood function of a data point containing signal as follows: <math display="block" id="M26"><mtable columnalign="left"><mtr><mtd><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mn>1</mn><mrow><mn>2</mn><mi>μ</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mtext>erf</mtext><mrow><mo>[</mo><mrow><mfrac><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub></mrow><mrow><msqrt><mn>2</mn></msqrt><mi>σ</mi></mrow></mfrac><mo>−</mo><mfrac><mi>σ</mi><mrow><msqrt><mn>2</mn></msqrt><mi>μ</mi></mrow></mfrac></mrow><mo>]</mo></mrow></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mtext>exp</mtext><mrow><mo>[</mo><mrow><mo>−</mo><mfrac><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub></mrow><mi>μ</mi></mfrac><mo>+</mo><mfrac><mrow><msup><mi>σ</mi><mn>2</mn></msup></mrow><mrow><mn>2</mn><msup><mi>μ</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow><mn>.</mn></mtd></mtr></mtable></math> Therefore, the likelihood function in Eq. (11) can be rewritten as follows: <math display="block" id="M27"><mtable><mtr><mtd><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>β</mi><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mi>β</mi><mtext> prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mi>β</mi></mrow><mo>)</mo></mrow><mtext>​</mtext><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mtext>B</mtext><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>=</mo><mi>β</mi><mfrac><mrow><mtext>exp</mtext><mrow><mo>[</mo><mrow><mo>−</mo><mfrac><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub></mrow><mi>μ</mi></mfrac><mo>+</mo><mfrac><mrow><msup><mi>σ</mi><mn>2</mn></msup></mrow><mrow><mn>2</mn><msup><mi>μ</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow></mrow><mrow><mn>2</mn><mi>μ</mi></mrow></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mtext>erf</mtext><mrow><mo>[</mo><mrow><mfrac><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub></mrow><mrow><msqrt><mn>2</mn></msqrt><mi>σ</mi></mrow></mfrac><mo>−</mo><mfrac><mi>σ</mi><mrow><msqrt><mn>2</mn></msqrt><mi>μ</mi></mrow></mfrac></mrow><mo>]</mo></mrow></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>+</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mfrac><mrow><mtext>exp</mtext><mrow><mo>[</mo><mrow><mo>−</mo><mfrac><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow><mrow><mn>2</mn><msup><mi>σ</mi><mn>2</mn></msup></mrow></mfrac></mrow><mo>]</mo></mrow></mrow><mrow><msqrt><mrow><mn>2</mn><mi>π</mi></mrow></msqrt><mi>σ</mi></mrow></mfrac><mo>.</mo></mtd></mtr></mtable></math> Finally, by substituting the likelihood expressed by Eqs. (19) and the prior probability given by Eq. (12) into Eqs. (13), we obtain the posterior probability for a data point i containing no signal: <math display="block" id="M28"><mtable columnalign="left"><mtr><mtd><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mtext>B</mtext><mi>i</mi></msub><mo>|</mo><msub><mi>y</mi><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>β</mi><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mtext>    </mtext><mo>=</mo><mfrac><mrow><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mtext>B</mtext><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>β</mi><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo></mrow><mrow><mi>p</mi><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>β</mi><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mrow></mfrac></mtd></mtr><mtr><mtd><mtext>    </mtext><mfrac><mrow><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mtext> prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mtext>B</mtext><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mrow><mrow><mi>β</mi><mtext>​</mtext><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mrow><mover accent="true"><mtext>B</mtext><mo stretchy="true">¯</mo></mover></mrow><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>+</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mtext> prob</mtext><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>|</mo><msub><mtext>B</mtext><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow></mrow></mfrac></mtd></mtr><mtr><mtd><mtext>    </mtext><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mi>f</mi><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>β</mi><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mo stretchy="false">)</mo></mrow></mfrac><mo>,</mo></mtd></mtr></mtable></math> where <math display="block" id="M29"><mtable columnalign="left"><mtr><mtd><mi>f</mi><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mn>,</mn><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>β</mi><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mi>β</mi><mrow><mn>1</mn><mo>−</mo><mi>β</mi></mrow></mfrac><msqrt><mrow><mfrac><mi>π</mi><mn>2</mn></mfrac></mrow></msqrt><mfrac><mi>σ</mi><mi>μ</mi></mfrac></mtd></mtr><mtr><mtd><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mtext>erf</mtext><mrow><mo>[</mo><mrow><mfrac><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub></mrow><mrow><msqrt><mn>2</mn></msqrt><mi>σ</mi></mrow></mfrac><mo>−</mo><mfrac><mi>σ</mi><mrow><msqrt><mn>2</mn></msqrt><mi>μ</mi></mrow></mfrac></mrow><mo>]</mo></mrow></mrow><mo>)</mo></mrow><mtext>exp</mtext><mrow><mo>[</mo><mrow><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub></mrow><mrow><msqrt><mn>2</mn></msqrt><mi>σ</mi></mrow></mfrac><mo>−</mo><mfrac><mi>σ</mi><mrow><msqrt><mn>2</mn></msqrt><msub><mi>μ</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow><mn>.</mn></mtd></mtr></mtable></math> Bayesian reweighted penalized least squares algorithm 2 Thus far, we have obtained the posterior probability of a data point containing no signal, expressed by Eq. (20). In this study, this probability value is used as the weight for each point in Eq. (5). Thus, the weight of each datum i can be expressed in a unified format: <math display="block" id="M30"><mtable columnalign="left"><mtr><mtd><msub><mi>w</mi><mi>i</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mfrac><mi>β</mi><mrow><mn>1</mn><mo>−</mo><mi>β</mi></mrow></mfrac><msqrt><mrow><mfrac><mi>π</mi><mn>2</mn></mfrac></mrow></msqrt><mi>F</mi><mo stretchy="false">(</mo><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mo stretchy="false">)</mo></mrow></mfrac></mtd></mtr><mtr><mtd><mtext>where   </mtext><mi>F</mi><mo stretchy="false">(</mo><mi>d</mi><mn>,</mn><mi>σ</mi><mn>,</mn><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mi>σ</mi><mi>μ</mi></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mtext>erf</mtext><mrow><mo>[</mo><mrow><mfrac><mi>d</mi><mrow><msqrt><mn>2</mn></msqrt><mi>σ</mi></mrow></mfrac><mo>−</mo><mfrac><mi>σ</mi><mrow><msqrt><mn>2</mn></msqrt><mi>μ</mi></mrow></mfrac></mrow><mo>]</mo></mrow></mrow><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mtext>exp</mtext><mrow><mo>[</mo><mrow><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mi>d</mi><mrow><msqrt><mn>2</mn></msqrt><mi>σ</mi></mrow></mfrac><mo>−</mo><mfrac><mi>σ</mi><mrow><msqrt><mn>2</mn></msqrt><mi>μ</mi></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow></mtd></mtr></mtable></math> In practice, nuisance parameters such as the signal chance β , signal scale μ , and Gaussian noise standard deviation σ for each data point must be properly assigned to use the algorithm. Similar to existing algorithms for PLS-based baseline estimation, we also divide the values of subtraction d = y - z into a positive part d + and negative part d - in the iterations of our algorithm. Then, an approximation of σ and μ can be obtained using different parts of the subtraction: the standard deviation σ can be estimated from the quadratic mean of d - , whereas the signal scale μ can be estimated from the mean value of d + . As mentioned in the previous section, β reflects the probability of signal contribution in the data. According to Eq. (21), if the signal-to-noise ratio (SNR) μ 2 / σ 2 [ 20 20 ] and signal scale μ are constant, the calculated weight w i can be expressed as a function of the subtraction d i for different values of β , as shown in Fig. 1 Fig. 1 . The weight functions w i corresponding to different β values are clearly different. Thus, a rational value for β is required. Given that the proportion and distribution of signals in the data are unknown, it is convenient to use the average probability of no signal appearance in the total data set as a uniform β factor assigned to each data point, as suggested by the concept of good-and-bad model [ 18 18 ]. In this study, the average probability of no signals in the total data set can be defined as follows: <math display="block" id="M32"><mrow><mo>〈</mo><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mtext>B</mtext><mi>i</mi></msub><mo>|</mo><msub><mi>y</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mo>〉</mo><mo>≡</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover></mstyle><mtext>prob</mtext><mrow><mo>(</mo><mrow><msub><mtext>B</mtext><mi>i</mi></msub><mo>|</mo><msub><mi>y</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi></mrow><mo>)</mo></mrow><mn>.</mn></mrow></math> Note that β is updated in each iteration and the iteration terminates when β approaches <math id="M33"><mrow><mn>1</mn><mo>−</mo><mo>〈</mo><mtext>prob</mtext><mo stretchy="false">(</mo><msub><mi>B</mi><mi>i</mi></msub><mo>|</mo><msub><mi>y</mi><mi>i</mi></msub><mn>,</mn><mi>σ</mi><mn>,</mn><mi mathvariant="script">I</mi><mo stretchy="false">)</mo><mo>〉</mo></mrow></math> . All in all, the algorithm flow of the Bayesian reweighted penalized least squares (BrPLS) method is shown in Algorithm 1 Algorithm 1 . 5 1 Algorithm 1 BrPLS method 5 1 Input: measured spectra y, smoothness parameter λ , termination criteria ratio , maximum iteration maxIter 5 1 Output: estimated baseline z 5 1 procedure BrPLS <math id="M34"><mrow><mi>y</mi><mn>,</mn><mi>λ</mi><mn>,</mn><mi>r</mi><mi>a</mi><mi>t</mi><mi>i</mi><mi>o</mi></mrow></math> 1 1   4 1 <math id="M35"><mrow><mi>z</mi><mo>←</mo><mi>y</mi></mrow></math> , <math id="M36"><mrow><mi>β</mi><mo>←</mo><mn>0.5</mn></mrow></math> 1 1   4 1 <math id="M37"><mrow><mi>H</mi><mo>=</mo><mi>λ</mi><msup><mstyle mathvariant="bold"><mtext>D</mtext></mstyle><mi>T</mi></msup><mstyle mathvariant="bold"><mtext>D</mtext></mstyle></mrow></math> 1 1   4 1 with D in Eq. (2) 1 1   4 1 w =[1, 1, …, 1] 1 1   4 1 for i =0 to maxIter do 1 1   1 1   3 1 for i =0 to maxIter do 1 1   1 1   1 1   2 1 make W a diagonal matrix as W i,i = w i <math id="M38"><mrow><mi>z</mi><mo>'</mo><mo>=</mo><mo stretchy="false">(</mo><mstyle mathvariant="bold"><mtext>W</mtext></mstyle><mo>+</mo><mi>H</mi><msup><mo stretchy="false">)</mo><mrow><mo>−</mo><mn>1</mn></mrow></msup><mstyle mathvariant="bold"><mtext>W</mtext></mstyle><mi>y</mi></mrow></math> 1 1   1 1   1 1   2 1 d = y - z’ 1 1   1 1   1 1   2 1 μ = mean of positive part of d 1 1   1 1   1 1   2 1 σ = quadratic mean of negative part of d <math id="M39"><mrow><mi>F</mi><mo>=</mo><mfrac><mi>σ</mi><mi>μ</mi></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mtext>erf</mtext><mrow><mo>(</mo><mrow><mfrac><mi>d</mi><mrow><msqrt><mn>2</mn></msqrt><mi>σ</mi></mrow></mfrac><mo>−</mo><mfrac><mi>σ</mi><mrow><msqrt><mn>2</mn></msqrt><mi>μ</mi></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><msup><mtext>e</mtext><mrow><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mi>d</mi><mrow><msqrt><mn>2</mn></msqrt><mi>σ</mi></mrow></mfrac><mo>−</mo><mfrac><mi>σ</mi><mrow><msqrt><mn>2</mn></msqrt><mi>μ</mi></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></msup></mrow></math> 1 1   1 1   1 1   2 1 <math id="M40"><mrow><mi>w</mi><mo>=</mo><mn>1</mn><mo>/</mo><mo stretchy="false">(</mo><mn>1</mn><mo>+</mo><mfrac><mi>β</mi><mrow><mn>1</mn><mo>−</mo><mi>β</mi></mrow></mfrac><msqrt><mrow><mfrac><mi>π</mi><mn>2</mn></mfrac></mrow></msqrt><mi>F</mi><mo stretchy="false">)</mo></mrow></math> 1 1   1 1   1 1   2 1 if <math id="M41"><mrow><mo>|</mo><mi>z</mi><mo>−</mo><mi>z</mi><mo>'</mo><mo>|</mo><mo>/</mo><mo>|</mo><mi>z</mi><mo>'</mo><mo>|</mo><mo><</mo><mi>r</mi><mi>a</mi><mi>t</mi><mi>i</mi><mi>o</mi></mrow></math> then 1 1   1 1   1 1   1 1   1 1 break 1 1   1 1   1 1   2 1 else 1 1   1 1   1 1   1 1   1 1 <math id="M42"><mrow><mi>z</mi><mo>←</mo><mi>z</mi><mo>'</mo></mrow></math> 1 1   1 1   1 1   2 1 end if 1 1   1 1   3 1 end for 1 1   1 1   3 1 if <math id="M43"><mrow><mo>|</mo><mi>β</mi><mo>+</mo></mrow></math> mean of <math id="M44"><mrow><mi>w</mi><mo>−</mo><mn>1</mn><mo>|</mo><mo><</mo><mi>r</mi><mi>a</mi><mi>t</mi><mi>i</mi><mi>o</mi></mrow></math> then 1 1   1 1   1 1   2 1 Break 1 1   1 1   3 1 else 1 1   1 1   1 1   2 1 <math id="M45"><mrow><mi>β</mi><mo>←</mo><mn>1</mn><mo>−</mo></mrow></math> mean of w 1 1   1 1   3 1 end if 1 1   4 1 end for 5 1 end procedure Algorithm performance 1 In this section, the performance of the BrPLS method is evaluated along with previous PLS-based methods. The advantages and limitations of these algorithms were revealed by their application to the analysis of synthetic spectra with known baselines and real-world Schottky spectra. All performance tests were performed with codes programmed in Python using the Numpy and Scipy libraries. Reliability tests 2 We used synthetic spectra to test the reliability of different PLS-based baseline-estimation algorithms. Synthetic spectra were generated by adding signal components and random noise to a pre-defined baseline as follows: <math display="block" id="M46"><mrow><mi>d</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><mi>s</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>+</mo><mi>b</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>+</mo><mi>n</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mn>,</mn></mrow></math> where s ( x ), b ( x ), n ( x ), and d ( x ) denote the signal, baseline, noise, and resulting synthetic data, respectively. The synthetic signals were eight Gaussian peaks whose intensities, positions, and widths are listed in Table 1 Table 1 . Typical baseline patterns in PIXE [ 21 21 - 23 23 21-23 ] and Schottky spectra [ 24 24 ] contain linear, concave, or convex regions. Therefore, we chose two types of baseline for the reliability tests: linear and sinusoidal. Random noise was generated by using a random number from a Gaussian distribution. The SNR in terms of energy was defined in the absence of a baseline, according to the following equation [ 20 20 ]: <math display="block" id="M47"><mrow><mi>S</mi><mi>N</mi><mi>R</mi><mo>=</mo><mn>10</mn><mtext> </mtext><msub><mrow><mtext>log</mtext></mrow><mrow><mn>10</mn></mrow></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mtext>s</mtext></msub><mo>/</mo><msub><mi>E</mi><mtext>n</mtext></msub><mo stretchy="false">)</mo><mn>,</mn></mrow></math> where E s refers to the expected energy value of the signal and E n refers to the expected noise energy value in the spectrum. The SNR of the low-noise spectrum was set to 40 dB and that of the high-noise spectrum was set to 20 dB. 1 1 Peak No. 1 1 Intensity 1 1 Position 1 1 Width 1 1 1 1 1 10 1 1 9 1 1 0.7 1 1 2 1 1 2 1 1 20 1 1 0.3 1 1 3 1 1 5 1 1 22 1 1 0.1 1 1 4 1 1 15 1 1 40 1 1 0.2 1 1 5 1 1 3 1 1 49 1 1 0.2 1 1 6 1 1 2 1 1 52 1 1 0.1 1 1 7 1 1 20 1 1 60 1 1 0.6 1 1 8 1 1 14 1 1 70 1 1 0.5 The maximum number of iterations was set to 50 for the AsSL, airPLS, arPLS, and BrPLS methods. The early termination criterion for AsSL, arPLS, and BrPLS was a standard deviation of the difference between two estimated baselines in nearby iterations less than 10 -6 . For airPLS, Eq. (8) was used as the early termination criterion. To start the tests of these four methods, an appropriate smoothness parameter λ in Eq. 3 should be determined. Recall that a larger λ leads to a smoother estimated baseline; if λ is too large or too small, the estimated baseline will not match the curvature of the baseline or will be highly distorted by high peaks. In addition, the sensitivities of the four methods to the same λ are different. To obtain the best λ for each method, the parameter space of λ was scanned from 10 2 to 10 12 in 10-base logarithmic steps of 0.1. By the end of the iterations, the algorithm reaches a final set of z i ; then, the root mean square error (RMSE) can be calculated as follows: <math display="block" id="M48"><mrow><mtext>RMSE</mtext><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>=</mo><msqrt><mrow><mfrac><mrow><msup><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><mrow><mo>(</mo><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow></mstyle></mrow><mn>2</mn></msup></mrow><mi>N</mi></mfrac></mrow></msqrt><mn>.</mn></mrow></math> This RMSE expression reflects the quadratic mean of the estimated signal. If the baseline shape is linear, as λ increases, the RMSE curve first increases rapidly and then converges, as shown in Figs. 2a, c Figs. 2a, c . In this case, λ values can be chosen at any location where the RMSE is stable. For example, the stable RMSE( λ ) region shown in Fig. 2a Fig. 2a corresponds to values of λ ranging from 10 8.1 to 10 12 and that shown in Fig. 2c Fig. 2c corresponds to values of λ ranging from 10 8 to 10 12 . If the baseline contains depression or protrusion regions, the RMSE as a function of λ can become more complicated, as shown in Figs. 2b, d Figs. 2b, d . A suitable λ value is always chosen at the stable plateau of RMSE( λ ), where the estimated baseline is close to the actual baseline. For the specific case of a sinusoidal baseline, as in Fig. 2b Fig. 2b , the optimal values of λ that can be selected range from 10 6.5 to 10 8 for the AsLS and airPLS methods, from 10 6.5 to 10 10.2 for the arPLS method, and from 10 6.5 to 10 9.9 for the BrPLS method. Similarly, for the case shown in Fig. 2d Fig. 2d , the optimal λ values that can be selected range from 10 7 to 10 8.6 for the AsLS and airPLS methods, and from 10 4.5 to 10 8.6 for the arPLS and BrPLS methods. Generally, the expected baseline must be as smooth as possible. Therefore, the maximum optional λ was used in the tests for each method, as listed in Tab. 2 Tab. 2 . 2 1 SNR 3 1 20 dB 3 1 40 dB 1 1 Baseline 1 1 Method 1 1 log( λ ) 1 1 R 2 1 1 RMSE 1 1 log( λ ) 1 1 R 2 1 1 RMSE 1 1 Linear 1 1 AsSL 1 1 12.0 1 1 -10.8717 1 1 3.0775 1 1 12.0 1 1 0.9418 1 1 2.8609 1 1   1 1 airPLS 1 1 12.0 1 1 -8.7724 1 1 3.0539 1 1 12.0 1 1 0.9827 1 1 2.8520 1 1   1 1 arPLS 1 1 12.0 1 1 0.9890 1 1 2.8746 1 1 12.0 1 1 0.9998 1 1 2.8521 1 1   1 1 BrPLS 1 1 12.0 1 1 0.9931 1 1 2.8731 1 1 12.0 1 1 0.9998 1 1 2.8511 1 1 Sinusoidal 1 1 AsSL 1 1 8.0 1 1 -11.6166 1 1 3.0141 1 1 8.6 1 1 0.6222 1 1 2.8135 1 1   1 1 airPLS 1 1 8.0 1 1 -9.0078 1 1 3.0225 1 1 8.6 1 1 0.8452 1 1 2.8374 1 1   1 1 arPLS 1 1 10.2 1 1 0.9791 1 1 2.8710 1 1 8.6 1 1 0.9985 1 1 2.8540 1 1   1 1 BrPLS 1 1 9.9 1 1 0.9815 1 1 2.8706 1 1 8.6 1 1 0.9990 1 1 2.8518 To quantitatively assess the performance of each method, we introduced the goodness of the baseline estimation factor R 2 [ 25 25 ]: <math display="block" id="M49"><mrow><msup><mi>R</mi><mn>2</mn></msup><mo>=</mo><mn>1</mn><mo>−</mo><mfrac><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>b</mi><mi>i</mi></msub><mo>−</mo><msub><mi>z</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></mstyle></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msup><mrow><mrow><mo>(</mo><mrow><mo>〈</mo><mi>b</mi><mo>〉</mo><mo>−</mo><msub><mi>b</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></mstyle></mrow></mfrac><mn>,</mn></mrow></math> where b i and <math id="M50"><mrow><mo>〈</mo><mi>b</mi><mo>〉</mo></mrow></math> denote the synthetic baseline and mean of all the baseline values, respectively. Ideally, the estimated baseline values match the synthetic values, i.e., R 2 = 1. The results are presented in Fig. 3 Fig. 3 and Table 2 Table 2 . The BrPLS and arPLS methods have R 2 values close to 1, but the R 2 values of the other two methods deviate from 1. The results suggest that the AsSL and airPLS methods are inadequate for spectral-baseline estimation. The differences in R 2 between the BrPLS and arPLS methods are less than 0.5%, indicating no significant differences between both methods in these analysis tests. In Figs. 3(c) and 3(d) Figs. 3(c) and 3(d) , the estimated baseline near the data boundary is slightly distorted within the 1- σ error interval of the noise owing to the limitations of the PLS-based methods, where only the 2nd-order difference calculation ( m = 2) of the baseline was used in Eqs. (1) and (2). In any case, the differences between the estimated and synthetic baselines of both the arPLS and BrPLS methods are within the 1- σ interval of the synthetic Gaussian noise. This means that the estimated baselines of both methods are reliable within the 1- σ error interval. Note also that the BrPLS method has a relatively better ability to handle the signal peaks in the convex regions of the baseline, as indicated by the arrows in Fig. 3(d) Fig. 3(d) . This will be further evaluated in the next section. Comparison between the arPLS and BrPLS methods 2 According to the results presented in the previous section, the BrPLS and arPLS methods outperform the other two methods regardless of the baseline shape. Moreover, the BrPLS method has a slight advantage over the arPLS method when dealing with cases in which the signal peaks are located in the convex regions of the baseline. In this section, we discuss and compare the advantages and limitations of both algorithms by applying them to the analysis of other synthetic spectra. The difference between the BrPLS and arPLS methods increases in regions where the baseline has protruding parts. To estimate such a baseline, the most challenging locations are typically those at which the signal peaks lie near the highest point of the convex part and at the inflection points of the rising and falling edges. Thus, a typical convoluted Landau-Gaussian function was used to generate the synthetic baseline in the synthetic spectrum, and three Gaussian signal peaks were placed in the aforementioned parts of the baselines. The intensities and widths of the signal peaks are listed in Table 3 Table 3 . The SNR was set to 40 dB. The maximum number of iterations and early termination criteria for both methods were the same as those in the previous tests. Similarly, the λ values were selected in the same manner as in the previous tests. 1 1 Picture label 2 1 Peak 1 2 1 Peak 2 2 1 Peak 3 3 1 arPLS 3 1 BrPLS 1 1   1 1 Intensity 1 1 Width 1 1 Intensity 1 1 Width 1 1 Intensity 1 1 Width 1 1 log( λ ) 1 1 R 2 1 1 RMSE 1 1 log( λ ) 1 1 R 2 1 1 RMSE 1 1 a 1 1 3 1 1 0.1 1 1 5 1 1 0.2 1 1 10 1 1 0.1 1 1 7.5 1 1 0.9994 1 1 0.5309 1 1 7.5 1 1 0.9997 1 1 0.5308 1 1 b 1 1 3 1 1 0.3 1 1 5 1 1 0.4 1 1 10 1 1 0.2 1 1 7.5 1 1 0.9986 1 1 0.7609 1 1 7.5 1 1 0.9987 1 1 0.7604 1 1 c 1 1 3 1 1 0.6 1 1 5 1 1 0.8 1 1 10 1 1 0.4 1 1 7.0 1 1 0.9904 1 1 01.0795 1 1 7.0 1 1 0.9910 1 1 1.0745 1 1 d 1 1 6 1 1 0.1 1 1 10 1 1 0.2 1 1 20 1 1 0.1 1 1 7.7 1 1 0.9990 1 1 1.0617 1 1 7.7 1 1 0.9993 1 1 1.0615 1 1 e 1 1 6 1 1 0.3 1 1 10 1 1 0.4 1 1 20 1 1 0.2 1 1 7.8 1 1 0.9954 1 1 1.5221 1 1 7.8 1 1 0.9964 1 1 1.5209 1 1 f 1 1 6 1 1 0.6 1 1 10 1 1 0.8 1 1 20 1 1 0.4 1 1 7.5 1 1 0.9721 1 1 2.1592 1 1 7.5 1 1 0.9771 1 1 2.1526 1 1 g 1 1 12 1 1 0.1 1 1 20 1 1 0.2 1 1 40 1 1 0.1 1 1 8.1 1 1 0.9938 1 1 2.1235 1 1 8.1 1 1 0.9963 1 1 2.1230 1 1 h 1 1 12 1 1 0.3 1 1 20 1 1 0.4 1 1 40 1 1 0.2 1 1 8.1 1 1 0.9878 1 1 3.0452 1 1 8.1 1 1 0.9903 1 1 3.0430 1 1 i 1 1 12 1 1 0.6 1 1 20 1 1 0.8 1 1 40 1 1 0.4 1 1 8.1 1 1 -1.3622 1 1 4.3993 1 1 8.1 1 1 0.9282 1 1 4.3130 The results are presented in Fig. 4 Fig. 4 and Table 3 Table 3 . Both methods produce reasonable baseline results. However, when the signal peaks are located at the rising edge of the baseline, the two methods yield unsatisfactory estimations. This is another inherent limitation of PLS-based methods as they strive to balance the smoothness of the baseline and fidelity to the original spectrum. Signal peaks located at the rapidly changing parts of the baseline represent the case in which the two factors are conflicting. As the peak width and height increase, the baseline obtained by the arPLS method further deviates from the actual value when the signal peak is located at the baseline protrusion regions. In contrast, the BrPLS method yields significantly better results, as shown in Fig. 4(i) Fig. 4(i) . This is because the BrPLS method considers the heights of different signal peaks and assigns different weight values to the parameter μ in Eq. (21). Conversely, the weight function in Eq. (9) for the arPLS method is a symmetric sigmoid curve that does not consider signal peak characteristics. Application to real-world Schottky spectra 2 Schottky spectroscopy [ 26 26 ] is an important method for beam diagnostics in heavy-ion accelerator facilities. It is routinely used in precision mass and lifetime measurement experiments in nuclear physics studies [ 27 27 ]. The corresponding Schottky spectrum reflects the power density of signals induced by circulating ions in a storage ring. The periodic motion of an ion is represented by the resonance peaks at each harmonic of the revolution frequency in the measured spectrum [ 28 28 ]. Accurate determination of the peak position and area in the power density spectrum is essential for data analysis. The revolution frequency, which is proportional to the mass-to-charge ratio of the ion, can be derived from the peak position. The peak area can be derived from the baseline-corrected spectrum. It is proportional to the number of corresponding ions and can be used to determine the lifetime of the ions in the storage ring. In the Experimental Storage Ring at Darmstadt [ 29 29 ] and experimental Cooler Storage Ring (CSRe) at Lanzhou [ 30 30 ], ions circulate with frequencies ranging from several hundred kHz to a few MHz. The bandwidth of the measured Schottky spectrum is typically less than 2 MHz [ 31 31 - 33 33 31-33 ]. The baseline of this type of spectrum can be approximated as a straight line and, in some cases, can be corrected using simple 1st-order polynomial fittings [ 34 34 , 35 35 ]. However, the baseline becomes complicated for the full spectrum, which contains several harmonics of the full momentum acceptance of the storage ring. Owing to complex signal processing in hardware [ 24 24 ], the baseline appearing in this spectrum does not always have a theoretical Lorentzian shape [ 36 36 ]. Different types of baseline distortions may occur during the actual signal processing. Thus, a method that can effectively estimate and correct Lorentzian-type baselines is required. The newly developed BrPLS method offers a new option in this direction. The power spectral density of Schottky noise is multiplicative [ 37 37 ] and our Bayesian assumption is based on the additive components in the data, as expressed in Eq. (10). Therefore, a logarithmic transformation of the power densities in the original Schottky spectra was applied first. The processed spectra were then adjusted to non-negative values by subtracting the minimum values. Experimental data were obtained from the measurements of 86 Kr projectile fragments [ 38 38 - 40 40 38-40 ] at CSRe in 2018. Each frame in the Schottky spectrum had a data acquisition time of 2.73 s. The statistical uncertainty of the amplitude at a given frequency point was approximately 1.4% [ 41 41 ]. Particle identification in the spectrum was performed according to [ 26 26 ]. During baseline estimation, steps similar to those of the previously described synthetic data treatment were applied. The parameters and other settings of the arPLS and BrPLS methods were the same as those described in the previous section. The results of a simple baseline fit with a local 1st-order polynomial (fitting frequency range between ± 30 kHz of the signal peak center) are also presented for comparison. The peak position and area under the frequency peak were determined by fitting Gaussian functions to the peaks after baseline corrections. The baseline estimation results for a single frame are shown in Fig. 5 Fig. 5 . Consistent with the previous section, the BrPLS method could effectively handle the presence of peaks in the baseline protrusion regions. In contrast, the arPLS method underestimated the baseline to the extent of introducing more spurious signal peaks, as indicated by the blue line in Fig. 6 Fig. 6 . In Fig. 5 Fig. 5 , five signal peaks are marked for further discussion. Peaks 1 and 4 belong to the same ion species 82 Se 34+ but differ in the harmonic number of the ion revolution frequency. Peak 1 is located in the flat baseline part whereas peak 4 is located in the middle of the falling edge of the baseline. Peaks 2 and 5 belong to the same ion species 84 Br 35+ and are located in the relatively flat part of the baseline. The small ion peak 3 belongs to 74 Ga 31+ ions and is located near the top of the baseline-protrusion region. The estimated revolution frequencies were consistent among the three methods within the 1 σ error range, as shown in Fig. 7 Fig. 7 . The central values vary slightly as a function of the peak position and method, particularly when the peak is located in the convex part of the baseline. The peak area as a function of the elapsed time was used to determine the ion storage lifetime after ion injection into the ring. More spectra shown in the top graph of Fig. 5 Fig. 5 were analyzed for this task, and the differences in the performance of the three methods became more apparent. Ideally, the ion lifetimes estimated from the ion peaks at different harmonic numbers were expected to be the same. As shown in Fig. 8 Fig. 8 , the three methods achieved this goal. However, arPLS incurred a relatively larger error for ion peak 4, which is located in the baseline protrusion region. In this region, the ion peak areas obtained by the arPLS method were somewhat overfitted, resulting in an underestimation of the baseline. Some blue points obtained using this method and depicted in Fig. 8(b) Fig. 8(b) deviate from the overall trend, leading to larger uncertainty in the corresponding lifetime, as shown in Fig. 8(c) Fig. 8(c) . The overall difference between the results of the polynomial fitting and BrPLS methods is negligible. It is expedient to use a local linear fit to the baseline to estimate the area of individual peaks when the baseline shape of the full spectrum is very complex. However, this method relies heavily on the degree of linearity of the baseline within the selected local range. Additionally, to estimate the neutral atomic lifetime of the corresponding nuclide of an ion from its storage lifetime in the ring, it is necessary to use other stored ions in the same spectrum as a reference. A local linear fit to individual peaks is challenging to guarantee the consistency of the baseline estimate for all the peaks located in different parts of the spectrum. In contrast, the results obtained by the newly developed BrPLS method are generally more stable and reliable than those of the other two methods. Summary and outlooks 1 This paper reports the BrPLS method for baseline estimation, which provides a unified and reliable weight assignment based on the Bayesian theorem. This method is an upgrade to the existing PLS method and is implemented as an iteration algorithm. The new method was applied to the treatment of both synthetic and experimental spectra. It outperforms existing PLS-based baseline estimation methods and yields relatively better results, especially when the signal peaks lie on the convex part of the baseline in the presence of noise. An open-source program implementing the BrPLS method can be downloaded from Github [ 43 43 ], providing an alternative baseline correction method for spectrum analysis. Further developments may include the following: (1) using Poisson noise instead of Gaussian noise ( ϵ i ) in Eq. (10); (2) improving the reliability of the baseline estimations using datum-specific scale factors ( μ i ) for each bin instead of an average μ in Eq. (17); (3) similarly, assigning datum-specific ( β i ) instead of an average β in Eq. (11) to the probability that the datum contains signal components. References F. Gan , G. H. Ruan , J. Y. Mo , Baseline correction by improved iterative polynomial fitting with automatic threshold . Chemom. Intell. Lab. Syst. 82 ( 1–2 ), 59–65 ( 2006 ). doi: 10.1016/j.chemolab.2005.08.009 http://doi.org/10.1016/j.chemolab.2005.08.009 Baseline correction by improved iterative polynomial fitting with automatic threshold Z.M. Zhang , S. Chen , Y.Z. Liang et al. , An intelligent background-correction algorithm for highly fluorescent samples in Raman spectroscopy: Background-correction algorithm for highly fluorescent samples in Raman spectroscopy . J. Raman Spectrosc. 41 ( 6 ), 659–669 ( 2009 ). doi: 10.1002/jrs.2500 http://doi.org/10.1002/jrs.2500 An intelligent background-correction algorithm for highly fluorescent samples in Raman spectroscopy: Background-correction algorithm for highly fluorescent samples in Raman spectroscopy A. Jirasek , G. Schulze , M.M.L. Yu et al. , Accuracy and Precision of Manual Baseline Determination . Appl. Spectrosc. 58 ( 12 ), 1488 - 1499 ( 2004 ). doi: 10.1366/0003702042641236 http://doi.org/10.1366/0003702042641236 Accuracy and Precision of Manual Baseline Determination W. von der Linden , V. Dose , L. Padayachee et al. , Signal and background separation . Phys. Rev. E. 59 ( 6 ), 6527–6534 ( 1999 ). doi: 10.1103/physreve.59.6527 http://doi.org/10.1103/physreve.59.6527 Signal and background separation M. A. Kneen , H. J. Annegarn , Algorithm for fitting XRF, SEM and PIXE X-ray spectra backgrounds . Nucl. Instr. Meth. B 109/110 , 209–213 ( 1996 ). doi: 10.1016/0168-583X(95)00908-6 http://doi.org/10.1016/0168-583X(95)00908-6 Algorithm for fitting XRF, SEM and PIXE X-ray spectra backgrounds X.G. Ma , Z.X. Zhang , Application of wavelet transform to background correction in inductively coupled plasma atomic emission spectrometry . Anal. Chim. Acta 485 ( 2 ), 233 - 239 ( 2003 ). doi: 10.1016/S0003-2670(03)00395-7 http://doi.org/10.1016/S0003-2670(03)00395-7 Application of wavelet transform to background correction in inductively coupled plasma atomic emission spectrometry M. Morháč , V. Matoušek , Peak Clipping Algorithms for Background Estimation in Spectroscopic Data . Appl. Spectrosc. 62 ( 1 ), 91–206 ( 2008 ). doi: 10.1366/000370208783412762 http://doi.org/10.1366/000370208783412762 Peak Clipping Algorithms for Background Estimation in Spectroscopic Data C. G. Ryan , E. Clayton , W. L. Griffin et al. , SNIP, a statistics-sensitive background treatment for the quantitative analysis of PIXE spectra in geoscience applications . Nucl. Instr. Meth. B 34 , 396–402 ( 2013 ). doi: 10.1016/0168-583X(88)90063-8 http://doi.org/10.1016/0168-583X(88)90063-8 SNIP, a statistics-sensitive background treatment for the quantitative analysis of PIXE spectra in geoscience applications J. C. Lindon , G. E. Tranter , D. W. Koppenaal , Encyclopedia of spectroscopy and spectrometry , 2nd edn . ( Elservier Academic Press , Oxford , 2010 ). Encyclopedia of spectroscopy and spectrometry E. T. Whittaker , On a new method of graduation . P. Edinburgh Math. Soc. 41 , 63–75 ( 1922 ). doi: 10.1017/S0013091500077853 http://doi.org/10.1017/S0013091500077853 On a new method of graduation R. J. Hodrick , E. C. Prescott , Postwar U.S. business cycles: An empirical investigation . J. Money, Credit Banking 29 ( 1 ), 1–16 ( 1997 ). doi: 10.2307/2953682 http://doi.org/10.2307/2953682 Postwar U.S. business cycles: An empirical investigation P. H. C. Eilers , A perfect smoother . Anal. Chem. 75 ( 14 ), 3631 - 3636 ( 2003 ). doi: 10.1021/ac034173t http://doi.org/10.1021/ac034173t A perfect smoother P. H. C. Eilers , V. Pesendorfer , R. Bonifacio , in Proc. 2017 9th International Workshop on the Analysis of Multitemporal Remote Sensing Images (MultiTemp) , Bruges , June 2017 . Proc. 2017 9th International Workshop on the Analysis of Multitemporal Remote Sensing Images (MultiTemp) P. H. C. Eilers , Parametric time warping . Anal. Chem. 76 ( 2 ), 404–411 ( 2004 ). doi: 10.1021/ac034800e http://doi.org/10.1021/ac034800e Parametric time warping Z. M. Zhang , S. Chen , Y. Z. Liang , Baseline correction using adaptive iteratively reweighted penalized least squares . Analyst 135 ( 5 ), 1138–1146 ( 2010 ). doi: 10.1039/b922045c http://doi.org/10.1039/b922045c Baseline correction using adaptive iteratively reweighted penalized least squares S. J. Baek , A. Park , Y. J. Anh , J. Choo , Baseline correction using asymmetrically reweighted penalized least squares smoothing . Analyst 140 ( 1 ), 250–257 ( 2015 ). doi: 10.1039/C4AN01061B http://doi.org/10.1039/C4AN01061B Baseline correction using asymmetrically reweighted penalized least squares smoothing W. K. Newey , J. L. Powell , Asymmetric Least Squares Estimation and Testing . Econometrica 55 ( 4 ), 819–847 ( 1987 ). doi: 10.2307/1911031 http://doi.org/10.2307/1911031 Asymmetric Least Squares Estimation and Testing D. S. Sivia , in MAXENT96 - Proceedings of the Maximum Entropy Conference , Port Elizabeth , 1996 , ed. by M. Sears , V. Nedeljkovic , N. E. Pendock , and S. Sibisi , p. 131 MAXENT96 - Proceedings of the Maximum Entropy Conference D. S. Sivia , J. Skilling , Data analysis: a Bayesian tutorial , 2nd edn . ( Oxford University Press , Oxford , 2006 ), pp. 5 – 13 Data analysis: a Bayesian tutorial B. Carlson , P. B. Crilly , Communication systems: An introduction to signals and noise in electrical communication , 5th edn . ( McGraw-Hill Higher Education , Boston , 2010 ), pp. 423 – 426 Communication systems: An introduction to signals and noise in electrical communication G. C.-Y. Chan , in Encyclopedia of Analytical Science , ( 3rd edn ) by P. Worsfold , C. Poole , A. Townshend , M. Miró ( Elsevier , Amsterdam Oxford Cambridge , 2019 ), pp. 194 – 299 doi: 10.1016/B978-0-12-409547-2.14487-7 http://doi.org/10.1016/B978-0-12-409547-2.14487-7 Encyclopedia of Analytical Science H.W. Wang , G.T. Fan , L.X. Liu et al. , Commissioning of laser electron gamma beamline SLEGS at SSRF . Nucl. Sci, Tech. 33 , 87 ( 2022 ). doi: 10.1007/s41365-022-01076-0 http://doi.org/10.1007/s41365-022-01076-0 Commissioning of laser electron gamma beamline SLEGS at SSRF R.M.J. Li , S.K. Liu , S. T. Lin et al. , Identification of anomalous fast bulk events in a p-type point-contact germanium detector . Nucl. Sci, Tech. 33 , 57 ( 2022 ). doi: 10.1007/s41365-022-01041-x http://doi.org/10.1007/s41365-022-01041-x Identification of anomalous fast bulk events in a p-type point-contact germanium detector C. Trageser , C. Brandau , C. Kozhuharov et al. , A new data acquisition system for Schottky signals in atomic physics experiments at GSI’s and FAIR’s storage rings . Phys. Scr. T166 , 014062 ( 2015 ). doi: 10.1088/0031-8949/2015/T166/014062 http://doi.org/10.1088/0031-8949/2015/T166/014062 A new data acquisition system for Schottky signals in atomic physics experiments at GSI’s and FAIR’s storage rings D. Chicco , M. J. Warrens , G. Jurman , The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation . Peer. J. Comput Sci. ( 2021 ). doi: 10.7717/peerj-cs.623 http://doi.org/10.7717/peerj-cs.623 The coefficient of determination R-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation Yu. A. Litvinov , H. Geissel , T. Radon et al. , Mass measurement of cooled neutron-deficient bismuth projectile fragments with time-resolved Schottky mass spectrometry at the FRS-ESR facility . Nucl. Phys. A 756 ( 1–2 ), 3–38 ( 2005 ). doi: 10.1016/j.nuclphysa.2005.03.015 http://doi.org/10.1016/j.nuclphysa.2005.03.015 Mass measurement of cooled neutron-deficient bismuth projectile fragments with time-resolved Schottky mass spectrometry at the FRS-ESR facility Yu. A. Litvinov . S. Bishop , K. Blaum et al. , Nuclear physics experiments with ion storage rings . Nucl. Instr. Meth. B 317 , 603–616 ( 2013 ). doi: 10.1016/j.nimb.2013.07.025 http://doi.org/10.1016/j.nimb.2013.07.025 Nuclear physics experiments with ion storage rings S. Chattopadhyay , Some fundamental aspects of fluctuations and coherence in charged-particle beams in storage rings , AIP Conf. Proc 127 , 467 ( 1985 ). doi: 10.1063/1.35175 http://doi.org/10.1063/1.35175 Some fundamental aspects of fluctuations and coherence in charged-particle beams in storage rings B. Franzke , The heavy ion storage and cooler ring project ESR at GSI , Nucl. Instr. Meth. B 24–25 , 18–25 ( 1987 ). doi: 10.1016/0168-583X(87)90583-0 http://doi.org/10.1016/0168-583X(87)90583-0 The heavy ion storage and cooler ring project ESR at GSI W.L. Zhan , J.W. Xia , B.W. Wei et al. HIRFL-CSR project , AIP Conf. Proc 600 , 175 ( 2001 ). doi: 10.1063/1.1435229 http://doi.org/10.1063/1.1435229 HIRFL-CSR project B. Franzke , H. Geissel , G. Münzenberg , Mass and lifetime measurements of exotic nuclei in storage rings , Mass Spec. Rev. 27 , 428–469 ( 2008 ). doi: 10.1002/mas.20173 http://doi.org/10.1002/mas.20173 Mass and lifetime measurements of exotic nuclei in storage rings B. Bosch , Yu. A. Litvinov , T. Stöhlker , Nuclear physics with unstable ions at storage rings , Prog. Part. Nucl. Phys. 73 , 84–140 ( 2013 ). doi: 10.1016/j.ppnp.2013.07.002 http://doi.org/10.1016/j.ppnp.2013.07.002 Nuclear physics with unstable ions at storage rings Yu. A. Litvinov , F. Bosch , Beta decay of highly charged ions , Rep. Prog. Phys 74 ( 1 ), 016301 ( 2011 ). doi: 10.1088/0034-4885/74/1/016301 http://doi.org/10.1088/0034-4885/74/1/016301 Beta decay of highly charged ions Yu. A. Litvinov , Ph.D. thesis, der Justus-Liebig Universität Gieflen , 2003 . X. L. Tu , X. C. Chen , J. T. Zhang et al. First application of combined isochronous and Schottky mass spectrometry: Half-lives of fully ionized 49 CR 24+ and 53 Fe 26+ atoms , Phys. Rev. C 97 , 014321 ( 2018 ). doi: 10.1103/PhysRevC.97.014321 http://doi.org/10.1103/PhysRevC.97.014321 First application of combined isochronous and Schottky mass spectrometry: Half-lives of fully ionized 49CR24+ and 53Fe26+ atoms X. C. Chen , Ph.D. thesis, Ruprecht-Karls-Universität Heidelberg , 2015 . X. C. Chen , Yu. A. Litvinov , M. Wang et al. , Denoising scheme based on singular-value decomposition for one-dimensional spectra and its application in precision storage-ring mass spectrometry . Phys. Rev. E 99 ( 6 ), 063320 ( 2019 ). doi: 10.1103/PhysRevE.99.063320 http://doi.org/10.1103/PhysRevE.99.063320 Denoising scheme based on singular-value decomposition for one-dimensional spectra and its application in precision storage-ring mass spectrometry C.W. Ma , D. Peng , H.L. Wei et al. , A Bayesian-neural-network prediction for fragment production in proton induced spallation reaction , Chin. Phys. C 44 ( 12 ), 124107 ( 2020 ). doi: 10.1088/1674-1137/abb657 http://doi.org/10.1088/1674-1137/abb657 A Bayesian-neural-network prediction for fragment production in proton induced spallation reaction C.W. Ma , H.L. Wei , X.Q. Liu et al. , Nuclear fragments in projectile fragmentation reactions , Prog. Part. Nucl. Phys. 121 , 103911 ( 2021 ). doi: 10.1016/j.ppnp.2021.103911 http://doi.org/10.1016/j.ppnp.2021.103911 Nuclear fragments in projectile fragmentation reactions C.W. Ma , H.L. Wei , X.X. Chen et al. , Precise machine learning models for fragment production in projectile fragmentation reactions using Bayesian neural networks , Chin. Phys. C 46 ( 7 ), 074104 ( 2022 ). doi: 10.1088/1674-1137/ac5efb http://doi.org/10.1088/1674-1137/ac5efb Precise machine learning models for fragment production in projectile fragmentation reactions using Bayesian neural networks B. Schlitt , Ph.D. thesis, Ruprecht-Karls-Universität Heidelberg , 1997 . F.G. Kondev , M. Wang , W.J. Huang et al. , The NUBASE2020 evaluation of nuclear physics properties . Chinese Phys. C 45 ( 3 ), 030001 ( 2021 ). doi: 10.1088/1674-1137/abddae http://doi.org/10.1088/1674-1137/abddae The NUBASE2020 evaluation of nuclear physics properties Q. Wang , baseline estimation program, Github ( 2021 ). https://github.com/NanaVan/baseline-estimatehttps://github.com/NanaVan/baseline-estimate https://github.com/NanaVan/baseline-estimatehttps://github.com/NanaVan/baseline-estimate Accessed 8 Nov 2021 baseline estimation program, Github All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Qian Wang and Xiang-Cheng Chen. The first draft of the manuscript was written by Qian Wang and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-01132-9 This work was partially supported by the National Key R&D Program of China (No. 2018YFA0404401), CAS Project for Young Scientists in Basic Research (No. YSBR-002), and Strategic Priority Research Program of the Chinese Academy of Sciences. (No. XDB34000000). 2022-08-12 2022-09-29 2022-10-02 2022-11-18 2022-12-16 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2022 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-01132-9 http://dx.doi.org/10.1007/s41365-022-01132-9 . Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2022 Nuclear Science and Techniques 11 33