# Source: content/notes/ml/probabilistic-and-instance-models.md
# Independent CPU example; use the curriculum environment.
# See /notes/ml/#example-environment or /notes/deep-learning/#example-environment.

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF
from sklearn.metrics import mean_squared_error

rng = np.random.default_rng(20)
X_train = np.linspace(-3, 3, 45)[:, None]
noise_std = 0.12
y_train = np.sin(X_train[:, 0]) + rng.normal(0, noise_std, len(X_train))
X_test = np.linspace(-2.9, 2.9, 70)[:, None]
y_test = np.sin(X_test[:, 0]) + rng.normal(0, noise_std, len(X_test))
kernel = ConstantKernel(1.0, constant_value_bounds="fixed") * RBF(1.0, length_scale_bounds="fixed")
model = GaussianProcessRegressor(kernel=kernel, alpha=noise_std ** 2,
                                 optimizer=None, normalize_y=False, random_state=20)
model.fit(X_train, y_train)
mean, latent_std = model.predict(X_test, return_std=True)
observation_std = np.sqrt(latent_std ** 2 + noise_std ** 2)
coverage = np.mean(np.abs(y_test - mean) <= 1.96 * observation_std)
assert mean.shape == y_test.shape and np.isfinite(mean).all()
assert np.all(observation_std >= latent_std)
assert np.sqrt(mean_squared_error(y_test, mean)) < 0.4
print("RMSE:", np.sqrt(mean_squared_error(y_test, mean)))
print("Observed 95% interval coverage:", coverage)
print("Log marginal likelihood:", model.log_marginal_likelihood_value_)