0
$\begingroup$

I am try to measure the error in the slope of the black dots. To do this I performed a linear regression to the data using from scipy.optimize import curve_fit on python and plotted it as shown by the dashed red lines. The legend shows the equation of the line of this fit along with the $R^{2}$ value.

I have also plotted some measurements of the same data but using a different technique which is the shown by the blue cross. The reason for this was to determine the systematic error in the measurement technique used to get the black dots. Similarly I have also performed a linear regression of this data as shown in the legend by the blue dash. So my question is, what is the total error in gradient of the data shown by the black dots.

plot

$\endgroup$

1 Answer 1

1
$\begingroup$

Without knowing the true slope there is no unique way of determining the error of the slope. So, all you can do is to select a method to determine the slope and then calculating the associated uncertainty. E.g. using the least square fit, the uncertainty of the slope estimate $\hat \beta$ is usually taken to be the standard deviation $$ Sd[\hat\beta_j] = \hat\sigma_\epsilon \sqrt{[(X^T X)^{-1}]_{jj}} $$ where $X$ is the so called design matrix. I use the index $jj$ to indicate that we only take the diagonal elements. The point estimate of the "random error" $\hat\sigma_\epsilon$ is often taken to be $\sqrt{\sum_{i=1}^N (y - y_{fit})^2/(N-p)}$, where $N$ is the number of data points and $p$ is the number of factors (including the intercept).

Often, I find it helpful to have an example to check my calculations. Therefore, here an example:

## Sample data:
x  = 1:10
y  = c(1.77954908388273, 2.44621066683337, -2.2018205040772, -0.764061711605929, -2.72481609322593, -3.51410342713565, -7.86973095833333, -8.78328718047611, -11.5217536084918, -6.80619245104001)
df = data.frame(x,y)

which look like this

enter image description here

library(MASS) # for inverse ginv()

## design matrix
X      = c(rep(1, N), x) 
dim(X) = c(N,2) # reshape the design matrix, so that its dimension is (N,2)

## Fit data (formulas not included above):
betaHat = ginv(t(X) %*% X) %*% t(X) %*% y 
yFit    = X %*% betaHat

## Estimate of the uncertainty:
res     = y - yFit   # residuals
sigma   = sqrt( sum(res^2)  / (N-2) ) # use p=2, because we fit slope and intercept
betaSD  = sigma * sqrt( diag(ginv(t(X) %*% X)) )

> Coefficients:
>            Estimate    Std. Error    t value   Pr(>|t|)    
> (Intercept)   3.6727     1.3757      2.670     0.028372 *  
> beta         -1.3943     0.2217     -6.289     0.000236 ***

Looking at your data, there are two regions where the data does not appear to be linear. This is at approx $x=15$ and $x = 115$. Thus, you should decide whether or not these regions are important. If they are not important you should consider omitting these regions, so that your fit is improved. Note that the data points at the "$x$-edge" carry more leverage. Hence, they are particularly dangerous to distort your fit.

$\endgroup$
1
  • $\begingroup$ Very descriptive and insightful answer! Thank you! I will find the sd of the slopes. I dont suppose there is much of a difference in your method and using from scipy.optimize import curve_fit as I have done to achieve the slopes above. $\endgroup$ Commented Dec 11, 2020 at 22:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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