相關向量機RVM算法介紹繼承sklearn-SVM-API實現回歸預測算例

一、相關向量機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 npfrom scipy.optimize import minimize
from scipy.special import expitfrom 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 paramsdef 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 phidef _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] = Trueif self.bias_used:keep_alpha[-1] = True# 如果模型中使用了偏置項if self.bias_used:# 如果最后一個基礎函數（偏置項）不在保留列表中，則設置bias_used為Falseif 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_ = Xself.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 = Nonereturn self  # 返回擬合后的模型實例

2、定義RVM-回歸Regression類——RVR

繼承BaseRVM基類功能，實現RVM的回歸預測功能。

class RVR(BaseRVM, RegressorMixin):"""相關向量機RVM回歸問題Regression——RVRImplementation 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_sself.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分類問題——簡稱RVCImplementation 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_posteriorsuper(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_posteriorreturn paramsdef _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, jacobiandef _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.xself.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]] = 1return super(RVC, self).fit(X, self.t)else:self.multi_ = Noneself.multi_ = OneVsOneClassifier(self)self.multi_.fit(X, y)return selfdef 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 reselse: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這條直線：