1
$\begingroup$

I have 5 data points with errors associated to them $y_i\pm dy_i$ and the corresponding $x_i$ values (which don't have uncertainties associated to them). I need to calculate the difference between the first of these points, $y_1$ and the rest, and fit a straight line to it (basically the plot will be $\Delta y$ vs $x$). For the other 4 points the error associated with them is just $d(\Delta y_i)=\sqrt{(dy_1)^2+(dy_i)^2}$ for $i$ from 2 to 5. About the first point itself, at that value of $x_1$, $\Delta y_1$ value should be zero. Also, given that this is the reference point, the error associated to that should be zero, too (right?). Now, when I want to make a least square fit, I need to weight the difference between the model and the data by $1/(d(\Delta y_i))$. However that is infinity in the case of the first point (which I guess it makes sense, as I am sure that the line should pass through that point). However I am not sure how to make it work numerically i.e. using a fitting program (Python for example). Should I just replace 0 by something like $10^{-15}$? Does anyone have any advice on how should I handle this infinity? Thank you!

$\endgroup$
  • $\begingroup$ Would you please post a minimal working example with the minimum amount of data that will reproduce the problem? $\endgroup$ – JJacquelin Oct 21 '19 at 6:28
  • $\begingroup$ @JJacquelin the OP is not describing a code problem, rather asks for advice on technique. Please see my answer. $\endgroup$ – James Phillips Oct 21 '19 at 16:49
0
$\begingroup$

Here is a graphical Python fitter with an example of making the first data point's uncertainty to be tiny - that is, the value is very certain - effectively forcing the straight line fit to pass through that point. For comparison the example includes a straight line fit where this is not done.

plot

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

# turn off "no covariance" warnings for simple example
import warnings
warnings.filterwarnings("ignore")

# note the single very small uncertainty for the first point
uncertainties = numpy.array([1.0E-10, 1.0, 1.0, 1.0, 1.0])


# make the first point far from the other Y values
xData = numpy.array([1.1, 2.2, 3.3, 4.4, 5.5])
yData = numpy.array([25.0, 20.2, 30.3, 40.4, 50.5])


def func(x, a, b):
    return (a * x) + b  


# these are the same as the scipy defaults
initialParameters = numpy.array([1.0, 1.0])

# curve fit the test data, first without uncertainties to
# get us closer to initial starting parameters
ssqParameters, pcov = curve_fit(func, xData, yData, p0 = initialParameters)

# now that we have better starting parameters, use uncertainties
fittedParameters, pcov = curve_fit(func, xData, yData, p0 = ssqParameters, sigma=uncertainties, absolute_sigma=True)

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print('Parameters:', fittedParameters)
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot with uncertanties
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # here is the curve without using uncertanties
    yModelNoUnc = func(xModel, *ssqParameters)

    # now the models as line plots
    axes.plot(xModel, yModel, label = 'Force first point')
    axes.plot(xModel, yModelNoUnc, label = 'No Uncertainties')
    plt.legend()

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for this! This is what I needed basically. My only concern was how to pick that very small value. If I pick e-10 or e-5 or e-15 would the result change significantly in general? I assume that the parameters of the fit and the value of chi square would be approximately the same. But would the uncertainty on the parameters/ confidence intervals still be the same? $\endgroup$ – user260669 Oct 21 '19 at 13:26
  • $\begingroup$ I chose a small uncertainty value, but you can make this 1.0E-20 and see that the fit still - in effect - passes through this point. $\endgroup$ – James Phillips Oct 21 '19 at 14:22
  • $\begingroup$ Whether that single data point's uncertainty value us 1.0E-10, 1.0E-15, or 1.0E-20 you get the same coefficient values with this example code. Even a value of 1.0E-9 would be "one billionth" so these values all work fine so long as the other uncertainty values are all 1.0. $\endgroup$ – James Phillips Oct 21 '19 at 14:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.