0
# PyMC Gaussian Processes
1

2
PyMC provides a comprehensive framework for Gaussian process (GP) modeling, offering flexible covariance functions, mean functions, and efficient implementations for various GP applications including regression, classification, and time series modeling.
3

4
## Core GP Implementations
5

6
### Marginal Gaussian Processes
7

8
The standard GP implementation for regression and interpolation:
9

10
```python { .api }
11
from pymc.gp import Marginal
12
import pymc as pm
13

14
class Marginal:
15
    """
16
    Marginal Gaussian Process implementation.
17
    
18
    Parameters:
19
    - cov_func: Covariance function
20
    - mean_func: Mean function (default: Zero)
21
    
22
    Methods:
23
    - marginal_likelihood: Add GP likelihood to model
24
    - conditional: Compute conditional distribution
25
    - predict: Make predictions at new locations
26
    """
27
    
28
    def __init__(self, cov_func, mean_func=None):
29
        pass
30
    
31
    def marginal_likelihood(self, name, X, y, noise=None, is_observed=True, **kwargs):
32
        """
33
        Add marginal likelihood to the model.
34
        
35
        Parameters:
36
        - name (str): Name for the GP likelihood
37
        - X (array): Input locations for training data
38
        - y (array): Observed output values  
39
        - noise (float/array): Observation noise variance
40
        - is_observed (bool): Whether y contains observed data
41
        
42
        Returns:
43
        - gp_likelihood: GP likelihood random variable
44
        """
45
        pass
46
    
47
    def conditional(self, name, Xnew, given=None, **kwargs):
48
        """
49
        Compute conditional GP distribution at new points.
50
        
51
        Parameters:
52
        - name (str): Name for conditional distribution
53
        - Xnew (array): New input locations
54
        - given (dict): Conditioning data {'X': X_obs, 'y': y_obs, 'noise': noise}
55
        
56
        Returns:
57
        - conditional_gp: Conditional GP distribution
58
        """
59
        pass
60

61
# Basic GP regression example
62
with pm.Model() as gp_model:
63
    # GP hyperparameters
64
    ls = pm.Gamma('ls', alpha=2, beta=1)  # Length scale
65
    eta = pm.HalfNormal('eta', sigma=2)   # Signal variance
66
    sigma = pm.HalfNormal('sigma', sigma=1)  # Noise variance
67
    
68
    # Covariance function
69
    cov_func = eta**2 * pm.gp.cov.ExpQuad(1, ls)
70
    
71
    # GP likelihood
72
    gp = pm.gp.Marginal(cov_func=cov_func)
73
    y_obs = gp.marginal_likelihood('y_obs', X=X_train, y=y_train, noise=sigma)
74
    
75
    # Sample posterior
76
    trace = pm.sample()
77
```
78

79
### Marginal Approximations
80

81
Efficient approximations for large datasets:
82

83
```python { .api }
84
from pymc.gp import MarginalApprox
85

86
class MarginalApprox:
87
    """
88
    Marginal GP with inducing point approximation.
89
    
90
    Parameters:
91
    - cov_func: Covariance function
92
    - approx (str): Approximation method ('FITC', 'VFE', 'DTC')
93
    
94
    Methods:
95
    - marginal_likelihood: Add approximate GP likelihood
96
    """
97

98
# Sparse GP with inducing points
99
with pm.Model() as sparse_gp:
100
    # Select inducing points
101
    Xu = X_train[::10]  # Subset of training points
102
    
103
    # GP with sparse approximation
104
    gp_sparse = pm.gp.MarginalApprox(cov_func=cov_func, approx='FITC')
105
    y_sparse = gp_sparse.marginal_likelihood(
106
        'y_sparse', X=X_train, Xu=Xu, y=y_train, noise=sigma
107
    )
108
```
109

110
### Sparse GP Implementations
111

112
```python { .api }
113
from pymc.gp import MarginalSparse
114

115
# Variational sparse GP
116
with pm.Model() as variational_sparse:
117
    # Inducing point locations and values
118
    Xu = pm.Data('Xu', inducing_inputs)
119
    f_u = pm.MvNormal('f_u', mu=np.zeros(n_inducing), cov=K_uu)
120
    
121
    gp_sparse = pm.gp.MarginalSparse(cov_func=cov_func)
122
    y_sparse = gp_sparse.marginal_likelihood(
123
        'y_sparse', X=X_train, Xu=Xu, fu=f_u, y=y_train, noise=sigma
124
    )
125
```
126

127
### Kronecker Structure GPs
128

129
Efficient GPs for grid-structured data:
130

131
```python { .api }
132
from pymc.gp import MarginalKron
133

134
class MarginalKron:
135
    """
136
    Kronecker-structured marginal GP for multidimensional grids.
137
    
138
    Parameters:
139
    - cov_funcs (list): List of 1D covariance functions
140
    - mean_func: Mean function
141
    """
142

143
# Kronecker GP for 2D grid data
144
with pm.Model() as kronecker_gp:
145
    # Separate covariance for each dimension
146
    cov_x = eta1**2 * pm.gp.cov.Matern52(1, ls1)
147
    cov_y = eta2**2 * pm.gp.cov.Matern52(1, ls2)
148
    
149
    # Kronecker GP
150
    gp_kron = pm.gp.MarginalKron(cov_funcs=[cov_x, cov_y])
151
    z_obs = gp_kron.marginal_likelihood('z_obs', X=grid_coords, y=grid_values, 
152
                                        noise=sigma)
153
```
154

155
## Latent Gaussian Processes
156

157
For classification and non-Gaussian likelihoods:
158

159
```python { .api }
160
from pymc.gp import Latent
161

162
class Latent:
163
    """
164
    Latent Gaussian Process for non-Gaussian likelihoods.
165
    
166
    Parameters:
167
    - cov_func: Covariance function  
168
    - mean_func: Mean function
169
    
170
    Methods:
171
    - prior: Define GP prior distribution
172
    - conditional: Conditional distribution at new points
173
    """
174
    
175
    def prior(self, name, X, **kwargs):
176
        """
177
        Define GP prior at input locations.
178
        
179
        Parameters:
180
        - name (str): Name for GP prior variable
181
        - X (array): Input locations
182
        
183
        Returns:
184
        - f_prior: GP prior distribution
185
        """
186
        pass
187

188
# GP classification example
189
with pm.Model() as gp_classification:
190
    # GP hyperparameters
191
    ls = pm.Gamma('ls', alpha=2, beta=1)
192
    eta = pm.HalfNormal('eta', sigma=2)
193
    
194
    # Covariance function
195
    cov_func = eta**2 * pm.gp.cov.Matern32(1, ls)
196
    
197
    # Latent GP
198
    gp = pm.gp.Latent(cov_func=cov_func)
199
    f = gp.prior('f', X=X_train)
200
    
201
    # Probit link for classification
202
    p = pm.invprobit(f)
203
    y_obs = pm.Bernoulli('y_obs', p=p, observed=y_binary)
204
```
205

206
### Latent Kronecker GPs
207

208
```python { .api }
209
from pymc.gp import LatentKron
210

211
# Latent Kronecker GP for spatial classification
212
with pm.Model() as spatial_classification:
213
    gp_kron = pm.gp.LatentKron(cov_funcs=[cov_x, cov_y])
214
    f_latent = gp_kron.prior('f_latent', X=spatial_coords)
215
    
216
    # Spatial logistic regression
217
    p_spatial = pm.invlogit(f_latent)
218
    y_spatial = pm.Bernoulli('y_spatial', p=p_spatial, observed=spatial_outcomes)
219
```
220

221
## Student-t Processes
222

223
Robust alternative to Gaussian processes:
224

225
```python { .api }
226
from pymc.gp import TP
227

228
class TP:
229
    """
230
    Student-t Process implementation.
231
    
232
    Parameters:
233
    - cov_func: Covariance function
234
    - nu: Degrees of freedom parameter
235
    """
236

237
# Robust regression with t-process
238
with pm.Model() as tp_model:
239
    # t-process parameters
240
    nu = pm.Gamma('nu', alpha=2, beta=0.1)  # Degrees of freedom
241
    
242
    # t-process
243
    tp = pm.gp.TP(cov_func=cov_func, nu=nu)
244
    y_tp = tp.marginal_likelihood('y_tp', X=X_train, y=y_train, noise=sigma)
245
```
246

247
## Hilbert Space Gaussian Processes
248

249
Efficient approximation for large-scale GPs:
250

251
```python { .api }
252
from pymc.gp import HSGP, HSGPPeriodic
253

254
class HSGP:
255
    """
256
    Hilbert Space Gaussian Process approximation.
257
    
258
    Parameters:
259
    - m (list): Number of basis functions per dimension
260
    - L (list): Boundary values per dimension  
261
    - cov_func: Covariance function to approximate
262
    """
263
    
264
    def __init__(self, m, L, cov_func):
265
        pass
266
    
267
    def prior_linearized(self, Xs):
268
        """
269
        Compute linearized prior representation.
270
        
271
        Parameters:
272
        - Xs (array): Input locations
273
        
274
        Returns:
275
        - phi: Basis function matrix
276
        - sqrt_psd: Square root of power spectral density
277
        """
278
        pass
279

280
# HSGP for large datasets
281
with pm.Model() as hsgp_model:
282
    # HSGP parameters
283
    m = [50]  # Number of basis functions
284
    L = [2.0]  # Boundary condition
285
    
286
    # HSGP approximation
287
    hsgp = pm.gp.HSGP(m=m, L=L, cov_func=cov_func)
288
    phi, sqrt_psd = hsgp.prior_linearized(X_train)
289
    
290
    # Linear model with HSGP features
291
    beta = pm.Normal('beta', 0, 1, shape=m[0])
292
    f_hsgp = phi @ (sqrt_psd * beta)
293
    
294
    # Likelihood
295
    y_hsgp = pm.Normal('y_hsgp', mu=f_hsgp, sigma=sigma, observed=y_train)
296
```
297

298
### Periodic HSGP
299

300
```python { .api }
301
class HSGPPeriodic:
302
    """
303
    HSGP for periodic functions.
304
    
305
    Parameters:
306
    - m (int): Number of basis functions
307
    - period (float): Period of the function
308
    - cov_func: Periodic covariance function
309
    """
310

311
# Periodic GP approximation
312
with pm.Model() as periodic_hsgp:
313
    # Periodic parameters
314
    period = 12.0  # Annual cycle
315
    m_periodic = 20
316
    
317
    hsgp_periodic = pm.gp.HSGPPeriodic(
318
        m=m_periodic, 
319
        period=period,
320
        cov_func=periodic_cov
321
    )
322
```
323

324
## Covariance Functions
325

326
### Stationary Covariance Functions
327

328
```python { .api }
329
import pymc.gp.cov as cov
330

331
# Exponentiated Quadratic (RBF/Squared Exponential)
332
class ExpQuad:
333
    """
334
    Exponentiated quadratic covariance function.
335
    
336
    Parameters:
337
    - input_dim (int): Number of input dimensions
338
    - ls (float/array): Length scale(s)
339
    - active_dims (list): Active input dimensions
340
    """
341
    
342
    def __init__(self, input_dim, ls, active_dims=None):
343
        pass
344

345
eq_cov = cov.ExpQuad(input_dim=2, ls=[1.0, 2.0])
346

347
# Matérn Covariance Functions
348
matern32 = cov.Matern32(input_dim=1, ls=1.5)  # Matérn 3/2
349
matern52 = cov.Matern52(input_dim=1, ls=1.5)  # Matérn 5/2
350

351
# Rational Quadratic
352
rat_quad = cov.RatQuad(input_dim=1, ls=1.0, alpha=2.0)
353

354
# Exponential (Matérn 1/2)
355
exponential = cov.Exponential(input_dim=1, ls=1.0)
356
```
357

358
### Non-Stationary Covariance Functions
359

360
```python { .api }
361
# Linear covariance
362
linear_cov = cov.Linear(input_dim=1, c=0.5)
363

364
# Polynomial covariance
365
poly_cov = cov.Polynomial(input_dim=1, c=1.0, d=3, offset=0.0)
366

367
# Gibbs (non-stationary)
368
def length_scale_func(x):
369
    return 0.1 + 0.9 * pm.math.sigmoid(x)
370

371
gibbs_cov = cov.Gibbs(input_dim=1, lengthscale_func=length_scale_func)
372
```
373

374
### Periodic Covariance Functions
375

376
```python { .api }
377
# Periodic covariance
378
periodic_cov = cov.Periodic(input_dim=1, period=12, ls=1.0)
379

380
# Cosine covariance
381
cosine_cov = cov.Cosine(input_dim=1, ls=1.0)
382
```
383

384
### Utility Covariance Functions
385

386
```python { .api }
387
# Constant covariance (bias term)
388
constant_cov = cov.Constant(c=2.0)
389

390
# White noise
391
white_noise = cov.WhiteNoise(sigma=0.1)
392

393
# Scaled covariance
394
scaled_cov = cov.ScaledCov(input_dim=1, scaling_func=scaling_function)
395

396
# Warped input covariance
397
warp_func = lambda x: pm.math.tanh(x)
398
warped_cov = cov.WarpedInput(input_dim=1, cov_func=base_cov, warp_func=warp_func)
399
```
400

401
## Covariance Function Operations
402

403
### Combining Covariance Functions
404

405
```python { .api }
406
# Addition (sum of independent processes)
407
combined_cov = cov1 + cov2 + cov3
408

409
# Multiplication (product of covariances)
410
product_cov = cov1 * cov2
411

412
# Example: Trend + Periodic + Noise
413
trend_cov = cov.Linear(input_dim=1, c=1.0)
414
periodic_cov = cov.Periodic(input_dim=1, period=12, ls=1.0) 
415
noise_cov = cov.WhiteNoise(sigma=0.1)
416

417
full_cov = trend_cov + periodic_cov + noise_cov
418
```
419

420
### Active Dimensions
421

422
```python { .api }
423
# Covariance functions for specific dimensions
424
cov_dim1 = cov.Matern52(input_dim=3, ls=1.0, active_dims=[0])      # Only 1st dim
425
cov_dim23 = cov.ExpQuad(input_dim=3, ls=[0.5, 2.0], active_dims=[1, 2])  # 2nd,3rd dims
426

427
# Combined multidimensional covariance
428
multi_cov = cov_dim1 + cov_dim23
429
```
430

431
## Mean Functions
432

433
### Standard Mean Functions
434

435
```python { .api }
436
import pymc.gp.mean as mean
437

438
# Zero mean (default)
439
class Zero:
440
    """Zero mean function."""
441
    pass
442

443
zero_mean = mean.Zero()
444

445
# Constant mean  
446
class Constant:
447
    """
448
    Constant mean function.
449
    
450
    Parameters:
451
    - c (float): Constant value
452
    """
453

454
constant_mean = mean.Constant(c=2.5)
455

456
# Linear mean
457
class Linear:
458
    """
459
    Linear mean function.
460
    
461
    Parameters:
462
    - coeffs (array): Linear coefficients
463
    - intercept (float): Intercept term
464
    """
465

466
linear_mean = mean.Linear(coeffs=[0.5, -1.2], intercept=1.0)
467
```
468

469
## Advanced GP Patterns
470

471
### Multi-Output GPs
472

473
```python { .api }
474
# Coregionalization model for multiple outputs
475
with pm.Model() as multioutput_gp:
476
    # Shared latent functions
477
    n_latent = 2
478
    n_outputs = 3
479
    
480
    # Mixing matrix
481
    W = pm.Normal('W', 0, 1, shape=(n_outputs, n_latent))
482
    
483
    # Latent GPs
484
    gps = []
485
    for i in range(n_latent):
486
        gp_i = pm.gp.Latent(cov_func=shared_cov)
487
        f_i = gp_i.prior(f'f_{i}', X=X_train)
488
        gps.append(f_i)
489
    
490
    # Linear combinations for each output
491
    F_latent = pm.math.stack(gps, axis=1)  # Shape: (n_obs, n_latent)
492
    F_outputs = pm.math.dot(F_latent, W.T)  # Shape: (n_obs, n_outputs)
493
    
494
    # Output likelihoods
495
    for j in range(n_outputs):
496
        y_j = pm.Normal(f'y_{j}', mu=F_outputs[:, j], sigma=sigma_j[j], 
497
                        observed=Y_train[:, j])
498
```
499

500
### Hierarchical GPs
501

502
```python { .api }
503
# Hierarchical GP with group-specific parameters
504
with pm.Model() as hierarchical_gp:
505
    # Global hyperparameters
506
    mu_ls = pm.Gamma('mu_ls', alpha=2, beta=1)
507
    sigma_ls = pm.HalfNormal('sigma_ls', sigma=0.5)
508
    
509
    # Group-specific length scales
510
    ls_groups = pm.Lognormal('ls_groups', mu=pm.math.log(mu_ls), 
511
                             sigma=sigma_ls, shape=n_groups)
512
    
513
    # Group-specific GPs
514
    for g in range(n_groups):
515
        cov_g = eta**2 * cov.Matern52(1, ls_groups[g])
516
        gp_g = pm.gp.Marginal(cov_func=cov_g)
517
        
518
        # Group data
519
        X_g = X_train[group_idx == g]
520
        y_g = y_train[group_idx == g] 
521
        
522
        y_obs_g = gp_g.marginal_likelihood(f'y_obs_{g}', X=X_g, y=y_g, noise=sigma)
523
```
524

525
### GP Regression with Heteroscedastic Noise
526

527
```python { .api }
528
# GP with input-dependent noise
529
with pm.Model() as heteroscedastic_gp:
530
    # Main regression GP
531
    gp_mean = pm.gp.Marginal(cov_func=cov_func)
532
    f_mean = gp_mean.prior('f_mean', X=X_train)
533
    
534
    # Noise variance GP (log scale)
535
    gp_noise = pm.gp.Marginal(cov_func=noise_cov_func)
536
    f_noise = gp_noise.prior('f_noise', X=X_train)
537
    log_sigma = pm.Deterministic('log_sigma', f_noise)
538
    
539
    # Heteroscedastic likelihood
540
    y_hetero = pm.Normal('y_hetero', mu=f_mean, sigma=pm.math.exp(log_sigma), 
541
                         observed=y_train)
542
```
543

544
## GP Utilities and Workflows
545

546
### Model Comparison
547

548
```python { .api }
549
# Compare different covariance functions
550
covariance_functions = {
551
    'matern32': cov.Matern32(1, ls),
552
    'matern52': cov.Matern52(1, ls), 
553
    'exp_quad': cov.ExpQuad(1, ls),
554
    'rational_quad': cov.RatQuad(1, ls, alpha)
555
}
556

557
models = {}
558
traces = {}
559

560
for name, cov_func in covariance_functions.items():
561
    with pm.Model() as model:
562
        gp = pm.gp.Marginal(cov_func=cov_func)
563
        y_obs = gp.marginal_likelihood('y_obs', X=X_train, y=y_train, noise=sigma)
564
        trace = pm.sample()
565
        
566
        models[name] = model
567
        traces[name] = trace
568

569
# Compare models
570
comparison = pm.compare(traces)
571
```
572

573
### Posterior Prediction Workflow
574

575
```python { .api }
576
# Complete GP prediction workflow
577
with pm.Model() as gp_model:
578
    # Model definition...
579
    trace = pm.sample()
580
    
581
# Predict at new locations
582
with gp_model:
583
    # Conditional distribution for predictions
584
    f_pred = gp.conditional('f_pred', Xnew=X_test)
585
    
586
    # Sample predictions
587
    pred_samples = pm.sample_posterior_predictive(
588
        trace, 
589
        var_names=['f_pred'],
590
        predictions=True
591
    )
592

593
# Extract prediction statistics
594
pred_mean = pred_samples.predictions['f_pred'].mean(dim=['chain', 'draw'])
595
pred_std = pred_samples.predictions['f_pred'].std(dim=['chain', 'draw'])
596
```
597

598
PyMC's Gaussian process framework provides a flexible and efficient foundation for non-parametric Bayesian modeling, supporting everything from simple regression to complex multi-output and spatiotemporal models.