• Python机器学习16——相关向量机(RVM)

    本系列基本不讲数学原理，只从代码角度去让读者们利用最简洁的Python代码实现机器学习方法。

    背景介绍

    学机器学习的应该都知道支持向量机（SVM），这个方法在深度学习兴起之前算是很热门的分类方法，在机器学习里面，分类算法属SVM效果比较好，回归算法属随机森林(RF)的效果比较好。

    虽然目前在深度学习神经网络算法的面前，它们的效果都已经黯然失色。但是学术界还是不少人使用这些传统的算法，因为数学理论强，很受老师们喜欢。

    相关向量机（RVM）也是，它是一种基于贝叶斯框架的方法，核心思想是先验后验概率找最大似然，在科研领域中，关于使用RVM的文章不在少数。

    但是由于我们最常用的sklearn库没有rvm的接口......所以没办法直接调用，这篇博客就是补充这个空白的。rvm使用Python实现相关向量机，并且做成我们最为熟悉的sklearn库接口，方便使用。

     

    算法原理

    我这里就不多介绍原理了，感兴趣同学直接看这个pdf，讲的还是不错的。

    https://files-cdn.cnblogs.com/files/axlute/RVMExplained.pdf

     

    代码实现

    定义RVM的类，回归和分类都有。

    1. """Relevance Vector Machine classes for regression and classification."""
    2. import numpy as np
    3.  
    4. from scipy.optimize import minimize
    5. from scipy.special import expit
    6.  
    7. from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin
    8. from sklearn.metrics.pairwise import (
    9.     linear_kernel,
    10.     rbf_kernel,
    11.     polynomial_kernel
    12. )
    13. from sklearn.multiclass import OneVsOneClassifier
    14. from sklearn.utils.validation import check_X_y
    15.  
    16.  
    17. class BaseRVM(BaseEstimator):
    18.  
    19.     """Base Relevance Vector Machine class.
    21.     Implementation of Mike Tipping's Relevance Vector Machine using the
    22.     scikit-learn API. Add a posterior over weights method and a predict
    23.     in subclass to use for classification or regression.
    24.     """
    25.  
    26.     def __init__(
    27.         self,
    28.         kernel='rbf',
    29.         degree=3,
    30.         coef1=None,
    31.         coef0=0.0,
    32.         n_iter=3000,
    33.         tol=1e-3,
    34.         alpha=1e-6,
    35.         threshold_alpha=1e9,
    36.         beta=1.e-6,
    37.         beta_fixed=False,
    38.         bias_used=True,
    39.         verbose=False
    40.     ):
    41.         """Copy params to object properties, no validation."""
    42.         self.kernel = kernel
    43.         self.degree = degree
    44.         self.coef1 = coef1
    45.         self.coef0 = coef0
    46.         self.n_iter = n_iter
    47.         self.tol = tol
    48.         self.alpha = alpha
    49.         self.threshold_alpha = threshold_alpha
    50.         self.beta = beta
    51.         self.beta_fixed = beta_fixed
    52.         self.bias_used = bias_used
    53.         self.verbose = verbose
    54.  
    55.     def get_params(self, deep=True):
    56.         """Return parameters as a dictionary."""
    57.         params = {
    58.             'kernel': self.kernel,
    59.             'degree': self.degree,
    60.             'coef1': self.coef1,
    61.             'coef0': self.coef0,
    62.             'n_iter': self.n_iter,
    63.             'tol': self.tol,
    64.             'alpha': self.alpha,
    65.             'threshold_alpha': self.threshold_alpha,
    66.             'beta': self.beta,
    67.             'beta_fixed': self.beta_fixed,
    68.             'bias_used': self.bias_used,
    69.             'verbose': self.verbose
    70.         }
    71.         return params
    72.  
    73.     def set_params(self, **parameters):
    74.         """Set parameters using kwargs."""
    75.         for parameter, value in parameters.items():
    76.             setattr(self, parameter, value)
    77.         return self
    78.  
    79.     def _apply_kernel(self, x, y):
    80.         """Apply the selected kernel function to the data."""
    81.         if self.kernel == 'linear':
    82.             phi = linear_kernel(x, y)
    83.         elif self.kernel == 'rbf':
    84.             phi = rbf_kernel(x, y, self.coef1)
    85.         elif self.kernel == 'poly':
    86.             phi = polynomial_kernel(x, y, self.degree, self.coef1, self.coef0)
    87.         elif callable(self.kernel):
    88.             phi = self.kernel(x, y)
    89.             if len(phi.shape) != 2:
    90.                 raise ValueError(
    91.                     "Custom kernel function did not return 2D matrix"
    92.                 )
    93.             if phi.shape[0] != x.shape[0]:
    94.                 raise ValueError(
    95.                     "Custom kernel function did not return matrix with rows"
    96.                     " equal to number of data points."""
    97.                 )
    98.         else:
    99.             raise ValueError("Kernel selection is invalid.")
    100.  
    101.         if self.bias_used:
    102.             phi = np.append(phi, np.ones((phi.shape[0], 1)), axis=1)
    103.  
    104.         return phi
    105.  
    106.     def _prune(self):
    107.         """Remove basis functions based on alpha values."""
    108.         keep_alpha = self.alpha_ < self.threshold_alpha
    109.  
    110.         if not np.any(keep_alpha):
    111.             keep_alpha[0] = True
    112.             if self.bias_used:
    113.                 keep_alpha[-1] = True
    114.  
    115.         if self.bias_used:
    116.             if not keep_alpha[-1]:
    117.                 self.bias_used = False
    118.             self.relevance_ = self.relevance_[keep_alpha[:-1]]
    119.         else:
    120.             self.relevance_ = self.relevance_[keep_alpha]
    121.  
    122.         self.alpha_ = self.alpha_[keep_alpha]
    123.         self.alpha_old = self.alpha_old[keep_alpha]
    124.         self.gamma = self.gamma[keep_alpha]
    125.         self.phi = self.phi[:, keep_alpha]
    126.         self.sigma_ = self.sigma_[np.ix_(keep_alpha, keep_alpha)]
    127.         self.m_ = self.m_[keep_alpha]
    128.  
    129.     def fit(self, X, y):
    130.         """Fit the RVR to the training data."""
    131.         X, y = check_X_y(X, y)
    132.  
    133.         n_samples, n_features = X.shape
    134.  
    135.         self.phi = self._apply_kernel(X, X)
    136.  
    137.         n_basis_functions = self.phi.shape[1]
    138.  
    139.         self.relevance_ = X
    140.         self.y = y
    141.  
    142.         self.alpha_ = self.alpha * np.ones(n_basis_functions)
    143.         self.beta_ = self.beta
    144.  
    145.         self.m_ = np.zeros(n_basis_functions)
    146.  
    147.         self.alpha_old = self.alpha_
    148.  
    149.         for i in range(self.n_iter):
    150.             self._posterior()
    151.  
    152.             self.gamma = 1 - self.alpha_*np.diag(self.sigma_)
    153.             self.alpha_ = self.gamma/(self.m_ ** 2)
    154.  
    155.             if not self.beta_fixed:
    156.                 self.beta_ = (n_samples - np.sum(self.gamma))/(
    157.                     np.sum((y - np.dot(self.phi, self.m_)) ** 2))
    158.  
    159.             self._prune()
    160.  
    161.             if self.verbose:
    162.                 print("Iteration: {}".format(i))
    163.                 print("Alpha: {}".format(self.alpha_))
    164.                 print("Beta: {}".format(self.beta_))
    165.                 print("Gamma: {}".format(self.gamma))
    166.                 print("m: {}".format(self.m_))
    167.                 print("Relevance Vectors: {}".format(self.relevance_.shape[0]))
    168.                 print()
    169.  
    170.             delta = np.amax(np.absolute(self.alpha_ - self.alpha_old))
    171.  
    172.             if delta < self.tol and i > 1:
    173.                 break
    174.  
    175.             self.alpha_old = self.alpha_
    176.  
    177.         if self.bias_used:
    178.             self.bias = self.m_[-1]
    179.         else:
    180.             self.bias = None
    181.  
    182.         return self
    183.  
    184.  
    185. class RVR(BaseRVM, RegressorMixin):
    186.  
    187.     """Relevance Vector Machine Regression.
    189.     Implementation of Mike Tipping's Relevance Vector Machine for regression
    190.     using the scikit-learn API.
    191.     """
    192.  
    193.     def _posterior(self):
    194.         """Compute the posterior distriubtion over weights."""
    195.         i_s = np.diag(self.alpha_) + self.beta_ * np.dot(self.phi.T, self.phi)
    196.         self.sigma_ = np.linalg.inv(i_s)
    197.         self.m_ = self.beta_ * np.dot(self.sigma_, np.dot(self.phi.T, self.y))
    198.  
    199.     def predict(self, X, eval_MSE=False):
    200.         """Evaluate the RVR model at x."""
    201.         phi = self._apply_kernel(X, self.relevance_)
    202.  
    203.         y = np.dot(phi, self.m_)
    204.  
    205.         if eval_MSE:
    206.             MSE = (1/self.beta_) + np.dot(phi, np.dot(self.sigma_, phi.T))
    207.             return y, MSE[:, 0]
    208.         else:
    209.             return y
    210.  
    211.  
    212. class RVC(BaseRVM, ClassifierMixin):
    213.  
    214.     """Relevance Vector Machine Classification.
    216.     Implementation of Mike Tipping's Relevance Vector Machine for
    217.     classification using the scikit-learn API.
    218.     """
    219.  
    220.     def __init__(self, n_iter_posterior=50, **kwargs):
    221.         """Copy params to object properties, no validation."""
    222.         self.n_iter_posterior = n_iter_posterior
    223.         super(RVC, self).__init__(**kwargs)
    224.  
    225.     def get_params(self, deep=True):
    226.         """Return parameters as a dictionary."""
    227.         params = super(RVC, self).get_params(deep=deep)
    228.         params['n_iter_posterior'] = self.n_iter_posterior
    229.         return params
    230.  
    231.     def _classify(self, m, phi):
    232.         return expit(np.dot(phi, m))
    233.  
    234.     def _log_posterior(self, m, alpha, phi, t):
    235.  
    236.         y = self._classify(m, phi)
    237.  
    238.         log_p = -1 * (np.sum(np.log(y[t == 1]), 0) +
    239.                       np.sum(np.log(1-y[t == 0]), 0))
    240.         log_p = log_p + 0.5*np.dot(m.T, np.dot(np.diag(alpha), m))
    241.  
    242.         jacobian = np.dot(np.diag(alpha), m) - np.dot(phi.T, (t-y))
    243.  
    244.         return log_p, jacobian
    245.  
    246.     def _hessian(self, m, alpha, phi, t):
    247.         y = self._classify(m, phi)
    248.         B = np.diag(y*(1-y))
    249.         return np.diag(alpha) + np.dot(phi.T, np.dot(B, phi))
    250.  
    251.     def _posterior(self):
    252.         result = minimize(
    253.             fun=self._log_posterior,
    254.             hess=self._hessian,
    255.             x0=self.m_,
    256.             args=(self.alpha_, self.phi, self.t),
    257.             method='Newton-CG',
    258.             jac=True,
    259.             options={
    260.                 'maxiter': self.n_iter_posterior
    261.             }
    262.         )
    263.  
    264.         self.m_ = result.x
    265.         self.sigma_ = np.linalg.inv(
    266.             self._hessian(self.m_, self.alpha_, self.phi, self.t)
    267.         )
    268.  
    269.     def fit(self, X, y):
    270.         """Check target values and fit model."""
    271.         self.classes_ = np.unique(y)
    272.         n_classes = len(self.classes_)
    273.  
    274.         if n_classes < 2:
    275.             raise ValueError("Need 2 or more classes.")
    276.         elif n_classes == 2:
    277.             self.t = np.zeros(y.shape)
    278.             self.t[y == self.classes_[1]] = 1
    279.             return super(RVC, self).fit(X, self.t)
    280.         else:
    281.             self.multi_ = None
    282.             self.multi_ = OneVsOneClassifier(self)
    283.             self.multi_.fit(X, y)
    284.             return self
    285.  
    286.     def predict_proba(self, X):
    287.         """Return an array of class probabilities."""
    288.         phi = self._apply_kernel(X, self.relevance_)
    289.         y = self._classify(self.m_, phi)
    290.         return np.column_stack((1-y, y))
    291.  
    292.     def predict(self, X):
    293.         """Return an array of classes for each input."""
    294.         if len(self.classes_) == 2:
    295.             y = self.predict_proba(X)
    296.             res = np.empty(y.shape[0], dtype=self.classes_.dtype)
    297.             res[y[:, 1] <= 0.5] = self.classes_[0]
    298.             res[y[:, 1] >= 0.5] = self.classes_[1]
    299.             return res
    300.         else:
    301.             return self.multi_.predict(X)

    好了，下面就可以像别的sklearn库里面的包一样使用了。

     

    代码测试

    我们对分类问题和回归问题都测试一下，并且和支持向量机做对比。

    1. import numpy as np
    2. import pandas as pd
    3. import matplotlib.pyplot as plt
    4. import seaborn as sns
    5. from sklearn.preprocessing import StandardScaler
    6. from sklearn.model_selection import train_test_split
    7. from sklearn.model_selection import KFold, StratifiedKFold
    8. from sklearn.model_selection import GridSearchCV
    9. from sklearn.metrics import plot_confusion_matrix
    10.  
    11. from sklearn.svm import SVC
    12. from sklearn.svm import SVR
    13. from sklearn.datasets import load_boston
    14. from sklearn.datasets import load_breast_cancer

    分类测试

    分类我们使用经典的鸢尾花数据集

    1. iris = load_breast_cancer() #加载数据
    2. X = iris.data
    3. y = iris.target
    4.  
    5. X_train, X_test, y_train, y_test =  train_test_split(X, y, test_size=0.2, stratify=y, random_state=0)
    6.  
    7. scaler = StandardScaler()
    8. scaler.fit(X_train)
    9. X_train_s = scaler.transform(X_train)
    10. X_test_s = scaler.transform(X_test)

    支持向量机，不同核函数的效果：

    1. #线性核函数
    2. model = SVC(kernel="linear", random_state=123)
    3. model.fit(X_train_s, y_train)
    4. print(model.score(X_test_s, y_test))
    5. #二次多项式核
    6. model = SVC(kernel="poly", degree=2, random_state=123)
    7. model.fit(X_train_s, y_train)
    8. print(model.score(X_test_s, y_test))
    9. #三次多项式
    10. model = SVC(kernel="poly", degree=3, random_state=123)
    11. model.fit(X_train_s, y_train)
    12. print(model.score(X_test_s, y_test))
    13. #径向核
    14. model = SVC(kernel="rbf", random_state=123)
    15. model.fit(X_train_s, y_train)
    16. print(model.score(X_test_s, y_test))
    17. #S核
    18. model = SVC(kernel="sigmoid",random_state=123)
    19. model.fit(X_train_s, y_train)
    20. print(model.score(X_test_s, y_test))

    相关向量机（RVM）效果：

    1. model = RVC(kernel="linear")
    2. model.fit(X_train_s, y_train)
    3. print(model.score(X_test_s, y_test))
    4.  
    5. model = RVC(kernel="rbf")
    6. model.fit(X_train_s, y_train)
    7. print(model.score(X_test_s, y_test))
    8.  
    9. model = RVC(kernel="poly")
    10. model.fit(X_train_s, y_train)
    11. print(model.score(X_test_s, y_test))

     

    效果差不多。 

    回归测试

    回归使用波士顿数据集

    1. # Support Vector Regression with Boston Housing Data
    2. X, y = load_boston(return_X_y=True)
    3. X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)
    4.  
    5. scaler = StandardScaler()
    6. scaler.fit(X_train)
    7. X_train_s = scaler.transform(X_train)
    8. X_test_s = scaler.transform(X_test)

    支持向量机效果：（不同核函数）

    1.  #线性核函数
    2. model = SVR(kernel="linear")
    3. model.fit(X_train_s, y_train)
    4. print(model.score(X_test_s, y_test))
    5. #二次多项式核
    6. model = SVR(kernel="poly", degree=2)
    7. model.fit(X_train_s, y_train)
    8. print(model.score(X_test_s, y_test))
    9. #三次多项式
    10. model = SVR(kernel="poly", degree=3)
    11. model.fit(X_train_s, y_train)
    12. print(model.score(X_test_s, y_test))
    13. #径向核
    14. model = SVR(kernel="rbf")
    15. model.fit(X_train_s, y_train)
    16. print(model.score(X_test_s, y_test))
    17. #S核
    18. model = SVR(kernel="sigmoid")
    19. model.fit(X_train_s, y_train)
    20. print(model.score(X_test_s, y_test))

     相关向量机（RVM）效果：

    1. model = RVR(kernel="linear")
    2. model.fit(X_train_s, y_train)
    3. print(model.score(X_test_s, y_test))
    4.  
    5. model = RVR(kernel="rbf")
    6. model.fit(X_train_s, y_train)
    7. print(model.score(X_test_s, y_test))
    8.  
    9. model = RVR(kernel="poly")
    10. model.fit(X_train_s, y_train)
    11. print(model.score(X_test_s, y_test))

    可以看到，在回归问题上，相关向量机比支持向量机的效果要好。

    结论：

    分类用SVM，回归用RVM 

    当然我这里只用了两个sklearn自带的数据集测试，结论肯定有点武断，有兴趣的同学可以用于别的数据集，然后做多次K折交叉验证，进一步对比他们的效果。

  • 相关阅读:
    【allegro 17.4软件操作保姆级教程八】布线操作基础之三
    【python】爬取杭州市二手房销售数据做数据分析【附源码】
    企业知识库有什么价值？
    Java 21 新特性：虚拟线程（Virtual Threads）
    可观测｜时序数据降采样在Prometheus实践复盘
    蓝牙信标的优势及应用场景
    装机指南1.0
    Oracle共享内存不释放
    linux常见命令以及jdk，tomcat环境搭建
    GDPU 数据结构 天码行空9
  • 原文地址：https://blog.csdn.net/weixin_46277779/article/details/127927306