MOGPTK (Multi-Output Gaussian Process Toolkit)

MOGPTK (Multi-Output Gaussian Process Toolkit) is a PyTorch-based library for interpretable probabilistic modelling of multichannel time series. It is designed to model signals by utilizing correlation information among different channels.

Core Capabilities

• Multi-Output GP Regression: Uses kernels like the Multi Output Spectral Mixture (MOSM) kernel to correlate data points across different channels.
• GPU Acceleration: Leverages PyTorch for significantly faster model training.
• Scalability: Includes sparse-variational GP regression functionality to handle larger datasets and reduce computation time.
• Signal Reconstruction: Particularly effective for reconstructing signals in cases of sporadic data loss due to its ability to learn cross-channel correlations.

Supported Inference Models

• Exact Gaussian process (maximum likelihood)
• Snelson
• OpperArchambeau
• Titsias
• Hensman

The following tutorials are available to guide you through using the toolkit:

• 00 - Quick Start: Basic usage of the toolkit in a notebook.
• 01 - Data Loading: Loading CSVs and DataFrames with date formatters.
• 02 - Data Preparation: Handling data, simulating sensor failure, and applying transformations.
• 03 - Parameter Initialization: Methods for initializing parameters for single-output (SM kernel) and multi-output (MOSM kernel) cases.
• 04 - Model Training: Training models while keeping specific parameters fixed.
• 05 - Error Metrics: Comparing models using various metrics.
• 06 - Custom Kernels and Mean Functions: Creating and training custom kernels and mean functions.
• 07 - Sparse Multi Input: Using sparse GPs with multiple input dimensions.
• 08 - Multi Likelihood Classification: Using different likelihoods (e.g., Bernoulli for classification and Student-T for regression) for different channels.

To set up a sparse Gaussian Process model, you need to define three components:

1. Kernel: Define the covariance function using mogptk.gpr classes (e.g., SquaredExponentialKernel) and specify the number of input_dims.
2. Inference: Choose a sparse inference method. For example, mogptk.Titsias allows you to specify the number of inducing_points and how to initialize them (e.g., init_inducing_points='random').
3. Model: Combine the data, kernel, and inference method into a mogptk.Model object.

You can manually initialize kernel parameters (like lengthscale) using the .assign() method to assist training convergence.

# Choose kernel
kernel = mogptk.gpr.SquaredExponentialKernel(input_dims=7)

# Choose inference
inference = mogptk.Titsias(inducing_points=100, init_inducing_points='random')

# Set up model
model = mogptk.Model(data, kernel, inference=inference)

# Initialize parameters
model.gpr.kernel.lengthscale.assign(np.random.rand(7))

To create a custom mean function that can be optimized by model.train, you must create a class that inherits from mogptk.gpr.Mean.

Key requirements:

1. Call super().__init__() in the __init__ method.
2. Define trainable parameters by assigning instances of mogptk.gpr.Parameter to class properties.
3. Implement the __call__(self, X) method, where X has shape (data points, input dims). The method must return Y values of shape (datapoints,).
4. Crucial: When using parameters inside __call__, you must 'call' them as methods (e.g., coefs = self.coefficients()) to retrieve their current values.

Note: mogptk.LoadFunction can be used to generate a Data object from a function, which is useful for creating ground-truth signals for testing.

import mogptk
import torch
import numpy as np

# Create ground truth data
f = lambda x: np.sin(x*2.0*np.pi) + 2*x - 0.2*x**2
data = mogptk.LoadFunction(f, start=0.0, end=10.0, n=100, var=0.5)

# Define custom mean class
class Mean(mogptk.gpr.Mean):
    def __init__(self):
        super(Mean, self).__init__()        
        self.coefficients = mogptk.gpr.Parameter([0.0, 0.0, 0.0])

    def __call__(self, X):
        # Must call the parameter as a method to get values
        coefs = self.coefficients()
        return coefs[0] + coefs[1]*X[:,1] + coefs[2]*X[:,1]**2

When using custom kernels, random initialization may lead to poor optimization. It is recommended to manually assign reasonable starting values to kernel parameters using the .assign() method.

To control the training process, you can toggle the .train attribute on the mean object, the kernel object, or individual Parameter objects. This allows for staged training (e.g., training the mean function first, then the kernel).

To train the model, use model.train() with arguments such as method (e.g., 'Adam'), lr (learning rate), iters (iterations), and error (e.g., 'MAE').

# Setup
mean = Mean()
kernel = mogptk.gpr.PeriodicKernel(input_dims=1)
mo_kernel = mogptk.gpr.IndependentMultiOutputKernel(kernel)
model = mogptk.Model(data, mo_kernel, mean=mean, name="Periodic")

# Manually initialize kernel parameters
kernel.lengthscale.assign(1.0)
kernel.period.assign(1.0)

# Stage 1: Train only the mean
mean.train = True
kernel.train = False
model.train(method='Adam', lr=0.05, iters=250, plot=True, error='MAE')

# Stage 2: Train only the kernel
mean.train = False
kernel.train = True
model.train(method='Adam', lr=0.05, iters=250, plot=True, error='MAE')

# Stage 3: Train both
mean.train = True
kernel.train = True
model.train(method='Adam', lr=0.05, iters=250, plot=True, error='MAE')

To use MOGPTK for Gaussian Process Regression, you can load data from a pandas DataFrame using the LoadDataFrame class. This class maps specific columns to your independent variable ($x$) and dependent variables ($y$).

Before loading, it is common practice to:

1. Replace invalid sensor readings (e.g., -200 or 200.0) with np.nan.
2. Resample/downsample the time series to a consistent frequency.
3. Normalize or detrend the data using TransformStandard() or TransformDetrend().

Example workflow:

import pandas as pd
import numpy as np
from mogptk import LoadDataFrame, TransformStandard

# ... (assume df is a prepared pandas DataFrame with a 'timestamp' column)

# Load data
data = LoadDataFrame(
    df, 
    x_col=['timestamp'],
    y_col=['NOx(GT)', 'CO(GT)']
)

# Select a single output and apply a transformation
data_nox = data['NOx(GT)']
data_nox.transform(TransformStandard())

```python
data = LoadDataFrame(
    df, x_col=['timestamp'],
    y_col=cols_meas
)

data_nox = data['NOx(GT)']
data_nox.remove(pct=0.25) # ratio of test data
data_nox.transform(TransformStandard())
```埋

To train a multi-output dataset where each channel requires a different likelihood (e.g., combining classification and regression), use mogptk.gpr.MultiOutputLikelihood. This allows you to pass specific likelihood instances for each output dimension.

In this pattern, you typically:

1. Create a mogptk.DataSet containing multiple mogptk.Data objects.
2. Define a mogptk.gpr.MultiOutputLikelihood containing the desired likelihoods in order.
3. Use an inference engine capable of handling non-Gaussian likelihoods, such as mogptk.Hensman.
4. Define a multi-output kernel like mogptk.gpr.LinearModelOfCoregionalizationKernel (LMC) to model cross-channel correlations.
5. Initialize a mogptk.Model and call .train().

import numpy as np
import torch
import mogptk

# 1. Prepare data
data = mogptk.DataSet(
    mogptk.Data(X1, y1), # e.g. Classification data
    mogptk.Data(X2, y2), # e.g. Regression data
)

# 2. Define multi-output likelihood
likelihood = mogptk.gpr.MultiOutputLikelihood(
    mogptk.gpr.BernoulliLikelihood(),
    mogptk.gpr.StudentTLikelihood(),
)

# 3. Choose inference (Hensman supports non-Gaussian likelihoods)
inference = mogptk.Hensman(likelihood=likelihood)

# 4. Choose kernel (LMC kernel for cross-channel features)
kernel = mogptk.gpr.LinearModelOfCoregionalizationKernel(
    mogptk.gpr.SquaredExponentialKernel(),
    mogptk.gpr.SquaredExponentialKernel(),
    output_dims=2,
)

# 5. Set up and train model
model = mogptk.Model(data, kernel, inference=inference)
model.train(iters=500, verbose=True, error='mape', lr=0.02)

# 6. Visualize
model.plot_prediction(n=100000)