# Weights for equidistant samples in power law fitting

I am working on fitting analytical curves to experimental data obtained in real viscoelastic tests (in fact, static creep tests). The setting of the problem is:

• the experimental data I have is a set of data points $$(t_{i},y_{i})$$ sampled equidistantly in time (i.e., $$t_{i+1}=t_{i} + \Delta t$$, in which $$\Delta t$$ is fixed and defined by the sampling frequency of the test
• In this field of knowledge, one possible analytical curve known to fit such experimental data is a power law, which can be simplified to $$Y(t)=A\cdot t^{B}$$, with $$B$$ being a positive number smaller than 1. This is the type of curve I am trying to fit to the experimental data I have.

Given that the nature of the curve to be fitted, I wondered if the fact that the experimental data is equidistantly sampled in time may negatively affect the fitting process, as I would have an uneven distribution of data points in different time ranges.

To circumvent that, I devised that I could take two approaches:

• resample the experimental data so it has more or less the same number of data points when analysed on a logarithmic scale
• apply weights during the fitting problem to account for this unbalancing

I would like some insights into these two approaches (and even if my initial worrying is valid). By the way, the numerical fitting is done in Python with SciPy's curve_fit method.

• Welcome to Scicomp. For your given experimental setup, how many data points are we talking about? If your storage and computation resources can handle it, more data is always better. If you have dynamical control over your measuring process, you can also adapt your timesteps as you go, e.g., halve your timestep if the last creep difference was over a certain threshold. Also, do you have a sample plot? Aug 7, 2023 at 12:12
• You could adjust your timesteps like: $\Delta t = min\{\Delta_{min}, 1/\Delta_{ y_{i-1}}\}$ Aug 7, 2023 at 14:03
• Hello @MPIchael! Thanks for the welcoming. Your suggestion to control the measuring process is good, and, in fact, on a real experient, it is typical to have data sampled like at 10 Hz at the first hour, 1 Hz in the first 3 days, and 0.01 Hz onwards, for a total for anywhere between 1 month and 1 year. I was worried on how to objectively define the sampling frequency (usually it is just based on what has been done before). I did not fully get your last suggestion on how to adust time steps. What would be $\Delta _{min}$ and $1/\Delta _{y_{i-1}}$?
– ren1
Aug 9, 2023 at 9:50
• In your measurement i'd call $\Delta_{y-1}$ the 'creep since the last timestep'. If the creep was large, you'd want to make smaller timesteps, if it was small you may take larger ones. The $\Delta_{min}$ might be a minimal fall-back timestep (0.01Hz) in your case. Aug 9, 2023 at 10:57

Given that the nature of the curve to be fitted is exponential

$$A\cdot t^B$$ is not exponential.

But I think it's fine to worry a little bit here, a curve fit could have an unwanted bias towards a region of low interest if we just happened to obtain most samples there. Whether that matters at all depends on how well the analytical function describes the data over the entire range.

It's not clear to me that there is any areas of low interest here, so I don't see anything immediately wrong with using your experimental data as-is. If anything I would expect the asymptotic behavior as $$t$$ gets large to be more important.

We also have the least squares, which ensures the fit doesn't get to far off track

If the true real world behavior isn't quite the expected power law, e.g. $$Y_{\mathrm{true}}(t)=t^{0.5}+0.1 t^{0.4}$$, then you'd had to decide what time range is most important to you and fit that.

But looking at a concrete example ($$A=1$$, $$B=0.5$$):

import scipy.optimize
import numpy as np

t_lin = np.linspace(0, 100, 100)
t_log = np.logspace(-100, 2, 100)
y_lin = t_lin**0.5 + 0.01*(np.random.rand(100)-0.5)
y_log = t_log**0.5 + 0.01*(np.random.rand(100)-0.5)

data_lin, _ = scipy.optimize.curve_fit(lambda t, A, B: A*t**B,  t_lin, y_lin, p0=(1,0.5))
data_log, _ = scipy.optimize.curve_fit(lambda t, A, B: A*t**B,  t_log, y_log, p0=(1,0.5))


Is my data_quad a better match because I tried to compensate my measurements so that y_quad has close to constant $$\Delta$$? No, they are all more or less equally good (just varies with the random noise applied):

>>> data_lin
array([0.99972186, 0.50007081])
array([1.00059216, 0.49983157])
>>> data_log
array([0.99722404, 0.50071082])


thanks to least squares keeping things in check.

I don't think resampling isn't a good strategy here; interpolating data (especially real world potentially noisy and scattered data) would introduce another source of errors you don't need. You could specify the optional sigma parameter to curve_fit, which determines the uncertainty in the data.

sigma : None or M-length sequence or MxM array, optional Determines the uncertainty in ydata. If we define residuals as r = ydata - f(xdata, *popt), then the interpretation of sigma depends on its number of dimensions:

• A 1-D sigma should contain values of standard deviations of errors in ydata. In this case, the optimized function is chisq = sum((r / sigma) ** 2).

which is a way to specify the weights you wanted (though intended for dealing with heteroscedasticity). Using a vector with suitably scaled values

data_lin2, _ = scipy.optimize.curve_fit(lambda t, A, B: A*t**B, t_lin, y_lin, p0=(1,0.5), sigma=np.sqrt(t_lin)+1)


but I doubt this makes much difference at all, and I'm less sure it is in any way an improvement.

Though if you have an estimate for the intended use of sigma here, you should of course definitely use it.

Testing the options here, plotting them all and comparing $$A$$ and $$B$$ is a no-brainer.

• Using synthetic data to verify the sampling is a great idea. Aug 8, 2023 at 1:16
• Thanks a lot, @Mikael! Your take on the matter was exactly what I needed to be sure for the way I was heading. I was afraid of missing something due to inexperience on the matter of curve fitting. And I did mix the things: it is not exponential, but power law.
– ren1
Aug 9, 2023 at 9:45