# 1. 1D example

$\begin{align}
    f_{simulator}\colon\mathbb{R}^{N\times D} &\to\mathbb{R}^{N} \\
    X &\mapsto \mathbf{y}
\end{align}$


```python
import numpy as np
import torch
```


```python
# set logger and enforce reproducibility
from GPErks.log.logger import get_logger
from GPErks.utils.random import set_seed
log = get_logger()
seed = 8
set_seed(seed)  # reproducible sampling
```


<br/>

**1D function example**: Forrester et al. (2008)

$f(x) = (6x - 2)^2 \sin(12x - 4)$

<br/>



```python
# function to learn (normally a high-dimensional, expensive deterministic model)
from GPErks.utils.test_functions import forrester
f = lambda x: forrester(x)
D = 1
```


```python
# build dataset
from GPErks.gp.data.dataset import Dataset
dataset = Dataset.build_from_function(
    f,
    D,
    n_train_samples=10,
    n_test_samples=10,
    design="srs",
    seed=seed,
    l_bounds=[0],
    u_bounds=[1]  # can put None if, as in this case, parameters range in [0, 1]
)
```

<br/>

**Gaussian process emulator (GPE):**
<br/>

$f(\mathbf{x}) = h(\mathbf{x}) + g(\mathbf{x})$

<br/>

**deterministic part:**
<br/>

$h(\mathbf{x}) := \beta_0 + \beta_1 x_1 + \dots + \beta_{D} x_{D}$

<br/>

**stochastic part:**
<br/>

$\begin{align}
    &g(\mathbf{x})\sim\mathcal{GP}(\mathbf{0},\,k_{\text{SE}}(d(\mathbf{x},\,\mathbf{x}'))) \\
    &k_{\text{SE}}(d(\mathbf{x},\,\mathbf{x}')) := \sigma_f^2\, e^{-\frac{1}{2}\,d(\mathbf{x},\,\mathbf{x}')} \\
    &d(\mathbf{x},\,\mathbf{x}') := (\mathbf{x}-\mathbf{x}')^\mathsf{T}\,\Lambda\,(\mathbf{x}-\mathbf{x}')
\end{align}$

<br/>

**likelihood:**
<br/>

$y=f(\mathbf{x}) + \varepsilon,\quad \varepsilon\sim\mathcal{N}(0,\,\sigma_n^2)$

<br/>



```python
# choose likelihood
from gpytorch.likelihoods import GaussianLikelihood
likelihood = GaussianLikelihood()
```


```python
# choose mean function
from gpytorch.means import LinearMean
mean_function = LinearMean(input_size=dataset.input_size)
```


```python
# choose covariance function (kernel)
from gpytorch.kernels import RBFKernel, ScaleKernel
kernel = ScaleKernel(RBFKernel(ard_num_dims=dataset.input_size))
```


```python
# choose metrics
from torchmetrics import MeanSquaredError, R2Score
metrics = [MeanSquaredError(), R2Score()]
```


```python
# define experiment
from GPErks.gp.experiment import GPExperiment
experiment = GPExperiment(
    dataset,
    likelihood,
    mean_function,
    kernel,
    n_restarts=3,
    metrics=metrics,
    seed=seed,  # reproducible training
    learn_noise=True  # y = f(x) + e, e ~ N(0, sigma^2I)
)
```


```python
# choose training options: device + optimizer
device = "cuda" if torch.cuda.is_available() else "cpu"
optimizer = torch.optim.Adam(experiment.model.parameters(), lr=0.1)
```


```python
# train model
from GPErks.train.emulator import GPEmulator
emulator = GPEmulator(experiment, device)
emulator.train(optimizer)
```


```python
# inference on stored test set
x_test = dataset.X_test
y_test = dataset.y_test

y_mean, y_std = emulator.predict(x_test)

for metric in metrics:
    print( metric(
        torch.from_numpy(y_mean), torch.from_numpy(y_test)
        ).item()
    )
```


```python
# perk n.1: automatic inference
from GPErks.perks.inference import Inference
inference = Inference(emulator)
inference.summary()  # can be retrieved from inference.scores_dct
print( inference.scores_dct )
```


```python
# nice plotting
x_train = dataset.X_train
y_train = dataset.y_train

xx = np.linspace(dataset.l_bounds[0], dataset.u_bounds[0], 1000)
yy_mean, yy_std = emulator.predict(xx)
yy_true = f(xx)

import matplotlib.pyplot as plt
height = 9.36111
width = 5.91667
fig, axis = plt.subplots(1, 1, figsize=(4*width/3, height/2))

axis.plot(xx, yy_true, c="C0", ls="--", label="true function")

CI = 2
axis.plot(xx, yy_mean, c="C0", label="predicted mean")
axis.fill_between(
    xx, yy_mean - CI * yy_std, yy_mean + CI * yy_std, color="C0", alpha=0.15, label="~95% CI"
)
axis.scatter(x_train, y_train, fc="C0", ec="C0", label="training data")

axis.scatter(x_test, y_test, fc="none", ec="C0", label="testing data")

axis.legend(loc="best")
fig.tight_layout()
plt.show()
```


```python
# check testing points
inference.plot()
```


```python
# draw samples from the posterior distribution
y_mean, y_std = emulator.predict(x_test)
print(y_mean.shape)
print(y_std.shape)
```


```python
y_samples = emulator.sample(x_test, n_draws=5)
print(y_samples.shape)
```


```python
y_samples = emulator.sample(xx, n_draws=5)

fig, axis = plt.subplots(1, 1, figsize=(4*width/3, height/2))

for i, ys in enumerate(y_samples):
    axis.plot(xx, ys, lw=0.8, label=f"posterior sample #{i+1}", zorder=1)
    
axis.plot(xx, yy_mean, c="k", lw=2, ls="--", label="posterior mean", zorder=2)
axis.scatter(x_train, y_train, fc="k", ec="k", label="training data", zorder=2)

axis.legend(loc="best")
fig.tight_layout()
plt.show()
```

96431b650b5c8b24907934bd94164771d9650536
ipynb
Jupyter Notebook
notebooks/example_1.ipynb
stelong/GPErks
7e8e0e4561c10ad21fba2079619418e416a167b6
MIT
notebooks/example_1.ipynb
stelong/GPErks
7e8e0e4561c10ad21fba2079619418e416a167b6
MIT
2021-12-10T14:16:51.000Z
2022-03-25T16:26:50.000Z
notebooks/example_1.ipynb
stelong/GPErks
7e8e0e4561c10ad21fba2079619418e416a167b6
MIT
2022-01-28T11:12:33.000Z
2022-01-28T11:12:33.000Z
Qwen/Qwen-72B
1. YES
2. YES