离散数据结构上的高斯过程

简介

高斯过程是机器学习中用于回归和分类任务的一种常用工具。然而，它们通常要求数据为固定长度的特征向量形式，这在某些应用中可能会受到限制。在本实验中，我们将探索如何通过定义一个直接对这些结构进行操作的核函数，将高斯过程应用于可变长度序列，如基因序列。我们将使用 scikit-learn 来实现我们的高斯过程模型。

虚拟机使用提示

虚拟机启动完成后，点击左上角切换到笔记本标签页，以访问 Jupyter Notebook 进行练习。

有时，你可能需要等待几秒钟让 Jupyter Notebook 完成加载。由于 Jupyter Notebook 的限制，操作的验证无法自动化。

如果你在学习过程中遇到问题，随时向 Labby 提问。课程结束后提供反馈，我们会及时为你解决问题。

定义序列核

我们定义一个 SequenceKernel 类，它继承自 scikit-learn 的 Kernel 和 GenericKernelMixin 类。这个核定义了可变长度序列之间的相似性度量。它通过 R 卷积来实现，这涉及到在一对字符串中的所有字母对之间对二元字母级核进行积分。然后，我们可以使用这个核在序列上执行回归和分类任务。

import numpy as np
from sklearn.gaussian_process.kernels import Kernel, Hyperparameter
from sklearn.gaussian_process.kernels import GenericKernelMixin
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.base import clone


class SequenceKernel(GenericKernelMixin, Kernel):
    """
    一个用于可变长度序列的最小（但有效）卷积核。
    """

    def __init__(self, baseline_similarity=0.5, baseline_similarity_bounds=(1e-5, 1)):
        self.baseline_similarity = baseline_similarity
        self.baseline_similarity_bounds = baseline_similarity_bounds

    @property
    def hyperparameter_baseline_similarity(self):
        return Hyperparameter(
            "baseline_similarity", "numeric", self.baseline_similarity_bounds
        )

    def _f(self, s1, s2):
        """
        一对序列之间的核值
        """
        return sum(
            [1.0 if c1 == c2 else self.baseline_similarity for c1 in s1 for c2 in s2]
        )

    def _g(self, s1, s2):
        """
        一对序列之间的核导数
        """
        return sum([0.0 if c1 == c2 else 1.0 for c1 in s1 for c2 in s2])

    def __call__(self, X, Y=None, eval_gradient=False):
        if Y is None:
            Y = X

        if eval_gradient:
            return (
                np.array([[self._f(x, y) for y in Y] for x in X]),
                np.array([[[self._g(x, y)] for y in Y] for x in X]),
            )
        else:
            return np.array([[self._f(x, y) for y in Y] for x in X])

    def diag(self, X):
        return np.array([self._f(x, x) for x in X])

    def is_stationary(self):
        return False

    def clone_with_theta(self, theta):
        cloned = clone(self)
        cloned.theta = theta
        return cloned

可视化序列相似性矩阵

我们可以使用我们的 SequenceKernel 来计算序列之间的相似性矩阵。我们将使用颜色映射来绘制这个矩阵，其中较亮的颜色表示较高的相似性。

import matplotlib.pyplot as plt

X = np.array(["AGCT", "AGC", "AACT", "TAA", "AAA", "GAACA"])

kernel = SequenceKernel()
K = kernel(X)
D = kernel.diag(X)

plt.figure(figsize=(8, 5))
plt.imshow(np.diag(D**-0.5).dot(K).dot(np.diag(D**-0.5)))
plt.xticks(np.arange(len(X)), X)
plt.yticks(np.arange(len(X)), X)
plt.title("Sequence similarity under the kernel")
plt.show()

对序列进行回归

我们可以使用我们的 SequenceKernel 对序列进行回归。我们将使用一个包含 6 个序列的数据集，并使用第 1、2、4 和 5 个序列作为训练集，对第 3 和第 6 个序列进行预测。

X = np.array(["AGCT", "AGC", "AACT", "TAA", "AAA", "GAACA"])
Y = np.array([1.0, 1.0, 2.0, 2.0, 3.0, 3.0])

training_idx = [0, 1, 3, 4]
gp = GaussianProcessRegressor(kernel=kernel)
gp.fit(X[training_idx], Y[training_idx])

plt.figure(figsize=(8, 5))
plt.bar(np.arange(len(X)), gp.predict(X), color="b", label="prediction")
plt.bar(training_idx, Y[training_idx], width=0.2, color="r", alpha=1, label="training")
plt.xticks(np.arange(len(X)), X)
plt.title("Regression on sequences")
plt.legend()
plt.show()

对序列进行分类

我们可以使用我们的 SequenceKernel 对序列进行分类。我们将使用一个包含 6 个序列的数据集，并根据序列中是否存在 'A' 进行训练。然后，我们将对另外 5 个序列进行预测，实际情况就是序列中是否至少有一个 'A'。在这里，我们做出了四个正确的分类，有一个分类错误。

X_train = np.array(["AGCT", "CGA", "TAAC", "TCG", "CTTT", "TGCT"])
## 序列中是否存在 'A'
Y_train = np.array([True, True, True, False, False, False])

gp = GaussianProcessClassifier(kernel)
gp.fit(X_train, Y_train)

X_test = ["AAA", "ATAG", "CTC", "CT", "C"]
Y_test = [True, True, False, False, False]

plt.figure(figsize=(8, 5))
plt.scatter(
    np.arange(len(X_train)),
    [1.0 if c else -1.0 for c in Y_train],
    s=100,
    marker="o",
    edgecolor="none",
    facecolor=(1, 0.75, 0),
    label="training",
)
plt.scatter(
    len(X_train) + np.arange(len(X_test)),
    [1.0 if c else -1.0 for c in Y_test],
    s=100,
    marker="o",
    edgecolor="none",
    facecolor="r",
    label="truth",
)
plt.scatter(
    len(X_train) + np.arange(len(X_test)),
    [1.0 if c else -1.0 for c in gp.predict(X_test)],
    s=100,
    marker="x",
    facecolor="b",
    linewidth=2,
    label="prediction",
)
plt.xticks(np.arange(len(X_train) + len(X_test)), np.concatenate((X_train, X_test)))
plt.yticks([-1, 1], [False, True])
plt.title("Classification on sequences")
plt.legend()
plt.show()

总结

在这个实验中，我们展示了如何使用自定义核函数在可变长度序列上使用高斯过程。我们展示了如何使用 scikit-learn 对序列数据执行回归和分类任务。