# Calculating error on a double natural log fit

I am writing residual gas analysis mass spectrometry data reduction software in Python. The evolution of gas intensity $$y$$ over time $$t$$ in the mass spec is roughly a double natural logarithmic function, or:

$$y = a \ln{(pt)} - b \ln{(qt)}$$

where the term $$a \ln{(pt)}$$ represents ingrowth from memory, and the term $$-b \ln{(qt)}$$ represents consumption from ionization.

The goal of fitting these data is to determine the theoretical maximum intensity $$y_0$$. This term $$y_0$$ can't be measured directly because, at the time of theoretical maximum gas intensity, $$t=0$$, the gas is still equilibrating in the mass spec.

Because our system doesn't "start the clock" for $$t$$ and start measuring $$y$$ until after the gas is fully equilibrated, which takes 30 seconds, I offset all timestamps by 30 seconds ($$t_{shifted} = t_{orig} + 30$$) to account for the equilibration period. This makes determining $$y_0$$ as easy as using curve_fit from scipy to get the parameters and evaluating the function at $$t=30$$. This time offset also avoids the need to evaluate $$\ln{(0)}$$ and pushes off the more extreme aspects of the natural log function outside of the fitting window, which better represents the raw data anyway.

Extrapolating the fit to $$y_0$$ at $$t=30$$ is referred to as the "intercept" due to the fact that $$t_{shifted}=30$$ is also $$t_{orig} = 0$$, but should not be confused with the actual intercept at $$\ln(0)$$, which obviously does not exist.

The problem comes when I try to estimate the error on $$y_0$$. I am way out of my depth, but I managed to make it this far:

$$y_0 = a \ln(30p) - b \ln(30q)$$

The partial derivatives of $$y_0$$ are:

$$\frac{\partial y_0}{\partial a} = \ln(30p)$$

$$\frac{\partial y_0}{\partial b} = -\ln(30q)$$

$$\frac{\partial y_0}{\partial p} = a/p$$

$$\frac{\partial y_0}{\partial q} = -b/q$$

After fitting the model using curve_fit, I obtain a covariance matrix pcov which captures the covariance (hence uncertainty) between the estimated parameters.

The variance of $$y_0$$ is:

$$\text{var} = \frac{\partial y_0}{\partial a} \verb|pcov[0,0]| + \frac{\partial y_0}{\partial b} \verb|pcov[1,1]| + \frac{\partial y_0}{\partial p} \verb|pcov[2,2]| + \frac{\partial y_0}{\partial q} \verb|pcov[3,3]|$$

And the error is then $$\sqrt{\text{var}}$$.

Sounds reasonable, right? Unfortunately, curve_fit sometimes fails to return a valid pcov due to highly scattered data and a poor fit, meaning I can't rely on this. When this technique does return a valid pcov, the results range from reasonable (within a few orders of magnitude of the linear model's intercept error) to totally unreasonable (>150% for a relatively decent fit).

I am, as I have said, a bit out of my depth here and would appreciate some guidance on how to estimate the error on $$y_0$$.

Here are the Python functions. Note that I also scale $$y$$ to around ~1, because raw intensity data are very small numbers that can cause problems with curve_fit (see https://stackoverflow.com/questions/78703318/scipy-returning-absurd-curve-fitting-results-using-double-nat-log-function/78703384#78703384).

def double_ln_func(t, a, b, p, q):
return a * np.log(p * t) - b * np.log(q * t)

def double_ln_fit(t, y):
# scale y to make the fitting process more stable
y_oom   = math.floor(math.log10(np.abs(y.max())))
y_scaled = y * 10 ** -y_oom

# offset the t values by 30 seconds to account for gas equilibration time
t_shifted = t + 30

# fit the data to the double nat log function
popt, pcov = curve_fit(double_ln_func, t_shifted, y_scaled, maxfev=10000)
fitted_t   = np.linspace(min(t_shifted), max(t_shifted), 100)
fitted_y   = double_ln_func(fitted_t, *popt)

# calculate intercept by evaluating the function at t = 30
intercept = double_ln_func(30, *popt)

# correct the intercept for the scaling factor
intercept = intercept * 10 ** y_oom

# partial derivatives
a, b, p, q = popt
partials = np.array([
np.log(p * 30),  # ∂intercept/∂a
-np.log(q * 30), # ∂intercept/∂b
a / p,           # ∂intercept/∂p
-b / q           # ∂intercept/∂q
])

# calculate the variance and error of the intercept and correct for the scaling factor
intercept_variance = np.dot(partials, np.dot(pcov, partials))
intercept_error    = np.sqrt(intercept_variance) * 10 ** y_oom

# re-scale the fitted y and t back to the original order of magnitude/scale
fitted_y = fitted_y * 10 ** y_oom
fitted_t = fitted_t - 30

return [fitted_t, fitted_y], intercept, intercept_error


And here's a set of raw data to play with which returns both nan (4 amu, 3 amu), plus one unreasonable error estimate of 1,458% (2 amu) and a reasonable error estimate of 0.82% (40 amu).


He 3798 Q1520_  2024-06-25T09:49:42
time_sec    2.22    3.25    4.24    40.27
2.691   1.918947E-9 1.353086E-9 1.083636E-9 2.424798E-11
8.393   1.924420E-9 1.352441E-9 1.081825E-9 2.451031E-11
13.995  1.924468E-9 1.350765E-9 1.082597E-9 2.493171E-11
19.595  1.925282E-9 1.349962E-9 1.081301E-9 2.494261E-11
25.195  1.927887E-9 1.349370E-9 1.080510E-9 2.491214E-11
30.895  1.930570E-9 1.348243E-9 1.079783E-9 2.517592E-11
36.495  1.932480E-9 1.347466E-9 1.079611E-9 2.510602E-11
42.095  1.931847E-9 1.348131E-9 1.079237E-9 2.513002E-11
47.795  1.934760E-9 1.346782E-9 1.078437E-9 2.512803E-11
53.395  1.929853E-9 1.345735E-9 1.078367E-9 2.498039E-11
58.996  1.934834E-9 1.345541E-9 1.077280E-9 2.520110E-11
64.596  1.932654E-9 1.344552E-9 1.077504E-9 2.550219E-11


And a screenshot showing the fits with the errors on the left:

Since writing this post, I've also attempted to use bootstrapping to estimate the error on the intercept:

def bootstrap_double_ln_fit(t, y, n=2000):
results = []
for i in range(n):
idx = np.random.choice(len(t), len(t))
t_sample = t[idx]
y_sample = y[idx]

popt, pcov = curve_fit(double_ln_func, t_sample, y_sample, maxfev=10000)
results.append(popt)

results = np.array(results)
std_params  = np.std(results, axis=0)

return std_params

def double_ln_fit(t, y):
# scale y to make the fitting process more stable
y_oom    = math.floor(math.log10(np.abs(y.max())))
y_scf    = 10 ** y_oom # y scaling factor for re-scaling results
y_scaled = y * 10 ** -y_oom # scale y to ~1

# offset the t values by 30 seconds to account for gas equilibration time
t_shifted = t + 30

# fit the data to the double nat log function
popt, pcov = curve_fit(double_ln_func, t_shifted, y_scaled, maxfev=10000)
fitted_t   = np.linspace(min(t_shifted), max(t_shifted), 100)
fitted_y   = double_ln_func(fitted_t, *popt)

intercept = double_ln_func(30, *popt) * y_scf

# re-scale fitted_t and fitted_y
fitted_t = fitted_t - 30
fitted_y = fitted_y * y_scf

# bootstrap the fit to get the standard error of the intercept
std_params = bootstrap_double_ln_fit(t_shifted, y_scaled)
intercept_error = std_params[0] * y_scf

return [fitted_t, fitted_y], intercept, intercept_error


This returns consistent errors, but in the thousands of percent range. It also makes the program consistently slower, which I'd prefer to avoid.

What are some other ways of calculating errors on $$y_0$$ that don't rely on pcov? Or perhaps I'm doing something wrong to begin with?

• You lost me by the time you referred to the "$t=0$ intercept," because it doesn't exist. Moreover, your model is overparameterized because $$a+b\log(pt)-c\log(qt)=(a+b\log(p)-c\log(q))+(b-c)\log(t)$$ demonstrates (i) $a$ is not defined and (ii) there are really only two parameters rather than five. Something clearly is very wrong about your formulation.
– whuber
Commented Jul 4 at 15:44
• Can you explain in more detail what you mean? Commented Jul 4 at 15:50
• The "$t=0$ intercept" is admittedly a confusing concept here. The timestamps are offset by 30 seconds so that "$t=0$" is really just $t=30$. For other models like linear, that timestamp offset is unnecessary, but for a natural log function it's required to a) be able to solve for the "intercept" which is again, now $t=30$, and b) push the more extreme parts of the natural log function as it approaches zero out of the range of the fit. Commented Jul 4 at 15:52
• Your model is $y = \alpha + \beta z$ where $z=\log(t),$ $\beta=b-c,$ and $\alpha=a+b\log(p)-c\log(q).$ This is an ordinary least squares model. It has two parameters, $\alpha$ and $\beta.$ It is mathematically impossible to estimate four parameters uniquely from the estimates of two parameters.
– whuber
Commented Jul 4 at 19:23
• can you plot the standard deviation calculation of your bootstrap test statistic as n bootstraps increase? you can save the re-calculated standard deviation every 100 iterations. then plot the value with the n iterations on x axis and estimate for the std.dev on y axis. that is a better way to check convergence. Commented Jul 5 at 16:50

It is surprising that you call "double log fit" because this is a single log fit : $$y=a\ln(pt)+b\ln(qt)$$ $$y=\ln(p^a t^a)+\ln(q^b t^b)=\ln\left((p^a t^a)(q^b t^b) \right)=\ln(p^aq^bt^{a+b})=\ln(p^aq^b)+\ln(t^{a+b})$$ $$y=A+B\:\ln(t)\quad\quad\begin{cases} A=\ln(p^aq^b) \\B=a+b \end{cases}$$ Regression calculus gives $$A$$ and $$B$$ but cannot discriminate $$a$$ and $$b$$ from the sum nor discriminate $$p$$ and $$q$$ from the product. If a software gives some numerical values $$a,b,p,q$$ this is not a unique solution but an hazardous particular solution.

If the above model equation was convenient the graph of $$y$$ function of $$\ln(t)$$ should be roughly linear. This isn't the case for example with 3amu (below).

The shift $$t_0=30$$ improves a lot the linearity.

The results of fitting given below were computed this way :

First : Linear regression for $$A$$ and $$B$$ of $$\quad y=A+B\:X\quad$$ with $$\quad X=\ln(t+30)$$

Second : Drawing of the $$y$$ computed $$Y=A+B\ln(t+30)$$.

The extrapolated $$y$$ at $$t=0$$ is : $$y_0=A+B\ln(30)$$.

RLMSE$$=\sqrt{\frac{1}{n}\Sigma_{k=1}^{k=n}(Y_k-y_k)^2}$$