model_pred, model_std = gpr.predict(prediction_data, return_std=True)

t = 1.96 #2 sided 90% Confidence interval z-infinity score
plt.figure(figsize=(8,4))
plt.subplot(2,1,1)
plt.plot(xl,f(xl),'b-',label='Source function')
plt.plot(xl,y_measured,'.',label='Measured points')
plt.plot(xl,model_pred,'--',label='GPR model',color='green')
plt.fill_between(xl,model_pred-t*model_std,model_pred+t*model_std,alpha=0.3,color='green')
plt.scatter([opt_val[0]],[opt_val[1]],label='GPR model minimized',color='red',marker='s',s=100,zorder=3)
plt.scatter([opt_val2[0]],[opt_val2[1]],label='GPR uncertainty minimized',color='black',marker='o',s=100,zorder=3)
plt.scatter([opt_val3[0]],[opt_val3[1]],label='Combined optimization',color='blue',marker='^',s=100,zorder=3)
plt.legend(loc='center left',bbox_to_anchor=(1, 0.5),fontsize=10,frameon=False)
plt.subplot(2,1,2)
plt.plot(xl,model_std,'k--',label='GPR uncertainty')
plt.scatter([opt_val2[0]],[opt_val2[2]],label='GPR uncertainty minimized',color='black',marker='o',s=100,zorder=3)
plt.legend(loc='center left',bbox_to_anchor=(1, 0.5),fontsize=10,frameon=False)
plt.tight_layout()
plt.show()