Fitting methods
For all options (except for the LIN, QDR and the EXP OR LIN),
the minpack.lm::nlsLM function with the LM (Levenberg-Marquardt algorithm)
algorithm is used. Note: For historical reasons for the Monte Carlo
simulations partly the function nls using the port algorithm.
The solution is found by transforming the function or using stats::uniroot.
LIN: fits a linear function to the data using
lm: $$y = mx + n$$
QDR: fits a linear function with a quadratic term to the data using
lm: $$y = a + bx + cx^2$$
EXP: tries to fit a function of the form
$$y = a(1 - exp(-\frac{(x+c)}{b}))$$
Parameters b and c are approximated by a linear fit using lm. Note: b = D0
EXP OR LIN: works for some cases where an EXP fit fails.
If the EXP fit fails, a LIN fit is done instead.
EXP+LIN: tries to fit an exponential plus linear function of the
form:
$$y = a(1-exp(-\frac{x+c}{b}) + (gx))$$
The \(D_e\) is calculated by iteration.
Note: In the context of luminescence dating, this
function has no physical meaning. Therefore, no D0 value is returned.
EXP+EXP: tries to fit a double exponential function of the form
$$y = (a_1 (1-exp(-\frac{x}{b_1}))) + (a_2 (1 - exp(-\frac{x}{b_2})))$$
This fitting procedure is not robust against wrong start parameters and
should be further improved.
GOK: tries to fit the general-order kinetics function after
Guralnik et al. (2015) of the form of
$$y = a (d - (1 + (\frac{1}{b}) x c)^{(-1/c)})$$
where c > 0 is a kinetic order modifier
(not to be confused with c in EXP or EXP+LIN!).
LambertW: tries to fit a dose-response curve based on the Lambert W function
according to Pagonis et al. (2020). The function has the form
$$y ~ (1 + (W((R - 1) * exp(R - 1 - ((x + D_{int}) / D_{c}))) / (1 - R))) * N$$
with \(W\) the Lambert W function, calculated using the package lamW::lambertW0,
\(R\) the dimensionless retrapping ratio, \(N\) the total concentration
of trappings states in cm^-3 and \(D_{c} = N/R\) a constant. \(D_{int}\) is
the offset on the x-axis. Please note that finding the root in mode = "extrapolation"
is a non-easy task due to the shape of the function and the results might be
unexpected.
Fit weighting
If the option fit.weights = TRUE is chosen, weights are calculated using
provided signal errors (Lx/Tx error):
$$fit.weights = \frac{\frac{1}{error}}{\Sigma{\frac{1}{error}}}$$
Error estimation using Monte Carlo simulation
Error estimation is done using a parametric bootstrapping approach. A set of
Lx/Tx values is constructed by randomly drawing curve data sampled from normal
distributions. The normal distribution is defined by the input values (mean = value, sd = value.error). Then, a dose-response curve fit is attempted for each
dataset resulting in a new distribution of single De values. The standard
deviation of this distribution becomes then the error of the De. With increasing
iterations, the error value becomes more stable. However, naturally the error
will not decrease with more MC runs.
Alternatively, the function returns highest probability density interval
estimates as output, users may find more useful under certain circumstances.
Note: It may take some calculation time with increasing MC runs,
especially for the composed functions (EXP+LIN and EXP+EXP).
Each error estimation is done with the function of the chosen fitting method.