"
]
}
],
"prompt_number": 276
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"c) Modify your model by including two covariance functions. Intitialize a covariance function with an exponentiated quadratic part, a Matern 3/2 part and a bias covariance. Set the initial lengthscale of the exponentiated quadratic to 80 years, set the initial length scale of the Matern 3/2 to 10 years. Optimize the new model and plot the fit again. How does it compare with the previous model? "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Exercise 5 c) answer"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 277
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"d) Repeat part c) but now initialize both of the covariance functions' lengthscales to 20 years. Check the model parameters, what happens now? "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Exercise 5 d) answer"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 278
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"e) Now model Model the data with a product of an exponentiated quadratic covariance function and a linear covariance function. Fit the covariance function parameters. Why are the variance parameters of the linear part so small? How could this be fixed?"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Exercise 5 e) answer"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 279
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## More Advanced: Uncertainty Propagation\n",
"\n",
"Let $x$ be a random variable defined over the real numbers, $\\Re$, and $f(\\cdot)$ be a function mapping between the real numbers $\\Re \\rightarrow \\Re$. Uncertainty\n",
"propagation is the study of the distribution of the random variable $f ( x )$.\n",
"\n",
"We will see in this section the advantage of using a model when only a few observations of $f$ are available. We consider here the 2-dimensional Branin test function\n",
"defined over [\u22125, 10] \u00d7 [0, 15] and a set of 25 observations as seen in Figure 3."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Definition of the Branin test function\n",
"def branin(X):\n",
" y = (X[:,1]-5.1/(4*np.pi**2)*X[:,0]**2+5*X[:,0]/np.pi-6)**2\n",
" y += 10*(1-1/(8*np.pi))*np.cos(X[:,0])+10\n",
" return(y)\n",
"\n",
"# Training set defined as a 5*5 grid:\n",
"xg1 = np.linspace(-5,10,5)\n",
"xg2 = np.linspace(0,15,5)\n",
"X = np.zeros((xg1.size * xg2.size,2))\n",
"for i,x1 in enumerate(xg1):\n",
" for j,x2 in enumerate(xg2):\n",
" X[i+xg1.size*j,:] = [x1,x2]\n",
"\n",
"Y = branin(X)[:,None]"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 23
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We assume here that we are interested in the distribution of $f (U )$ where $U$ is a\n",
"random variable with uniform distribution over the input space of $f$. We will focus on\n",
"the computation of two quantities: $E[ f (U )]$ and $P( f (U ) > 200)$.\n",
"\n",
"## Computation of $E[f(U)]$\n",
"\n",
"The expectation of $f (U )$ is given by $\\int_x f ( x )\\text{d}x$. A basic approach to approximate this\n",
"integral is to compute the mean of the 25 observations: `np.mean(Y)`. Since the points\n",
"are distributed on a grid, this can be seen as the approximation of the integral by a\n",
"rough Riemann sum. The result can be compared with the actual mean of the Branin\n",
"function which is 54.31.\n",
"\n",
"Alternatively, we can fit a GP model and compute the integral of the best predictor\n",
"by Monte Carlo sampling:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Fit a GP\n",
"# Create an exponentiated quadratic plus bias covariance function\n",
"kg = GPy.kern.RBF(input_dim=2, ARD = True)\n",
"kb = GPy.kern.Bias(input_dim=2)\n",
"k = kg + kb\n",
"\n",
"# Build a GP model\n",
"m = GPy.models.GPRegression(X,Y,k)\n",
"m.Gaussian_noise.variance = 1e-5\n",
"display(m)\n",
"m.Gaussian_noise.variance.fix()\n",
"display(m)\n",
"\n",
"# constrain parameters to be bounded\n",
"#m.constrain_bounded('.*rbf_var',1e-3,1e5)\n",
"#m.constrain_bounded('.*bias_var',1e-3,1e5)\n",
"#m.constrain_bounded('.*rbf_len',.1,200.)\n",
"\n",
"# fix the noise variance\n",
"\n",
"#m.constrain_fixed('.*noise',1e-5)\n",
"# Randomize the model and optimize\n",
"m.randomize()\n",
"m.optimize()\n",
"\n",
"# Plot the resulting approximation to Brainin\n",
"# Here you get a two-d plot becaue the function is two dimensional.\n",
"m.plot()\n",
"display(m.add.rbf.lengthscale)\n",
"\n",
"# Compute the mean of model prediction on 1e5 Monte Carlo samples\n",
"Xp = np.random.uniform(size=(1e5,2))\n",
"Xp[:,0] = Xp[:,0]*15-5\n",
"Xp[:,1] = Xp[:,1]*15\n",
"Yp = m.predict(Xp)[0]\n",
"np.mean(Yp)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"creating /var/folders/22/6ls22g994bdfdpwx4f9gcmsw0000gn/T/scipy-neil-mcRzlh/python27_intermediate/compiler_08f57918657bc3b6af9ef49d104a532a\n"
]
},
{
"html": [
"\n",
"\n",
"\n",
"Model: GP regression
\n",
"Log-likelihood: -74285.8774977
\n",
"Number of Parameters: 5
\n",
"
\n",
"\n",
"\n",
"\n",
" | GP_regression. | \n",
" Value | \n",
" Constraint | \n",
" Prior | \n",
" Tied to | \n",
"
\n",
"| add.rbf.variance | 1.0 | +ve | | |
\n",
"| add.rbf.lengthscale | (2,) | +ve | | |
\n",
"| add.bias.variance | 1.0 | +ve | | |
\n",
"| Gaussian_noise.variance | 1e-05 | +ve | | |
\n",
"
"
],
"metadata": {},
"output_type": "display_data",
"text": [
""
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"WARNING: reconstraining parameters GP_regression.Gaussian_noise.variance\n"
]
},
{
"html": [
"\n",
"\n",
"\n",
"Model: GP regression
\n",
"Log-likelihood: -74285.8774977
\n",
"Number of Parameters: 5
\n",
"
\n",
"\n",
"\n",
"\n",
" | GP_regression. | \n",
" Value | \n",
" Constraint | \n",
" Prior | \n",
" Tied to | \n",
"
\n",
"| add.rbf.variance | 1.0 | +ve | | |
\n",
"| add.rbf.lengthscale | (2,) | +ve | | |
\n",
"| add.bias.variance | 1.0 | +ve | | |
\n",
"| Gaussian_noise.variance | 1e-05 | fixed | | |
\n",
"
"
],
"metadata": {},
"output_type": "display_data",
"text": [
""
]
},
{
"html": [
"\n",
"\n",
"\n",
"\n",
" | Index | \n",
" GP_regression.add.rbf.lengthscale | \n",
" Constraint | \n",
" Prior | \n",
" Tied to | \n",
"
\n",
"| [0] | 3.00660818886e-08 | +ve | | N/A |
\n",
"| [1] | 0.511855770259 | +ve | | N/A |
\n",
"
"
],
"metadata": {},
"output_type": "display_data",
"text": [
"\u001b[1mGP_regression.add.rbf.lengthscale\u001b[0;0m:\n",
"Param([ 3.00660819e-08, 5.11855770e-01])"
]
},
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 24,
"text": [
"66.628745898828541"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD7CAYAAABOi672AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VHW+//HXmZpkkklPCBAIHQIoRZoIRAXFXtG1rLu6\nKorrqut1vSq7ovfa1nZ37auCZX/iKgroqigooS5SRAHpNQk1fVImmXZ+f5xEkpBITDI5OSef5+MR\nJ/lmhHeG5J3vfOd7zgEhhBBCCCGEEEIIIYQQQgghhBBCCCGEEEJ0EEq4/uBJkyapy5YtC9cfL4QQ\nZrUMyGrsE5aw/Y3LlqGqqmHfHn74Yd0zGPlNHj957OTxa9kbMKmpXg1bYQshhGhbUthCCGEQUthN\nyMrK0juCocnj13Ly2LWOmR+/sL3oCKg16zFCCCGaSVEUaKKbZYYthBAGIYUthBAGIYUthBAGIYUt\nhBAGIYUthBAGIYUthBAGIYUthBAGIYUthBAGIYUthBAGIYUthBAGIYUthBAGIYUthBAGIYUthDCF\ngqDKDcV+Li/ykxc054nnpLCFEIY33xvklPwqyiKPYnEdYkRBFXMqg5jtjKFyelUhhGEVBFXu9AT4\nLlRJ7/gc0i0enITYpropKEmnWyiGf8TZ6GYNZ9W1LTm96i9UEQB/SO8UQoifUzurrow8QlrCDk61\nHCCVnaSyhxHKPpLiduGMOcjwfPPMtqWwG/i2ELp9AmMWQ06F3mmEEI3xqiq/LQkwMsrL3ogjnK14\nKaKARDYSxTq8lDFJ8bDLkc+Z7jJuLwlwxASTMCnsGqoKL++Ci1bAW6Ph2p4wegl8eVjvZEKIhiIV\nhU3JDqp80ThKerNUjSWCXvzAKXzPqTjoxhI1Ebcng/0VcWxMtpNmoGWRpsgaNtoSyPT1sLkEXhgC\nD3wIPRPghsnwuw0wvQ/MzASL8f+9hTAVVVX5R2WImeXVTIg/ylH7UQKodAsms6IwjfuinPwx2opN\nMc4P78+tYXf6wt7hgStWwWkJcHk03Poe3JMFu/Jh9T549Xp4aDfE2OHdMZDo1DuxEKKhAwGVm0v8\n5CqVWBSICUbxVpyNQXbjLSJIYTdhRb5W1o8NBf9BeGwxvHcDfPAxDOoFUelw/yfw3m/hSy98lAfL\nz4L0KL2TCyEaUlWVf1WF8KtwTaTFULPqulqzS2Q2cBTYXGdsFpAHbKx5m9rqhDpZfARu7QO39IEF\nm+G1q2FSP3h5Hrz6Edw0Du6cCF/vgGeGwSA3bCjSO7UQojFVwCa/yjq/SmXHniu2mO0kn58DvAC8\nU2dMBZ6reTM8R51fWQ4rLFyqvX+kAL7dBA4bVPpPvK8ZVAfhrX2w9Jj2rGFGX+gVrXeq8KnCwy6+\nooQcEuhFP87BgUvvWOFTthP2vgpVhyFlCvT8NVjseqcKizW+ENeXVJAcdQCHxU9mfl9mx7o4J8Jc\nP7QnK+wVQEYj48Z8rnESsxfAv+YCTigsgdNvgKtvgR7peidre6oKF66AJUePj/1jL6w6C4bE6Zcr\nXKrwsIj/poJjAOSwmr0s5VyeMGdpF62DZVkQrNQ+zn0fDn8Cpy/QNVZb86oqf/EEeafaQ//E9URb\ny7HhJyOigJuKTuX8qniecdtwm2THQEt//dwJrAHuB2LaLo5+fH6Yt7j+WCgEX67Sys1svjpSv6wB\nPH54bJs+ecJtF1/+VNa1PBxiD9/olCjMts46Xta1Di2EgtW6xAmHNb4Qw/N9fG/Lp2fyFvpa99CT\nA/QglwGWnXRL3Eqe8whD86v5qsoEm7A5+Qy7Ma8AjwJu4GlgOvBMY3ecNWvWT+9nZWWRlZXVgr+u\nfVR4IRg8cbzIA/5A++cJt82ljY9vKmnfHO2lhANNjOe0c5J2Urq58XHPZkg6vX2zhIFXVTmr0M/k\n6Ar2ufI4nY34cNCf3QD8yCCGKxtZFenidNXJ1KJYDqY6OuRe7OzsbLKzs5t135YUdu00pRR4CXiZ\nZhR2R+eKBJsNGnZzYhw4TLjsN6yJZY/h8e2bo73E04sc1jQyntH+YdpD3DDw5p44Hjus/bOEQaSi\nsDLRzm9KXCSRzsbo4QxStrKFTFQsKKhsVIfT1duVtWVuliZ23ANnGk5mH3nkkSbv25IlkbSaWxtw\nLfB5C/6MDsdhh1812O9itcLU8frkCbezU+G8tPpj8Q54KFOfPOHWj3OIoUu9sVh60IezdEoUZpmz\nwNbgFeTu0yBxjC5xwmGEw8KGZAdnqEnsyx/K3kAf9qi92UcGO0P9yCkYwgB/KpuSnUxymuPFx5PN\nsOcCk4AkIBd4GMgChgE+YDnaEokp/PpC+PV4OPcuSIqHxS/Colwo8eqdrO0pCiw8A+YegG+OQY8o\nbYtjd5PuMXcSw1SeYjdLanaJ9KYPZ2EnUu9o4RE/AqZsgr2va7tEUs+B9Kv0TtXmHIrC/7htXBoR\ny/VFY4mIOIxNCZJf2Z3341ymKepaJyvsaxoZmx2OIHqprLMGUumHS2uW95LiYdhA+HiP9rGq1r+v\nGdgtcEMv7a0zcOAik0v0jtF+XL1g6ON6p2gXIx0Wvk+O5vny3lSpcF+yFZdJdobU1ZI1bNO4oCtc\nthJCwLQRcMtcKCiHx2dAajI88AnM+RbeuAGuXA2HvDA6Ue/UQojG2AG3RcESUnGYr6uBTn5oOkB+\nFdy+AbZ54C994K8LwOsHXxCGdoWLJ8GDW+GGDHh0CDiteicWQjS0O6ByU4kfn72MCCVESZWbt+Pt\nnCrnEmk2QxQ2aMsd7+fAXRvh1t6QUantGFnih2+L4O0xcHqS3imFEA2FVJUXKoL8b4WPCQmHqbIV\n4AQIJrC6KI07IiN4INqK3UDnFZHCbqZDXrhlnXZSKLsFrusJT5wCrk69cCREx1Q7q662e4h05zJc\nKcHHPiwoRJDOt2o8Snl3Sr1xhpptS2H/AqoKy/LBaYFxMqsWokPyqiopR3xkRVewP2YPV1FBLvsZ\nyDoA9jKWVFL5gHgyK3vycUksOR30wJmGpLCFEKbzeVWQ6SV+hsUW4Ik4xDilgGoOoaLiJJX1aipu\nXxqri5P5P7eDX0VaasuwQ5PCFkKYUklI5Z7SAMuDlQxMyCPOUooVOKJGc7C4BwNw8UqsjVQDzKxr\nSWELIUzts9rZtrsIm6Lyn9IEQ82q65LCFkKYXklI5X5PgEoVnnbb6GKgWXVdP1fYsv9BCGEKsQpM\ndlooVyHFGBtCfjGZYQshDO9IUOX20gC7FO3AGXfAzZtxdnrZjDfLbs01HYUQosNSVZX3KoMMK/AS\niM6jW9wO0mJ3ER2bw5jCSl6pCBIy0cRRZthCCEOqO6tOicthsFIE7MMCWOjBd2oiFZ503P5YQ822\nZYYthDCVkpDK4GM+go5ibHF7OF05hItD9OZ7evId8RxjrHIA3AdIdBUyNN/HvoDxJ5BS2EIIw4lV\nYGaMlTUVLjKCsWxX46jEgp1TiWQ4RwlyiBS6hqJZWRbNvS4rPUxw4jZZEhFCGNaP/hDXFfuJjyom\n4MpjvFKEQpDVaiKR3m7sLUvi3Xg7YxzGmZvKPmwhhGn5VJWHy4K8U+1leEIufiXEzqJ0ptqieMZt\nM9yFDKSwhRCmt7I6xJXFfjwqzIu3cX6EMddApLCFEJ1CsKZzrAY7HL0uOdJRCNEpGLmom8M4K/FC\nCNHJSWELIYRBSGELIYRBSGE3IhCCkLxeKoToYKSwG1h0GHr+G077CraU6J1GCCGOk8KuUeaH6etg\n+np4ezTM6AtZS+GpbRAM6Z1OCCFkHzYAy47BjWshKwWeHwaxDm18f4U2Xh2Et8dAvxh9cwohzE8O\nnDmJAZ9D/xhYeAY0PIq13A8jvoKRCTB3nD75hBCdh5xe9SSyz4RSP1y4Agqrj4/v8MDYJXB6Erw5\nSr98QggBMsP+iT8ED26CeXnw4enacsiMDfDYULi5N5j8ACohRAchSyK/wEe5cOVqSHHC5xO1pRAh\nhGgvUti/UFVQuzXoyb6EEAYmJ39qpkAI/rYO5uyAGCv8YRhcM1jvVOFTWApP/j9YuhHSU+Deq+GM\nU/ROFT7bq+DJfNhcBcMj4IEU6OPUO1X4bCOPJfxACRUMJp2pjCAK837Bbx3T3vwq/CoJ7uhy4iYC\no5MZdh2/Xgj/PAjEoj0y+fD8OLh7tM7BwsAfgBG/gy37jo/ZrPD18zBxmH65wmW/D4bvgpLg8bFk\nK/zQH9Ls+uUKl80c4O98hsrxn8EMUniQK7GE9cdeH/+TC3/JrT82owu81FufPK0hu0SaYVMB/LMQ\ncAF5wGEgER7aou3DNptPVtYva4BAEJ5+X5884fZyYf2yBsgPwhtF+uQJty/4rl5ZA+znGNvIbeL/\nMK7qEDx76MTx149Cvr/984STFDbwfTGctxLwAweBIFAN5EGlAmcthcNeXSO2uX2Hf9m40e3zNT6+\nv4lxoyvA84vGjawoAKWNTKr8Khw02b+vFDZw2So4Ixliyhp8IgTDVSgPwB0bdIkWNhNP/WXjRjfR\n9cvGja4fXX/RuJF1sUP/iBPHk20wKLL984STFDbw2QTYVArDhoBS5xFxRQBpkBoBr52mW7ywGJ0J\nt11Sf6x3V5h5gz55wu3mBDgjqv7YlGi4Jk6fPOF2KWOIJ7re2LkMpyvm26eqKPBib4is87NrU+CF\n3uA0WcPJi441yvxwyzrYXAIXxkC1Ch8Uwy194C+ZYDXZP3yt1Zvhm++gRypMOxMizbuJgKAKn3pg\nSxUMj4TzYsy3i6CuKnysZzclVDKYdHqRqneksDrkgw8KIKDCFYnQq5FZtxHIPuxmUlV4cRf8YSO4\nbNoRj+el6Z1KCNGZSGH/QgU15xNJMvFsUwjRMUlhCyGEQcg+bCGEMAEpbCGEMIiTFfZs4Ciwuc5Y\nDLAQyAEWQIO9Q0IIIcLiZIU9B5jaYOx2tLLuh3YQ921hyCWEEKKBkxX2CqC4wdho4E20g7dnA2PC\nkEsIIX6RveVw8Qo4eyn8WKp3mvBoyRr2KGB7zfvb0QpcCCF04QvCE1th9GLtcn5XdIespfDAJqgM\n6J2ubbWksE18bJgQwkhW5sPwr2BlAcwZAh9/Bi99CK/11y7zN2QRfGGiE5q15AIG64BBwMaa23VN\n3XHWrFk/vZ+VlUVWVlYL/johhKivKgi/3wCLjsBzw2DVWrjlS3jqYoiJgLveh8kD4PmJcOd3MDwO\nXh8FcQ69k58oOzub7OzsZt23ObPlDOBTYGjNx38C0mtunwH21dw2JAfOCCHC4uujcNd3sHoy7DkC\nl74B3/0JEmvOvlhWBRP+Bk9eBJP6w9TlcGtvuC5D19jN0poDZ+YCq4H+QC5wI/AK0APYAXQDXm2r\noEII0RyqCl0iwW0HFa2oE+ucKjcmAtLc2ucibZAeBWaYPp5sSeSaJsYvaWJcCCFEmMiRjkIIYRBS\n2EIIYRBS2EIIYRBS2EIIYRBS2EIIYRBS2EIIQ6p7mEdjh3yc7PNGJIUthDCcXtGwuRQe3gIpbjhS\nBvfOB68P/EF49AvYkAsZifD8DvjyCAx065269eQSYUIIQzpYCbeuh4NeeH4I/GMxLN2lHSAzojs8\ncBE8uAMsCsweBX1j9E7cPHJNRyGEKakqvLMf7vsBZvSFMVZtbLsFHt8GfxkMv++nlbZRSGELIUyt\ndra91aOt83aLMtasui4pbCGE6akqLDgI/hBcmW6sWXVdP1fYLTm9qhBCdDj7K+Dl3VAegEw3DInT\nO1Hbk10iQghDC6nw0i4YtRjO7QI39YIzs+HxrRAI6Z2ubUlhCyEMa185nJ0N/zwAH4+G1Svhw0Xw\n8UjIPgbjlpjr+o5S2EIIw1HrzKrPT4PbYuDKl2BoV5g6CC57GS4Eftdbu77jEyaZbcsathDCcFYV\nwFPbYcXZQCVMfBu+mgHD07XPXzAYpr4Cc66D9VPgguXQPwauSNc1dqvJDFsIYThVQa2AB7nB64f0\n+ONlDTAgFQamap/r6YJh8eAN6pe3rUhhCyGEQUhhCyGEQUhh11GOn2fYyoVkcwUrmMMegqa4dGfj\ndu6FS6dD7CkwZCr8c4HeicLr6+9g/J3gvhAm3Q2rtuidKLxe3wODPoe4j2HaKm2fslmpKjz0KnS7\nGFLOhz88B0ETLIE0JC861vFnNrGeopqP/MxmL9WEuI1+uuYKh4pKOPM6OHRU+/jHnfDrP0JUBFw+\nVd9s4bBpD5z/APj82sfLN8GU++CH16Ffd32zhcPb+7RDtWvNy4MNxbD9PHBY9csVLnnH4PGFxz9+\nYR6kjgHO0i1SWMgMu8YBKuqU9XELyCOACfYDNfDRouNlXdeL77Z/lvbw2r+Pl3UtbzXM/kKfPOH2\nwq4Tx/ZVwGeH2z9LuHmrodBz4vjRIigoaf884SSFXaMUf6PjFQQImHBZpLC48fGCE39nmUJjP9AA\nBSY6qKKuQl/j4wXV7ZujPXirQG1iTlVW2b5Zwk0Ku8Yg3MTjOGF8BPFEYL7nkOef2fj4BU2MG90F\nY5oYH9u+OdrLBWknjlkVmNrIuNElxGpLeQ1FOKFnl/bPE05S2DXsWJjJEFx1lvXTiOBeBumYKnwG\n9IZnHgRbnVcxJoyCB2folymcrj0brp9cf2z6RXDJeH3yhNujQ2BUwvGPHRZ4cQSkR+mXKZx6dIGu\nScc/TnDDqX3AYrKGkxcd6xhNIh8zgXUUEoGVkSRgM/HvtHtvhmsughXroEdXGDdC70ThY7XCuw/C\n/dfAln0wvC8M6KF3qvBJcMLaKbAiHw57ISsFUhqZhZpFpBN++AgWrwV/AM4ZA1fM0TtV25PCbiAK\nG5NI1TtGu+maCldfqHeK9jOkl/bWWUxI1jtB+3HY4QKTPmOqZd7poxBCmIwUthBCGIQUthDCcMy3\n0bZ5pLCFEIZk0Es2tooUthBCGIQUthBCGIQUthDCtFSTLXZLYQshTEsx2UK3FLYQQhiEFLYQQhiE\nFLYQwjQqq8BTrneK8JHCFkIYTlkAXDVnQvpmJ/RL1s6LPXkGjLsJij3a2Dc7tfu4rOBp/JT3hiKF\nLYQwnL3l0NsFpV7469cw8xy49s8Q9EK8FS79L7h7IsxeA3nF0Ctau+KO0YXzNVRVNdueGiFEh3DH\nBhjkhqM7tUJ2HYJla+HYSpWICJWkURb69oWMsVDshcmnwwe5MM8AZ/NTtK0tjXazzLCFEIaztxzi\nFHh5BXTxwqKVcHCZit22E2/lVnKXq2zZAeW7YP4PYPVp/4/RyfmwhRCGs78CVm6FMzNg9jzwbYTM\n/odwucqwWv0cPryPnYt7U+GHcy+CT9bDfqfeqVtPClsIYTh2C8RFQXkE+PwwZjRYVTcVFaWkpNhQ\nlDickbC3CmJc4IwCexMX6jWS1hT2fsADBAE/MLotAgkhxMnEOWBQGny4Cq6aDFF2ePepaMrK/Ph8\nXpKS+nHrQ/BjHvxwFG4dCl/k6Z269VpT2CqQBRS1TRQhhGieODvExkC5D6ZdCtf/Ce7+o8KypaeQ\nkqKSkqLw4Qp44X644gPolQJxR/VO3XqtfdHRZEfqCyGMIM4OpX7I6gt5VTByECT1h6+/VnjvPQvj\npkKsC4IuOK0HVKoQ79A7deu1prBV4BtgAXBx28TR37Ggym+Kfdzv8VMl2xKF6JCibVAegBHdYcth\nuG6qdsX0Y8cgPx+yN8K1U2HrERjeHSrqHGhjZK35EsYDh4FBwKfAWuBIW4TSy8feIHd7i+kRt5Od\nIScji/rxdkw0pzlk96MQHUntWfhsVu0UqhEO7TYpSRtX0cbKAVvNj68ZlgNaU9iHa263AZ8AFwGv\n173DrFmzfno/KyuLrKysVvx14VMcUplRWs0mZx7D439koLIDv2Jjc0I+V5YN5oaqLsyMseEw27ka\nhRC6y87OJjs7u1n3bWlhRwFWoAxIBs4Fnm94p7qF3VF9URXk9opS0uO2M8GylUxlG/3YS5ViJYVi\nusYU8lUgk0+L+vKO28VQu8y2hRBtp+Fk9pFHHmnyvi0t7FRgfs37hcCzQG4L/yzd/LUswCuhfHon\nbGWwsoXBbCGTbQziR4K4iKEcl+IhwlbOpoQKzi7MZE50LBdEWPWOLkSn1/AVps7wilNLC3sfMKwt\ng+hhlMPC82VRWIgmQBw+EvGSQgmF+LFRRjyVJOJT4nEHY4kMORlkkxm2EB1BwwXKzrBgaYLXTVvu\nTKeF94jh6mNDSU5ys84aQykxHMKJioOdDGEnmUT5BrO9uAvfJDrpbesM3xZCiI6oUxc2aKX9YVwE\nV+b34qZEN9/ao8kjGgUb+xhMpncICzwxZCc66CVlLYTQkTy/ByY5LXwc72BOYRLn+cawltEcZBzj\nKk5joSeGpVLWQogOQAq7xgSnhWkRFt4rjeQVMnlSHcD/eew8GG0jQ8paiA6llwte3g0hF3ywEQ76\nYMX38M/P4dPlsHAZVNrh9dUQ5Yant0PvaL1Tt55cwKDGU2UB3qgM8pvkY6y1bKMaO1MDQ3iiwM3s\nOBsXys4QIToMVYW5OXDPRjg7FlYug5GpsOUbKCuHrEshOxfOOQu+8MCjQ2B6H7AYYO71cxcwkMIG\nHisL8I7Xx+jkvViVbfRXd1KBjV1KJsnBTD7M78rrsQ4uiTRnaed7wB0JTrveScIvGITCQu2IOEsn\neH5Z5YMyLyTH6p0kPI544bYNsMsDAzzwzXfgD8LEoZCfBrEOeHM0ZLj0Ttp8csWZn/FoWYB/+svI\nSP6BGP7DqOASTgstYkxoCWNC3+CzrGFsynZuK69kvjeod9w2tWI7DL0fUm6D1Nth1jy9E4XX7NmQ\n3gNSu0BGL5g7V+9E4RMKwf2zIfkaSLkGRtwJ63bqnartdYmE+ePh7gz4vBpKS6CyAL5W4MokWJxl\nrLI+mU49w36+PMAbwSK6ubcyUPmRzNB6MtnCQHYQUG1sZRDblVP50TKSbWomuwoH8lZ0LFMijP97\nrrAMMu6C8qr646/9Dm49W59M4bR8OUzKqj9mscDab2HkSF0ihdUzH8F9b9YfS4iBA29BdKQukcJq\n2lUwbzFwIxAJvAETBsLyZToHawGZYTdhSbXKyOgyxishRgaD9K/aSVr1UZxVfpxqNbFBD119uQxS\nnYxWFM6M9rLMZ4LLVgDz1p5Y1gBzDPgN3hxvvXXiWCgE777b7lHaxZzFJ44VlcEna9o/S7iVlsL8\n+UAikIZ2HHYSrFgBe/bom62tdfp92ApgQcGCilL7pmrjiqqiqEEsaCdOMcDrFc3mDzQ+7mti3Oia\n/Hp97ZujvXSmf99KPwSvAM4ClqJdA+v3wDrtiulm0qln2PVYYkGJOnHcmgCK+R6my0aBvZHXUK8e\n1/5Z2sPVVzUxfnX75mgvV088cSzKCReNaf8s4bShCM5ZDwmnAcuAAOADVkNMOly1H7KP6RqxTZmv\niVpKUcA5on5pW1xg76VfpjDqlgBz74Rkt/ax1QI3ToI/nq9vrnC58EJ4ZBZE1qzfulzw9F9h0iRd\nY4XNzGvgmknHd8J0iYcPH4REt7652kp1EGZuhvOXw3AglA9J5cAu4ADEHQN7FZyuwvVr4M4NUO7X\nOXQb6NQvOl5Q6CcpLoe+1hwyyKE7OaSre0lRjxJQdrFXGUkumRykLwfpTV5VD3r6Uvhft3lWkqr9\nsCUXusZDWrzeacKvuBj27oV+/cBtkvL6OQcL4GgJDM0Au3m+bXloE6wphPNt8PYquHcUzHwOnrsD\nHDb479dhxrXwwQGYOAC2xWiXCHvTAJcK/7kXHU30T9hGFEBxA51gUzLa3uuRvfVO0X7i4825K6Qp\n3ZK0N7Mp9cNl3aE6F6YMALcFTsuEaTXPEBd9B7YgXDwUjpXBtUPg/Rx9M7cFWRIRQgiDkMIWQgiD\n6PSFHVIhiEqQUM1/A6hUouIjhJ8gfgI1o+bYgS2EOVTUbFGsqNmaWeHVzjFS+z5ARXX9+xpdp17D\nvjjCwrPlMRyKtVGoVOGlAAuHiWUrfmIoIcghfGzHz3bVwo5yFx/EdPrfcULo7oru8Kv/wCUp8OV2\nKKuEg0Vw9gztaunbcsDeGz77Hi4/D+79Hv4xSu/UrdepC3u6y4q3PIEXi07BkeDEr0Tix0UZcfiw\nsZtB7GAgxaFBbC/ozTvuKCY4pbCF0NuZqfDDudqJn1yjwOuBsr4wNgkIwe5EOKRC9wmw3QfrpkAv\nOb3qz+rw2/pqvVAe5NnKas5OysFn2UoXthPCzl4G0zMwmPcKUnkv3sFkKWshOpR6p1mNg8WLweOF\nqy6ARWXGOq1qLTm9ajO8VBHkr+V+7kgqZL51J3asXO4fwmOFLv4Vb+csKWshOqza06xuKNJKrW8M\nvDnKmLNqKexmerEiyBsVQRanQHlIZfhRWJBgJ0vKWogOT1XhxV1QGYT7BhprVl2XHDjTTJdGWHii\nLEAyToKqilPxSVkLYQDVQXj0R3hjL9gtkFMJT50C0SY7/k3aqA6V47/WLDUfCyE6tvVFMPIrWJcP\nffIgeRfkeuCUL8114ieQwq4npB4vbAVk37UQHVh1UDunyAXLYQSwYQkMcMAIF6xZDKeHtBM//d4k\nJ34CWRKpR619U9WaW50DCSEatbUUrloN6RHQKxd2+WCkD9YuAVck9IuCI/shvhoOxWqz7bnjYEyi\n3slbR2bYdaRaIdmiMCLfz8RCP1NNcCkwIczo5d1wXhpMDkG6C+4dDmXFsOFdWP0mJETA5WkwsSf0\nLYEZfeGZ7Xqnbj2ZYdcRqSisSbYzuzKEW4FrI6WwheiIVLSL61YVQY94beaZlggRTu3z6anabUai\ndra+DJd2Olajk8JuwKko3O5q5FIsQgihM5lCCiGEQUhhCyEMqXZPQEjVrvAXrLOtKxjSxkI1t2bZ\nPyCFLYQwHJsCgRAkREFRJfRMg32Hjn9+70HISIPCSkiM0u5rN0HbmeBLEEJ0Ng4L+ELQxQ1HPDAw\nA3blQqDmvNdb98HgPtrnuri1+zpM0HYm+BKEEJ2N01pT2DFaKUdFQNck2JMHJWXgqdB2ipitsGWX\niBDCcGpn2GmxcNijjWX20mbWhaXa+4qifa6LG3ZXSmELIYQubAp4VW0Nu7hSG8voCjlHwOfX1rQB\nSrwQHwmiPHJFAAAJc0lEQVSBCrAa9Ox9dZngd07bqgipTC/x80RZgJAcmy5Eh/TTOX/qlHDt6VRV\n6p9atfY+JuhrKey6dgZCjC3wU6HCZ9UhLi4KUBSS0hZCdAxS2DU+8gYZX+DnTpeVd+NsLE2009+m\nMDLfx3qfnLdPCKE/WcMGni4P8FJFkC8S7Ax3KLzMftKJ5LnYLpzuUDivyM/bcTbOj5BD1oUwErOt\nanb6ws4LqjxZFuSHFAcua4Dp/IdeLOQQLtYxjT9HDibOYue2Ej9TnBbsihlWwoToPMz0E9vpC/vR\nsgC3uKzkW8v4E59zHvOYXLIUr93JZ65D3Mw0nnJm0dum8GZliNvkxFBCCJ106sLeEQgxvyrIoylH\neZ2P+a36EeN3rCXqX35ww/VXv09q2hHuVw5xlfsiHi6M4YZIC1FGvbqnEMLQOnVhP+TxcUbSJo5Z\n5jHdv5DBK3ZjeTvEgc/AGQmpORVcfMMiupx6lIX2IwyKu4y/VfThgRiTXdlTCGEInXaXSGVIZaW1\ngO7WVZzBWvrn5GLNhvLlFhYVwqrDVvwr7TiXKowo3MIo1pLmXMs7AY/e0YUQTTDbi4wNtaawJwLb\ngF3AnW0Tp/0EAZslhE0J4UDBEoxAwUYopKACVVYrqiUCJWgFFRxYiUAlaLJL83pVldW+EDkBk3+n\n18gvhFVroahY7yTtY99B+M8mqPbpnST8FAUOHYV9B/ROEj6tWRL5GzAdOAB8CcwFCtoilF5UVUUN\nBgEI+nz4KgI4VRvmep35uLmVQWaUBiipuVr8NZEW5sTZcJh0J8xDT8Izr4LPBxERMPMueOguvVOF\nh7cKrvszzM/WPk6KgzdmwiWTdI0VNgfy4L35kH8UCEFMEvy4Q+9Uba+lM+zYmtvlaIX9FTCmTRLp\nRYWg319vKBQMEfQHdAoUXjkBlRtKtLIG7XDe97whni0P6porXBZ8AY//XStrgKoqmPkUfL1C31zh\n8ugbx8saoKAEfvUQ5Jv0mcWN92jPnmqVlcNV0/XLEy4tLexRQN1rEG8FxrY+jn5CwWCjl6VQQyHz\nXK6ijgVVQRr7VfRhlbmWfGp9+O/Gx+d91r452suHX584VlUN/zbpL6ilq04c27rz+PmxzSKsu0Rm\nzZr10/tZWVlkZWWF869rHXOuAjQpsollj4h2ztFeIpv4wmqvsm02kU18XU09DkZns9HoBMQIq3vZ\n2dlkZ2c3674tLex1wNN1Ph4MLGp4p7qF3dFZrFYUi0LD6bTFajFlmV8RaeE+D5Q2ePZws0kPDLrx\napj9fv1dBFYr/GaafpnC6XcXwz3P1x9LjoeLJuiTJ9yuvADez64/duZ42GmAb+eGk9lHHnmkyfu2\ndEmktOZ2IpABTAG+beGf1WFY7PX3V1vtthPGzCLBorAo0c5Iu/bbKNECT8RYuSnKAN/hLTB+NLz7\nAvTsrn3cJwM+eA2GDdE1VtjcdQ385WaIi9E+HjsUvnwBXJH65gqXfzwNA/uC1QaKBbqlwb9e1TtV\n22vNksjdwGuAHfg7Bt8hAtpE2mK1AkFsTif2SAcK5t0PNdZhYX2yA09IxaWA1QjPH1vhusvh2svA\nUwbuGGM8XW4pRYFHpsOffweVVeCO1jtReMVEw7lZcNM0SImHL/4DyYl6p2p7rSnsZcCgtgrS3uwK\nBEI2qlQnFYpC0FGJavVjdSg4gajqaiyKH9UJQWsEXkJ4sBOpmu/gUHcnOtReUSDWrXeK9mOzmb+s\n67LbwGnS1yWgEx+aHqEoTCOZvOpzWRRRQUSPCsZe+h1RngBXfwY4wX5JCM8lkSxNOIOv1anklE3i\nnohO9NMuhOhQOm1hA8yMtpF5rAevp9zCG9Y0jp76EZPvW0riqWWoDsg7L5nPYs9jBVcypfpMvq1y\n8OvkTns0vxBCZ526sFOtCne4rHzqieHZ+F9xP2kc7tKFydctoYoI/m29CA/TeE09lTNKQzwWYzX9\nOq8QouPq1IUNcG+0lX7HfNzhs/O642wepysvWQfgw0kWF3M93ZnjDRGhwCURMrsWQuin0xd2rEXh\nlVgb5xX5+bvbxsNRg9lBL2woZKgRPFIW5PXKIP9OtKPI7FoIoaNOX9gAV0Ra6WNTuLLIz2p/iGfd\nkXhCcH6JnyoVNiQ76GKVshZC6Eue49cYZtf2JOcGoetRHwPzfQyzW/g60S5lLYToEGSGXUecRWF+\nvI35VSFiFIUpsmYthOhApLAbUBSFyyPNeXi2EGYRbYOFB+FwH3DY4P0N2mH3H2dDlBO6p8IXW6HE\nC5UqfHEYRsbrnbr1wvlcX1XNfr0eIYQufEF4Yhu8tBtu6wYffgUDkqFnFfh94EmCVfvhxvPh1cNw\ndQ94bChEGWCKWrO5odFulsIWQhjW98Xw27WQ5oQ+pfDGMqgOwI3joawbfF8Kc0bDGcl6J20+KWwh\nhGnVnW3fkwExNnh8j7Fm1XVJYQshTO/7YrhuDZT64f1xxppV1yWFLYToFGorx8jHuP1cYRvsyYIQ\nQjTNyEXdHLLRWAghDEIKWwghDEIKWwghDEIKWwghDEIKWwghDEIKWwghDEIKWwghDEIKWwghDEIK\nuwnZ2dl6RzA0efxaTh671jHz4yeF3QQz/6O3B3n8Wk4eu9Yx8+MnhS2EEAYhhS2EEAYRzlOlZAOT\nwvjnCyGEGS0DsvQOIYQQQgghhBBCCCEMYRaQB2yseZuqaxpjmAhsA3YBd+qcxYj2A5vQvt/W6hvF\nEGYDR4HNdcZigIVADrAAiNYhl9DBw8Af9Q5hMBvRSrsnsB1I0jeO4ewDEvQOYSATgOHUL+w/AS8A\nTuBF4L90yBUWsq3v5Ex+0aE2FVtzuxw4AHwFjNEvjmHJ91zzrQCKG4yNBt4EqtFm4Kb5HpTCPrk7\ngTXA/WhPtUTTRqHNqmttBcbqlMWoVOAbtKfyF+ucxajqfh9uRytwU5CL8MJioEsj4w8BrwCPAm7g\naWA68Ez7RROd0HjgMDAI+BRtHfuIromMR56hCE4FVukdooOLRVvDrvUCcIFOWczgOeAWvUMYQAb1\n17A/QlvXBhgJzGvvQOEiSyI/L63m1gZcC3yuYxYjKK25nYj2QzQF+Fa3NMYTxfFlt2TgXGCRfnEM\n61vgJiCy5naNvnFEe3kHbYvVerTZjrx6f3KT0Lb17Qb+oHMWo+kFfF/z9jVa2YifNxc4hPYCYy5w\nI7KtTwghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQghj+f9noJhkIwn0TAAAAABJRU5ErkJg\ngg==\n",
"text": [
""
]
}
],
"prompt_number": 24
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6\n",
"\n",
"a) Has the approximation of the mean been improved by using the GP model?"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"# Exercise 6 a) answer "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"b) One particular feature of GPs we have not use for now is their prediction variance. Can you use it to define some confidence intervals around the previous result?"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Exercise 6 b) answer"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computation of $P( f (U ) > 200)$\n",
"\n",
"In various cases it is interesting to look at the probability that $f$ is greater than a given\n",
"threshold. For example, assume that $f$ is the response of a physical model representing\n",
"the maximum constraint in a structure depending on some parameters of the system\n",
"such as Young\u2019s modulus of the material (say $Y$) and the force applied on the structure\n",
"(say $F$). If the later are uncertain, the probability of failure of the structure is given by\n",
"$P( f (Y, F ) > \\text{f_max} )$ where $f_\\text{max}$ is the maximum acceptable constraint.\n",
"\n",
"### Exercise 7\n",
"\n",
"a) As previously, use the 25 observations to compute a rough estimate of the probability that $f (U ) > 200$."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Exercise 7 a) answer"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"b) Compute the probability that the best predictor is greater than the threshold."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Exercise 7 b) answer"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"c) Compute some confidence intervals for the previous result"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Exercise 7 c) answer"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These two values can be compared with the actual value {$P( f (U ) > 200) = 1.23\\times 10^{\u22122}$ .\n",
"\n",
"We now assume that we have an extra budget of 10 evaluations of f and we want to\n",
"use these new evaluations to improve the accuracy of the previous result.\n",
"\n",
"### Exercise 8\n",
"\n",
"a) Given the previous GP model, where is it interesting to add the new observations if we want to improve the accuracy of the estimator and reduce its variance?"
]
},
{
"cell_type": "raw",
"metadata": {},
"source": [
"# Exercise 8 a) answer"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"b) Can you think about (and implement!) a procedure that updates sequentially the model with new points in order to improve the estimation of $P( f (U ) > 200)$?"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Exercise 8 b) answer"
],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}![\"Figure](\"http://ml.dcs.shef.ac.uk/gpss/gpws14/figa.png\")
\n",
"![\"Figure](\"http://ml.dcs.shef.ac.uk/gpss/gpws14/figb.png\")
\n",
"![\"Figure](\"http://ml.dcs.shef.ac.uk/gpss/gpws14/figc.png\")
\n",
"![\"Figure](\"http://ml.dcs.shef.ac.uk/gpss/gpws14/figd.png\")
\n",
"![\"Figure](\"http://ml.dcs.shef.ac.uk/gpss/gpws14/fige.png\")
\n",
"![\"Figure](\"http://ml.dcs.shef.ac.uk/gpss/gpws14/figf.png\")
\n",
"