![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\" "\\"Open"){kind=icon}
![\"Download\"](\"https://img.shields.io/badge/Github-Download-blue.svg\" "\\"Download"){kind=icon}
"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[11.3 The Kennedy-O’Hagan Predictive Model](https://ndcbe.github.io/cbe67701-uncertainty-quantification/11.03-Contributed-Example.html#11.3-The-Kennedy-O’Hagan-Predictive-Model)",
"section": "11.3 The Kennedy-O’Hagan Predictive Model"
}
},
"source": [
"# 11.3 The Kennedy-O’Hagan Predictive Model"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 1,
"link": "[11.3 The Kennedy-O’Hagan Predictive Model](https://ndcbe.github.io/cbe67701-uncertainty-quantification/11.03-Contributed-Example.html#11.3-The-Kennedy-O’Hagan-Predictive-Model)",
"section": "11.3 The Kennedy-O’Hagan Predictive Model"
}
},
"source": [
"Created by Jiale Shi (jshi1@nd.edu)\n",
"These examples and codes were adapted from:\n",
"* McClarren, Ryan G (2018). *Uncertainty Quantification and Predictive Computational Science: A Foundation for Physical Scientists and Engineers*, *Chapter 11 : Predictive Models Informed by Simulation, Measurement, and Surrogates*, Springer, https://link.springer.com/chapter/10.1007/978-3-319-99525-0_11\n",
" and its supplementary material **Calibration Simple with Discrep.ipynb**\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbpages": {
"level": 1,
"link": "[11.3 The Kennedy-O’Hagan Predictive Model](https://ndcbe.github.io/cbe67701-uncertainty-quantification/11.03-Contributed-Example.html#11.3-The-Kennedy-O’Hagan-Predictive-Model)",
"section": "11.3 The Kennedy-O’Hagan Predictive Model"
}
},
"outputs": [],
"source": [
"## import all needed Python libraries here\n",
"from pyDOE import *\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"%matplotlib inline\n",
"\n",
"import math\n",
"from scipy.stats import multivariate_normal\n",
"from scipy.stats import norm\n",
"import numpy as np\n",
"np.random.seed(1000)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbpages": {
"level": 1,
"link": "[11.3 The Kennedy-O’Hagan Predictive Model](https://ndcbe.github.io/cbe67701-uncertainty-quantification/11.03-Contributed-Example.html#11.3-The-Kennedy-O’Hagan-Predictive-Model)",
"section": "11.3 The Kennedy-O’Hagan Predictive Model"
}
},
"outputs": [],
"source": [
"#covariance function\n",
"def cov(x,y,beta,l,alpha):\n",
" exponent = np.sum(beta*np.abs(x-y)**alpha)\n",
" return 1/l * math.exp(-exponent)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"nbpages": {
"level": 1,
"link": "[11.3 The Kennedy-O’Hagan Predictive Model](https://ndcbe.github.io/cbe67701-uncertainty-quantification/11.03-Contributed-Example.html#11.3-The-Kennedy-O’Hagan-Predictive-Model)",
"section": "11.3 The Kennedy-O’Hagan Predictive Model"
}
},
"outputs": [],
"source": [
"#likelihood function\n",
"def likelihood(z,x,beta,lam,alpha, beta_t, lam_t, alpha_t, meas_cov, N,M):\n",
" Sig_z = np.zeros((N+M,N+M))\n",
" #fill in matrix with sim covariance\n",
" for i in range(N+M):\n",
" for j in range(i+1):\n",
" tmp = cov(x[i,:],x[j,:],beta,lam,alpha)\n",
" if (i < N):\n",
" tmp += cov(x[i,0], x[j,0], beta_t, lam_t, alpha_t)\n",
" Sig_z[i,j] = tmp\n",
" Sig_z[j,i] = tmp\n",
" #add in measurement error cov\n",
" Sig_z[0:N,0:N] += meas_cov\n",
" #print(Sig_z)\n",
" likelihood = multivariate_normal.logpdf(z,mean=0*z, cov=Sig_z,allow_singular=True)\n",
" return likelihood"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbpages": {
"level": 2,
"link": "[11.3.1 Toy example of KOH model](https://ndcbe.github.io/cbe67701-uncertainty-quantification/11.03-Contributed-Example.html#11.3.1-Toy-example-of-KOH-model)",
"section": "11.3.1 Toy example of KOH model"
}
},
"source": [
"## 11.3.1 Toy example of KOH model\n",
"To demonstrate the behavior of the predictive model, we consider a simple simulation code\n",
"\n",
"$$\\eta (x,t) = \\sin xt$$\n",
"\n",
"We also consider measurement generated from the function\n",
"\n",
"$$y(x) = \\sin 1.2 x + 0.1 x + \\epsilon =\\eta(x,1.2)+0.1 x + \\epsilon$$\n",
"\n",
"where $\\epsilon$ is a measurement error that is normally distributed with mean 0 and standard deviation 0.005.\n",
"We will use the Kennedy-O'Hagan model to estimate the calibration parameter, which in this case has a true value of $t=1.2$, and fit **a discrepancy function**. We know that the true discrepancy function is linear, and we can compare our estimate to the true function.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"nbpages": {
"level": 2,
"link": "[11.3.1 Toy example of KOH model](https://ndcbe.github.io/cbe67701-uncertainty-quantification/11.03-Contributed-Example.html#11.3.1-Toy-example-of-KOH-model)",
"section": "11.3.1 Toy example of KOH model"
},
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFaNJREFUeJzt3X2MHPV9x/HPp4fBpybi0uAk+ICatugaFydxdSIPRkIkpOdEET5cSIA+kJYIRQlKkNprfaJKSZTqHFlqURXUxilpHhSeRI6zG9xeIXZApUnK0SMYG65yKYS7Q8GBHG3SdWM73/5xu7C37D3s7ezM7M77JZ12d3Y881slfHb2+3sYR4QAAMXyC1k3AACQPsIfAAqI8AeAAiL8AaCACH8AKCDCHwAKiPAHgAIi/IEqtl9j+2nbV1dte63tH9i+vGbfM22/YPvCmu1fs/21tNoMrIaZ5AUsZPu3JH1d0saIOGr7byS9MSK219n3DyTtkPTWiDhW/rdflfQbEfFCqg0HGkD4A3XY/rKk0yR9QdI3JJ0fEc8tsu8/SZqU9FlJj0saioi7U2oqsCqEP1CH7ddJOixpjebD/O+X2PccSY9K+hdJP4uIyxfbF8gLwh9YhO37Jb1L0pkR8dIy+35S0l9I+tWI+GEa7QOaQYcvUIft35W0QdL9kj5X3naO7Z9U/mr+ySFJPyL40S5OyboBQN7YfoOkv5L0QUlPSjpk+7aIeFDSazJtHJAQrvyBV/u8pLGIOFDu5P0TSV+0fVrG7QISQ/gDVWwPSrpQ0lBlW0T8naRpSZ/Kql1A0ujwBYAC4sofAAqI8AeAAiL8AaCACH8AKKDcjvM/44wzYsOGDVk3AwDayiOPPPKjiFi33H65Df8NGzZoYmIi62YAQFux/cxK9qPsAwAFRPgDQAER/gBQQIQ/ABRQIuFv+0u2n7f9+CLv2/Zf2z5i+zHbv5nEeQEAq5PUaJ8va34lxK8u8v77JJ1X/nu7pL8pPwJowNjkjHaNT2l2rqT1Pd0aGujT4OberJuFNpTIlX95nfMXl9hlm6SvxrzvSuqxfWYS5waKYmxyRsOjBzUzV1JImpkraXj0oMYmZ7JuGtpQWjX/XknPVr2eLm9bwPZ1tidsTxw9ejSlpgHtYdf4lErHTy7YVjp+UrvGpzJqEdpZWuHvOttetZZ0ROyOiP6I6F+3btkJakChzM6VGtoOLCWt8J+WdHbV67MkzaZ0bqAjrO/pbmg7sJS0wn+vpN8vj/p5h6SXyrfHA7BCQwN96l7TtWBb95ouDQ30SQdGkj1Z0sdD7iQ11PN2Sd+R1Gd72va1tj9q+6PlXfZJekrSEUlflPSxJM4LFMng5l6NbN+k3p5uWVJvT7dGtm+aH+3zwM6VH2glwd7I8dCWEhnqGRFXLfN+SPp4EucCimxwc+/KhnYeGJEuHq7/3gM7F38PhZHbVT0BLOPAyMIr9JtOn3+8aMfqAn6p4/Fl0XFyewP3/v7+YElnYIVuOl266aXFX9cGe8ViwV7779E2bD8SEf3L7ceVP9Aplrtyr4Q8wQ4R/kBnSDrgL9qRTLuQW6zqCXSCRmryKwl2avwdj/AHOtFSAZ9msDNfILcIf6AT5eXKnfkCuUX4A0AB0eELIFnMF2gLhD+AZDGstC1Q9gHSQucncoTwB9JSxM5P5gvkFuEPoHWo8ecWNX+glRbp/Hyy72P69asoAyE7XPkDrXTxsMa2HdabT94pSdpw7DZtOHabLjt8ETdeR6YIf6DFuPF6dsYmZ7Rl536du+Nebdm5v32+cFMYHED4Ay1WucH6zSe2192O1hibnNHw6EHNzJUUkmbmShoePdgeXwApDA4g/IEWq9xg/eYTl9fdjtbgF9fSCH+gxZa88TpaZrFfVrn9xXVgpDwprjwjuvK8RSUgRvsALVa55+6u8SnNzpW0vqdbQwN9K7sXL1ZtfU+3ZuoEfW5/caU8M5rwB1Kw4huvIzFDA30aHj24oPTDL65XEP4AOlJb/+JKYWY0N3AH8uzAyHwpoPIILGOlN3CnwxfIs8qQvyKuC4SWIvwBoICo+QN5s9jNULgpChLElT+QNxcPzw/zqwz1q3686aXOD37ue5AKwh9oFUJsdejfSAXhD7RKEiFWGfLHTVGQMGr+QJ5VSjxFKPVw0/dUMc4faFb1GPzaEKsgxFaOm743ZaXj/LnyB5r1wM6FV+gprs8CrBY1fwD5Qv9GKrjyB1ZjJTVqQmx1KI+lgpo/0CzKO8gR1vYBACyK8AeaRXkHbYjwB5pFjRptKJHwt73V9pTtI7ZfdRlk+8O2j9p+tPz3kSTOCwBYnaZH+9juknSLpPdKmpb0sO29EXG4Ztc7I+L6Zs8HAGheElf+F0g6EhFPRcTPJN0haVsCxwXyjYXb0MaSCP9eSc9WvZ4ub6v127Yfs3237bPrHcj2dbYnbE8cPXo0gaYBLcTqk2hjSYS/62yrnTzwD5I2RMRbJN0v6Sv1DhQRuyOiPyL6161bl0DTAAD1JDHDd1pS9ZX8WZJmq3eIiBeqXn5R0ucSOC+QPlafRIdIIvwflnSe7XMlzUi6UtLV1TvYPjMiniu/vFTSEwmcF0hfqxZuq14ZFEhB02WfiDgh6XpJ45oP9bsi4pDtz9i+tLzbJ2wfsv19SZ+Q9OFmzwt0FPoPkLJEFnaLiH2S9tVs+1TV82FJXNagszCzF22MVT2B1Wq2TEP/ATJE+ANZ4cYvyBBr+wBAARH+QB7Qf4CUEf5AHlDjR8oIf2AlWMcHHYbwBySNTc5oy879OnfHvdqyc7/GJmcW7sA4fHQYRvug8MYmZzQ8elCl4yclSTNzJQ2PHpQkDW6ut0Yh0P648kfh7Rqfejn4JemGU+5W6fhJvfDNT5eHYJbH31eeUwJCByD8UXizc6UFr284ZVSS9Nmfbpsfe18Zf195TucsOgDhj8Jb39Pd0HagExD+KLyhgT798amjenrt1Xp67fyCtE+vvVoPHbvslRIP4/DRYejwReENbu7VmD6jLeO/o9m5kv5r7dUa23Z4YWcvpR50GMIf0PwXwMthfxOjfND5KPsAtSjxoAAIf6AWJR4UAOEPAAVE+ANAARH+ADN2UUCEP8CibSggwh8ACohx/igmbp6OgiP8UUzcPB0FR9kHAAqI8Efx1I7uYUYvCojwR/HUju6hxo8CIvwBoIDo8EUxMLoHWIDwRzEwugdYgLIPABQQ4Y/iYXQPQNkHnWdscka7xqc0O1fSn/3iHr3+A3/OLRmBGoQ/OsrY5IyGRw+qdPykJOnak3fqzaPbJXFrRqAaZR90lF3jUy8Hf0Xp+EntGp/KqEVAPnHlj44yO1fSDafcrRtOGX1529Nrr5aOSTrAsE6ggvBHR1nf0y39RNpw7DZJ88G/4dht6u3p1kMXvzvj1gH5QdkHHWVooG/BVb8kda/p0tBAX0YtAvIpkfC3vdX2lO0jtl81js72abbvLL//PdsbkjgvUKvSqdvb0y1LurXrQxrZvonOXqBG02Uf212SbpH0XknTkh62vTciDlftdq2kH0fEr9m+UtLnJH2o2XMDL6tZvuGhY5dJayVduEMi+IFXSaLmf4GkIxHxlCTZvkPSNknV4b9N0k3l53dL+rxtR0QkcH6A5RuABiVR9umV9GzV6+nytrr7RMQJSS9Jen0C5wbm1a7RD2BJSYS/62yrvaJfyT6yfZ3tCdsTR48eTaBpKIzqFTtZvgFYVhLhPy3p7KrXZ0maXWwf26dIOl3Si7UHiojdEdEfEf3r1q1LoGnodGOTM9qyc78kacvO/RqbnGEsP7ACSYT/w5LOs32u7VMlXSlpb80+eyVdU35+uaT91PvRrCdvH9bgno3znbua7+Qd3LNRT95O+APLabrDNyJO2L5e0rikLklfiohDtj8jaSIi9kq6VdLXbB/R/BX/lc2eF7j2mfdq5tiFkl6ZzCVJvc9066EsGwa0gURm+EbEPkn7arZ9qur5MUlXJHEuQJJ0YESzc2+p+9bsXCnlxgDthxm+aE8P7JxfyqHs5hPbX35evR1AfYQ/2tbQQJ+613RJkm4+cbkklnIAVoqF3dAeKuP4q4Z0Du7ZqMEu6dZTP6TP/nSb1vd0a2igj6UcgBUg/NEeHtg5P2u3zizea8t/AFaOsg8AFBBX/sivmsXadNPp848X7WAWL9Akwh/5xWJtQMtQ9gGAAiL80R4o8wCJIvzRHlisDUgU4Q8ABUT4I1+4KQuQCsIf+VI9tBNAyxD+AFBAjPNH9paazEVHL9AShD+yx2QuIHWUfQCggAh/pGu50TxM5gJSQfgjXUuN5jkwQo0fSAnhj/xgmCeQGjp80XqM5gFyxxGRdRvq6u/vj4mJiaybgaTVjuap/WKo4IsBWBXbj0RE/3L7ceWPbDHME8gENX+ki9E8QC4Q/kjXUqUcvhiA1BD+yA9q/EBqCH8AKCDCHwAKiPAHgAIi/AGggAh/pIPbMwK5QvgjHazbA+QK4Q8ABUT4I3mVEs+BkfKSDeWF3CrPKQEBmWNhNySv3ho9rNsDpGKlC7tx5Q8ABcSqnkjGcmv2s24PkCuUfZA8SjxAZlJZz9/2L0m6U9IGSU9L+mBE/LjOficlHSy//EFEXNrMedFexiZntGt8SrNzJa3v6dbQQJ8GN/dm3Syg0Jqt+e+Q9K2IOE/St8qv6ylFxNvKfwR/p6sq8YxNzmh49KBm5koKSTNzJQ2PHtTY5Ex27QPQdPhvk/SV8vOvSBps8njoBFVLM+8an1Lp+MkFb5eOn9Su8am0WwWgSrPh/8aIeE6Syo9vWGS/tbYnbH/X9qJfELavK+83cfTo0SabhjyYnSs1tB1AOpat+du+X9Kb6rx1YwPnOSciZm3/iqT9tg9GxH/W7hQRuyXtluY7fBs4PnJqfU+3ZuoE/fqe7gxaA6Bi2Sv/iLgkIs6v87dH0g9tnylJ5cfnFznGbPnxKUnflrQ5sU+AXBsa6FP3mq4F27rXdGlooC+jFgGQmi/77JV0Tfn5NZL21O5g+3W2Tys/P0PSFkmHmzwv8maRJRsGN/dqZPsm9fZ0y5J6e7o1sn0To32AjDU1zt/26yXdJekcST+QdEVEvGi7X9JHI+Ijtt8l6QuSfq75L5ubI+LW5Y7NOP82w9h+IBdSGecfES9Iek+d7ROSPlJ+/q+SNjVzHgBAsljeAau33JIOAHKL8MfqXTz8SshT9gHaCqt6Fhnr6gOFRfgXWZK3VmTVTqCtEP5IBjV+oK1Q8y8aOmkBiPAvnoQ7aVmuGWhPhD9WrbJcc2XVzspyzZL4AgByjpp/kTXZSctyzUD7IvyLrMkaP8s1A+2L8MeqLbYsM8s1A/lH+KO+pSaAld9juWagfRH+qG+pCWDl91iuGWhfjPZBUwY39xL2QBsi/PGKpSaASUwOAzpIUzdzaSVu5pKxpSaAsYInkFsrvZkLNf+8anTFTVboBNAAwj+vGl1xM8kVOqWlJ4CxgifQ9gh/1LdUHZ8aP9D26PDNk0ZX3GSFTgCrRIdvXjXaqUonLADR4QsAWALhn1eNdqrSCQugAYR/XjVas6fGD6ABhD8AFBDhX1RMCgMKjfAvqqQnhQFoK4Q/ABQQ4b8a7VoyOTBSng9QngxWed6unwfAqhH+q9FOJZPqYL94eH4iWGUyWOU5I4WAwiH8O107fVEBSA1r+6xUp62jw6QwoNBY22c18r6OTu0XVUW7flEBWLGVru3DlX8WDoy0NoQvHn7l+Hn/ogKQCWr+q9FsyYQ6PICMEf6r0U6lE2r7AOro3PDP29j1rMbYt9MXFYDUdG6Hb55r3XluG4C2lsrNXGxfYfuQ7Z/bXvRktrfanrJ9xHZ71yHy9osCAFah2bLP45K2S3pwsR1sd0m6RdL7JG2UdJXtjU2et740SitJdNZShweQsaaGekbEE5Jke6ndLpB0JCKeKu97h6Rtkg43c+662mWII3V4ABlLY5x/r6Rnq15PS3p7vR1tXyfpOkk655xzWt+yleq02b0ACm/Z8Ld9v6Q31XnrxojYs4Jz1PtZULeXOSJ2S9otzXf4ruDYi0uytNIuvygAYIWWDf+IuKTJc0xLOrvq9VmSZps85vK4IgeARaUxzv9hSefZPtf2qZKulLQ3hfO2Bp21ADpAs0M9L7M9Lemdku61PV7evt72PkmKiBOSrpc0LukJSXdFxKHmmp0hflEA6ADNjva5R9I9dbbPSnp/1et9kvY1cy4AQHI6d3kHAMCiCH8AKCDW88+hsckZ7Rqf0uxcSet7ujU00KfBzb1ZNwtAByH8E5REaI9Nzmh49KBKx09KkmbmShoePShJfAEASAxln4RUQntmrqTQK6E9NjnT0HF2jU+9HPwVpeMntWt8KsHWAig6wj8hSYX27Fypoe0AsBqEf0KSCu31Pd0NbQeA1SD8E5JUaA8N9Kl7TdeCbd1rujQ00LfqtgFALcI/IUmF9uDmXo1s36Tenm5ZUm9Pt0a2b6KzF0CiGO2TkEo4JzFEc3BzL2EPoKUI/wQR2gDaBWUfACggwh8ACojwB4ACIvwBoIAIfwAoIEc0d5/0VrF9VNIzTRziDEk/Sqg57YbPXjxF/dxScT/7Yp/7lyNi3XL/OLfh3yzbExHRn3U7ssBnL95nL+rnlor72Zv93JR9AKCACH8AKKBODv/dWTcgQ3z24inq55aK+9mb+twdW/MHACyuk6/8AQCLIPwBoIA6Ovxt77L9pO3HbN9juyfrNqXF9hW2D9n+ue2OHwZne6vtKdtHbO/Iuj1psf0l28/bfjzrtqTJ9tm2D9h+ovz/809m3aa02F5r+99sf7/82T+9muN0dPhLuk/S+RHxFkn/IWk44/ak6XFJ2yU9mHVDWs12l6RbJL1P0kZJV9nemG2rUvNlSVuzbkQGTkj6o4h4s6R3SPp4gf43/z9J746It0p6m6Sttt/R6EE6Ovwj4p8j4kT55XclnZVle9IUEU9ERGN3j29fF0g6EhFPRcTPJN0haVvGbUpFRDwo6cWs25G2iHguIv69/Px/JD0hqRA304h5Pym/XFP+a3jkTkeHf40/lPSPWTcCLdEr6dmq19MqSBBAsr1B0mZJ38u2Jemx3WX7UUnPS7ovIhr+7G1/Jy/b90t6U523boyIPeV9btT8z8Svp9m2VlvJZy8I19nGGOYCsP0aSd+QdENE/HfW7UlLRJyU9LZyP+Y9ts+PiIb6fdo+/CPikqXet32NpA9Iek902KSG5T57gUxLOrvq9VmSZjNqC1Jie43mg//rETGadXuyEBFztr+t+X6fhsK/o8s+trdK+lNJl0bE/2bdHrTMw5LOs32u7VMlXSlpb8ZtQgvZtqRbJT0REX+ZdXvSZHtdZeSi7W5Jl0h6stHjdHT4S/q8pNdKus/2o7b/NusGpcX2ZbanJb1T0r22x7NuU6uUO/WvlzSu+Y6/uyLiULatSoft2yV9R1Kf7Wnb12bdppRskfR7kt5d/m/7Udvvz7pRKTlT0gHbj2n+wue+iPhmowdheQcAKKBOv/IHANRB+ANAARH+AFBAhD8AFBDhDwAFRPgDQAER/gBQQP8P72dNaBisarUAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}