I am using scipy curve_fit to fit some data and would like to get confidence intervals for my fitted line. The end goal will be to predict a y value, its upper and lower bounds for any given value of x.
My data looks at how number of properties change with temperature. It is saved in a dataframe. The shape of temp (x) vs property (y) sigmoidal (S-shaped) and I am fitting using a logistic function:
def logistic(x, L=1, k=1, x_0=0, b=0):
y = L / (1 + np.exp(-k*(x - x_0)))+ b
return y
My code is as follows:
#make temporary df
df_temp=df[["Temp (C)", prop]].dropna()
#make temporary means, standard devs and standard errors df
means_stds = df_temp.groupby("Temp (C)").agg(["mean", "std", "count"]).dropna()
means_stds[(prop, "stderr")] = means_stds[(prop, "std")]/np.sqrt(means_stds[(prop, "count")])
means_stds=means_stds[(means_stds !=0).all(1)]
#initial guesses
p0 = [max(df_temp[prop]) - min(df_temp[prop]), initial_guesses_dic[prop], np.median(df_temp['Temp (C)']), min(df_temp[prop])]
#x, y and sigma data for input to curve fit
x = means_stds.index
y = means_stds[(prop, "mean")]
sig=means_stds[(prop, "std")]
#fit data
popt_dic[prop], _pcov = curve_fit(logistic, x, y, p0=p0, sigma=sig, absolute_sigma=True, method='lm')
perr=np.sqrt(np.diag(_pcov))
#plot
plt.plot(x, y, marker=".", linewidth=0, label='all data')
plt.plot(x, logistic(x, *popt_dic[prop]), linestyle="--", color="r", label='fit of all data')
plt.fill_between(x, logistic(x, *(popt_dic[prop]+2*perr)), logistic(x, *(popt_dic[prop]-2*perr)), alpha=0.25, label='confidence interval')
plt.ylabel(prop)
plt.legend()
plt.show()
This all works seems to work fine except that my confidence intervals don't look right: example fit of data and example fit of data 2
I have also tried to calculate the confidence intervals by just changing the L and b values (height and minimum y-value respectively) of the fit:
#fit data
popt_dic[prop], _pcov = curve_fit(logistic, x, y, p0=p0, method='lm', sigma=sig)
#find lower and upper bounds (using L and b only - keep steepness and midpoint same)
perr=np.sqrt(np.diag(_pcov))
popt_lwr_dic[prop]=popt_dic[prop] - 2 *perr
popt_lwr_dic[prop][1:3]=popt_dic[prop][1:3]
popt_upr_dic[prop]=popt_dic[prop] + 2 *perr
popt_upr_dic[prop][1:3]=popt_dic[prop][1:3]
#plot
plt.plot(x, y, marker=".", linewidth=0, label='all data')
plt.plot(x, logistic(x, *popt_dic[prop]), linestyle="--", color="r", label='fit of all data')
plt.fill_between(x, logistic(x, *popt_lwr_dic[prop]), logistic(x, *popt_upr_dic[prop]), alpha=0.25, label='confidence interval')
plt.ylabel(prop)
plt.legend()
plt.show()
This seems to produce good plots but I'm not sure the maths behind it is sound.
What is the correct way to find the confidence intervals for this fitted line?