扩展Scikit-Learn -- EuroPython 2014 笔记

Posted 2021-06-28 有数可据

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了扩展Scikit-Learn -- EuroPython 2014 笔记相关的知识，希望对你有一定的参考价值。

本文是EuroPython 2014 嘉宾 Florian Wilhelm 的讲座《Extending Scikit-Learn with your own Regressor》的笔记。我写笔记是因为全部翻译比较费时间，而且很多英文句子翻译过来就很怪了。

Theil Sen 回归

嘉宾首先介绍了Scikit-learn和最小二乘法线性回归。

这里，最小二乘法有个问题，它对异常值(Outliers)比较敏感。

在这里插入图片描述

解决的方法是使用Theil Sen算法。Theil Sen算法的思想是。用斜率的中位数替换斜率的平均值。我们都知道使用平均值统计是会有偏差的，往往和我们的感受相差较大，这是因为马云一个人，和我们很多人平均，大家都成了百万富翁。所以平均值意义不大。中位数才可靠。

Theil Sen的算法非常简单。

在这里插入图片描述

所有的点对，计算斜率
计算斜率的中位数

如果要使用Scikit-Learn编写Theil Sen，就要先了解Scikit-Learn的结构。Florian直接帮我们画了UML：

在这里插入图片描述

还给出了示例代码：

在这里插入图片描述

我的实现

我找来了Theil Sen的Scikit-Leran示例，并作图如下：

在这里插入图片描述

然后，几行代码，就实现了Theil Sen，当然，我这里只管一元一次。

class MyTheilSenRegressor(LinearModel, RegressorMixin):
    def __init__(self):
        ...
        
    def fit(self, X, y):
        n_samples, n_features = X.shape
        
        # get the slop using for loops
        slops = []
        for i in range(n_samples):
            for j in range(i+1, n_samples):
                slop = (y[j] - y[i]) / (X[j][0] - X[i][0])
                slops.append(slop)
        final_slop = np.median(slops)
        x_bar = np.mean(X)
        y_bar = np.mean(y)
        intercept = y_bar - final_slop * x_bar
        self.intercept_ = intercept
        self.coef_ = np.array([final_slop])
        return self

然后，把我的算法放到前面的示例里：

在这里插入图片描述

我发现，我的代码居然还比Scikit-Learn快。其实还可以更快的，我决定用numba。

import numba

def get_slop(X, y):
    slops = []
    for i in range(n_samples):
        for j in range(i+1, n_samples):
            slop = (y[j] - y[i]) / (X[j][0] - X[i][0])
            slops.append(slop)
    final_slop = np.median(slops)
    return final_slop

这里测试结果是26毫秒。

然后，我直接在这个函数上面加上了装饰器@numba.njit()，结果报错。

TypingError: Failed in nopython mode pipeline (step: nopython frontend)
final_slop = np.median(slops)

好吧，我可以把numpy函数拿到numba函数的外面。如下：

@numba.njit()
def get_slops(X, y):
    slops = []
    for i in range(n_samples):
        for j in range(i+1, n_samples):
            slop = (y[j] - y[i]) / (X[j][0] - X[i][0])
            slops.append(slop)
    #final_slop = np.median(slops)
    return slops

def get_slop2(X, y):
    slops = get_slops(X, y)
    return np.median(slops)

这里测试结果是949微妙，比纯python版快了几十倍。

重写算法：

class FastTheilSenRegressor(LinearModel, RegressorMixin):
    def __init__(self):
        ...
        
    def fit(self, X, y):
        n_samples, n_features = X.shape
        
        # get the slop using for loops
        final_slop = get_slop2(X, y)
        x_bar = np.mean(X)
        y_bar = np.mean(y)
        intercept = y_bar - final_slop * x_bar
        self.intercept_ = intercept
        self.coef_ = np.array([final_slop])
        return self