fast.ai 机器学习笔记(二)(3)https://developer.aliyun.com/article/1482666
给定变量找到最佳分割点
让我们尝试编写找到分割点的函数。因此,我们需要实现find_better_split
。它将接受一个变量的索引,并找出最佳的分割点,确定它是否比我们目前为止的任何分割更好,对于第一个变量,答案总是肯定的,因为到目前为止最好的分割点是没有分割,这是无穷糟糕的。
所以让我们首先确保我们有东西可以进行比较。我们要进行比较的是 scikit-learn 的随机森林。我们需要确保 scikit-learn 的随机森林获得与我们完全相同的数据,因此我们首先创建集成,从中提取一棵树,然后找出这棵树使用了哪个特定的随机样本x
和y
,然后将它们存储起来,以便我们可以将它们传递给 scikit-learn(这样我们就有完全相同的信息)。
ens = TreeEnsemble(x_sub, y_train, 1, 1000) tree = ens.trees[0] x_samp,y_samp = tree.x, tree.y
所以让我们继续使用 scikit-learn 创建一个随机森林。一个树(n_estimators
),一个决策(max_depth
),没有自助采样,所以整个数据集。所以这应该与我们即将创建的东西完全相同。让我们试试看。
m = RandomForestRegressor( n_estimators=1, max_depth=1, bootstrap=False ) m.fit(x_samp, y_samp) draw_tree(m.estimators_[0], x_samp, precision=2)
我们需要定义find_better_split
函数,并且它需要一个变量。让我们定义我们的x
(即自变量),说好,它是树中的所有内容,但只有在这个节点中的那些索引,而在树的顶部是所有内容,只有这一个变量(var_idx
)。然后对于我们的y
,它就是在这个节点中的索引处的因变量。所以这就是我们的x
和y
。
让我们现在逐个检查我们独立变量中的每个值。我会告诉你接下来会发生什么。假设我们的独立变量是 YearMade,它不会按顺序排列。所以我们要去到第一行,然后说好,这里的 YearMade 是 3。那么我要做的是尝试计算如果我们决定以数字 3 为分支时的得分。我需要知道哪些行大于 3,哪些行小于或等于 3,它们将成为我的左侧和右侧。然后我们需要一个得分。我们可以使用很多得分,所以在随机森林中,我们称之为信息增益。信息增益就像我们的得分因为我们将数据分成这两组而变得更好了多少。我们可以用很多方法来计算它:基尼系数、交叉熵、均方根误差等等。如果你考虑一下,有一个均方根误差的替代公式,数学上与一个约束尺度内相同,但稍微容易处理一些,那就是我们要找到一个分割点,使得这两组数据的标准差尽可能低。所以我想找到一个分割点,把所有的猫放在这边,所有的狗放在那边。所以如果这些都是猫,那些都是狗,那么这边的标准差为零,那边的标准差也为零。否则,这是一群完全随机混合的猫和狗,那是一群完全随机混合的猫和狗,它们的标准差会高得多。明白了吗?事实证明,如果找到一个最小化这两组标准差或者具体来说是两个标准差的加权平均的分割点,数学上与最小化均方根误差是相同的。如果你想的话,课后你可以自己证明这一点。
首先,我们需要找到,将其分成两组[37:29]。那么所有大于三的东西在哪里?4、6 和 4。所以我们需要它们价格的标准差。
接下来是标准差小于或等于 3,我们只需取这两者的加权平均值,这就是我们的得分。如果我们在 3 上分割,那就是我们的得分。
然后下一步是尝试在 4 上分割,尝试在 1 上分割,尝试在 6 上分割,多余地再次尝试在 4 上分割,然后再次在 1 上分割,找出哪个效果最好。这就是我们的代码:
def find_better_split(self, var_idx): x,y = self.x.values[self.idxs,var_idx], self.y[self.idxs] for i in range(1,self.n-1): lhs = x<=x[i] rhs = x>x[i] if rhs.sum()==0: continue lhs_std = y[lhs].std() rhs_std = y[rhs].std() curr_score = lhs_std*lhs.sum() + rhs_std*rhs.sum() if curr_score<self.score: self.var_idx,self.score,self.split = var_idx,curr_score,x[i]
我们将逐行进行,假设左侧是x
中小于或等于特定值的任何值。右侧是x
中大于此特定值的每个值。
lhs
和 rhs
中将包含什么数据类型?它们实际上会包含什么?它们将是布尔数组,我们可以将其视为零和一。因此,lhs
将是一个数组,每次它不小于或等于时为 false;否则为 true,而 rhs
将是相反的布尔数组。现在我们不能对空集合取标准差,所以如果没有任何大于这个数字 (x[i]
) 的值,那么 rhs
将全部为 false,这意味着总和为零。在这种情况下,让我们不再继续这一步,因为没有什么可以取标准差,显然这不是一个有用的分割。
假设我们已经走到这一步,现在我们可以计算左侧和右侧的标准差,然后取加权平均值或求和,这两者对于一个标量来说是相同的,因此这就是我们的得分。然后我们可以检查这个得分是否比迄今为止的最佳得分更好,我们迄今为止的最佳得分,最初将其初始化为无穷大,因此最初这是更好的。如果更好,让我们存储所有我们需要的信息:哪个变量找到了这个更好的分割,我们找到的得分是多少,以及我们分割的值是多少。就是这样。如果我们运行这个,并且我正在使用%timeit
,它会看这个命令运行需要多长时间,并试图给出一个统计上有效的度量,这样你就可以看到,它已经运行了 10 次以获得平均值,然后又运行了 7 次以获得运行间的平均值和标准差,所以它花了我 76 毫秒加减 11。
%timeit find_better_split(tree,1) tree ''' 76.6 ms ± 11.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) n: 1000; val:10.079014121552744; score:681.0184057251435; split:3744.0; var:MachineHoursCurrentMeter '''
所以让我们来检查这是否有效。find_better_split(tree, 0)
,0 代表YearMade
,1 代表MachineHoursCurrentMeter
,所以当我们用 1 时,我们得到了MachineHoursCurrentMeter
,得分为 681.0184057251435,然后我们再次用零运行,得到了更好的分数(658),并分割了 1974。
find_better_split(tree,0); tree ''' n: 1000; val:10.079014121552744; score:658.5510186055949; split:1974.0; var:YearMade '''
所以 1974 年,让我们与上面的 scikit-learn 的随机森林进行比较,是的,这棵树也是这样做的。所以我们确认了这种方法给出了与 sklearn 的随机森林相同的结果。你还可以在这里看到值 10.08 与 sklearn 的根节点的值匹配。所以我们有了一个可以处理一个分割的东西。
问题:为什么我们不在 x 上放一个 unique?因为我还没有尝试优化性能。你可以在 Excel 中看到我检查了这个1
两次,4
两次,这是不必要的。
好的,Yannet 已经在考虑性能,这是好事。那告诉我这段代码的计算复杂度是多少?
!
O(n²) 是因为有一个循环和 x<=x[i]
,我们必须检查每个值,看它是否小于 x[i]
。了解如何快速计算计算复杂度是很有用的。我可以保证你做的大多数面试都会要求你即兴计算计算复杂度。而且当你编码时,你希望它成为第二天性。这种技巧基本上是“有循环吗?”如果有,那么显然我们要做这个 n
次,所以涉及到一个 n
。循环里面还有循环吗?如果有,那么你需要将它们两个相乘。在这种情况下,没有。循环里面有任何不是常数时间操作的东西吗?所以你可能会看到一个排序在里面,你只需要知道排序是 nlog(n)
—— 这应该是第二天性的。如果你看到一个矩阵相乘,你需要知道那是什么。在这种情况下,有一些东西在进行逐元素数组操作,所以要留意任何地方,numpy 在对数组的每个值做一些操作。在这种情况下,它正在检查每个 x
的值是否小于一个常数,所以它将不得不这样做 n
次。所以要将这个扩展成一个计算复杂度,你只需要将循环中的事物数量乘以循环内部的最高计算复杂度,n
次 n
是 n²
。
问题:在这种情况下,我们不能只是预先对列表进行排序,然后进行一次nlog(n)
的计算吗?有很多事情可以做来加快速度,所以在这个阶段我们只关心的是计算复杂度。但绝对可以。当然,它肯定不是最好的。所以接下来我们要做的就是这个。就像好吧,n²不太好,所以让我们试着让它变得更好。
所以这是我尝试改进的地方。首先,标准差的方程式是什么?
实际上,我们通常不使用那个公式,因为它要求我们多次计算x
减去平均值。有谁知道只需要 x 和 x²的公式吗?
www.wolframalpha.com/input/?i=standard+deviation
这是一个非常好的知识点,因为现在你可以计算任何东西的方差或标准偏差。你只需要首先抓取列本身。列的平方。只要你把它们存储在某个地方,你就可以立即计算标准偏差。
所以这对我们有用的原因是,如果我们首先对我们的数据进行排序。然后如果你考虑一下,当我们一步一步地向下走时,每一组都与左边的前一组完全相同,只是多了一件东西,右边则少了一件东西。因此,我们只需要跟踪 x 的总和和 x²的总和,我们只需在左边添加一个东西,x²再添加一个东西,在右边移除一个东西。因此,我们不必每次都遍历整个数据集,因此我们可以将其转化为 O(n)算法。这就是我在这里所做的一切:
tree = TreeEnsemble(x_sub, y_train, 1, 1000).trees[0] def std_agg(cnt, s1, s2): return math.sqrt((s2/cnt) - (s1/cnt)**2) def find_better_split_foo(self, var_idx): x,y = self.x.values[self.idxs,var_idx], self.y[self.idxs] sort_idx = np.argsort(x) sort_y,sort_x = y[sort_idx], x[sort_idx] rhs_cnt,rhs_sum,rhs_sum2 = self.n, sort_y.sum(), (sort_y**2).sum() lhs_cnt,lhs_sum,lhs_sum2 = 0,0.,0. for i in range(0,self.n-self.min_leaf-1): xi,yi = sort_x[i],sort_y[i] lhs_cnt += 1; rhs_cnt -= 1 lhs_sum += yi; rhs_sum -= yi lhs_sum2 += yi**2; rhs_sum2 -= yi**2 if i<self.min_leaf or xi==sort_x[i+1]: continue lhs_std = std_agg(lhs_cnt, lhs_sum, lhs_sum2) rhs_std = std_agg(rhs_cnt, rhs_sum, rhs_sum2) curr_score = lhs_std*lhs_cnt + rhs_std*rhs_cnt if curr_score<self.score: self.var_idx,self.score,self.split = var_idx,curr_score,xi
我对数据进行排序,然后我会跟踪右侧的事物数量(rhs_cnt
),右侧事物的总和(rhs_sum
)和右侧的平方和(rhs_sum2
)。最初所有事物都在右侧。因此最初n
是计数,y.sum()
是右侧的总和,y²(y**2
)的总和是右侧的平方和。然后最初左侧没有任何事物,因此为零。然后我们只需要循环遍历每个观察值:
- 左手计数加一,右手计数减一。
- 将值加到左手总和,从右手总和减去。
- 将值的平方加到左手,从右手减去。
现在我们需要小心,因为如果我们说小于或等于一,例如,我们不会停在第一行,而是必须将该组中的所有内容都包括在内。
所以我要做的另一件事是确保下一个值不同于这个值。如果是的话,我会跳过它。所以我只是要再次检查这个值和下一个值不相同(**if** xi==sort_x[i+1]:
)。只要它们不相同,我就可以继续前进,通过传入计数、总和和平方和来计算我的标准偏差。那个公式就在那里:
现在我们可以对右侧进行同样的操作,这样我们就可以像之前一样计算加权平均分数,下面的所有行都是一样的。
所以我们把 O(n²)的算法转换成了 O(n)的算法。一般来说,像这样的东西会给你带来比将某些东西推送到 Spark 集群或者更快的 RAM 或者在 CPU 中使用更多核心等更多价值。这是你想要改进你的代码的方式。具体来说,编写代码时不要过多考虑性能。运行它,看看它是否对你的需求足够快。如果是,那么你就完成了。如果不是,进行性能分析。在 Jupyter 中,你可以使用%prun
,它会告诉你算法中时间花在哪里。然后你可以去看看实际花费时间的部分,思考它在算法上是否尽可能高效。在这种情况下,我们运行它,从 76 毫秒降到不到 2 毫秒。现在一些新手可能会认为“哦,太好了,我节省了 60 多毫秒”,但关键是这将被运行数千万次。所以 76 毫秒版本太慢了,对于实际使用的任何随机森林来说都是不切实际的。而另一方面,我找到的 1 毫秒版本实际上是相当可接受的。
%timeit find_better_split_foo(tree,1) tree ''' 2.2 ms ± 148 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) n: 1000; val:10.079014121552744; score:658.5510186055565; split:1974.0; var:YearMade '''
然后检查,数字应该与之前完全相同,而且确实如此。
find_better_split_foo(tree,0); tree ''' n: 1000; val:10.079014121552744; score:658.5510186055565; split:1974.0; var:YearMade '''
现在我们有了一个函数find_better_split
,它可以做我们想要的事情,我想把它插入到我的DecisionTree
类中。这是一个非常酷的 Python 技巧。Python 可以动态执行所有操作,因此我们实际上可以说DecisionTree
中名为find_better_split
的方法就是我刚刚创建的那个函数。
DecisionTree.find_better_split = find_better_split_foo
将其放在该类中。现在我告诉你这件事稍微令人困惑的地方是,左边的find_better_split
和右边的find_better_split
实际上并没有任何关系。它们只是恰好以相同的顺序拥有相同的字母。所以我可以将其命名为find_better_split_foo
,然后我可以调用它。现在我的函数实际上被称为find_better_split_foo
,但我期望调用的方法是名为 DecisionTree.find_better_split 的东西。所以在这里,我可以说:
DecisionTree.find_better_split = find_better_split_foo
在使用的每种语言中,了解命名空间的工作原理是很重要的。其中最重要的一点是了解它是如何确定名称所指的内容的。所以这里(DecisionTree.find_varsplit
)意味着DecisionTree
类内部定义的find_better_split
,而不是其他地方。右边的这个意味着全局命名空间中的find_better_split_foo
。许多语言没有全局命名空间,但 Python 有。因此,即使它们恰好以相同的顺序拥有相同的字母,它们也不以任何方式指向相同的内容。就像这边的家庭可能有一个叫 Jeremy 的人,而我的家庭也有一个叫 Jeremy 的人。我们的名字恰好相同,但我们并不是同一个人。
现在我们已经将find_better_split
方法放入了具有这个新定义的DecisionTree
中,当我现在调用TreeEnsemble
构造函数时,决策树集合构造函数会调用create_tree
,create_tree
实例化DecisionTree
,DecisionTree
调用find_varsplit
,它会遍历每一列以查看是否可以找到更好的分割点,我们现在已经定义了find_better_split
,因此当我们创建TreeEnsemble
时,它已经执行了这个分割点。
tree = TreeEnsemble(x_sub, y_train, 1, 1000).trees[0]; tree ''' n: 1000; val:10.079014121552744; score:658.5510186055565; split:1974.0; var:YearMade '''
好的。这很不错,对吧?我们一次只做一点点,测试每一步。当你们实现随机森林解释技术时,你们可能想尝试以这种方式编程,检查每一步是否与 scikit-learn 所做的匹配,或者与你们构建的测试匹配。
完整的单棵树[55:13]
在这一点上,我们应该尝试更深入地探究。现在让我们将 max_depth 设置为 2。这是 scikit-learn 所做的。在 YearMade 74 处中断后,它接着在 MachineHoursCurrentMeter 2956 处中断。
m = RandomForestRegressor( n_estimators=1, max_depth=2, bootstrap=False ) m.fit(x_samp, y_samp) draw_tree(m.estimators_[0], x_samp, precision=2)
所以我们有一个叫做find_varsplit
的东西,它只是遍历每一列,尝试看看是否有更好的分割点。但实际上,我们需要再进一步。我们不仅需要遍历每一列,看看这个节点是否有更好的分割点,而且还需要看看我们刚刚创建的左侧和右侧是否有更好的分割点。换句话说,左侧和右侧应该成为决策树本身。
在右侧创建这棵树和在左侧创建这棵树之间没有任何区别,除了左侧包含 159 个样本,右侧包含一千个。
因此,第一行代码与之前完全相同。然后我们检查它是否是叶节点。如果是叶节点,那么我们就没有更多的事情要做了。这意味着我们就在底部,没有进行分割,所以我们不需要做任何进一步的操作。另一方面,如果它不是叶节点,那么我们需要将其分割成左侧和右侧。现在,早些时候,我们创建了一个左侧和右侧的布尔数组。最好是有一个索引数组,因为我们不想在每个节点中都有一个完整的布尔数组。因为请记住,尽管在这个大小的树中看起来似乎没有很多节点,但当它完全展开时,底层(即如果最小叶大小为 1)包含与整个数据集相同数量的节点。因此,如果每个节点都包含整个数据集大小的完整布尔数组,那么内存需求会增加。另一方面,如果我们只存储此节点中所有内容的索引,那么它将变得越来越小。
def find_varsplit(self): for i in range(self.c): self.find_better_split(i) if self.is_leaf: return x = self.split_col lhs = np.nonzero(x<=self.split)[0] rhs = np.nonzero(x>self.split)[0] self.lhs = DecisionTree(self.x, self.y, self.idxs[lhs]) self.rhs = DecisionTree(self.x, self.y, self.idxs[rhs])
np.nonzero 与x<=self.split
完全相同,它得到布尔数组,但将其转换为true
的索引[58:07]。因此,这个lhs
现在是左侧和右侧的索引列表。现在我们有了左侧和右侧的索引,我们现在可以继续创建一个决策树。所以self.lhs
是我们左侧的决策树,self.rhs
是我们右侧的决策树。我们不需要做其他事情。我们已经写好了这些。我们已经有一个可以创建决策树的构造函数。所以当你真正思考这在做什么时,会有点让人头疼,对吧?因为find_varsplit
被调用的原因是因为决策树构造函数调用了它。但是find_varsplit
本身又调用了决策树构造函数。所以我们实际上有循环递归。我并不聪明到足以能够思考递归,所以我选择不去想。我只是写出我的意思,然后不再考虑。我想要什么?找到一个变量分割。我必须遍历每一列,看看是否有更好的东西,如果成功进行了分割,找出左侧和右侧,然后将它们转换为决策树。现在尝试思考这两种方法如何相互调用会让我发疯,但我不需要这样做。我知道我有一个有效的决策树构造函数,我知道我有一个有效的find_varsplit
,所以就这样。这就是我进行递归编程的方式,就是假装我没有。我只是忽略它。这是我的建议。你们中很多人可能足够聪明,能够比我更好地思考这个问题,那就好。如果你能的话。
DecisionTree.find_varsplit = find_varsplit
所以现在我已经写好了,我可以将其打补丁到 DecisionTree 类中,一旦我这样做了,TreeEnsemble 构造函数将会使用它,因为 Python 是动态的。
tree = TreeEnsemble(x_sub, y_train, 1, 1000).trees[0]; tree ''' n: 1000; val:10.079014121552744; score:658.5510186055565; split:1974.0; var:YearMade '''
现在我可以检查1:00:31。我的左手边应该有 159 个样本和值为 9.66。
tree.lhs ''' n: 159; val:9.660892662981706; score:76.82696888346362; split:2800.0; var:MachineHoursCurrentMeter '''
右手边,841 个样本和 10.15。
tree.rhs ''' n: 841; val:10.158064432982941; score:571.4803525045031; split:2005.0; var:YearMade '''
左手边的左手边,150 个样本和 9.62。
tree.lhs.lhs ''' n: 150; val:9.619280538108496; score:71.15906938383463; split:1000.0; var:YearMade '''
所以你可以看到,因为我并不聪明到足以编写机器学习算法,不仅我第一次写不正确,通常每一行我写的都是错误的。所以我总是从这样的假设开始,我刚刚输入的代码几乎肯定是错误的。我只需要看看为什么以及如何。所以我只是确保。最终我会到达这样一个点,让我很惊讶的是,它不再出错了。所以在这里,我可以感觉到好吧,如果所有这些事情碰巧与 scikit-learn 完全相同,那将是令人惊讶的。所以看起来还不错。
tree.lhs.rhs ''' n: 9; val:10.354428077535193 '''
预测 [1:01:43]
现在我们有了一个可以构建整个树的东西,我们想要有一个可以计算预测的东西。所以提醒一下,我们已经有了一个可以为TreeEnsemble
计算预测的东西(通过调用tree.predict(x)
),但在DecisionTree
中没有叫做tree.predict
的东西,所以我们需要写一个。
为了让这更有趣,让我们开始增加我们使用的列数。
cols = [ 'MachineID', 'YearMade', 'MachineHoursCurrentMeter', 'ProductSize', 'Enclosure','Coupler_System', 'saleYear' ]
让我们再次创建我们的TreeEnsemble
。
%time tree = TreeEnsemble(X_train[cols], y_train, 1, 1000).trees[0] x_samp,y_samp = tree.x, tree.y ''' CPU times: user 288 ms, sys: 12 ms, total: 300 ms Wall time: 297 ms '''
这一次,让我们将最大深度设为 3。
m = RandomForestRegressor( n_estimators=1, max_depth=3, bootstrap=False ) m.fit(x_samp, y_samp) draw_tree(m.estimators_[0], x_samp, precision=2, ratio=0.9, size=7)
所以现在我们的树变得更加有趣。现在让我们定义如何为树创建一组预测。因此,树的一组预测就是每一行的预测。就是这样。这就是我们的预测。因此,树的预测是数组中每一行的预测。所以再次,我们跳过思考,思考很困难。所以让我们继续推迟。这个**for** xi **in** x
很方便,对吧?请注意,无论数组的秩如何,您都可以使用 numpy 数组中的for
blah。无论数组中的轴数是多少。它的作用是遍历主轴。
这些概念在我们进入越来越多的神经网络时将变得非常重要,因为我们将一直在进行张量计算。因此,向量的主轴是向量本身。矩阵的主轴是行。三维张量的主轴是表示切片的矩阵等等。在这种情况下,因为 x 是一个矩阵,这将循环遍历行。如果您以这种方式编写您的张量代码,那么它将很好地推广到更高的维度。在这个 x
中有多少维度并不重要。这将循环遍历每个主轴。因此,我们现在可以称之为 DecisionTree.predict
。
def predict(self, x): return np.array([self.predict_row(xi) for xi in x])
所以我需要做的就是编写predict_row
。我一直在拖延思考,这很好,实际上我需要做工作的地方,现在基本上是微不足道的。如果我们在叶节点,那么预测值就等于我们在原始树构造函数中计算的那个值(即y
的平均值)。如果不是叶节点,那么我们必须弄清楚是沿左路径还是右路径进行预测。因此,如果这一行中的变量(xi[self.var_idx]
)小于或等于我们决定拆分的值,则我们沿左路径前进;否则我们沿右路径前进。然后,确定我们想要的路径/树之后,我们只需在其上调用predict_row
。再次,我们无意中创建了递归的东西。如果是叶节点,则返回该值;否则根据需要返回左侧或右侧的预测值。
def predict_row(self, xi): if self.is_leaf: return self.val t = self.lhs if xi[self.var_idx]<=self.split else self.rhs return t.predict_row(xi) DecisionTree.predict_row = predict_row
注意这里的self.lhs **if** xi[self.var_idx]<=self.split **else** self.rhs
,这个 if 与上面的 if 没有任何关系:
if something: x= do1() else: x= do2()
上面的这个 if 是一个控制流语句,告诉 Python 沿着这条路径或那条路径进行一些计算。下面的这个 if 是一个返回值的运算符。
x = do1() if something else do2()
所以你们做过 C 或 C++的人会认出它与这个是完全相同的(即三元运算符):
x = something ? do1() : do2()
基本上我们要做的是,我们要得到一个值,如果something
为真,我们会说这个值是(do1()
),否则是另一个值(do2()
)。你可以用冗长的方式来写,但那将需要写 4 行代码来做一件事,而且还需要你编写的代码,如果你自己或向别人阅读时,表达方式并不自然。我想说“我要走的树是左边,如果变量小于分割值,否则是右边。所以我想按照我思考或说代码的方式来编写我的代码。因此,这种三元运算符对此非常有帮助。
所以现在我已经对一行进行了预测,我可以将其放入我的类中:
DecisionTree.predict = predict
现在我可以计算预测。
%time preds = tree.predict(X_valid[cols].values) ''' CPU times: user 156 ms, sys: 4 ms, total: 160 ms Wall time: 162 ms '''
现在我可以将我的实际数据与我的预测数据进行对比。当你做散点图时,通常会有很多点重叠在一起,所以一个好的技巧是使用 alpha。Alpha 表示透明度,不仅在 matplotlib 中,在世界上几乎所有的图形包中都是如此。因此,如果将 alpha 设置为小于 1,那么这意味着你需要将 20 个点叠加在一起才能完全显示为蓝色。这是一个很好的方法来看看有多少点重叠在一起 - 散点图的一个好技巧。
plt.scatter(preds, y_valid, alpha=0.05)
这是我的 R²。
metrics.r2_score(preds, y_valid) ''' 0.50371522136882341 '''
那么现在让我们继续进行一个没有最大分裂次数的随机森林,我们的树集合也没有最大分裂次数,我们可以将我们的 R²与他们的 R²进行比较。
m = RandomForestRegressor( n_estimators=1, min_samples_leaf=5, bootstrap=False ) %time m.fit(x_samp, y_samp) preds = m.predict(X_valid[cols].values) plt.scatter(preds, y_valid, alpha=0.05)
metrics.r2_score(preds, y_valid) ''' 0.47541053100694797 '''
它们并不相同,但实际上我们的稍微好一点。我不知道我们做了什么不同,但我们会接受它😊 所以现在我们有了一个对于一个只有一棵树的森林,在使用一个真实的实际数据集(推土机的蓝皮书)进行验证时,与 scikit-learn 相比提供了同样好的准确性。
把它放在一起
让我们继续完善这个。现在我想要做的是创建一个包含这段代码的包。我通过创建一个方法,再创建一个方法,然后将它们拼接在一起来创建这个包。现在我回到笔记本中,收集了所有实现方法的单元格,然后将它们全部粘贴在一起。
class TreeEnsemble(): def __init__(self, x, y, n_trees, sample_sz, min_leaf=5): np.random.seed(42) self.x,self.y,self.sample_sz,self.min_leaf = x,y,sample_sz,min_leaf self.trees = [self.create_tree() for i in range(n_trees)] def create_tree(self): idxs = np.random.permutation(len(self.y))[:self.sample_sz] return DecisionTree( self.x.iloc[idxs], self.y[idxs], idxs=np.array(range(self.sample_sz)), min_leaf=self.min_leaf ) def predict(self, x): return np.mean([t.predict(x) for t in self.trees], axis=0)def std_agg(cnt, s1, s2): return math.sqrt((s2/cnt) - (s1/cnt)**2)class DecisionTree(): def __init__(self, x, y, idxs, min_leaf=5): self.x,self.y,self.idxs,self.min_leaf = x,y,idxs,min_leaf self.n,self.c = len(idxs), x.shape[1] self.val = np.mean(y[idxs]) self.score = float('inf') self.find_varsplit() def find_varsplit(self): for i in range(self.c): self.find_better_split(i) if self.score == float('inf'): return x = self.split_col lhs = np.nonzero(x<=self.split)[0] rhs = np.nonzero(x>self.split)[0] self.lhs = DecisionTree(self.x, self.y, self.idxs[lhs]) self.rhs = DecisionTree(self.x, self.y, self.idxs[rhs]) def find_better_split(self, var_idx): x,y = self.x.values[self.idxs,var_idx], self.y[self.idxs] sort_idx = np.argsort(x) sort_y,sort_x = y[sort_idx], x[sort_idx] rhs_cnt,rhs_sum,rhs_sum2 = self.n,sort_y.sum(),(sort_y**2).sum() lhs_cnt,lhs_sum,lhs_sum2 = 0,0.,0. for i in range(0,self.n-self.min_leaf-1): xi,yi = sort_x[i],sort_y[i] lhs_cnt += 1; rhs_cnt -= 1 lhs_sum += yi; rhs_sum -= yi lhs_sum2 += yi**2; rhs_sum2 -= yi**2 if i<self.min_leaf or xi==sort_x[i+1]: continue lhs_std = std_agg(lhs_cnt, lhs_sum, lhs_sum2) rhs_std = std_agg(rhs_cnt, rhs_sum, rhs_sum2) curr_score = lhs_std*lhs_cnt + rhs_std*rhs_cnt if curr_score<self.score: self.var_idx,self.score,self.split = var_idx,curr_score,xi @property def split_name(self): return self.x.columns[self.var_idx] @property def split_col(self): return self.x.values[self.idxs,self.var_idx] @property def is_leaf(self): return self.score == float('inf') def __repr__(self): s = f'n: **{self.n}**; val:**{self.val}**' if not self.is_leaf: s += f'; score:**{self.score}**; split:**{self.split}**; var: **{self.split_name}**' return s def predict(self, x): return np.array([self.predict_row(xi) for xi in x]) def predict_row(self, xi): if self.is_leaf: return self.val t = self.lhs if xi[self.var_idx]<=self.split else self.rhs return t.predict_row(xi)
就是这样。这就是我们一起编写的代码。
ens = TreeEnsemble(X_train[cols], y_train, 5, 1000) preds = ens.predict(X_valid[cols].values) plt.scatter(y_valid, preds, alpha=0.1, s=6);
metrics.r2_score(y_valid, preds) ''' 0.71011741571071241 '''
这里我们有一个蓝色推土机的模型,使用了我们完全从头开始编写的随机森林,R²为 71。这很酷。
性能和 Cython
当我尝试比较这个与 scikit-learn 的性能时,这个要慢得多,原因是虽然很多工作是由 numpy 完成的,numpy 是优化良好的 C 代码,但想想树的最底层。如果我们有一百万个数据点,树的底层有大约 500,000 个决策点,底下有一百万个叶子。这就像调用了 500,000 个分割方法,其中包含多次调用 numpy,而 numpy 只有一个要计算的项目。这是非常低效的。这是 Python 在性能方面特别不擅长的事情(即多次调用大量函数)。我们可以看到它并不差。对于 15 年前被认为是相当大的随机森林来说,这被认为是相当不错的性能。但是现在,这至少比应该的速度慢了几百倍。
scikit-learn 的开发人员为了避免这个问题所做的是,他们使用了一种叫做 Cython 的东西来实现。Cython 是 Python 的一个超集。所以你写的任何 Python 代码基本上都可以作为 Cython 来使用。但是 Cython 运行方式有所不同。它不是直接传递给 Python 解释器,而是将其转换为 C 语言,编译,然后运行该 C 代码。这意味着,第一次运行时会花费一些时间,因为需要进行翻译和编译,但之后运行会快得多。所以我想快速向你展示一下这是什么样子,因为你肯定会遇到 Cython 可以帮助你工作的情况,而你大部分一起工作的人可能从未使用过它(甚至可能不知道它的存在),所以拥有这种超能力是非常棒的。
在笔记本中使用 Cython,你可以这样说:
%load_ext Cython
这里是一个 Python 函数fib1
:
def fib1(n): a, b = 0, 1 while b < n: a, b = b, a + b
这里是一个 Cython 函数。它与顶部的%%cython
完全相同。实际上,它的运行速度大约是fib1
的两倍,因为它进行了编译。
%%cython def fib2(n): a, b = 0, 1 while b < n: a, b = b, a + b
这里是同样的版本,我使用了一个特殊的 Cython 扩展叫做cdef
,它定义了返回值和每个变量的 C 数据类型。基本上这就是你可以用来开始加快运行速度的技巧。在那一点上,现在它知道它不只是一个名为 T 的 Python 对象。所以 fib3,它和之前完全一样,但我们说我们传递给它的东西的数据类型是什么,然后定义每个变量的数据类型。
%%cython def fib3(int n): cdef int b = 1 cdef int a = 0 cdef int t = 0 while b < n: t = a a = b b = a + b
所以如果我们这样做,现在我们有了一个快 10 倍的东西。
%timeit fib1(50) ''' 705 ns ± 62.5 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) ''' %timeit fib2(50) ''' 362 ns ± 26.5 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) ''' %timeit fib3(50) ''' 70.7 ns ± 4.07 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each) '''
这并不需要太多额外的工作,只是用一点 Python 和一些标记,所以知道它的存在是很好的,因为如果有一些定制的东西你想要做,实际上要去 C 语言编译并链接回来是很痛苦的。而在这里做起来相当容易。
问题:当你在使用 Cython 版本时,对于 numpy 数组,是否有特定的 C 类型[1:17:16]?是的,有很多特定的内容用于将 Cython 与 numpy 集成,有一个完整的页面介绍了这些内容。所以我们不用担心过多细节,但你可以阅读那个页面,基本上可以了解基本思想。
使用 NumPy 与 Cython — Cython 0.29a0 文档
有这个cimport
,基本上是将某种类型的 Python 库导入到代码的 C 部分,然后你可以在 Cython 中使用它。这很简单直接。
所以你现在的任务是实现:
- 基于树方差的置信度
- 特征重要性
- 部分依赖
- 树解释器
对于那个随机森林。删除冗余特征根本不使用随机森林,所以你不必担心这个。外推不是一种解释技术,所以你也不必担心这个。所以只是其他的。所以基于树方差的置信度,我们已经编写了那段代码,所以我怀疑我们在笔记本中拥有的完全相同的代码应该继续工作。所以你可以尝试确保让它工作。特征重要性是通过变量洗牌技术实现的,一旦你让它工作,偏依赖只是几行代码之遥,因为你不是洗牌一列,而是用一个常数值替换它。几乎是相同的代码。
然后树解释器,这将需要您编写一些代码并思考。一旦您编写了树解释器,如果您愿意,您就非常接近创建特征重要性的第二种方法 - 即在所有行中累加重要性的方法。这意味着,您将非常接近进行交互重要性。事实证明,xgboost 实际上有一个非常好的交互重要性库,但似乎没有一个适用于随机森林,因此您可以从使其在我们的版本上运行开始(如果您想进行交互重要性),然后您可以使其在原始的 scikit-learn 版本上运行,这将是一个很酷的贡献。有时,针对自己的实现编写代码更好,因为您可以清楚地看到发生了什么。
如果在任何时候遇到困难,请在论坛上提问。关于如何寻求帮助的整个页面都在维基上有。当你在 Slack 上向同事寻求帮助,当你在 Github 或 Discourse 上向技术社区的人寻求帮助时,正确地寻求帮助将有助于让人们愿意帮助你并能够帮助你。
- 搜索你遇到的错误,看看是否已经有人问过。
- 你已经尝试过如何修复它了吗?
- 你认为出了什么问题?
- 你用的是什么样的电脑?它是如何设置的?软件版本是什么?
- 你确切地输入了什么,确切地发生了什么?
你可以通过截屏来做到这一点,所以确保你有一些非常容易使用的截屏软件。所以如果我要截屏,我只需按下一个按钮,选择区域,复制到剪贴板,转到论坛,粘贴进去,然后就完成了(你甚至可以将图像缩小!)。
更好的做法是,如果有几行代码和错误消息需要查看,可以创建一个 Gist。Gist 是一个很方便的 Github 工具,基本上可以让你分享代码。如果我想要创建一个这样的 Gist,我实际上有一个扩展:
点击那个,给它起个名字,然后点击“公开”。这样就可以将我的 Jupyter 笔记本公开分享。然后我可以复制那个 URL,复制链接位置,然后粘贴到我的论坛帖子中。然后当人们点击它时,他们会立即看到我的笔记本。
现在,这个特定的按钮是一个扩展,所以在 Jupyter 上,您需要点击 Nbextensions,然后点击 Gist-it。当您在那里时,您还应该点击 Collapsible Headings,这是我使用的一个方便的功能,让我可以折叠和展开内容。
如果您打开 Jupyter 时没有看到这个 Nbextensions 按钮,那么只需搜索 Jupyter Nbextensions — 它会告诉您如何使用 pip 安装并设置它。
神经网络广义定义 [1:23:20]
除了作业之外,我们已经完成了随机森林,直到下一个课程,当你看到 GBMs 时,我们已经完成了决策树集成。我们将转向广义的神经网络。神经网络将使我们能够超越随机森林的最近邻方法。所有随机森林能做的就是对已经看到的数据进行平均。它不能外推或计算。线性回归可以计算和外推,但只能以非常有限的方式。神经网络给我们带来了两全其美的好处。
我们将从将它们应用于非结构化数据开始。非结构化数据指的是像素、声波振幅或单词 - 数据中所有列中的所有内容都是相同类型,而不是数据库表中有收入、成本、邮政编码和州名(结构化数据)。我们也将用它来处理结构化数据,但稍后再做。非结构化数据稍微容易一些,也是更多人长期以来一直在应用深度学习的领域。
如果您也在学习深度学习课程,您会发现我们将从两个不同的方向接近相同的结论。因此,深度学习课程从解决复杂的卷积神经网络开始,使用复杂的优化方案,我们将逐渐深入了解它们的工作原理。而机器学习课程则更多地从随机梯度下降的实际工作原理开始,我们可以用单层来创建逻辑回归等内容。当我们添加正则化时,它如何给我们提供岭回归、弹性网络套索等内容。当我们添加更多层时,它如何让我们处理更复杂的问题。在这个机器学习课程中,我们只会看到全连接层,我认为下个学期与 Yannet 一起,您可能会看到一些更复杂的方法。因此,这个机器学习课程中,我们将更多地关注矩阵的实际运算过程,而深度学习则更多地关注如何以世界级水平解决真实世界的深度学习问题的最佳实践。
下周,我们将研究经典的 MNIST 问题,即如何识别数字。如果你感兴趣,你可以提前尝试使用随机森林来解决这个问题,你会发现效果不错。考虑到随机森林基本上是一种最近邻的类型(它在树空间中找到最近的邻居),那么随机森林绝对可以识别出这个 9,这些像素与我们在其他图像中看到的像素相似,而且平均来说,它们也是 9。因此,它绝对可以使用随机森林解决这类问题。但我们最终会受到数据限制,因为每次我们增加一个决策点,我们的数据大致减半,所以这就限制了我们可以进行的计算量。而神经网络,我们将能够使用大量的参数,通过我们将学习的正则化技巧,我们将能够进行大量的计算,实际上我们几乎没有什么限制可以计算的结果。
祝你在随机森林解释方面好运,下次再见。