{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sequential Design\n",
"This notebook shows a sequential design example in `hetGPy`. It combines several aspects of previous examples, such as acquisition functions like **Expected Improvement (EI)** or **Integreated Mean-Squared Prediction Error (IMPSE)**.\n",
"\n",
"Note that these examples represent advanced use of the package."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sequential Design for Expected Improvement\n",
"\n",
"Consider the following example from Figure 10 of [Jones, Schonlau, and Welch (1998)](https://link.springer.com/article/10.1023/A:1008306431147) and from [(Schonlau 1997)](https://uwspace.uwaterloo.ca/handle/10012/190). The code and narrative text for this example were adpated from chapter 7.2 of [Gramacy (2020)](https://bookdown.org/rbg/surrogates/chap7.html#chap7ei).\n",
"\n",
"Note the large gap along the x-axis before the final input. This leads to a large prediction interval between those inputs. We can calculate the EI along a dense grid of inputs and notice the tradeoff between \"exploration\" (choosing a design in the high noise region) and \"exploitation\" (choosing a design near the current observed globabl minimum)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%config InlineBackend.figure_formats = ['svg']\n",
"%matplotlib inline\n",
"import numpy as np\n",
"from hetgpy import homGP\n",
"from hetgpy.optim import crit_EI\n",
"from scipy.stats import norm\n",
"import matplotlib.pyplot as plt\n",
"X = np.array([1, 2, 3, 4, 12]).reshape(-1,1)\n",
"Y = np.array([0, -1.75, -2, -0.5, 5])\n",
"\n",
"model = homGP()\n",
"model.mle(X=X,\n",
" Z=Y,\n",
" covtype='Gaussian',\n",
" init={'theta':np.array([10.0]),'g':1e-8}\n",
")\n",
"Xgrid = np.linspace(1,12,1000).reshape(-1,1)\n",
"preds = model.predict(x=Xgrid)\n",
"preds['upper'] = norm.ppf(0.95, loc = preds['mean'], scale = np.sqrt(preds['sd2'])).squeeze()\n",
"preds['lower'] = norm.ppf(0.05, loc = preds['mean'], scale = np.sqrt(preds['sd2'])).squeeze()\n",
"\n",
"fig, ax = plt.subplots(ncols=2,figsize=(9,6),sharex=True)\n",
"ax[0].scatter(X,Y,color='blue',label='Initial')\n",
"ax[0].plot(Xgrid,preds['mean'],color='blue')\n",
"ax[0].plot(Xgrid,preds['upper'],linestyle='dashed',color='blue')\n",
"ax[0].plot(Xgrid,preds['lower'],linestyle='dashed',color='blue')\n",
"\n",
"EI_grid = crit_EI(x=Xgrid,model=model)\n",
"ax[1].plot(Xgrid,EI_grid,color='teal',label='EI')\n",
"ax[1].set_ylim([0,0.135])\n",
"\n",
"ax[0].legend(edgecolor='black');ax[1].legend(edgecolor='black');\n",
"ax[0].set_xlabel('X');ax[1].set_xlabel('X');\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can then acquire a new design at the maximum EI location and update our model:\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"EI_grid = crit_EI(x=Xgrid,model=model)\n",
"design = Xgrid[EI_grid.argmax()]\n",
"Xnew = design.reshape(-1,1)\n",
"Znew = model.predict(Xnew)['mean']\n",
"model.update(Xnew=Xnew,Znew=Znew)\n",
"\n",
"ax[0].scatter(Xnew,Znew,label='Acquisition',s=100,color='green',zorder=2)\n",
"ax[0].text(Xnew,Znew,s='0',color='black',ha='center',va='center',zorder=3)\n",
"ax[0].legend(edgecolor='black');\n",
"EI_grid = crit_EI(x=Xgrid,model=model)\n",
"ax[1].clear()\n",
"ax[1].plot(Xgrid,EI_grid,color='teal',label='EI')\n",
"ax[1].legend()\n",
"ax[0].set_xlabel('X');ax[1].set_xlabel('X');\n",
"fig.tight_layout()\n",
"fig\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And then we can acquire another update to capitalize on the \"exploration\" component (and reduce our uncertainty on the righthand side of the domain):"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"EI_grid = crit_EI(x=Xgrid,model=model)\n",
"design = Xgrid[EI_grid.argmax()]\n",
"Xnew = design.reshape(-1,1)\n",
"Znew = model.predict(Xnew)['mean']\n",
"model.update(Xnew=Xnew,Znew=Znew)\n",
"\n",
"preds = model.predict(x=Xgrid)\n",
"preds['upper'] = norm.ppf(0.95, loc = preds['mean'], scale = np.sqrt(preds['sd2'])).squeeze()\n",
"preds['lower'] = norm.ppf(0.05, loc = preds['mean'], scale = np.sqrt(preds['sd2'])).squeeze()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ax[0].scatter(Xnew,Znew,s=100,color='green',zorder=2)\n",
"ax[0].text(Xnew,Znew,s='1',color='black',ha='center',va='center',zorder=3)\n",
"ax[0].plot(Xgrid,preds['lower'],color='green',linestyle='dashed')\n",
"ax[0].plot(Xgrid,preds['upper'],color='green',linestyle='dashed')\n",
"ax[0].legend(edgecolor='black');\n",
"EI_grid = crit_EI(x=Xgrid,model=model)\n",
"ax[1].clear()\n",
"ax[1].plot(Xgrid,EI_grid,color='teal',label='EI')\n",
"ax[1].legend(edgecolor='black')\n",
"ax[0].set_xlabel('X');ax[1].set_xlabel('X');\n",
"fig.tight_layout()\n",
"fig"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sequential Design for Integrated Mean-Squared Prediction Error\n",
"\n",
"\n",
"Instead of trying to find the global optimum of a function, a different goal is to efficiently learn the response surface, and thus pose the sequential design problem as one where we wish to minimize the (global) predictive variance.\n",
"\n",
"Consider the function $f(x)=2(\\exp(-30(x - 0.25)^2) + \\sin(x^2)) - 2$ (implemented as `f1d2` in `hetGPy`) from [Boukouvalas and Cornford (2009)](https://scholar.google.com/scholar?q=Boukouvalas+A,+Cornford+D+(2009).+%E2%80%9CLearning+Heteroscedastic+Gaussian+Processes+for+Complex+Datasets.%E2%80%9D&hl=en&as_sdt=0&as_vis=1&oi=scholart). The code below adapts the example from [Binois et. al 2021](https://www.jstatsoft.org/article/view/v098i13), Figure 7.\n",
"\n",
"We use a heteroskedastic noise function $f_n(x)=1/3\\exp(\\sin(2\\pi x))$ so the overall data-generating process is: \n",
"$$\n",
"\\begin{align*}\n",
"Y \\sim N(f(x),r(x))\n",
"\\end{align*}\n",
"$$\n",
"\n",
"\n",
"Since choosing a design location that is either a replicate (and thus from a discrete set of inputs) or a new one (which is continuous), we are mixing discrete and continuous optimization. This poses an issue for choosing replicates (since up to floating point error, we are unlikely to choose _exactly_ the same design location) We also use a replication-biased lookahead heuristic (governed by the $h$ parameter in `IMSPE_optim`) from [Binois et. al 2019](https://arxiv.org/abs/1710.03206) that tends to prefer replication. Details can be found in that paper or in [Gramacy 2020 Ch. 10.3.2](https://bookdown.org/rbg/surrogates/chap10.html#chap10lookahead).\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from hetgpy import hetGP\n",
"from hetgpy.IMSE import IMSPE_optim\n",
"\n",
"rand = np.random.default_rng(13)\n",
"def f1d2(x): \n",
" return 2*(np.exp(-30*(x - 0.25)**2) + np.sin(x**2)) - 2\n",
"def fn(x): \n",
" return 1/3 * (np.exp(np.sin(2 * np.pi * x)))\n",
"def fr(x):\n",
" x = x.squeeze().reshape(-1)\n",
" return f1d2(x) + rand.normal(size = max(x.shape[0],1), scale = fn(x), loc = 0)\n",
"\n",
"X = np.linspace(0, 1, 15).reshape(-1,1)\n",
"Y = fr(X).squeeze()\n",
"Xorig = X.copy()\n",
"mod = hetGP()\n",
"mod.mle(X = X, Z = Y, lower = np.array([0.0001]), upper = np.array([1.0]),known={},init={})\n",
"opt = IMSPE_optim(mod, h = 5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code cell wraps the routine in a loop to acquire 500 samples sequentially. Every 50 iterations, we re-train on the entire set of data acquired thus far (to prevent the log-likelihood from getting stuck in local optima).\n",
"\n",
"Note that this cell takes a few minutes to run (hence the progress bar)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 500/500 [05:16<00:00, 1.58it/s]\n"
]
}
],
"source": [
"import warnings\n",
"from tqdm import tqdm\n",
"warnings.simplefilter(\"ignore\") # ignore occasional RuntimeWarnings in ll optimization\n",
"for i in tqdm(range(500)):\n",
" opt = IMSPE_optim(mod, h = 5)\n",
" X = np.concatenate([X, opt['par']],axis=0)\n",
" Ynew = fr(opt['par'])\n",
" Y = np.concatenate([Y,Ynew])\n",
" mod.update(Xnew = opt['par'], Znew = Ynew, ginit = mod.g*1.01,maxit=100,known={})\n",
" if (i + 1) % 50 == 0:\n",
" mod2 = hetGP()\n",
" X, Y = X.copy(), Y.copy()\n",
" mod2.mle(X = X,\n",
" Z = Y, \n",
" lower = np.array([0.0001]),\n",
" upper = np.array([1.0])\n",
" )\n",
" if mod2.ll > mod.ll:\n",
" mod = mod2.copy()\n",
"grid = np.linspace(0,1,1000).reshape(-1,1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The resulting figure shows that we acquire many samples (and replicates) in the high noise region. We also see that we get a faithful estimate of the data-generating process (shown in red). Gray bars shown the original design locations (with no replicates)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"preds = mod.predict(grid)\n",
"preds['upper'] = norm.ppf(0.95, loc = preds['mean'], scale = np.sqrt(preds['sd2']+preds['nugs']).squeeze())\n",
"preds['lower'] = norm.ppf(0.05, loc = preds['mean'], scale = np.sqrt(preds['sd2']+preds['nugs']).squeeze())\n",
"fig, ax = plt.subplots()\n",
"ltrue = norm.ppf(0.05,loc=f1d2(grid),scale=fn(grid))\n",
"utrue = norm.ppf(0.95,loc=f1d2(grid),scale=fn(grid))\n",
"ax.vlines(Xorig,color='gray',ymin=-100,ymax=-4.5,label='Original',alpha=0.5)\n",
"# Estimate\n",
"ax.plot(grid,preds['mean'],label='Estimate',color='blue')\n",
"ax.plot(grid,preds['upper'],color='blue',linestyle='dashed')\n",
"ax.plot(grid,preds['lower'],color='blue',linestyle='dashed')\n",
"# Ground truth\n",
"ax.plot(grid,f1d2(grid),color='red',label='Truth')\n",
"ax.plot(grid,ltrue,color='red',linestyle='dashed')\n",
"ax.plot(grid,utrue,color='red',linestyle='dashed')\n",
"\n",
"ax.scatter(X,Y,color='green',label='Acquistions')\n",
"ax.set_ylim(-5, 5)\n",
"ax.legend(edgecolor='black',facecolor='white',framealpha=1.0)\n",
"ax.set_xlabel('X');ax.set_ylabel('f1d2(X)')\n",
"fig.tight_layout()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}