prediction_data = pd.DataFrame(np.array([xl]).T, columns=features)
model_pred, model_std = gpr.predict(prediction_data, return_std=True)

t = 1.96 #2 sided 90% Confidence interval z-infinity score

plt.plot(xl,f(xl),label='source function')
plt.plot(xl,y_measured,'.',label='measured points')
plt.plot(xl,model_pred,'--',label='trained model',color='green')
plt.fill_between(xl,model_pred-t*model_std,model_pred+t*model_std,alpha=0.3,color='green')
plt.legend()
plt.show()