[FIXED] Python for confidence limit

Issue

I try to simulate the confidence limit of the mean value in the figure above. But something wrong with my code. The result is far different than the above figure. I used confidence interval to obtain the slope and intercept. Could someone one give me hints ? Thanks

Here is my code.

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import t
from scipy import stats
x = np.array([3,7,11,15,18,27,29,30,30,31,31, 32,33,33,34,36,36,36,37,38,39,39,39,40,41,42,42,43,44,45,46,47,50])
y = np.array([5,11,21,16,16,28,27,25,35,30,40,32,34,32,34,37,38,34,36,38,37,36,45,39,41,40,44,37,44,46,46,49,51])

mean_x = np.mean(x)
n = len(x)
res = stats.linregress(x, y)
tinv = lambda p, df: abs(t.ppf(p/2, df))
ts = tinv(0.05, n-2)
bins = np.linspace(0,3,54)

plt.plot(x, y, 'o', label='Data Points')
plt.plot(x,res.intercept + res.slope*x, 'c', label='fitted line')
plt.plot(x,(res.intercept+ts*res.intercept_stderr)+(res.slope + ts*res.stderr)*x,'b', label='Upper Limit')
plt.plot(x,(res.intercept-ts*res.intercept_stderr)+(res.slope - ts*res.stderr)*x, 'g' ,label='Lower Limit')
plt.legend()
plt.show()

Solution

There are lots of ways to estimate the error of your model. As an example, here is a linear fit of the error:

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import t
from scipy import stats
x = np.array([3,7,11,15,18,27,29,30,30,31,31, 32,33,33,34,36,36,36,37,38,39,39,39,40,41,42,42,43,44,45,46,47,50])
y = np.array([5,11,21,16,16,28,27,25,35,30,40,32,34,32,34,37,38,34,36,38,37,36,45,39,41,40,44,37,44,46,46,49,51])

n = len(x)
res = stats.linregress(x, y)
tinv = lambda p, df: abs(t.ppf(p/2, df))
ts = tinv(0.05, n - 2)
bins = np.linspace(0, 3, 54)

pred = res.intercept + res.slope * x
errors = np.abs(y - pred)

plt.scatter(x, y, label='Data Points')
plt.plot(x, pred, 'c', label='fitted line')


error_res = stats.linregress(x, errors)
pred_err = error_res.intercept + error_res.slope * x
plt.plot(x, pred + 1 * pred_err, label='Upper confidence')
plt.plot(x, pred - 1 * pred_err, label='Lower confidence')

plt.legend()
plt.savefig('hey')

Answered By - kilojoules