From 9c7c2b826720ce03003252f9842d0e6261cf7f56 Mon Sep 17 00:00:00 2001 From: preethamak <167765989+preethamak@users.noreply.github.com> Date: Sat, 11 Oct 2025 23:56:55 +0530 Subject: [PATCH] Create 186. Gaussian_process_for_regression --- 186. Gaussian_process_for_regression | 174 +++++++++++++++++++++++++++ 1 file changed, 174 insertions(+) create mode 100644 186. Gaussian_process_for_regression diff --git a/186. Gaussian_process_for_regression b/186. Gaussian_process_for_regression new file mode 100644 index 00000000..2012e9dc --- /dev/null +++ b/186. Gaussian_process_for_regression @@ -0,0 +1,174 @@ +import numpy as np +import math + +# --- Kernels ------------------------------------------------------------- + +def linear_kernel(x, x_prime, sigma_b=0.0, sigma_v=1.0): + return sigma_b**2 + sigma_v**2 * np.dot(x, x_prime.T) + +def rbf_kernel(x, x_prime, sigma=1.0, length_scale=1.0): + sqdist = np.sum((x[:, None] - x_prime[None, :]) ** 2, axis=2) + return sigma**2 * np.exp(-0.5 * sqdist / length_scale**2) + +def matern_kernel(x, x_prime, length_scale=1.0, nu=1.5): + d = np.sqrt(np.sum((x[:, None] - x_prime[None, :]) ** 2, axis=2)) + if nu == 0.5: + k = np.exp(-d / length_scale) + elif nu == 1.5: + sqrt3 = np.sqrt(3) + k = (1 + sqrt3 * d / length_scale) * np.exp(-sqrt3 * d / length_scale) + elif nu == 2.5: + sqrt5 = np.sqrt(5) + k = (1 + sqrt5 * d / length_scale + 5 * d**2 / (3 * length_scale**2)) * np.exp(-sqrt5 * d / length_scale) + else: + raise NotImplementedError + return k + +def periodic_kernel(x, x_prime, sigma=1.0, length_scale=1.0, period=1.0): + diff = np.pi * np.abs(x[:, None] - x_prime[None, :]) / period + return sigma**2 * np.exp(-2 * (np.sin(diff) ** 2) / length_scale**2) + +def rational_quadratic_kernel(x, x_prime, sigma=1.0, length_scale=1.0, alpha=1.0): + sqdist = np.sum((x[:, None] - x_prime[None, :]) ** 2, axis=2) + return sigma**2 * (1 + sqdist / (2 * alpha * length_scale**2)) ** (-alpha) + +# --- Base Class ------------------------------------------------------------- + +class _GaussianProcessBase: + def __init__(self, kernel="rbf", noise=1e-5, kernel_params=None): + self.kernel_name = kernel + self.noise = noise + self.kernel_params = kernel_params or {} + + def _select_kernel(self, x1, x2): + if self.kernel_name == "linear": + return linear_kernel(x1, x2, **self.kernel_params) + elif self.kernel_name == "rbf": + return rbf_kernel(x1, x2, **self.kernel_params) + elif self.kernel_name == "matern": + return matern_kernel(x1, x2, **self.kernel_params) + elif self.kernel_name == "periodic": + return periodic_kernel(x1, x2, **self.kernel_params) + elif self.kernel_name == "rational_quadratic": + return rational_quadratic_kernel(x1, x2, **self.kernel_params) + else: + raise ValueError("Unknown kernel") + + def _compute_covariance(self, X1, X2): + return self._select_kernel(X1, X2) + +# --- Gaussian Process Regression -------------------------------------------- + +class GaussianProcessRegression(_GaussianProcessBase): + def fit(self, X, y): + self.X_train = np.array(X) + self.y_train = np.array(y) + K = self._compute_covariance(self.X_train, self.X_train) + self.K = K + self.noise * np.eye(len(X)) + self.K_inv = np.linalg.inv(self.K) + + def predict(self, X_test, return_std=False): + X_test = np.array(X_test) + K_s = self._compute_covariance(X_test, self.X_train) + K_ss = self._compute_covariance(X_test, X_test) + 1e-8 * np.eye(len(X_test)) + mu = K_s @ self.K_inv @ self.y_train + if return_std: + cov = K_ss - K_s @ self.K_inv @ K_s.T + std = np.sqrt(np.diag(cov)) + return mu, std + return mu + + def log_marginal_likelihood(self): + n = len(self.y_train) + sign, logdet = np.linalg.slogdet(self.K) + return -0.5 * self.y_train.T @ self.K_inv @ self.y_train - 0.5 * logdet - 0.5 * n * np.log(2 * np.pi) + + def optimize_hyperparameters(self): + pass # Not needed for this task +import numpy as np +import math + +# --- Kernels ------------------------------------------------------------- + +def linear_kernel(x, x_prime, sigma_b=0.0, sigma_v=1.0): + return sigma_b**2 + sigma_v**2 * np.dot(x, x_prime.T) + +def rbf_kernel(x, x_prime, sigma=1.0, length_scale=1.0): + sqdist = np.sum((x[:, None] - x_prime[None, :]) ** 2, axis=2) + return sigma**2 * np.exp(-0.5 * sqdist / length_scale**2) + +def matern_kernel(x, x_prime, length_scale=1.0, nu=1.5): + d = np.sqrt(np.sum((x[:, None] - x_prime[None, :]) ** 2, axis=2)) + if nu == 0.5: + k = np.exp(-d / length_scale) + elif nu == 1.5: + sqrt3 = np.sqrt(3) + k = (1 + sqrt3 * d / length_scale) * np.exp(-sqrt3 * d / length_scale) + elif nu == 2.5: + sqrt5 = np.sqrt(5) + k = (1 + sqrt5 * d / length_scale + 5 * d**2 / (3 * length_scale**2)) * np.exp(-sqrt5 * d / length_scale) + else: + raise NotImplementedError + return k + +def periodic_kernel(x, x_prime, sigma=1.0, length_scale=1.0, period=1.0): + diff = np.pi * np.abs(x[:, None] - x_prime[None, :]) / period + return sigma**2 * np.exp(-2 * (np.sin(diff) ** 2) / length_scale**2) + +def rational_quadratic_kernel(x, x_prime, sigma=1.0, length_scale=1.0, alpha=1.0): + sqdist = np.sum((x[:, None] - x_prime[None, :]) ** 2, axis=2) + return sigma**2 * (1 + sqdist / (2 * alpha * length_scale**2)) ** (-alpha) + +# --- Base Class ------------------------------------------------------------- + +class _GaussianProcessBase: + def __init__(self, kernel="rbf", noise=1e-5, kernel_params=None): + self.kernel_name = kernel + self.noise = noise + self.kernel_params = kernel_params or {} + + def _select_kernel(self, x1, x2): + if self.kernel_name == "linear": + return linear_kernel(x1, x2, **self.kernel_params) + elif self.kernel_name == "rbf": + return rbf_kernel(x1, x2, **self.kernel_params) + elif self.kernel_name == "matern": + return matern_kernel(x1, x2, **self.kernel_params) + elif self.kernel_name == "periodic": + return periodic_kernel(x1, x2, **self.kernel_params) + elif self.kernel_name == "rational_quadratic": + return rational_quadratic_kernel(x1, x2, **self.kernel_params) + else: + raise ValueError("Unknown kernel") + + def _compute_covariance(self, X1, X2): + return self._select_kernel(X1, X2) + +# --- Gaussian Process Regression -------------------------------------------- + +class GaussianProcessRegression(_GaussianProcessBase): + def fit(self, X, y): + self.X_train = np.array(X) + self.y_train = np.array(y) + K = self._compute_covariance(self.X_train, self.X_train) + self.K = K + self.noise * np.eye(len(X)) + self.K_inv = np.linalg.inv(self.K) + + def predict(self, X_test, return_std=False): + X_test = np.array(X_test) + K_s = self._compute_covariance(X_test, self.X_train) + K_ss = self._compute_covariance(X_test, X_test) + 1e-8 * np.eye(len(X_test)) + mu = K_s @ self.K_inv @ self.y_train + if return_std: + cov = K_ss - K_s @ self.K_inv @ K_s.T + std = np.sqrt(np.diag(cov)) + return mu, std + return mu + + def log_marginal_likelihood(self): + n = len(self.y_train) + sign, logdet = np.linalg.slogdet(self.K) + return -0.5 * self.y_train.T @ self.K_inv @ self.y_train - 0.5 * logdet - 0.5 * n * np.log(2 * np.pi) + + def optimize_hyperparameters(self): + pass # Not needed for this task