| decompose_OSLcurve {OSLdecomposition} | R Documentation |
Multi-exponential CW-OSL decomposition
Description
Function for determining the signal component amplitudes of a multi-exponential decay curve if the signal component decay parameters are already given. Thus, this function decomposes CW-OSL curves with known components of unknown intensity.
Usage
decompose_OSLcurve(
curve,
components,
background.fitting = FALSE,
algorithm = "det",
error.estimation = "empiric",
verbose = TRUE
)
Arguments
curve |
data.frame or matrix or RLum.Data.Curve (required):
CW-OSL curve x-Axis: |
components |
data.frame or numeric vector (required):
Either a vector containing the decay parameters of the CW-OSL components or a table (data.frame), usually the table returned by fit_OSLcurve.
In case of a vector: It is recommended to use less than 7 parameters. The parameters will be sorted in decreasing order.
In case of a data.frame, one column must be named |
background.fitting |
logical (with default):
if |
algorithm |
character string (with default):
Choice of curve decomposition approach. Either |
error.estimation |
character string (with default):
integral error estimation approach, either |
verbose |
logical (with default): enables console text output |
Details
The function assumes multiple exponentially decaying signal components with first-order kinetics:
I(t) = n_1 \lambda_1 exp(-\lambda_1 t) + n_2 \lambda_2 exp(-\lambda_2 t) + ... + n_K \lambda_K exp(-\lambda_K t)
with I(t) the CW-OSL signal, n the signal component intensity,
\lambda the signal component decay constant and K the number of signal components.
For the actual decomposition procedure, the integrated version of this formula is used, see Mittelstrass et al. (2021) for details.
Decomposition algorithm
The calculation procedure depends on the function argument algorithm.
This function includes two different decomposition algorithms: "det" for determinant solution
and "nls" for nonlinear least squares estimate
algorithm = "det" (default)
The function calculates the CW-OSL component intensities by building an equation system which is then solved by a determinant-based approach (Cramers rule). This purely analytical approach gives the algorithm a solution in all possible cases, even if the measurement consists just of noise or the wrong model is used. There are also no 'false minima' events. The statistical error is calculated by applying the propagation of uncertainty method on Cramers rule.
The precision of this algorithm as well as the propagation of eventual systematic errors of the decay rate values,
depend on the integration intervals, given by the columns $t.start, $t.end, $ch.start and $ch.end
of the data.frame used as input for the argument components.
In principle, these can be chosen freely. Reasonable integration intervals are defined by optimise_OSLintervals.
If not defined, the logarithmic mean values between life times (reciprocal decay rate) of subsequent components are
used as interval borders.
algorithm = "nls"
As alternative algorithm, Levenberg-Marquardt nonlinear regression is available, see minpack.lm::nlsLM for details.
The results are identical to that of the "det" algorithm in accuracy and precision. But there is the slight chance (< 1 %)
of fitting failure when using the "nls" algorithm. Also, the statistical errors are underestimated by 20-80 %
in most cases. As advantage, the "nls" algorithm is less sensitive against systematic errors
caused by uncorrected signal background.
algorithm = "det+nls"
Both algorithms can be combined. Then, "det" provides the startings values and the error estimations for
"nls" and returns replacement results, in case "nls" fails. "nls" compensates for potential systematic
errors in the fast and medium components intensity values due to uncorrected signal background. However, the
background signal will still affect slow component results. The slowest component will be overestimated while
the second slowest component will be underestimated. If these components are of particular interest,
it is recommended to set background.fitting = TRUE
All three methods were tested at 5x10^6 simulated CW-OSL curves by Mittelstrass (2019) for their performance (+++ reliable results in all cases; ++ reliable in >95% of cases: + reliable in most cases):
| Characteristics | det | nls | det+nls |
| Decomposition success rate | 100 % | >99 % | 100 % |
| Component intensity accuracy | +++ | +++ | +++ |
| Accuracy in case of uncorrected background | + | ++ | ++ |
| Error estimation accuracy | +++ | + | ++ |
In summary, algorithm = "det" is recommended for the most cases. If the signal background level is
significant (> 2 % of initial signal) but was not corrected, algorithm = "det+nls" is the better choice.
Setting background.fitting = TRUE is usually not recommended, only in case slow components shall
be investigated in measurements with uncorrected background.
Error estimation
In case of algorithm = "det" or "det+nls" the Propagation of Uncertainty method is used to
transform signal bin error values (column $bin.error) into component intensity error values (column $n.error). The signal bin error
calculation depends on the argument error.estimation, see below.
If algorithm = "nls" is used, the error values provided by minpack.lm::nlsLM are returned.
error.estimation = "empiric" (default)
The standard deviation of each signal bin (signal bin = signal value of an integrated time interval) is calculated from the corrected sample variance between the CW-OSL model and the actual CW-OSL curve for that interval. Thus, statistical errors are monitored accurately without any prior knowledge required. However, potential systematic errors are monitored insufficiently. Also, at least two (better more) data points per signal bin are needed to estimate its standard deviation. If a signal bin consists just of one data point, its square root value is taken as standard deviation, in accordance to the Poisson distribution.
error.estimation = "poisson" or numeric value
Alternatively the standard error can be calculated by approximating a Poisson distributed signal error, known as Shot noise. This is suitable if the lack of data points on the x-axis circumvents an empiric error estimation, like with spatially or spectrally resolved CCD measurements. Also the parameter can be set to a numeric value, which represents the detector noise in cts / s and is assumed to be normal distributed. The detector noise will be added on top of the Poisson distributed shot noise.
error.estimation = "only.bin.RSS"
The error estimation is omitted but the residual sum of squares (RSS) between input curve and combined
signal component curves is calculated. However, the RSS value is divided into sections according to
the signal bins (column $bin.RSS). The full RSS value can be calculated by summing over the complete column.
The RSS value is usually used a minimization target in fitting algorithms, like done in fit_OSLcurve.
The values of the $bin.RSS column allows for weighted fitting by applying pre-factors to the bin RSS values.
For further speed advance, the calculation of $components$n.residual and $components$initial.signal is
also omitted.
error.estimation = "none"
The error estimation is omitted. This option saves significant computing time, if the error estimation is
not of significance. For further speed advance, the calculation of $components$n.residual and
$components$initial.signal is also omitted.
Systematic errors
The ratio of the error values of both error estimation methods can be used to detect (but not quantify) systematic
errors. "poisson" error values are not affected by systematic errors, while "empiric" errors are.
If the detector noise is known and taken into account, the relation between both values for a given
signal bin should be about empiric / poisson = 1. In case of systematic errors, this ratio increases.
Value
The input table components data.frame will be returned with added or overwritten
columns: $n, $n.error, $n.residual, $bin, $bin.error, $bin.RSS, $initial.signal.
Which columns are written depends on the selected parameters. If an input data.frame contains already
one of the above columns but parameters are selected which do not re-calculate the values, the values
of the columns are set to NA.
Last updates
2022-07-25, DM: Extended algorithm for bin-wise RSS calculation and added error estimation option "only.bin.RSS"
Author(s)
Dirk Mittelstraß, dirk.mittelstrass@luminescence.de
Please cite the package the following way:
Mittelstraß, D., Schmidt, C., Beyer, J., Heitmann, J. and Straessner, A.: R package OSLdecomposition: Automated identification and separation of quartz CW-OSL signal components, in preparation.
References
Mittelstraß, D., 2019. Decomposition of weak optically stimulated luminescence signals and its application in retrospective dosimetry at quartz (Master thesis). TU Dresden, Dresden.
See Also
fit_OSLcurve, optimise_OSLintervals, RLum.OSL_decomposition, minpack.lm::nlsLM
Examples
# Set some arbitrary decay parameter for a dim CW-OSL measurement of quartz
components <- data.frame(name = c("fast", "medium", "slow"),
lambda = c(2, 0.5, 0.02),
n = c(1000, 1000, 10000))
# Simulate the CW-OSL curve and add some signal noise and some detection background
curve <- simulate_OSLcomponents(components, simulate.curve = TRUE,
add.poisson.noise = TRUE, add.background = 40)
# Decompose the simulated curve
components <- decompose_OSLcurve(curve, components)
# Display the component separation results
plot_OSLcurve(curve, components)
### Decomposition including signal background fitting:
# Define optimized integration intervals, including an interval for the background
components <- optimise_OSLintervals(components, curve, background.component = TRUE)
# Decompose again and view results
components <- decompose_OSLcurve(curve, components, background.fitting = TRUE)
plot_OSLcurve(curve, components)