等渗回归和PAVA算法

磐创AI

发布于 2021-04-21 10:46:35

3.3K1

发布于 2021-04-21 10:46:35

文章被收录于专栏：磐创AI技术团队的专栏磐创AI技术团队的专栏

磐创AI分享

来源 | analyticsvidhya

作者 | NOOBMASTER21

编译 | Flin

介绍

等渗回归是很少被谈论但肯定是最酷的回归技术之一。我之所以说“很少谈论”，是因为与线性回归不同，它不经常被讲授或使用。等渗回归做出一个更笼统的假设，即最能代表数据的函数是单调的，而不是线性的（是的，线性也是单调的，反之亦然）。

“Isotonic（等渗）”一词源自两个希腊词根：“ iso ”和“ tonos ”；“ iso ”的字面意思是相等，“ tonos ”的意思是伸展。

因此，等渗回归（也称为等长回归），对数据拟合一个分段常数非递减（阶梯状）函数，因此提供了线性回归的替代方法，线性回归本质上对数据拟合一条直线。

与线性回归相比，这是等渗回归的样子。

如果你像我一样，那么你可能想知道算法背后的数学知识，然后才能开始使用它们。对于任何希望使用该算法获得最佳结果的人来说，确切地知道是什么原因导致了算法的行为方式是一个好习惯。

所以这是我们要做的：

首先，我们将阐述等渗回归所解决的问题。
然后，我们将看一些数学并尝试了解解决方案。在这一部分中，我们将学习PAVA算法。
最后，我们将研究它在python中的实现及其主要应用程序。

问题表述

等渗回归解决了以下问题：给定的顺序 n 个数据点 y1,y2….yn ，我们如何才能通过β1,….βn 单调序列总结这一序列？

我们稍后会再讲。首先，我们需要形式化单调序列。这些序列可以用另一种方式来考虑。我们可以根据约束将点所在的空间划分为不同的部分。然后，我们可以对部分中存在的点进行线性回归。这将给我们一个等渗图，类似于上面的图。

要制定单调序列及其含义，我们必须做一个假设。

同方差的正态误差

像其他线性模型一样，我们可以假定这种回归类型中的误差是同方差的。换句话说，所有误差将具有有限的方差。由于误差不依赖于预测值 xi，我们可以制定一个可以有效拟合数据的线性模型。同样，我们可以假设因变量是随机的，并且它们服从正态分布。

和

以上条件是由于单调性约束。

处理重复的预测变量值

我们可以看到，当 xi = xj 意味着 μi = μj （相同的x值在y的分布中有相同的平均值）。

因此，对于所有重复的预测器或x值，我们可以有一个均值参数，换句话说，对于每个唯一的预测器或x值，我们在y中只能有一个均值参数。

它甚至意味着什么？为什么我们要首先考虑这一点？我的意思是我们怎么会有重复的x值？

我将在这里解释问题的一部分，另一部分需要更多解释。那么，重复或相同的x值意味着什么？

正如我们在定义本身中看到的，等渗回归以单调方式拟合数据。因此，在拟合数据时，如果算法发现违反此单调性约束的点，则该点将与相邻的x值合并在一起，以形成我们之前考虑的块或单调序列。

很酷的是，单调序列或块中的所有x值都将具有相同的y值。你可以在上面的等渗曲线中看到这一点，等渗曲线的平坦部分或拉伸部分（没错，等拉伸回归就是从这里得名）。

这些显示为绘图中的步骤。现在，我们要在这些所谓的块中找到y的值。你现在可能已经可以看到PAVA算法是如何解决的。现在，让我们确切地了解y值的计算方式。

令 z1,….,zk 为按排序顺序的唯一x值（以后将具有相同的y值）。

然后，我们可以将所有唯一x值的权重定义为：

因此，现在y值变为：

现在，当我们将y值除以它们各自的权重时，y值的分布将变为：

现在，此分布表示将要预测的唯一y值，其分量为 v1,……,vk。

这些组件将是单调的，因为这是我们必须遵循的约束。

还，

现在我们需要回到数学描述。

负对数似然

你一定听说过线性回归中的最大似然估计及其最终如何给出最佳拟合线。通常，我们尝试使似然函数最大化，但是如果我们取似然函数的对数并将整个表达式乘以-1，则得到的是负对数似然，它会最小化而不是因为-1而最大化。因此，基本上，我们通过最小化来最大化。

我们通过取正态分布（假设y值来自这个分布）的对数乘以-1来获得。

拉格朗日

如果你对约束优化略知一二，那么你很可能听说过拉格朗日函数。每当我们面临优化（在此处最小化）上述成本函数的任务时，在一些约束条件(这里是单调性)下，我们使用拉格朗日乘子。

Lagrange乘数：https://en.wikipedia.org/wiki/Lagrange_multiplier

上式称为拉格朗日方程。求解该方程式将为我们提供负对数似然函数的最小值，从而最终使可能性最大化，从而确保与数据的最佳拟合。

请注意，除了对数似然函数中两个已经存在的术语之外，又增加了一项。这个额外的术语是由于我们正在考虑的约束条件而产生的结果。因此，我们的问题是一个约束优化问题。

KKT条件

对于无约束凸优化问题，我们知道当一个点的梯度为零时可以达到最小值。KKT条件是我们现在面临的约束凸优化问题的全局最小值的等价条件。

这些条件是：

拉格朗日最小化
原问题的可行性
对偶问题的可行性
互补松弛性

对于第一个条件（拉格朗日极小化），

请注意，关于y值或预测的值，我们如何获取Lagrangian的偏导数。

需要注意的是，如果xi=zm ，∂µi/∂νm是等于1，否则等于0。因此，仅当xi = zm （其中z m代表唯一的x值或我们将拥有唯一y值的x值）时，第一行总和中的项才为非零。

可以看出，该公式仅对1<m <k<="" strong="" style="font-size: inherit;color: inherit;line-height: inherit;">有效。以使其适用于1 <= M <= K，我们定义 λ0 =λķ = 0。</k<="">

将上述方程式设置为零将为我们提供第一个 Karush-Kuhn-Tucker 条件。让我们把它设为零，乘以 σ2 ，

这里 κm = σ2λm已被引入。

可以看出，当 κm = σ2λm时，κm 为零、正或负，因此我们可以在其余的KKT条件下使用 κm 。

现在第二个条件（原问题的可行性）：

对于第三个条件（对偶问题的可行性）：

最后，（**互补松弛性 **）条件：

现在你可能想知道为什么我们要引入 κm = σ2λm 。原因是，当我们估计平均参数或y值时，换句话说，我们通常不考虑方差。为了实现不含方差的最大似然估计，我们引入了一个只直接依赖于拉格朗日乘子λm的新变量。

现在我们有了KKT条件，我们准备好计算y值。

定义块

首先，我们仅应用第一个（拉格朗日导数等于零）和第四个（互补松弛性）条件。

我们可以将y值的空间划分为等量连续的块，如果该块中的值不等于任一侧的平均值，则这些块的长度将为1。

此时，我们将唯一的y值称为 ν，它是特定块的平均值。

因此，让我们考虑一个块。

一个块只有一个唯一的y值，在该块中我们将其称为v （均值参数） 。

假设，

上标中的星号用于区分两个不同的块。因此，序列， νj∗ = · · · = νj∗∗−1 表示的块。

现在该边缘的情况下，我们定义 ν0 =-∞和ν K + 1 = +∞。

互补松弛性条件表明，

因此，现在使用第一个条件，

其中

解出上面的式子，

现在，这向我们揭示了一条非常重要且很酷的信息。仅通过应用第一和第四条件，我们发现，

等量块中的平均值是vj值的加权平均值，vj值是块的yi值的未加权平均值。

PAVA算法（The Pool Adjacent Violators Algorithm）

PAVA算法基本上可以执行其名称所暗示的功能。它检查这些点，如果发现违反约束的点，则将其值与相邻成员合并，最终形成一个块。

具体而言，PAVA会执行以下操作：

**[初始化步骤] **将v和κ设置为满足条件1,3和4的任何值。
[终止步骤]如果满足条件2（原始条件），则终止。
[**池相邻的违背者]选择任意的j，使得 **νj > νj+1.。然后“合并”包含j和j + 1的块，使其成为一个块（该合并块的nu或µ值将再次成为该合并块的值的加权平均值）。[此步骤保持条件1、3和4。]
再次执行步骤2。

现在，我们对这种算法的工作原理有了基本了解，让我们在python中实现它：

首先，我们需要导入所需的库：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.isotonic import IsotonicRegression
from matplotlib.collections import LineCollection
import random

然后让我们为其生成一些点和一些随机噪声：

n=100 #no. of points
x=np.arange(n)
y= 50*np.log(1+x) + np.random.randint(-50,50,size=(n,))

现在让我们线性等渗地拟合数据，

ir=IsotonicRegression()
lr=LinearRegression()
y_ir=ir.fit_transform(x,y)
lr.fit(x[:,np.newaxis],y)

生成一些线段并绘制结果：

lines=[[[i,y[i]],[i,y_ir[i]]] for i in range(n)]
lc=LineCollection(lines)
plt.figure(figsize=(16,4))
plt.plot(x,y,'r.',markersize=12)
plt.plot(x,y_ir,'b.-',markersize=12)
plt.plot(x,lr.predict(x[:, np.newaxis]), 'y-')
plt.gca().add_collection(lc)
plt.legend(('Data','Isotonic Fit','Linear Fit'))
plt.title("Isotonic Regression")

最后结果：

人们还可以看到该算法实时工作

应用领域

现在的问题是，等渗回归在哪里使用？

等渗回归通常用于分类器的概率校准。分类器的概率校准处理优化分类器输出的过程，以便可以将分类器模型的输出直接解释为置信度。例如：如果分类器以90％的概率将一定数量的样本归为某个特定类别，则这些样本中大约90％实际上应该属于该特定类别。

主要有两种方法：

Platt’s Scaling：将逻辑回归模型拟合到分类器模型的输出。
等渗回归：使等渗或阶梯状曲线适合模型的输出。

现在，我们将看到等渗回归如何用于概率校准：

注意：尽管等渗回归在校正单调失真方面比Platt Scaling更强大，但它很容易使数据过拟合。线性回归的灵活性带来了需要更多数据的代价。

让我们导入库：

from sklearn.datasets import make_classification
from sklearn.svm import SVC
from sklearn.calibration import CalibratedClassifierCV
from sklearn.model_selection import train_test_split
from sklearn.calibration import calibration_curve

生成数据并将其分为训练集和测试集：

X, y = make_classification(n_samples=10000, n_classes=2, random_state=42)
trainX, testX, trainy, testy = train_test_split(X, y, test_size=0.2, random_state=42)

拟合两种模型：一种不带校准，一种带校准：

model1=SVC()
model1.fit(trainX,trainy)
preds1=model1.decision_function(testX)

model2=SVC()
cal=CalibratedClassifierCV(model2,method='isotonic',cv=5)
cal.fit(trainX,trainy)
preds2=cal.predict_proba(testX)[:,1]
Finding the calibration curves:

查找校准曲线：

fop1, mpp1= calibration_curve(testy, preds1, n_bins=10, normalize=True)
fop2, mpp2= calibration_curve(testy, preds2, n_bins=10)

绘制曲线：

plt.figure(figsize=(12,8))
plt.plot(mpp1,fop1,marker='.')
plt.plot(mpp2,fop2,marker='.')
plt.plot([0,1],[0,1],linestyle='--',color='black')
plt.legend(("Support Vector Classfier","Calibrated SVC","Perfectly Calibrated"))

结果：