在机器学习领域，高斯过程是一种强大的非参数贝叶斯方法，通常用于回归和分类任务。然而，传统的高斯过程模型要求输入数据必须是固定长度的特征向量。本文将探讨如何通过核函数直接在离散数据结构上应用高斯过程，例如可变长度序列、树和图。
具体来说，本文中的输入变量是一些以可变长度字符串形式存储的基因序列，这些字符串由字母‘A’、‘T’、‘C’和‘G’组成。输出变量在回归任务中是浮点数，在分类任务中是True/False标签。通过定义一个基于R-卷积的核函数，可以在一对字符串的所有字母对之间整合一个二进制字母级核函数，从而计算基因序列之间的相似度。
本文将生成三个图表。第一个图表使用色图可视化核函数的值，即序列之间的相似度，颜色越亮表示相似度越高。第二个图表展示了一个包含6个序列的数据集上的回归结果。这里使用第1、2、4和5个序列作为训练集，对第3和第6个序列进行预测。第三个图表展示了一个分类模型，通过在6个序列上训练并对另外5个序列进行预测。这里的真值仅仅是序列中是否至少包含一个‘A’。模型在这里正确分类了四个，而在一个上失败了。
import numpy as np
from sklearn.base import clone
from sklearn.gaussian_process import GaussianProcessClassifier, GaussianProcessRegressor
from sklearn.gaussian_process.kernels import GenericKernelMixin, Hyperparameter, Kernel
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
kernel = SequenceKernel()
# 序列相似性矩阵在核函数下
X = np.array(["AGCT", "AGC", "AACT", "TAA", "AAA", "GAACA"])
K = kernel(X)
D = kernel.diag(X)
import matplotlib.pyplot as plt
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("核函数下的序列相似性")
plt.show()
# 回归
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="预测")
plt.bar(training_idx, Y[training_idx], width=0.2, color="r", alpha=1, label="训练")
plt.xticks(np.arange(len(X)), X)
plt.title("序列上的回归")
plt.legend()
plt.show()
# 分类
X_train = np.array(["AGCT", "CGA", "TAAC", "TCG", "CTTT", "TGCT"])
# 序列中是否包含 'A'
Y_train = np.array([True, True, True, False, False, False])
gp = GaussianProcessClassifier(kernel=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="训练")
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="真实")
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="预测")
plt.xticks(np.arange(len(X_train) + len(X_test)), np.concatenate((X_train, X_test)))
plt.yticks([-1, 1], [False, True])
plt.title("序列上的分类")
plt.legend()
plt.show()
在上述代码中，首先定义了一个名为SequenceKernel的类，它是一个用于可变长度序列的最小但有效的卷积核。然后，使用这个核函数来计算序列之间的相似度，并在核函数下生成序列相似性矩阵。接下来，展示了如何在序列数据上进行回归和分类任务。
在回归部分，使用前四个序列作为训练集，对剩下的两个序列进行预测。在分类部分，训练了一个模型来预测序列中是否至少包含一个‘A’，并在测试集上进行了预测。