一、相关向量机RVM与支持向量机SVM对比

1、相关向量机（RVM）

①定义与原理

相关向量机（Relevance Vector Machine, RVM）是一种基于概率模型的机器学习算法，主要用于分类和回归分析。基于稀疏贝叶斯学习框架，通过自动选择一小部分相关向量来进行回归或分类任务。RVM使用贝叶斯公式对权重进行更新和调整，通过最大化似然函数来更新权重，并利用高斯过程进行建模。这种方法可以在训练时自动调整权重，以适应数据的变化和不确定性。

常用的核函数有rbf（高斯）径向基核函数、liner线性核函数、poly多项式核函数。高斯核函数是径向基核函数的一种，但在RVM中rbf径向基核函数特指高斯核函数。

②结构与组成

RVM由输入层、隐含层和输出层组成。输入层负责接收外部数据，隐含层负责数据计算和加工，输出层负责输出结果。其中，隐含层中的权重和偏置可以通过训练得到，以实现最优的分类或回归效果。

③优势与应用

RVM高效、灵活、易于理解！与传统的支持向量机（SVM）相比，RVM需要的训练时间更短，且可以在处理高维数据时表现更优秀。此外，RVM还具有很好的可解释性，可以帮助用户理解模型的决策过程。

RVM可以应用于各种领域，如金融、医疗、图像处理等。例如，在金融领域中，RVM可以用于信用评级和风险评估；在医疗领域中，RVM可以用于疾病诊断和预测；在图像处理领域中，RVM可以用于图像分类和识别。

---

2、支持向量机（SVM）

①定义与原理

支持向量机（Support Vector Machine, SVM）是一种二分类模型，其目标是寻找一个超平面来对样本进行分割，分割的原则是间隔最大化。SVM通过引入核函数方法，克服了维数灾难和非线性可分的问题，向高维空间映射时没有增加计算的复杂性。SVM的最终决策函数只由少数的支持向量所确定，因此计算的复杂性取决于支持向量的数目，而不是样本空间的维数。

②优势与局限

SVM算法具有在小样本情况下表现优异、泛化能力强等优点。然而，SVM对于大规模训练样本难以实施，因为SVM算法借助二次规划求解支持向量，这其中会涉及m阶矩阵的计算，当矩阵阶数很大时将耗费大量的机器内存和运算时间。此外，经典SVM算法只给出了二分类的算法，对于多分类问题解决效果并不理想。

---

3、RVM与SVM的对比

• 优化目标：
  • RVM：基于贝叶斯框架，通过最大化边缘似然函数来选择相关向量，实现稀疏化。
  • SVM：寻找最大间隔的超平面进行样本分割。
• 计算复杂度：
  • RVM：在处理高维数据时表现更优秀，训练时间相对较短。
  • SVM：对于大规模训练样本，计算复杂度较高，特别是在高维空间中。
• 模型稀疏性：
  • RVM：通过自动选择相关向量，模型更加稀疏。
  • SVM：虽然决策函数由支持向量决定，但支持向量的数量可能较多，不一定总是稀疏的。
• 输出类型：
  • RVM：可以给出概率型的输出，具有更好的可解释性。
  • SVM：输出为分类标签或回归值，不直接提供概率信息。
• 适用场景：
  • RVM：由于其高效、灵活和稀疏性，适用于需要快速训练和解释的场景。
  • SVM：在小样本情况下表现优异，广泛应用于各种分类和回归问题。

---

二、Python算例——RVM回归预测&分类

1、定义RVM基类——BaseRVM

由于sklearn中没现成的相关向量机库，可以借助sklearn-SVM库的API（继承它一些功能），再根据需要编写一些独属RVM的类，即可实现RVM模型。由于自己学习笔记的特征，博主特地加了详细注释：

import numpy as np

from scipy.optimize import minimize
from scipy.special import expit

from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin
from sklearn.metrics.pairwise import (
    linear_kernel,
    rbf_kernel,
    polynomial_kernel
)
from sklearn.multiclass import OneVsOneClassifier
from sklearn.utils.validation import check_X_y

#sklearn中不存在标准的RVM需要的类，属于SVM的变种。需要自己定义RVM的类，回归和分类都有。(对象的三大属性：封装、继承、多态)
class BaseRVM(BaseEstimator):
    """提供RVM共通的方法或属性——无论分类还是回归都会用到
    继承了scikit-learn API中的BaseEstimator类。添加后验加权方法和预测,用于分类或回归的子类。
    """

    def __init__(
            self,
            kernel='rbf',
            degree=3,
            coef1=None,
            coef0=0.0,
            n_iter=3000,
            tol=1e-3,
            alpha=1e-6,
            threshold_alpha=1e9,
            beta=1.e-6,
            beta_fixed=False,
            bias_used=True,
            verbose=False
    ):   #定义属性变量时，已经设置好了默认值
        """类定义的一般格式。使用时通过对象.属性名的方式访问，如BaseRVM.kernel（类中的函数也是：对象.函数名方式访问使用）"""
        self.kernel = kernel  #基函数（核函数），一般rbf为径向基函数，poly为多项式核
        self.degree = degree  #当kernel为'poly'时，此参数指定多项式的度。对于其他类型的核，此参数将被忽略。
        self.coef1 = coef1    #当使用'rbf'或'poly'核时，此参数用于指定核函数的第一个系数。对于'rbf'，影响径向基函数的宽度或平滑度。对于'poly'，影响多项式核的系数。如果设置为None，则使用默认的核函数参数
        self.coef0 = coef0    #当使用'poly'（多项式）核时，此参数指定多项式核的常数项。对于其他类型的核，此参数将被忽略
        self.n_iter = n_iter  #迭代的最大次数。RVM通过迭代优化过程来找到最佳的相关性向量
        self.tol = tol        #收敛容差。当算法连续迭代中alpha值的变化小于此容差时，认为算法已经收敛，停止迭代
        self.alpha = alpha    #所有基础函数的初始alpha值（即每个基础函数的先验精度）。在RVM中，alpha值越大，对应的基础函数越有可能成为相关性向量。
        self.threshold_alpha = threshold_alpha  #用于裁剪（prune）过程的阈值。alpha值大于此阈值的基础函数将被保留为相关性向量
        self.beta = beta      #噪声精度的初始值。在RVM中，beta值控制了对观测噪声的估计
        self.beta_fixed = beta_fixed #指定beta值是否在整个训练过程中保持不变。如果为True，则beta值在训练过程中不会更新
        self.bias_used = bias_used   #指定是否在模型中包含偏置项（即常数项）。对于许多实际问题，包含偏置项可以提高模型的灵活性
        self.verbose = verbose       #控制是否打印训练过程中的额外信息。如果为True，则训练过程中会打印迭代次数、alpha值、beta值等信息.

    def get_params(self, deep=True):
        """将类的所有参数（即属性）以字典的形式返回。这是遵循scikit-learn API的一种常见做法，使得模型的参数可以很容易地被获取和检查。
        deep参数通常用于控制是否递归地获取嵌套对象（例如，如果模型的某些属性也是具有get_params方法的scikit-learn兼容对象）的参数。
        然而，在这个例子中，BaseRVM类并没有这样的嵌套对象，因此deep参数没有被用到。
        """
        params = {
            'kernel': self.kernel,
            'degree': self.degree,
            'coef1': self.coef1,
            'coef0': self.coef0,
            'n_iter': self.n_iter,
            'tol': self.tol,
            'alpha': self.alpha,
            'threshold_alpha': self.threshold_alpha,
            'beta': self.beta,
            'beta_fixed': self.beta_fixed,
            'bias_used': self.bias_used,
            'verbose': self.verbose
        }
        return params

    def set_params(self, **parameters):
        """set_params方法也是scikit-learn API中的一部分，用于通过关键字参数（kwargs）设置对象的参数（即属性）。
        set_params方法使得用户可以一次性地更新模型的多个参数，而不需要单独调用每个属性的setter方法。
        **parameters可变关键字参数允许传入任意数量的关键字参数。
    通过遍历parameters字典的项（即参数名和对应的值），使用setattr函数动态地将每个参数名（parameter）作为属性名，将对应的值（value）作为属性值，设置到当前对象（self）上。
        """
        for parameter, value in parameters.items():
            setattr(self, parameter, value)  #setattr函数设置对象的属性值。第一个参数是要设置属性的对象，第二个参数是属性名（作为字符串），第三个参数是要设置的新属性值。
        return self  #set_params方法返回当前对象（self），这允许链式调用，即可以在一个表达式中连续调用多个方法。例如有一个模型对象model想要一次性设置kernel和degree参数：model.set_params(kernel='rbf', degree=3)

    def _apply_kernel(self, x, y):
        """根据选择的核函数，进行对应的运算。线性、径向基、多项式、自定义核"""
        if self.kernel == 'linear':
            phi = linear_kernel(x, y)
        elif self.kernel == 'rbf':
            phi = rbf_kernel(x, y, self.coef1)
        elif self.kernel == 'poly':
            phi = polynomial_kernel(x, y, self.degree, self.coef1, self.coef0)
        elif callable(self.kernel):
            phi = self.kernel(x, y)
            if len(phi.shape) != 2:
                raise ValueError(
                    "Custom kernel function did not return 2D matrix"
                )
            if phi.shape[0] != x.shape[0]:
                raise ValueError(
                    "Custom kernel function did not return matrix with rows"
                    " equal to number of data points."""
                )
        else:
            raise ValueError("Kernel selection is invalid.")

        if self.bias_used: #如果self.bias_used为True，则在返回的核函数结果phi上添加一个偏置项。
            phi = np.append(phi, np.ones((phi.shape[0], 1)), axis=1)  #通过在phi矩阵的末尾添加一列全1来实现的，使用np.append函数沿着列方向（axis=1）添加
        return phi

    def _prune(self):
        """根据 alpha 值的阈值来移除基础函数"""
        # 根据alpha值是否小于阈值threshold_alpha来确定哪些基础函数需要保留
        keep_alpha = self.alpha_ < self.threshold_alpha

        # 如果没有任何alpha值小于阈值，则至少保留第一个基础函数
        # 如果模型中使用了偏置项（bias_used为True），则也保留最后一个基础函数（通常是偏置项）
        if not np.any(keep_alpha):
            keep_alpha[0] = True
            if self.bias_used:
                keep_alpha[-1] = True

        # 如果模型中使用了偏置项
        if self.bias_used:
            # 如果最后一个基础函数（偏置项）不在保留列表中，则设置bias_used为False
            if not keep_alpha[-1]:
                self.bias_used = False
            # 更新relevance_数组，移除不需要的基础函数的相关性信息
            # 注意这里使用keep_alpha[:-1]来排除偏置项（如果存在且不需要保留）
            self.relevance_ = self.relevance_[keep_alpha[:-1]]
        else:
            # 如果模型中没有使用偏置项，则直接根据keep_alpha来更新relevance_数组
            self.relevance_ = self.relevance_[keep_alpha]

        # 更新alpha_数组，只保留需要的基础函数的alpha值
        self.alpha_ = self.alpha_[keep_alpha]
        # 更新alpha_old数组，与alpha_相同，只保留需要的值
        self.alpha_old = self.alpha_old[keep_alpha]
        # 更新gamma数组，只保留与保留的基础函数相关的值
        self.gamma = self.gamma[keep_alpha]
        # 更新phi矩阵，只保留与保留的基础函数相关的列
        self.phi = self.phi[:, keep_alpha]
        # 更新sigma_矩阵，这是一个协方差矩阵，只保留与保留的基础函数相关的行和列
        self.sigma_ = self.sigma_[np.ix_(keep_alpha, keep_alpha)]
        # 更新m_数组，只保留需要的基础函数的m值（可能是某种均值或参数）
        self.m_ = self.m_[keep_alpha]
        """通过修改与这些基础函数相关的各种属性（如 alpha_, gamma, phi, sigma_, m_ 等）来实现对基础函数的剪枝（即移除不重要的基础函数）
        贝叶斯学习的特性:通过自动确定哪些输入特征（基础函数）对模型预测最为重要，从而实现模型的稀疏性"""

    def fit(self, X, y):
        """Fit the RVR to the training data."""
        # 检查输入数据X和y是否符合要求
        X, y = check_X_y(X, y)

        # 获取样本数量和特征数量
        n_samples, n_features = X.shape

        # 应用核函数计算phi矩阵，phi矩阵是输入数据X通过核函数转换后的结果
        # 在RVM中，phi矩阵通常用于表示基础函数的激活程度
        self.phi = self._apply_kernel(X, X)

        # 获取基础函数的数量，即phi矩阵的列数
        n_basis_functions = self.phi.shape[1]

        # 初始化relevance_为X（这里可能是个错误，通常relevance_应该用于存储与基础函数相关性的某种度量）
        # 在实际实现中，relevance_可能需要在后续步骤中根据alpha值或其他标准来更新
        self.relevance_ = X
        self.y = y  # 存储目标变量

        # 初始化alpha_（基础函数的权重或重要性度量）为alpha（一个超参数）乘以全1向量
        self.alpha_ = self.alpha * np.ones(n_basis_functions)
        self.beta_ = self.beta  # 初始化beta_（噪声精度或噪声方差的倒数）

        # 初始化m_（基础函数的均值或参数）为零向量
        self.m_ = np.zeros(n_basis_functions)

        # 保存旧的alpha_值以便后续比较
        self.alpha_old = self.alpha_

        # 执行迭代优化过程
        for i in range(self.n_iter):
            # 执行后验更新步骤，通常涉及计算m_和sigma_（协方差矩阵）
            self._posterior()

            # 更新gamma，它可能与alpha_和sigma_的对角线元素有关
            self.gamma = 1 - self.alpha_ * np.diag(self.sigma_)

            # 更新alpha_，通常基于gamma和m_的平方
            self.alpha_ = self.gamma / (self.m_ ** 2)

            # 如果beta_不是固定的，则更新beta_
            if not self.beta_fixed:
                # beta_的更新公式通常涉及样本数量、gamma的总和以及残差平方和
                self.beta_ = (n_samples - np.sum(self.gamma)) / (
                    np.sum((y - np.dot(self.phi, self.m_)) ** 2))

                # 执行剪枝步骤，移除不重要的基础函数
            self._prune()

            # 如果启用了详细输出，则打印当前迭代的信息
            if self.verbose:
                print("Iteration: {}".format(i))
                print("Alpha: {}".format(self.alpha_))
                print("Beta: {}".format(self.beta_))
                print("Gamma: {}".format(self.gamma))
                print("m: {}".format(self.m_))
                print("Relevance Vectors: {}".format(self.relevance_.shape[0]))
                print()

            # 检查alpha_的变化是否小于容忍度tol，并且迭代次数大于1，如果是，则提前终止迭代
            delta = np.amax(np.absolute(self.alpha_ - self.alpha_old))

            if delta < self.tol and i > 1:
                break

                # 更新alpha_old为当前alpha_的值，以便下一次迭代比较
            self.alpha_old = self.alpha_

            # 如果模型中使用了偏置项，则将偏置项的值存储在self.bias中
        if self.bias_used:
            self.bias = self.m_[-1]
        else:
            self.bias = None

        return self  # 返回拟合后的模型实例

2、定义RVM-回归Regression类——RVR

继承BaseRVM基类功能，实现RVM的回归预测功能。

class RVR(BaseRVM, RegressorMixin):
    """相关向量机RVM回归问题Regression——RVR
    Implementation of Mike Tipping's Relevance Vector Machine for regression
    借助scikit-learn API
    """
    def _posterior(self):
        """计算权重（或基础函数参数）的后验分布。通过求解后验均值（m_）和后验协方差（sigma_）来实现"""
        i_s = np.diag(self.alpha_) + self.beta_ * np.dot(self.phi.T, self.phi) #中间变量i_s
        self.sigma_ = np.linalg.inv(i_s)  #通过中间变量i_s计算后验协方差矩阵sigma_
        self.m_ = self.beta_ * np.dot(self.sigma_, np.dot(self.phi.T, self.y)) #计算后验均值向量m_ ，其中 self.y 是目标变量。

    def predict(self, X, eval_MSE=False):
        """对新的输入数据 X 进行预测，并可选地计算预测的均方误差（MSE）"""
        phi = self._apply_kernel(X, self.relevance_)  #计算新的输入数据 X 与当前相关性向量（self.relevance_）之间的核函数值，存储在 phi 中
        y = np.dot(phi, self.m_) #预测y值
        if eval_MSE:
            MSE = (1 / self.beta_) + np.dot(phi, np.dot(self.sigma_, phi.T))
            return y, MSE[:, 0] #返回预测值 y 和MSE的每一行的第一个元素（假设MSE是一个二维数组，每行对应一个输入样本的MSE）
        else:
            return y

3、定义RVM-分类Classification类——RVC

继承BaseRVM基类功能，实现RVM的分类功能。

class RVC(BaseRVM, ClassifierMixin):
    """相关向量机RVM分类问题——简称RVC
    Implementation of Mike Tipping's Relevance Vector Machine for
   借助 scikit-learn API.
    """

    def __init__(self, n_iter_posterior=50, **kwargs):
        """Copy params to object properties, no validation."""
        self.n_iter_posterior = n_iter_posterior
        super(RVC, self).__init__(**kwargs)

    def get_params(self, deep=True):
        """Return parameters as a dictionary."""
        params = super(RVC, self).get_params(deep=deep)
        params['n_iter_posterior'] = self.n_iter_posterior
        return params

    def _classify(self, m, phi):
        return expit(np.dot(phi, m))

    def _log_posterior(self, m, alpha, phi, t):

        y = self._classify(m, phi)

        log_p = -1 * (np.sum(np.log(y[t == 1]), 0) +
                      np.sum(np.log(1 - y[t == 0]), 0))
        log_p = log_p + 0.5 * np.dot(m.T, np.dot(np.diag(alpha), m))

        jacobian = np.dot(np.diag(alpha), m) - np.dot(phi.T, (t - y))

        return log_p, jacobian

    def _hessian(self, m, alpha, phi, t):
        y = self._classify(m, phi)
        B = np.diag(y * (1 - y))
        return np.diag(alpha) + np.dot(phi.T, np.dot(B, phi))

    def _posterior(self):
        result = minimize(
            fun=self._log_posterior,
            hess=self._hessian,
            x0=self.m_,
            args=(self.alpha_, self.phi, self.t),
            method='Newton-CG',
            jac=True,
            options={
                'maxiter': self.n_iter_posterior
            }
        )

        self.m_ = result.x
        self.sigma_ = np.linalg.inv(
            self._hessian(self.m_, self.alpha_, self.phi, self.t)
        )

    def fit(self, X, y):
        """Check target values and fit model."""
        self.classes_ = np.unique(y)
        n_classes = len(self.classes_)

        if n_classes < 2:
            raise ValueError("Need 2 or more classes.")
        elif n_classes == 2:
            self.t = np.zeros(y.shape)
            self.t[y == self.classes_[1]] = 1
            return super(RVC, self).fit(X, self.t)
        else:
            self.multi_ = None
            self.multi_ = OneVsOneClassifier(self)
            self.multi_.fit(X, y)
            return self

    def predict_proba(self, X):
        """Return an array of class probabilities."""
        phi = self._apply_kernel(X, self.relevance_)
        y = self._classify(self.m_, phi)
        return np.column_stack((1 - y, y))

    def predict(self, X):
        """Return an array of classes for each input."""
        if len(self.classes_) == 2:
            y = self.predict_proba(X)
            res = np.empty(y.shape[0], dtype=self.classes_.dtype)
            res[y[:, 1] <= 0.5] = self.classes_[0]
            res[y[:, 1] >= 0.5] = self.classes_[1]
            return res
        else:
            return self.multi_.predict(X)

4、RVR算例——boston房价数据回归预测

boston数据包含10+个特征变量X和一个目标变量y（房价）。代码只能输入多维特征X（二维矩阵），不能输入一个特征x。

对比了SVM和RVM分别使用不同核函数的回归预测效果，评价指标是R^2。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.model_selection import KFold, StratifiedKFold
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import plot_confusion_matrix
from sklearn.svm import SVC
from sklearn.svm import SVR
from sklearn.datasets import load_boston
from sklearn.datasets import load_breast_cancer
plt.rcParams['font.sans-serif'] = ['SimSun'] # 指定默认字体为宋体
plt.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
#回归测试——波士顿房价数据集（多个特征变量X构成一个二维矩阵，目前不支持单一的特征变量；一个目标变量y房价）
X, y = load_boston(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)
scaler = StandardScaler()
scaler.fit(X_train)
X_train_s = scaler.transform(X_train)
X_test_s = scaler.transform(X_test)

#支持向量机SVM效果：（不同核函数）
 #线性核函数
model = SVR(kernel="linear")
model.fit(X_train_s, y_train)
print(f"SVM-线性核函数R^2 score:{model.score(X_test_s, y_test)}")

#二次多项式核
model = SVR(kernel="poly", degree=2)
model.fit(X_train_s, y_train)
print(f"SVM-二次多项式核函数R^2 score:{model.score(X_test_s, y_test)}")

#三次多项式
model = SVR(kernel="poly", degree=3)
model.fit(X_train_s, y_train)
print(f"SVM-三次多项式核函数R^2 score:{model.score(X_test_s, y_test)}")
##（高斯）径向基核函数
model = SVR(kernel="rbf")
model.fit(X_train_s, y_train)
print(f"SVM-径向基核函数R^2 score:{model.score(X_test_s, y_test)}")
#S核
model = SVR(kernel="sigmoid")
model.fit(X_train_s, y_train)
print(f"SVM-sigmoid核函数R^2 score:{model.score(X_test_s, y_test)}")


# 相关向量机（RVM）效果：

#线性核函数
model = RVR(kernel="linear")
model.fit(X_train_s, y_train)
print(f"RVM-线性核函数R^2 score:{model.score(X_test_s, y_test)}")

#（高斯）径向基核函数
model = RVR(kernel="rbf")
model.fit(X_train_s, y_train)
print(f"RVM-径向基核函数R^2 score:{model.score(X_test_s, y_test)}")

#多项式核函数，上面基类中默认时3次多项式）
model = RVR(kernel="poly")
model.fit(X_train_s, y_train)
print(f"RVM-三次多项式核函数R^2 score:{model.score(X_test_s, y_test)}")

运行结果如下，可见RVM使用径向基核函数时0.92302表现最佳，选用RVM-rbf核函数做最终预测。

5、选定得分最高的RVR-rbf核函数作为最终预测

并加入可视化方式：实际值-预测值图&残差图。并输入新的X特征，预测得到房价。

#（高斯）径向基核函数
model = RVR(kernel="rbf")
model.fit(X_train_s, y_train)
print(f"RVM-径向基核函数R^2 score:{model.score(X_test_s, y_test)}")
#结果可视化：实际值-预测值散点图，越接近y=x直线说明越好。残差图：预测值与实际值的差异
y_pred = model.predict(X_test_s)
plt.scatter(y_test, y_pred)
plt.xlabel('Actual Price')
plt.ylabel('Predicted Price')
plt.title('Actual vs Predicted Prices')
plt.show()  #实际值-预测值散点图
plt.scatter(y_pred, y_test - y_pred)
plt.xlabel('Predicted Price')
plt.ylabel('Residuals')
plt.title('Residual Plot')
plt.axhline(y=0, color='black', lw=1)  # 添加水平线表示0残差
plt.show() #残差图
#利用模型预测新的房价y值
X_new = np.array([[
    0.00632,  # 假设的CRIM值
    18.0,     # 假设的ZN值
    2.31,     # 假设的INDUS值
    0,        # 假设的CHAS值（0表示不在查尔斯河边）
    0.538,    # 假设的NOX值
    6.575,    # 假设的RM值
    65.2,     # 假设的AGE值
    4.0900,   # 假设的DIS值
    1,        # 假设的RAD值
    296.0,    # 假设的TAX值
    15.3,     # 假设的PTRATIO值
    396.90,   # 假设的B值（根据公式计算或从某处获得）
    4.98      # 假设的LSTAT值
]])
X_new_s = scaler.transform(X_new)
predicted_price = model.predict(X_new_s)
print(f"Predicted Price: {predicted_price[0]}")

实际价格与预测价格关系图如下，理论上应该接近y=x这条直线：

预测价格与残差的关系图，理论上应该接近0：

输入X特征向量二维数组，利用拟合的模型预测的价格结果如下：