3

I am trying to fit a power law of the form a*x**b+c to some data points, using curve_fit from scripy.optimize. Here's the MWE:

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

def func_powerlaw(x, m, c, c0):
    return c0 + x**m * c

x = np.array([1.05, 1.0,  0.95, 0.9,  0.85, 0.8,  0.75, 0.7,  0.65, 0.6,  0.55])
y = np.array([1.26, 1.24, 1.2,  1.17, 1.1,  1.01, 0.95, 0.84, 0.75, 0.71, 0.63])
dy = np.array([0.078]*11)

fig, (a1) = plt.subplots(ncols=1,figsize=(10,10))
a1.errorbar(x, y, yerr = dy, ls = '', marker='o')

popt, pcov = curve_fit(func_powerlaw, x, y, sigma = dy, p0 = [0.3, 1, 1], bounds=[(0.1, -2, -2), (0.9, 10, 2)], absolute_sigma=False, maxfev=10000, method = 'trf')
perr=np.sqrt(np.diag(pcov))

xp = np.linspace(x[0],x[-1], 100)
a1.plot(xp, func_powerlaw(xp, *popt), lw=3, zorder = 1, c = 'b')
print(popt, perr)

Output: [0.35609897 3.24929422 -2.] [0.47034928 3.9030258 3.90965249]

Here's the Figure.

For all three parameters, the errors are larger than the value estimations themselves. Judging from experience this cannot be right, since the line fits the data points very well. Even if I don't set any bounds and/or initial guess, the values change, but the errors remain too high. The only boundary required is that 0.1<=m<=0.9. What am I doing wrong? Any help is greatly appreciated!

George
  • 451
  • 1
  • 6
  • 17
  • 2
    The uncertainty `dy` means that any value `y_i +/- dy` is compatible with the observed data `y`. Large errors on the observed data can lead to large errors on the fitted parameters. Try the following: `for y in np.random.uniform(y - dy, y + dy, size=(25, len(y))): ...` and then perform the fitting for each randomly sampled `y` (sampled within its uncertainty boundaries) but *without* indicating an error (`sigma`) to the fit procedure. Then record the resulting parameter estimates (`popt`) and take a look at `np.mean(popts), np.std(popts)`. The second value should match the one from your fit. – a_guest Sep 24 '21 at 16:02
  • Take my opinion with an ocean of salt, as I never do these kind of things, but could it be possible that it is overfitting? – chess_lover_6 Sep 24 '21 at 16:02
  • @a_guest Following this logic, shouldn't removing the sigmas altogether result in a better fit? Because the opposite seems to be true. – George Sep 24 '21 at 16:35
  • @chess_lover_6 I've also tried `lmfit`, which is a wrapper of curve_fit as far as I understood and which prints the χ^2 values. They are bad, not even close to 1. So I don't think it is because of overfitting. – George Sep 24 '21 at 17:02
  • @George First of all, using `sigma` with values that are all equal together with `absolute_sigma=False` does not have any effect since, as mentioned in the docs, *if `False` (default), only the relative magnitudes of the sigma values matter*. Nevertheless, it seems that this isn't relevant for your problem at hand, since when you remove the `sigma` altogether, this still results in a large `perr`. If I understood correctly, this is your actual question: The resulting curve fits the data points well (judging from the plot, or R2 score for example), however the reported `perr` is very large. – a_guest Sep 25 '21 at 08:15
  • If so, please update your question to focus on this part of the problem in order to make it a minimal example. – a_guest Sep 25 '21 at 08:16

1 Answers1

2

The size of the errors in your fit parameters is partly driven here by the size of the measurement errors (dy in your code). As you can see in the plot the sigmas are large relative to the dispersion of the points, so a wide range of curves could fit the data. The dy values are hard coded here, are they real values? Try making the smaller and see how that affects the curve fit errors. Also they look like absolute sigmas in the plot, so you should set the absolute_sigma flag to True.

As an added note, if you don't supply measurement errors in the sigma parameter the default is 1.0 (very large in your case) and not 0.0 for every y value. This explains why omitting sigma gives larger fit errors.

Conor
  • 1,028
  • 1
  • 8
  • 15
  • dy are indeed real values. If the sigma values were to blame for the large errors, shouldn't removing them altogether decrease the error values? Trying that increases the error values by a lot. So does setting the absolute_sigma value to true. – George Sep 24 '21 at 16:31
  • If you pass in None (the default) for sigma, the default measurement error is used for each y point, which is 1.0. In your case this is a large error and thus the error on the fitted parameters is large. The default value is not very helpful in this case! – Conor Sep 25 '21 at 15:24
  • Try setting the measurement error to 0.01 (for example) and see what the fit looks like. Remember that these are assumed to be 1 sigma errors, so if the measurements are normally distributed around the true values, the best fit line should pass through roughly 2/3 of the error bars. – Conor Sep 25 '21 at 15:26