背景

  最近在用fasttext进行文本分类，说到文本分类就不得不提起word2vec。Word2Vec是由Google的Mikolov等人提出的一个词向量计算模型，其重要意义在于将自然语言转换成了计算机能够理解的向量，能够更好地计算词与词之间的相似性。

One-Hot representation

  一般比较常见的表示词向量的方式就是独热码表示(One-Hot representation)，它首先将文章中不同的单词提取出来作为一个词集，然后把词集的大小作为词向量的维度，对于每个具体的词将对应的位置置1。例如，我们又下面5个词汇组成的词集[King,Queue,MAN,WOMAN,CHILD]，则词Queue的One-Hot representation表达形式为[0,1,0,0,0]，同理WOMAN的表达形式为[0,0,0,1,0]。

  One-Hot representation用来表示词向量简单粗暴，但是有很多问题：

1. 向量维度和词集大小一致，能够达到百万级，不利于内存计算
2. 两个不一样的词做内积永远为0，无法表达词之间的相似关系

Distributed representation

  Distributed representation被翻译为稠密表示(对于One-Hot representation只有在词汇对应位置上才置1，所以它是稀疏的，而Distributed representation不需要那么多维度来表示所有词汇，它可能在每个维度上都有值，所以在表现形式上更稠密），它通过训练的方式把每个词映射到一个较短的词向量上来，所有的词向量构成了向量空间，进而可以通过统计学方法来研究词与词之间的关系。

原理

  word2vec接收分词后的文本作为输入，通过神经网络学习输入的每个词汇前后N个单词可能出现的词（这里应该是预测前面所有和后面所有的单词，但是简化为N-gram模型，相关原理详见文章最后的参考链接），最后产生一个稠密向量表示每个词。因为网络能够学习到一个词前后可能会出现什么词，所以通过word2vec产生的词向量可以计算词之间的相似性。

网络结构

  word2vec是一个简单的三层神经网络，即一个输入层、一个隐藏层、一个输出层。输入层就是词的One-Hot representation，隐层的神经网络单元数量就是embedding size，即最终词向量的维度，隐藏层之后不需要使用激活函数，直接接到输出层。输出层是softmax(这里实际是Hierarchical Softmax，为了简化直接使用softmax)，得到每个预测结果的概率。其网络结构示意图如下所示。

  假设词集大小为10000，词向量的维度为300，即embedding_size=300

输入层

  一个词的one-hot表达形式，长度为10000，例如[0,0,1,0…,0,0]

隐藏层

  隐藏层的神经元数量就是词向量的长度，隐层参数是一个[10000,300]的矩阵。实际上，这个矩阵就是最后的词向量。其实输入和隐藏层的作用就是把词的one-hot表达形式转化为词向量的表达形式，下面的图更清楚一些

输出层

  输出层就是一个10000维的向量，每个值代表一个词的概率。

skip-gram

  skip-gram核心思想是通过中心词来预测周围词。假设中心词是cat，窗口长度为2，则根据cat预测左边两个词和右边两个词。例如，对于文本“the quick brown fox jumps over the lazy dog”，窗口长度为2时，有

1. 当中心词为the时，有(the,quick),(the,brown）
2. 当中心词为quick时，有(quick,the),(quick,brown),(quick,fox）
3. 当中心词为brown时，有(brown,the),(brown,quick),(brown,fox),(brown,jumps)
4. 当中心词为fox时，有(fox,quick),(fox,brown),(fox,jumps),(fox,over)

上面的窗口移动了4次，产生了13个样本。

CBOW

  CBOW(continuous bag of words)的核心思想与skip-gram相反，它通过中心词周围的词来预测中心词。例如,对于文本“the quick brown fox jumps over the lazy dog”，窗口长度为2时，有

1. 当中心词为the时，有([quick,brown],the）
2. 当中心词为quick时，有([the,brown,fox],quick）
3. 当中心词为brown时，有([the,quick,fox,jumps],brown)
4. 当中心词为fox时，有([quick,brown,jumps,over],fox)

上面的窗口移动了4次，但是只产生了4个样本。这时候input是4个词，label是一个词，经过隐藏层之后输入的4个词被映射成了4个EmbedingSize维度的向量，对4个向量求平均后才能作为下一层的输入。

   两个模型相比，skip-gram模型能产生更多的训练样本，抓住更多词与词之间语义上的细节，在语料多足够好的理想条件下，skip-gram模型优于CBOW模型。在语料较少的情况下，难以抓住足够词与词之间之间的细节，CBOW模型求平均的特性反而效果可能更好。

代码实现

  tensorflow官网的示例中给出了word2vec的代码实现，但是该实现太繁琐，现给出简单版实现，详见word2vec。

网络结构

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

# coding=utf-8

import os
import tensorflow as tf
from data import SkipGramDataSet
import numpy as np

dataset = SkipGramDataSet(os.path.join(os.path.curdir, 'test.txt'))

VOCAB_SIZE = dataset.vocab_size
print 'vocab_size:{}'.format(VOCAB_SIZE)
EMBEDDING_SIZE = 128
LEARNING_RATE = 0.01

TRAIN_STEPS = 10000

BATCH_SIZE = 32
WINDOW_SIZE = 2


class Word2Vec(object):

    def __init__(self):
        self.graph = tf.Graph()
        with self.graph.as_default():
            # 输入层，维度为词集大小
            with tf.name_scope('inputs'):
                self.x = tf.placeholder(shape=(None, VOCAB_SIZE), dtype=tf.float32)
                self.y = tf.placeholder(shape=(None, VOCAB_SIZE), dtype=tf.float32)
            # 隐藏层，w1就是词向量
            with tf.name_scope('layer1'):
                self.W1 = tf.Variable(tf.random_uniform([VOCAB_SIZE, EMBEDDING_SIZE], -1, 1), dtype=tf.float32,
                                      name='w1')
                self.b1 = tf.Variable(tf.random_normal([EMBEDDING_SIZE], dtype=tf.float32))
            # hidden是把输入的one-hot转化为词向量的结果
            hidden = tf.add(self.b1, tf.matmul(self.x, self.W1))
            with tf.name_scope('layer2'):
                self.W2 = tf.Variable(tf.random_uniform([EMBEDDING_SIZE, VOCAB_SIZE], -1, 1), dtype=tf.float32)
                self.b2 = tf.Variable(tf.random_normal([VOCAB_SIZE]), dtype=tf.float32)
            # 输出层是softmax求概率之后的结果
            self.prediction = tf.nn.softmax(tf.add(tf.matmul(hidden, self.W2), self.b2))
            # 损失函数是交叉熵
            log = self.y * tf.log(self.prediction)
            self.loss = tf.reduce_mean(-tf.reduce_sum(log, reduction_indices=[1], keep_dims=True))
            # 梯度下降
            self.opt = tf.train.GradientDescentOptimizer(LEARNING_RATE).minimize(self.loss)

    # 把词的下标值转化为one-hot表示
    def _one_hot_input(self, dataset):
        # features和labels记录了词在词集中的位置
        features, labels = dataset.generate_batch_inputs(BATCH_SIZE, WINDOW_SIZE)
        f, l = [], []
        for w in features:
            # 产生全0向量
            tmp = np.zeros([VOCAB_SIZE])
            # 下标位置置1
            tmp[w] = 1
            f.append(tmp)
        for w in labels:
            tmp = np.zeros([VOCAB_SIZE])
            tmp[w] = 1
            l.append(tmp)
        return f, l

    def train(self, dataset, n_iters, ):
        with tf.Session(graph=self.graph) as sess:
            sess.run(tf.global_variables_initializer())
            for i in range(n_iters):
                features, labels = self._one_hot_input(dataset)
                predi, loss, w1 = sess.run([self.prediction, self.loss],
                                           feed_dict={self.x: features, self.y: labels})
                print 'loss:{}'.format(loss)

skip-gram数据准备

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

class DataSet(object):
    def __init__(self, file):
        self.file = file
        self.data_index = 0
        self._build_dataset()

    def _build_dataset(self):
        if not os.path.exists(self.file):
            raise ValueError("file doesn't exists --> {}".format(self.file))
        f = open(self.file, 'r')
        # 保存词集
        self.data = tf.compat.as_str(f.read()).split()
        if f:
            f.close()
        c = collections.Counter(self.data).most_common()
        # 计算词集大小
        self.vocab_size = len(c)
        self.counter = c.insert(0, ('UNK', -1))
        self.vocab_size += 1
        # 词-下标字典
        self.word2id = dict()
        # 下标-词字典
        self.id2word = dict()
        for word, _ in c:
            self.word2id[word] = len(self.word2id)
            self.id2word[len(self.id2word)] = word

    def generate_batch_inputs(self, batch_size, window_size):
        raise NotImplementedError()


class SkipGramDataSet(DataSet):
    def generate_batch_inputs(self, batch_size, window_size):
        features = np.ndarray(shape=(batch_size,), dtype=np.int32)
        labels = np.ndarray(shape=(batch_size,), dtype=np.int32)
        i = 0
        while True:
            if self.data_index == len(self.data):
                self.data_index = 0
            # 窗口的左侧位置
            left = max(0, self.data_index - window_size)
            # 窗口的右侧位置
            right = min(len(self.data), self.data_index + window_size + 1)
            # 遍历窗口里的每个单词
            for k in range(left, right):
                if k != self.data_index:
                    # 输入是中心词
                    features[i] = self.word2id[self.data[self.data_index]]
                    # label值是中心词周围在窗口内的值
                    labels[i] = self.word2id[self.data[k]]
                    i += 1
                    if i == batch_size:
                        return features, labels
            self.data_index += 1

参考

1. Word2vec原理浅析及tensorflow实现
2. 自己动手实现word2vec（Skip-gram模型）

SoftMax简介

  SoftMax回归是Logistic回归模型在多分类问题上的推广，在多分类问题中，类标签$y$可以取两个以上的值。本人之前对softmax回归的名字由来一直有一个疑问，为什么叫softmax回归而不叫其他回归？

  softmax回归有人翻译成柔性回归。假设有这样一种场景，在进行minist数据集分类时，输出端10个神经元其中三个值是0、1、2，为了得到一个概率分布，我们可能会用0/3，1/3，2/3归一化，但是这样不符合要求。因为这样取该值概率为0，永远无法取得，但实际应该有一定概率可以取到改值。同样2比1大，我们就说有2/3的概率取得2对应的值，这样也不对。用softmax，加入指数后即使输出为0，也可以求得一个小的非0值，虽然很小但是仍有可能取得该值。对于2来说，加入指数后呈指数行增长，比1增长的快，所以使得2的概率比2/3大。

公式推导

  假如有$m$个训练样本${(x^{(1)},y{(1)}),(x^{(2)},y{(2)}),…,(x^{(m)},y{(m)})}$，输入特征$x^{(i)}\in \mathcal{R}^{n+1}$，类标记为$y_i \in {0,1,…,k}$。假设函数为对于每一个样本估计所属类别的概率$p(y=j|x)$，即

$$
h_{\theta}(x^{(i)})=\left[
\begin{split}
{p(y^{(i)}=1|x^{i};\theta)}& \\
{p(y^{(i)}=2|x^{i};\theta)}& \\
{\vdots}& \\
{p(y^{(i)}=k|x^{i};\theta)} &\\
\end{split}
\right]=\frac{1}{\sum_{j=1}^{k} e^{\theta^{T}_{j}x^{(i)}}}\left[\begin{split}
e^{\theta^{T}_{1}x^{(i)}}& \\
e^{\theta^{T}_{2}x^{(i)}}& \\
{\vdots}& \\
e^{\theta^{T}_{k}x^{(i)}} &\\
\end{split} \right]
$$

其中$\theta$表示向量，且$\theta_{i} \in \mathcal{R}^{n+1}$。对于每一个样本估计其所属类别的概率为
$$
p(y^{(i)}=j|x^{(i)};\theta)=\frac{e^{\theta_{j}^{T}x^{(i)}}}{\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}}}
$$

代价函数

  引入指示函数$I$，表示样本$i$是否属于第$j$类，所以对于softmax回归的代价函数为

$$
J(\theta)=-\frac{1}{m}[\sum_{i=1}^{m}\sum_{j=1}^{k}I({y^{(i)}}=j)log\frac{e^{\theta_{j}^{T}x^{(i)}}}{\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}}}]
$$

求解

  对于上述代价函数，使用梯度下降算法对其进行求解，首先对其进行求梯度

$$
\frac{\partial J(\theta)}{\partial \theta_j}=\frac{ \partial \lbrace-\frac{1}{m}[\sum_{i=1}^{m}\sum_{j=1}^{k}I({y^{(i)}}=j)log\frac{e^{\theta_{j}^{T}x^{(i)}}}{\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}}}] \rbrace}{\partial \theta_j}
$$

对于一个样本$i$只能属于一个类别$j$，

• 若$y^{(i)}=j$，则$I({y^{(i)}}=j)=1$
  $$
  \begin{split}
  \frac{\partial J(\theta)}{\partial \theta_j} &=-\frac{1}{m}\sum_{i=1}^{m}\frac{\partial log\frac{e^{\theta_{j}^{T}x^{(i)}}}{\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}}}}{\partial \theta_j} \\
  &=-\frac{1}{m}\sum_{i=1}^{m}[\frac{\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}}}{e^{\theta_{j}^{T}x^{(i)}}} \frac{e^{\theta^{T}_{l}x^{(i)}} \cdot x^{(i)} \cdot \sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}} - e^{\theta^{T}_{j}x^{(i)}} \cdot x^{(i)} \cdot e^{\theta^{T}_{j}x^{(i)}}}{(\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}})^2}]\\
  \end{split}
  $$

• 若$y^{(i)} \neq j$，$y^{(i)} \neq j^{‘}$，则$I({y^{(i)}}=j)=0$，$I({y^{(i)}}=j^{‘})=1$

$$
\begin{split}
\frac{\partial J(\theta)}{\partial \theta_j} &= -\frac{1}{m}\sum_{i=1}^{m}\frac{ \partial log\frac{e^{\theta_{j^{‘}}^{T}x^{(i)}}}{\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}}}}{\partial \theta_j} \\
&= -\frac{1}{m}\sum_{i=1}^{m} [\frac{\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}}}{e^{\theta_{j^{‘}}^{T}x^{(i)}}} \frac{-e^{\theta^{T}_{j^{‘}}x^{(i)}} \cdot x^{(i)} \cdot e^{\theta^{T}_{j}x^{(i)}}}{(\sum_{l=1}^{k} e^{\theta^{T}_{l}x^{(i)}})^2} ]\
&=-\frac{1}{m}\sum_{i=1}^{m} [-\frac{e^{\theta^{T}_{j}x^{(i)}}}{\sum_{l=1}^{k}e^{\theta^{T}_{l}x^{(i)}}} \cdot x^{(i)}]
\end{split}
$$

综上有

$$
\frac{\partial J(\theta)}{\partial \theta_j} = -\frac{1}{m}\sum_{i=1}^{m}[x^{(i)} \cdot (I \lbrace y^{(i)} =j \rbrace - p(y^{(i)}=j|x^{(i)};\theta))]
$$

接下来可以对$\theta_j$使用梯度下降。

与logistic的关系

softmax回归的参数特点

  softmax回归中存在参数冗余的问题，简单来讲就是参数中有些参数是没有用的，为了证明这点，假设从参数$\theta_j$中减去向量$\psi$，假设函数为
$$
\begin{split}
p(y^{(i)}=j|x^{(i)};\theta)&=\frac{e^{(\theta_{j}-\psi)^{T}x^{(i)}}}{\sum_{j=1}^{k} e^{(\theta_l-\psi)^{T}x^{(i)}}} \\
&= \frac{e^{\theta_j^T \cdot x^{(i)}} \cdot e^{-\psi^T \cdot x^{(i)}}}{\sum_{j=1}^{k} e^{\theta_l^{T} \cdot x^{(i)}}\cdot e^{-\psi^{T} \cdot x^{(i)}}}\\
&=\frac{e^{\theta_{j}^{T}x^{(i)}}}{\sum_{j=1}^{k} e^{\theta_l^{T}x^{(i)}}}
\end{split}
$$

从softmax推导出logtistic

  logistic算法是softmax回归的特殊情况，即$k=2$时，此时softmax回归有

$$
h_{\theta}^{x} = \frac{1}{e^{\theta_{1}^{T} \cdot x} + e^{\theta_{2}^{T} \cdot x}}\left[
\begin{aligned}
e^{\theta_{1}^{T} \cdot x} \
e^{\theta_{2}^{T} \cdot x}
\end{aligned}
\right]
$$

利用softmax回归参数冗余的特点，另$\psi=\theta_1$，从两个向量中减去该向量，得到

$$
\begin{split}
h_{\theta}^{x} &= \frac{1}{e^{(\theta_{1}-\psi)^{T} \cdot x} + e^{(\theta_{2}-\psi)^{T} \cdot x}}\left[
\begin{aligned}
e^{(\theta_{1}-\psi)^{T} \cdot x} \
e^{(\theta_{2}-\psi)^{T} \cdot x}
\end{aligned}
\right] \\
&=\left[
\begin{aligned}
\frac{1}{1+e^{(\theta_2-\theta_1)^T \cdot x}} \
\frac{e^{(\theta_2-\theta_1)^T \cdot x}}{1+e^{(\theta_2-\theta_1)^T \cdot x}}
\end{aligned}
\right]
\end{split}
$$
上述表达式和logistic是一致的

tensorflow实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

# coding=utf-8

import tensorflow as tf

from tensorflow.examples.tutorials.mnist import input_data

# 数据集
mnist = input_data.read_data_sets('/tmp/data/', one_hot=True)
# 超参数
learning_rate = 0.01
training_epochs = 10
batch_size = 100
display_step = 1
# 输入数据
x = tf.placeholder(tf.float32, [None, 784])
y = tf.placeholder(tf.float32, [None, 10])
# 参数
W = tf.Variable(tf.zeros([784, 10]))
b = tf.Variable(tf.zeros([10]))
# softmax函数
pred = tf.nn.softmax(tf.matmul(x, W))
# 损失函数，这里因为y用one_hot表示，所以可以直接用矩阵代替指示函数I
cost = tf.reduce_mean(-tf.reduce_sum(y * tf.log(pred), reduction_indices=1))
# W梯度下降
W_grad = -tf.matmul(tf.transpose(x), y - pred)
# b梯度下降
b_grad = -tf.reduce_mean(-tf.matmul(tf.transpose(x), y - pred), reduction_indices=0)
# W更新方式
new_W = W.assign(W - learning_rate * W_grad)
# b更新方式
new_b = b.assign(b - learning_rate * b_grad)

init = tf.global_variables_initializer()

with tf.Session() as sess:
    init.run()
    for epoch in range(training_epochs):
        avg_cost = 0
        total_batch = int(mnist.train.num_examples / batch_size)
        for i in range(total_batch):
            batch_xs, batch_ys = mnist.train.next_batch(batch_size)
            _, _, c = sess.run([new_W, new_b, cost], feed_dict={x: batch_xs, y: batch_ys})
            avg_cost += c / total_batch
        if (epoch + 1) % display_step == 0:
            print "Epoch:", '%04d' % (epoch + 1), "cost=", "{:.9f}".format(avg_cost)
    correct_prediction = tf.equal(tf.argmax(pred, 1), tf.argmax(y, 1))
    accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
    print "Accuracy:", accuracy.eval({x: mnist.test.images[:3000], y: mnist.test.labels[:3000]})

参考

1. Softmax回归
2. 利用-tf-gradients-在-TensorFlow-中实现梯度下降

  朴素贝叶斯(naive bayes)算法是基于贝叶斯定理与特征条件独立假设的分类方法。对于给定的训练数据集，首先基于特征条件独立假设学习输入、输出的联合概率分布，然后基于此模型，对给定的输入$x$，利用贝叶斯定理求出后验概率最大的输出$y$。朴素贝叶斯法实现简单，学习与预测的效率都很高，是一种常用方法。

数学原理

  设输入空间$\mathcal X \subseteq R^n$为$n$维向量的集合，输出空间为类标记集合$\mathcal Y={c_1,c_2,…,c_K}$。输入为特征向量$x \in \mathcal X$，输出为类标记$y \in \mathcal Y$。$X$是定义在输入空间$\mathcal X$上的随机变量，$Y$是定义在输出空间$\mathcal Y$上的随机变量。$P(X,Y)$是$X$和$Y$的联合概率分布。训练数据集
$$
T={(x_1,y_1),(x_2,y_2),\ldots,(x_N,y_N)}
$$
由$P(X,Y)$独立同分布产生。

  朴素贝叶斯法通过训练数据集学习联合概率分布$P(X,Y)$。具体地，学习以下先验概率分布及条件概率分布。先验概率分布，

$$P(Y=c_K),\quad k=1,2,\ldots,K$$

条件概率分布

$$P(X=x|Y=c_k)=P(X^{(1)}=x^{(1)},…,X^{(n)}=x^{(n)}),\quad k=1,2,\ldots,K$$

于是学习到联合概率分布$$P(X,Y)$$

  条件概率分布$P(X=x|Y=c_k)$有指数量级的参数，其估计实际是不可行的。事实上，假设$x^{(j)}$可取值有$S_j$个，$j=1,2,\ldots,n$，$Y$可取值有$K$个，那么参数个数为$K \prod_{j=1}^{n}S_j$。

  朴素贝叶斯算法对条件概率分布作了条件独立的假设，即用于分类的特征在类别确定的条件下都是条件独立的。
$$
\begin{split}
P(X=x|Y=c_k) &= P(X^{(1)}=x^{(1)},\ldots,X^{(n)}=x^{(n)}|Y=c_k)\\
&=\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)
\end{split}
\tag{1}
$$
这一假设使朴素贝叶斯算法变得简单，但是有时会牺牲一定的分类准确率。

  朴素贝叶斯法分类时，对给定的输入$x$，通过学习到的模型计算后验概率分布$P(Y=c_k|X=x)$，将后验概率最大的类作为$x$的类输出。后验概率计算根据贝叶斯定理进行：

$$
P(Y=c_k|X=x)=\frac{P(X=x|Y=c_k)P(Y=c_k)}{\sum_k P(X=x|Y=c_k) P(Y=c_k)}
\tag{2}
$$
把1式代入2式，有

$$
P(Y=c_k|X=x)=\frac{P(Y=c_k)\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)}{\sum_k P(Y=c_k)\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)}
$$

所以输出为
$$
y=f(x)=\mathop{\arg\max}_{c_k} \frac{P(Y=c_k)\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)}{\sum_k P(Y=c_k)\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)}
\tag{3}
$$

由于3式中分母对所有$c_k$都相同，所以
$$
y=f(x)=\mathop{\arg\max}_{c_k} P(Y=c_k)\prod_{j=1}^{n}P(X^{(j)}=x^{(j)}|Y=c_k)
\tag{4}
$$

参数估计

极大似然估计

  在朴素贝叶斯算法中，学习意味着估计$P(Y=c_k)$和$P(X^{(j)}=x^{(j)}|Y=c_k)$。可以使用极大似然估计法估计相应的概率。先验概率$P(Y=c_k)$的极大似然估计是
$$
P(Y=c_k)=\frac{\sum_{i=1}^{N}I(y_i=c_k)}{N},\quad k=1,2,\ldots,K
$$

设第$j$个特征$x^{(j)}$可能取值的集合为$\lbrace a_{j1},a_{},\ldots,a_{jS_j} \rbrace$，条件概率$P(X^{(j)}=a_{jl}|Y=c_k)$的极大似然估计是
$$
P(X^{(j)}=a_{jl}|Y=c_k)= \frac{\sum_{i=1}^{N}I(x^{(j)}_{i}=a_{jl},y_i=c_k)}{\sum_{i=1}^{N}I(y_i=c_k)}
$$
式中，$x_{i}^{(j)}$是第$i$个样本的第$j$个特征；$a_{jl}$是第$j$个特征可能取的第$l$个值；$I$为指示函数。

贝叶斯估计

  用极大似然估计可能出现所要估计的值为$0$的情况，这时会影响到后验概率的计算结果，使分类产生偏差。解决这一问题的方法是采用贝叶斯估计，条件概率的贝叶斯估计是
$$
P_{\lambda}(X^{(j)}=a_{jl}|Y=c_k)=\frac{\sum_{i=1}^{N}I(x_{i}^{(j)}=a_{jl},y_i=c_k)+\lambda}{\sum_{i=1}^{N}I(y_i=c_k)+S_j\lambda}
$$
式中$\lambda \geq 0$。等价于在随机变量各个取值的频数上赋予一个正数$\lambda>0$。当$\lambda=0$时就是极大似然估计。常取$\lambda=1$，这时称为拉普拉斯平滑(Laplace smoothing)。显然对于任何$l=1,2,\ldots,S_j,\quad k=1,2,\ldots,K$，有

$$
P_{\lambda}(X^{(j)}=a_{jl}|Y=c_k)>0\\
\sum_{i=1}^{S_{jl}}P(X^{(j)}=a_{jl}|Y=c_k)=1
$$
同样，先验概率的贝叶斯估计是
$$
P_{\lambda}(Y=c_k)=\frac{\sum_{i=1}^{N}I(y_i=c_k)+\lambda}{N+K\lambda}
$$

示例

  训练数据如下表，学习一个朴素贝叶斯分类器并确定$x=(2,S)^T$的类标记$y$，表中$X^{(1)}、X^{(2)}$为特征，取值的集合分别为$A_1={1,2,3},A_2={S,M,L}$，$Y$为类标记，$Y \in C=\lbrace 1,-1 \rbrace$

使用贝叶斯估计进行建模

$$P(Y=1)=\frac{9+1}{15+2}=\frac{10}{17},P(Y=-1)=\frac{6+1}{15+2}=\frac{7}{17}$$
$$
P(X^{(1)}=1|Y=1)=\frac{2+1}{9+3}=\frac{3}{12},P(X^{(1)}=2|Y=1)=\frac{3+1}{9+3}=\frac{4}{12},P(X^{(1)}=3|Y=1)=\frac{4+1}{9+3}=\frac{5}{12}
$$
$$
P(X^{(1)}=1|Y=-1)=\frac{3+1}{6+3}=\frac{4}{9},P(X^{(1)}=2|Y=-1)=\frac{2+1}{6+3}=\frac{3}{9},P(X^{(1)}=3|Y=-1)=\frac{1+1}{6+3}=\frac{2}{9}
$$
$$
P(X^{(2)}=S|Y=1)=\frac{1+1}{9+3}=\frac{2}{12},P(X^{(2)}=M|Y=1)=\frac{4+1}{9+3}=\frac{5}{12},P(X^{(2)}=L|Y=1)=\frac{4+1}{9+3}=\frac{5}{12}
$$
$$
P(X^{(2)}=S|Y=-1)=\frac{3+1}{6+3}=\frac{4}{9},P(X^{(2)}=M|Y=-1)=\frac{2+1}{6+3}=\frac{3}{9},P(X^{(2)}=L|Y=-1)=\frac{1+1}{6+3}=\frac{2}{9}
$$

对于给定的$x=(2,S)^T$计算，
$$P(Y=1)P(X^{(1)}=2|Y=1)P(X^{(2)}=S|Y=1)=\frac{10}{17} \times\frac{4}{12} \times \frac{2}{12}=0.0327$$
$$P(Y=-1)P(X^{(1)}=2|Y=-1)P(X^{(2)}=S|Y=-1)=\frac{7}{17} \times\frac{3}{9} \times \frac{4}{9}=0.0610$$
由于$P(Y=-1)P(X^{(1)}=2|Y=-1)P(X^{(2)}=S|Y=-1)$最大，所以$y=-1$

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

# coding=utf-8
from numpy import *
import os


def createVocabList(dataSet):
    vocabSet = set()
    for document in dataSet:
        vocabSet = vocabSet | set(document)
    return list(vocabSet)


# 制作词袋模型
def setOfWords2Vec(vocabList, inputSet):
    returnVec = [0] * len(vocabList)
    for word in inputSet:
        if word in vocabList:
            returnVec[vocabList.index(word)] = 1
        else:
            print 'the word {} is not in my Vocabulary!'.format(word)
    return returnVec


def trainNB0(trainMatrix, trainCategory):
    numTrainDocs = len(trainMatrix)
    numWords = len(trainMatrix[0])
    # 包含敏感词的此条出现的概率
    pAbusive = sum(trainCategory) / float(numTrainDocs)
    # 构造单词出现次数列表 [1,1,...,1],这里使用了拉普拉斯平滑的思想
    p0Num = ones(numWords)
    p1Num = ones(numWords)
    p0Denom = 2.0
    p1Denom = 2.0
    for i in range(numTrainDocs):
        # 是否是侮辱性文件
        if trainCategory[i] == 1:
            # 对侮辱性文件的向量进行加和,[0,1,1,....] + [0,1,1,....]->[0,2,2,...]
            p1Num += trainMatrix[i]
            # 计算所有侮辱性文件中出现的单词总数
            p1Denom += sum(trainMatrix[i])
        else:
            p0Num += trainMatrix[i]
            p0Denom += sum(trainMatrix[i])
    # 在1类别下，每个单词出现的概率,[1,2,3,5]/90->[1/90,...]
    p1Vect = log(p1Num / p1Denom)
    #  在0类别下，每个单词出现的概率
    p0Vect = log(p0Num / p0Denom)
    return p0Vect, p1Vect, pAbusive


def classifyNB(vec2Classify, p0Vec, p1Vec, pClass1):
    # p1Vec是每个词在类别1里出现的概率，vec2Classify表示这个词在该文章中是否出现过，乘起来表示文章中每个词出现的概率
    # 这里进行了log变换，所以写成了加法
    p1 = sum(vec2Classify * p1Vec) + log(pClass1)
    p0 = sum(vec2Classify * p0Vec) + log(1 - pClass1)
    if p1 > p0:
        return 1
    else:
        return 0


def testParse(bigString):
    import re
    listOfTokens = re.split('r\W*', bigString)
    return [tok.lower() for tok in listOfTokens if len(tok) > 2]


def spamTest():
    docList = []
    classList = []
    fullText = []
    for i in range(1, 26):
        # 读取非垃圾邮件文本
        wordList = testParse(open(os.getcwd() + '/email/spam/{}.txt'.format(i)).read())
        docList.append(wordList)
        fullText.append(wordList)
        classList.append(1)
        # 读取垃圾邮件文本
        wordList = testParse(open(os.getcwd() + '/email/ham/{}.txt'.format(i)).read())
        docList.append(wordList)
        fullText.append(wordList)
        classList.append(0)
    # 建立词集
    vocabList = createVocabList(docList)
    # 训练集语料的下标
    trainingSet = range(50)
    testSet = []
    # 选10个语料作为测试集
    for i in range(10):
        randIndex = int(random.uniform(0, len(trainingSet)))
        testSet.append(trainingSet[randIndex])
        # 删除后训练集变成40个
        del (trainingSet[randIndex])
    trainMat = []
    trainClasses = []
    for docIndex in trainingSet:
        trainMat.append(setOfWords2Vec(vocabList, docList[docIndex]))
        trainClasses.append(classList[docIndex])
    # 计算各种概率
    p0V, p1V, pSpam = trainNB0(array(trainMat), array(trainClasses))
    errorCount = 0
    # 进行测试
    for docIndex in testSet:
        # 构建词向量
        wordVector = setOfWords2Vec(vocabList, docList[docIndex])
        # 进行分类
        if classifyNB(array(wordVector), p0V, p1V, pSpam) != classList[docIndex]:
            errorCount += 1
    print 'the error rate is {}'.format(float(errorCount) / len(testSet))


def main():
    spamTest()


if __name__ == '__main__':
    main()

代码详见朴素贝叶斯代码实现

参考

1. 机器学习实战，人民邮电出版社
2. 统计学习方法，清华大学出版社
3. AI Learning

简介

  xgboost(Extreme Gradient Boosting)是一款常用的机器学习框架，其诞生于2014年12月，因具有良好的学习效果以及高效的训练速度得到广泛关注。仅在2015年，在Kaggle竞赛中获胜的29个算法中，有17个使用了XGBoost库，而作为对比，近年大热的深度神经网络方法，这一数据则是11个。在KDDCup 2015竞赛中，排名前十的队伍全部使用了XGBoost库。XGBoost不仅学习效果很好，而且速度也很快，相比梯度提升算法在另一个常用机器学习库scikit-learn中的实现，XGBoost的性能经常有十倍以上的提升。

  xgboost与另一种ensemble模型gbdt(Gradient boosting decision tree)有着紧密的联系，两者均与决策树、提升模型、梯度提升等概念相关。本文从决策树出发，逐步进行xgboost推导，并结合实例分析手动计算过程，最终给出代码实现，希望有助于大家理解xgboost相关原理。

决策树(Decision Tree)

  决策树是传统机器学习中常见的一种模型，其易于理解，可解释性强，计算速度快。

分类决策树

  构建分类决策树模型的过程中最关键的步骤在于如何选择结点分裂准则，常见的分裂准则有ID3、C4.5等，相关概念可参考决策树。

CART

  CART模型是决策树中比较特殊的一种模型，其既可以用于分类又可以用于回归预测，分类树采用基尼指数最小化准则，回归采用平方误差最小化准备，相关概念可参考CART-分类与回归树。

提升树(Boosting Tree)

提升方法(boosting)

  提升方法的思想：对于一个复杂任务来说，将多个专家系统的判断进行适当的综合所得出的判断，要比其中任何一个专家单独的判断好，简而言之就是“三个臭皮匠顶个诸葛亮”的道理。

  对于分类问题而言，给定一个训练样本集，求比较粗糙的分类规则（弱分类器）要求比精确的分类规则（强分类器）容易得多。提升方法就是从弱学习算法出发，反复学习，得到一系列弱分类器，然后组合这些弱分类器，构成一个强分类器。

加法模型

  加法模型表示如下

$$f(x)=\sum_{m=1}^{M}\beta_mb(x;\gamma_m)$$
其中$b(x;\gamma_m)$为基函数，$\gamma_m$为基函数的参数，$\beta_mb$为基函数的系数。

• 注：傅立叶变换可以看成一种加法模型

前向分布算法

  在给定训练数据和损失函数$L(y,f(x))$的情况下，学习加法模型$f(x)$成为经验风险极小化问题：

$$\mathop{min}_{\beta_m,\gamma_m}\sum_{i=1}^{N}L[y_i,\sum_{m=1}^{M}\beta_mb(x_i;\gamma_m)]$$

使用前向分布算法解决这一优化问题的思路是：因为学习的是加法模型，如何能够从前向后，每一步只学习一个基函数及其系数，逐步逼近优化目标函数，就可以简化优化的复杂度。具体地，只需优化如下损失函数：

$$\mathop{min}_{\beta,\gamma}\sum_{i=1}^{N}L(y_i,\beta b(x_i;\gamma))$$

提升树模型

  以决策树为基函数的提升方法称为提升树。提升树模型可以表示为决策树的加法模型

$$f_M(x)=\sum_{m=1}^{M}T(x,\Theta_m)$$

其中$T(x,\Theta_m)$表示决策树，$\Theta_m$表示决策树的参数， $M$表示决策树的个数。对于分类问题决策树采用二叉分类树，对于回归问题决策树采用二叉回归树。

  前向回归算法：

$$
\begin{split}
f_0(x) &= 0 \\
f_m(x) &= f_{m-1}(x)+T(x;\Theta_m),\quad m=1,2,…,M \\
f_M(x) &= \sum_{m=1}^{M}T(x;\Theta_m)
\end{split}
$$

给定当前模型$f_{m-1}(x)$，需求解，
$$
\hat \Theta_m=\mathop{\arg\min}_{\Theta_m}\sum_{i=1}^{N}L(y_i,f_{m-1}(x_i)+T(x;\Theta_m))
$$

得到$\hat \Theta_m$即，第$m$棵树的参数。

  当使用平方误差损失函数的时候，
$$L(y,f(x))=(y-f(x))^2$$

其损失变为

$$
\begin{split}
L(y,f_{m-1}(x)+T_m(x;\Theta))&=[y-f_{m-1}(x) - T_m(x;\Theta)]^2\\
&=[r-T_m(x;\Theta)]^2
\end{split}
$$

其中$r=y-f_{m-1}(x)$就是当前拟合数据的残差。所以对于回归问题（回归问题的损失函数一般为MSE）的提升树算法来说，只需要简单拟合当前模型的残差。

梯度提升决策树(GBDT)

  提升树利用加法模型与前向分布算法实现学习的优化过程，当损失函数是平方损失函数时，每一步优化比较简单。但对于一般损失函数而言，往往每一步优化并不那么容易，所以需要寻找其他方式求解。

  损失函数如下
$$
L(\Theta_m)=\sum_{i=1}^{N}L(y_i,f_{m-1}(x_i)+T(x;\Theta_m))
\tag{4.1}
$$

我们的目的是让4.1式最小，尝试用泰勒一阶展开有

$$
L(\Theta_m) \approx \sum_{i=1}^{N}[L(y_i,f_{m-1}(x_i)) + \frac{\partial{L(y_i,f_{m-1}(x_i))}}{\partial{L(f_{m-1}(x_i))}} \cdot T(x;\Theta_m)]
\tag{4.2}
$$

• 注： 泰勒一阶展开公式$f(x+\Delta x)=f(x) + \frac{\partial f(x)}{\partial x} \cdot \Delta x$

由于迭代后损失函数减小，所以
$$
L(\Theta_m)- \sum_{i=1}^{N}L(y_i,f_{m-1}(x_i))<0
$$
即
$$
\sum_{i=1}^{N}[\frac{\partial{L(y_i,f_{m-1}(x_i))}}{\partial{L(f_{m-1}(x_i))}} \cdot T(x;\Theta_m)]<0
\tag{4.3}
$$
如果让提升树$T(x;\Theta_m)$每次拟合$L(y_i,f_{m-1}(x_i))$的负梯度

$$
T(x;\Theta_m)=-\frac{\partial{L(y_i,f_{m-1}(x_i))}}{\partial{L(f_{m-1}(x_i))}}
\tag{4.4}
$$
可以保证4.3式始终小于0。

例子

  训练数据如下表所示，$x$的取值范围为区间$[0.5,10.5]$，$y$的取值范围为区间$[5.0,10.0]$，使用GBDT算法解决该回归问题。

  回归问题的损失函数使用平方误差损失函数，所以该问题退化为使用回归树拟合残差。

• 求解第一棵回归树

  由于只有$x$一个变量，因此最优切分变量为$x$。接下来假设9个切分点为$[1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5]$，计算每个切分点的输出值。如$s=1.5$时，$R_1={1},R_2={2,3,4,5,6,7,8,9,10}$，这两个区域的输出值分别为$c_1=5.56,c_2=\frac{1}{9}(5.7+5.91+6.4+6.8+7.05+8.9+8.7+9+9.05)=7.50$，所以有

接下来计算每个切分点的误差，如$s=1.5$时，$loss(s=1.5)=\frac{1}{2}(5.56-5.56)^2+\sum_{i=2}^{10}\frac{1}{2}(y_i-7.5)^2=0+15.72=15.72$，所以有

其中当$s=6.5$时，$loss(s)$最小，因此第一个划分变量是$j=x,s=6.5$。所以回归树为

$$
T_1(x) = \begin{cases}
6.24\quad x \lt 6.5 \\
8.91 \quad x \geq 6.5 \\
\end{cases}
$$
$$ f_1(x) = T_1(x)$$
残差为$r_{2i}= y_i - f_1(x)$，如下

$$L(y,f_1(x)) = \sum_{i=1}^{10}(y_i-f_1(x))^2=1.93$$

• 接下来对上一步残差进行拟合

其中当$s=3.5$时，$loss(s)$最小，因此第一个划分变量是$j=x,s=3.5$。
所以回归树为
$$
T_2(x) = \begin{cases}
-0.52\quad x \lt 3.5 \\
0.22 \quad x \geq 3.5 \\
\end{cases}
$$
$$ f_2(x) = f_1(x)+T_2(x) = \begin{cases}
5.72\quad x \lt 3.5 \\
6.46 \quad 3.5 \le x \lt 6.5 \\
9.13 \quad x \ge 6.5 \\
\end{cases} $$
用$f_2(x)$拟合训练数据的平方损失误差是
$$L(y,f_2(x)) = \sum_{i=1}^{10}(y_i-f_2(x))^2=0.79$$

  假设此时损失函数满足要求，则最终的提升树为
$$
f(x) = f_2(x)=
\begin{cases}
5.72\quad x \lt 3.5 \\
6.46 \quad 0.5 \le x \lt 6.5 \\
9.13 \quad x \ge 6.5 \\
\end{cases}
$$

xgboost

  xgboost的损失函数与gbdt有所不同，其形式如下
$$
L(\Theta_m)=\sum_{i=1}^{N}L(y_i,f_{m-1}(x_i)+T(x;\Theta_m)) + \Omega(T(x;\Theta_m))+ C
\tag{5.1}
$$
其中$i$是指地$i$个样本，$f_{m-1}(x_i)$表示第$m-1$个模型对样本$i$的预测值，$T(x;\Theta_m)$表示第$m$个树模型，$\Omega(T(x;\Theta_m)$表示正则项，$C$表示计算过程中的一些常数项。

对5.1式连加部分尝试用泰勒二阶展开有
$$
L(\Theta_m) \approx \sum_{i=1}^{N}[L(y_i,f_{m-1}(x_i)) + \frac{\partial{L(y_i,f_{m-1}(x_i))}}{\partial{L(f_{m-1}(x_i))}} \cdot T(x;\Theta_m) + \frac{1}{2} \cdot \frac{\partial^2{L(y_i,f_{m-1}(x_i))}}{\partial{L(f_{m-1}(x_i))^2}} \cdot T^2(x;\Theta_m)]+ \Omega(T(x;\Theta_m)) + C
\tag{5.2}
$$

• 注： 泰勒二阶展开式$f(x+\Delta x)=f(x) + \frac{\partial f(x)}{\partial x} \cdot \Delta x + \frac{\partial^2 f(x)}{\partial^2 x} \cdot \Delta^2 x$

此处定义样本$i$的一阶和二阶导数为
$$
g_i=\frac{\partial{L(y_i,f_{m-1}(x_i))}}{\partial{L(f_{m-1}(x_i))}} \\
h_i=\frac{\partial^2{L(y_i,f_{m-1}(x_i))}}{\partial{L(f_{m-1}(x_i))^2}}
$$
所以5.1式可以表示为
$$
L(\Theta_m) \approx \sum_{i=1}^{N}[L(y_i,f_{m-1}(x_i)) + g_i \cdot T(x;\Theta_m) + \frac{1}{2} \cdot h_i \cdot T^2(x;\Theta_m)]+ \Omega(T(x;\Theta_m)) + C
\tag{5.3}
$$
其中$\sum_{i=1}^{N}L(y_i,f_{m-1}(x_i))$是定值可以和$C$项合并，有
$$
L(\Theta_m) \approx \sum_{i=1}^{N}[g_i \cdot T(x;\Theta_m) + \frac{1}{2} \cdot h_i \cdot T^2(x;\Theta_m)]+ \Omega(T(x;\Theta_m)) + C
\tag{5.4}
$$
  决策树的正则项一般考虑的是叶子结点的个数和叶子权值，比较常见的正则项如下
$$
\Omega(T(x;\Theta_m)) = \gamma \cdot T_t+\lambda \cdot \frac{1}{2}\sum_{j=1}^{T}w^2_j
$$
5.4式变为
$$
L(\Theta_m) \approx \sum_{i=1}^{N}[g_i \cdot T(x;\Theta_m) + \frac{1}{2} \cdot h_i \cdot T^2(x;\Theta_m)]+ \gamma \cdot T_t+\lambda \cdot \frac{1}{2}\sum_{j=1}^{T}w^2_j + C
\tag{5.5}
$$
  为了进一步整合3式，对$T(x;\Theta_m)$的表达形式进行修改，$T(x;\Theta_m)$重新定义为
$$
T(x;\Theta_m)=w_{q(x)}
$$
其中$w$是叶子结点的权重值，$q$函数表示输入$x$和叶子结点的对应关系。所以5.5式变为
$$
L(\Theta_m) \approx \sum_{i=1}^{N}[g_i \cdot w_{q(x_i)} + \frac{1}{2} \cdot h_i \cdot w^2_{q(x_i)}]+ \gamma \cdot T_t+\lambda \cdot \frac{1}{2}\sum_{j=1}^{T}w^2_{q(x_i)} + C
\tag{5.6}
$$
对5.6式做进一步推导，$\sum_{i=1}^{N}[g_i \cdot w_{j}(x_i)]$表示第$m$棵树的输出的乘积，这是从样本角度来看的；从树的角度来看，每个样本一定且唯一被分到其中的一个叶子结点，一个叶子结点可以分到多个样本，可以将每个叶子结点的样本梯度相加，再乘以叶子结点的权值，最后对所有叶子结点求和，所以有
$$
\sum_{i=1}^{N}g_i \cdot w_{q(x_i)}=\sum_{j=1}^{T}(\sum_{i \in I_j}g_i)w_{j} \\
\sum_{i=1}^{N}h_i \cdot w^2_{q(x_i)}=\sum_{j=1}^{T}(\sum_{i \in I_j}h_i)w^2_{j}
$$

把位于同一叶子上的一阶偏导和记做$G$，同一叶子上的二阶偏导和记为$H$，所以最终1式变为
$$
\begin{split}
L(\Theta_m) &\approx \sum_{j=1}^{T}[(\sum_{i \in I_j}g_i)w_{j}+ \frac{1}{2} \cdot (\sum_{i \in I_j}h_i)w^2_{j}]+ \gamma \cdot T_t+\lambda*\frac{1}{2}\sum_{j=1}^{T}w^2_{j} + C \\
&\approx \sum_{j=1}^{T}[G_i \cdot w_{j}+ \frac{1}{2} \cdot (\lambda+ H_i )w^2_{j}]+ \gamma \cdot T_t + C
\end{split}
\tag{5.7}
$$
5.7式对$w_mj$求偏导令其为0，有
$$
\begin{split}
\frac{\partial L(\Theta_m)}{\partial w_{mj}}=0 =>\\
G_j+(H_j+\lambda)=0 =>\\
w_{j}=-\frac{G_j}{H_j+\lambda}
\end{split}
$$
然后带回公式7有
$$
L(\Theta_m) \approx -\frac{1}{2}\sum_{j=1}^{T}\frac{G^2_j}{H_j+\lambda} + \gamma \cdot T_t
\tag{5.8}
$$
公式5.8就是最终的损失函数。

  该损失函数怎么使用呢？其实它可以指导建树，建树的时候最关键的一步就是选择一个分裂的准则，也就如何评价分裂的质量。比如在上面GBDT的介绍里，选择MSE来评价分裂的质量。在xgboost里分裂准则直接与损失函数挂钩，具体来说是通过损失函数计算增益$Gain$
$$
Gain=\frac{1}{2}[\frac{G^2_L}{H_L+\lambda}+\frac{G^2_R}{H_R+\lambda}-\frac{(G_L+G_R)^2}{(H_L+H_R)+\lambda}] - \gamma
$$
其中$ \frac{(G_L+G_R)^2}{(H_L+H_R)+\lambda}$表示对于一个结点，如果不分裂的话此时的损失。如果在这个节点进行分裂，则分裂后的左右子结点的损失分别为$\frac{G^2_L}{H_L+\lambda}$，$\frac{G^2_R}{H_R+\lambda}$，找到一种分裂使$Gain$最大。

示例

  数据集如下，15条样本数据，2个特征：

参数设置：树的深度设置为3，棵树设置为2，学习率设置为0.1，两个正则参数$\lambda=1,\gamma=0$。分类问题，损失函数选择logloss。

  下面计算logloss的一阶导数和二阶导数，假设预测值为$\hat{y_i}$。
$$
L(y_i,\hat{y_i})=-y_ilog(y_{i,pred}) - (1-y_i)log(1-y_{i,pred})
$$
其中$y_{i,pred}=\frac{1}{1+e^{-\hat{y_i}}}$

一阶导数
$$
L’(y_i,\hat{y_i})=(\frac{y_i}{y_{i,pred}}-\frac{1-y_i}{1-y_{i,pred}})y^{‘}_{i,pred}=\frac{y_i-y_{i,pred}}{y_{i,pred}(1-y_{i,pred})}*y_{i,pred}(1-y_{i,pred})=y_i-y_{i,pred}
$$
二阶导数
$$
L’’(y_i,\hat{y_i})=(y_i-y_{i,pred})’=-y_{i,pred}(1-y_{i,pred})
$$

设置初始预测值$y_{i,pred}=0.5$，举例计算$ID_1$样本的一阶二阶导数。
$g_1=y_i-y_{i,pred}=0-0.5=-0.5$
$h_1=-y_{i,pred}(1-y_{i,pred})=-0.5*(1-0.5)=-0.25$
所以有

进行结点分裂。首先看特征$x_1$，其取值为$[1,2,3,6,7,8,9,10]$，需要依次计算每个取值进行分裂的增益。

以取值1为划分时

左结点样本$I_L=[]$，右结点样$I_L=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]$

左结点为空，所以$G_L=0$，$H_L=0$

右结点一阶导数和为$G_R=\sum_{i \in I_R}g_i=-1.5$，二阶导数和为$H_R=\sum_{i \in I_R} h_i=3.75$

所以增益$Gain=0$。

以取值2为划分时

左结点样本$I_L=[1,4]$，右结点样$I_L=[2,3,5,6,7,8,9,10,11,12,13,14,15]$

左结一阶导数和为$G_L=\sum_{i \in I_L}g_i=0.5-0.5=0$，二阶导数和为$H_L=\sum_{i \in I_L} h_i=0.25+0.25=0.5$

右结点一阶导数和为$G_R=\sum_{i \in I_R}g_i=0.5-0.5+…+0.5-0.5=-1.5$，二阶导数和为$H_R=\sum_{i \in I_R} h_i=0.25+…+0.25=3.25$

所以增益为
$$
\begin{split}
Gain&=[\frac{G^{2}_{L}}{H_{L}+\lambda}+\frac{G^{2}_{R}}{H_{R}+\lambda}-\frac{(G_{L}+G_{R})^2}{(H_{L}+H_{R})+\lambda}]\\
&=[0+\frac{1.5^2}{(-3.25)+1}-\frac{(0+1.5)^2}{(-0.5-3.25)+1}]\\
&=0.0557275541796
\end{split}
$$

所以有

特征$x1$以$x1<10$分裂时可以得到最大的增益0.615204678363。

看特征$x_2$，其取值为$[-5,-2,0,2,5]$，需要依次计算每个取值进行分裂的增益。

以特征$x_2$分裂，最大的增益是0.218623481781，小于以$x_1$分裂的增益，所以第一次分裂以$x_1<10$进行分裂。

由于设置的最大深度是3，此时只有1层，所以还需要继续往下分裂。

左子节点的样本集合$I_L=[1,2,3,4,5,6,7,8,9,10,11,12,14,15] $

右子节点的样本集合$I_R=[13]$

右子结点只剩下一个样本，不需要再进行分裂，可以计算结点值
$$
w=-\frac{G_R}{H_R+\lambda}=-\frac{0.5}{1+0.25}=-0.4
$$
下面对左子结点$I_L=[1,2,3,4,5,6,7,8,9,10,11,12,14,15]$进行分裂，分裂的时候把此时的结点看成根节点，同样也是需要遍历所有特征(x1,x2)的所有取值作为分裂点，选取增益最大的点。

遍历$x_1$

可以发现$x_1$在$x<3$时能取得最大增益0.253968253968。

遍历$x_2$

可以发现$x_2$在$x<2$时能取得最大增益0.4444。

综上本次分裂应该选择$x<2$作为分裂点，分裂后左右子结点集合如下

左子节点的样本集合$I_L=[1,3,5,6,8,10,11,15] $

右子节点的样本集合$I_R=[2,4,7,9,12,14]$

到此为止，第一棵树构建完毕，下面给出表达形式。根据加法模型，第一棵树的预测结果应该为
$$y=f_0(x_i)+f_1(x_i)$$
其中$f_0(x_i)$是初始值，可以看成第0棵树。在第1步中，我们假设初始值$y_{i,pred}$=0.5，由于使用logloss损失函数，所以其原始输出需要进行一个sigmoid的反变换$f_0=ln(\frac{0.5}{1-0.5})=0$，所以第一棵树的预测结果就是
$$y=f_0(x_i)+f_1(x_i)=0+w_{q(x_i)}=w_{q(x_i)}$$
第一棵树的结构如下图

  建立第2棵树，整个过程和第1棵树一样，但是需要更新$y_{i,pred}$，也就是说在第一棵树的基础上进行第二棵树的构建，这里的思想和GBDT一致。

  以ID=1样本为例，更新$y_{i,pred}$。对于ID=1经过第一棵树，落在-0.04结点，经过sigmoid映射后的值为$p_{1,pred}=\frac{1}{1+e^{0.04}}=0.490001$，所以有

根据上表可推导出第二棵树的结构如下图

• 注： 由于hexo排版问题，表格中的数据进行了截位处理，可能精度不够，可以使用下面的程序进行计算。

代码

  针对上面的例子，简单实现了xgboost的求解过程，代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212

# coding=utf-8
import math

import numpy as np
import pandas as pd
import pygraphviz as pgv
from TreeNode import TreeNode


# 选择最佳切分点
def chooseBestSplit(dataMat, hi, gi, r, T=0):
    # 获取特征数量
    featNum = dataMat.shape[1] - 1
    # 最大增益
    max_gain = float("-inf")
    # 切分特征下标
    split_feat = -1
    # 切分点的值
    split_value = -1
    # 切分后左子树的一阶导数和
    split_g_l = -1
    # 切分后右子树的一阶导数和
    split_g_r = -1
    # 切分后左子树的二阶导数和
    split_h_l = -1
    # 切分后右子树的二阶导数和
    split_h_r = -1
    # 遍历特征
    for feat in range(featNum):
        print '{}th feature'.format(feat)
        featList = dataMat[:, feat].T.tolist()[0]
        uniqueVals = sorted(set(featList))
        # 遍历特征值
        for value in uniqueVals:
            # 挑选比特征值小的样本，即左子树样本
            left_points = np.where(dataMat[:, feat] < value)
            # 挑选比特征值大的样本，即右子树样本
            right_points = np.where(dataMat[:, feat] >= value)
            # 左子树一阶导数和
            g_l = G(left_points, gi)
            # 右子树一阶导数和
            g_r = G(right_points, gi)
            # 左子树二阶导数和
            h_l = H(left_points, hi)
            # 右子树二阶导数和
            h_r = H(right_points, hi)
            # 计算分裂增益
            g = gain(g_l, h_l, g_r, h_r, r)
            print '{}-g_l:{}, h_l:{}, g_r:{}, h_r:{}-g:{}'.format(value, g_l, h_l, g_r, h_r, g)
            if g >= max_gain:
                max_gain = g
                split_feat = feat
                split_value = value
                split_g_l = g_l
                split_g_r = g_r
                split_h_l = h_l
                split_h_r = h_r
    return split_feat, split_value, split_g_l, split_g_r, split_h_l, split_h_r


# 进行分裂
def binSplitDataSet(dataMat, split_feat, split_value, hi, gi):
    left_points = np.where(dataMat[:, split_feat] < split_value)
    right_points = np.where(dataMat[:, split_feat] >= split_value)
    return dataMat[left_points[0]], dataMat[right_points[0]], hi[left_points[0]], hi[right_points[0]], gi[
        left_points[0]], gi[right_points[0]]


def createTree(dataMat, hi, gi, depth=0, max_depth=3, r=1, eta=0.1):
    # 选择最佳切分特征、最佳切分点、左右子树的一阶二阶导数和
    feat, val, g_l, g_r, h_l, h_r = chooseBestSplit(dataMat, hi, gi, r)
    root = TreeNode(feat, val)
    # 结点分裂，返回左右子结点每个样本的一阶导数和二阶导数值
    lSet, rSet, hi_l, hi_r, gi_l, gi_r = binSplitDataSet(dataMat, feat, val, hi, gi)
    # 如果数据集中样本个数大于1并且树的深度小于3层
    if len(set(lSet[:, -1].T.tolist()[0])) > 1 and depth + 1 < max_depth:
        root.left = createTree(lSet, hi_l, gi_l, depth + 1)
    else:
        leaf = TreeNode(-1, -1)
        leaf.weight = eta * cal_weight(g_l, h_l, r)
        leaf.isLeaf = True
        root.left = leaf

    if len(set(rSet[:, -1].T.tolist()[0])) > 1 and depth + 1 < max_depth:
        root.right = createTree(rSet, hi_r, gi_r, depth + 1)
    else:
        leaf = TreeNode(-1, -1)
        leaf.weight = eta * cal_weight(g_r, h_r, r)
        leaf.isLeaf = True
        root.right = leaf
    return root


# 计算叶子结点权重
def cal_weight(g, h, r):
    return -g / (h + r)


# 计算一阶导数和
def G(points, hi):
    return np.sum(hi[points])


# 计算二阶导数和
def H(points, gi):
    return np.sum(gi[points])


# 计算增益
def gain(g_l, h_l, g_r, h_r, r):
    left_gain = math.pow(g_l, 2) / (h_l + r)
    right_gain = math.pow(g_r, 2) / (h_r + r)
    all_gain = math.pow(g_l + g_r, 2) / (h_l + h_r + r)
    return left_gain + right_gain - all_gain


# 计算每个样本的一阶导数值
def g_i(y_pred, y_i):
    return y_pred - y_i


# 计算每个样本的二阶导数值
def h_i(y_pred):
    return y_pred * (1 - y_pred)


def load_data(path):
    data = pd.read_csv(path, dtype=np.float64, delimiter='\t', header=None)
    dataMat = np.mat(data.values)
    return dataMat


# 初始化一阶导数和二阶导数值
def init_base_score(trees, dataMat):
    label = dataMat[:, -1]
    if len(trees) == 0:
        base_score = np.zeros((dataMat.shape[0], 1))
        # 初始值设置为0.5，即base_score
        base_score += 0.5
        gi = g_i(base_score, label)
        hi = h_i(base_score)
    else:
        # 上一次预测值
        pred_res = predict(trees, dataMat)
        gi = g_i(pred_res, label)
        hi = h_i(pred_res)
    return hi, gi


# 预测函数
def predict(trees, dataMat):
    pred_res = np.zeros((dataMat.shape[0], 1), dtype=np.float64)
    for tree in trees:
        for i in range(dataMat.shape[0]):
            # 获取输入数据在每棵树上的输出
            weight = tree.get_weight(dataMat[i, :])
            # sigmoid变换
            pred_res[i, 0] += 1 / (1 + math.exp(-weight))
    return pred_res


# 画图
def draw_tree(root, i):
    A = pgv.AGraph(directed=True, strict=True)
    display(root, A)
    A.graph_attr['epsilon'] = '0.01'
    print A.string()  # print dot file to standard output
    A.write('tree_{}.dot'.format(i))
    A.layout('dot')  # layout with dot
    A.draw('tree_{}.png'.format(i))  # write to file


def display(root, A):
    if not root:
        return
    A.add_node(root.uid, label='x[{}]<{}'.format(root.split_feat, root.split_value))
    if root.left:
        if root.left.isLeaf:
            A.add_node(root.left.uid, label='leaf={}'.format(root.left.weight))
            A.add_edge(root.uid, root.left.uid, label='yes', color='red')
        else:
            A.add_node(root.left.uid, label='x[{}]<{}'.format(root.left.split_feat, root.left.split_value))
            A.add_edge(root.uid, root.left.uid, label='yes', color='red')
            display(root.left, A)

    if root.right:
        if root.right.isLeaf:
            A.add_node(root.right.uid, label='leaf={}'.format(root.right.weight))
            A.add_edge(root.uid, root.right.uid, label='no', color='blue')
        else:
            A.add_node(root.right.uid, label='x[{}]<{}'.format(root.right.split_feat, root.right.split_value))
            A.add_edge(root.uid, root.right.uid, label='no', color='blue')
            display(root.right, A)


def main():
    dataMat = load_data('./data.txt')
    root = None
    trees = []
    tree_num = 2
    for i in range(tree_num):
        print '{}th tree building'.format(i)
        hi, gi = init_base_score(trees, dataMat)
        root = createTree(dataMat, hi, gi)
        trees.append(root)
    for i in range(len(trees)):
        print trees[i]
        draw_tree(trees[i], i)


if __name__ == '__main__':
    main()

代码详见xgboost简单实现

参考

1. 机器学习实战，人民邮电出版社
2. 统计学习方法，清华大学出版社
3. xgboost原理分析以及实践
4. XGBoost Documentation

  分类与回归树(classification and regression tree, CART)模型是应用比较广泛的决策树学习方法之一，它既可以用于分类也可以用于回归。CART算法使用二元切分递归地处理每个特征，如果特征值大于给定值就走左子树，否则就走右子树。

最小二乘回归树生成算法

数学描述

  假设$X$和$Y$分别为输入和输出变量，并且$Y$是连续变量，给定训练数据集
$$D={(x_1,y_1),(x_2,y_2),…,(x_N,y_N)}$$
生成对应的回归树。

  一个回归树对应着输入空间(即特征空间)的一个划分以及在划分的单元上的输出值。假设已将输入空间划分为$M$个单元$R_1,R_2,…,R_m$，并且在每个单元$R_m$上有一个固定的输出值$C_m$，于是回归树模型可以表示为

$$f(x) = \sum_{m=1}^{M}c_mI(x \in R_m)$$

其中函数$I$为指示函数。

  当输入空间的划分确定时，可以用平方误差$\sum_{x_i \in R_m}(y_i-f(x_i))^2$来表示回归树对于训练数据的预测误差。用平方误差最小的准则求解每个单元上的最优输出值，可以发现单元$R_m$上的最优$\hat c_m$是$R_m$上所有示例$x_i$对应的输出$y_i$的均值，即
$$\hat c_m=avg(y_i|x_i \in R_m)$$

  上面内容讲述了在输入空间划分确定的情况下，如何衡量决策树的效果以及在效果最好的情况下反推每个区域输出值$c_m$。接下来我们来看如何划分输入空间。

  首先，CART树是一棵二叉树，所以我们需要选择一个特征$x^{(j)}$和它的取值$s$，作为整个数据集的切分量和切分点，并定义由其切分的两个区域
$$R_1(j,s)=\lbrace x|x^{(j)} \leq s \rbrace和R_2(j,s)=\lbrace x|x^{(j)} \gt s \rbrace $$
然后寻找最优缺切分变量$j$和最优切分点$s$，即求解

$$min_{j,s}[min_{c1}\sum_{x_i \in R_1(j,s)}(y_i-c_1)^2+min_{c2}\sum_{x_i \in R_2(j,s)}(y_i-c_2)^2]$$

• 注：这里是要找到使$R_1$数据的预测误差最小的$c_1$和使$R_2$数据的预测误差最小的$c_2$，之前说过使误差最小的取值就是分别取两个数据集的均值，然后还要保证$R_1$和$R_2$误差和最小。

对于固定的输入变量$j$可以找到最优的切分点$s$：

$$\hat c_1=avg(y_i|x_i \in R_1(j,s)) 和 \hat c_2=avg(y_i|x_i \in R_2(j,s))$$

• 注：简单讲，这里的$c_1$和$c_2$是每次按照切分点$s$分割成两波数据的均值，这里不明白的可以看后面的例子。

  遍历所有输入变量，找到最优的切分变量$j$，构成一对$(j,s)$，并依此将输入空间划分为两个区域。接着对每个区域重复上述划分过程，直到满足条件为止。这样就生成了一棵回归树，这样的回归树通常称为最小二乘回归树。

示例

  上述数学描述比较晦涩，举个简单的例子，输入数据见下表

• 第一次划分

  由于只有$x$一个变量，因此最优切分变量为$x$。接下来假设9个切分点为$[1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5]$，计算每个切分点的输出值。如$s=1.5$时，$R_1={1},R_2={2,3,4,5,6,7,8,9,10}$，这两个区域的输出值分别为$c_1=5.56,c_2=\frac{1}{9}(5.7+5.91+6.4+6.8+7.05+8.9+8.7+9+9.05)=7.50$，所以有

接下来计算每个切分点的误差，如$s=1.5$时，$loss(s=1.5)=\frac{1}{2}(5.56-5.56)^2+\sum_{i=2}^{10}\frac{1}{2}(y_i-7.5)^2=0+15.72=15.72$，所以有

其中当$s=6.5$时，$loss(s)$最小，因此第一个划分变量是$j=x,s=6.5$。

$$
T=\begin{cases}
6.24 \quad x \le 6.5 \\
8.91 \quad x \gt 6.5\\
\end{cases}
$$

• 第二次划分

    第一个划分变量将数据集划分成了两部分即$R_1={1,2,3,4,5,6},R_2={7,8,9,10}$，接下来分别在$R_1,R_2$上进行划分。

  •   对于$R_1$

  取切分点$[1.5,2.5,3.5,4.5,5.5]$，对应的输出值为

  误差为

所以$s=3.5$时，$loss(s)$最小

  假设在生成3个区域后停止划分，那么最终生成的回归树如下：

$$
T=\begin{cases}
5.72 \quad x \le 3.5 \\
6.75 \quad 3.5 \lt x \leq 6.5 \\
8.91 \quad x \gt 6.5\\
\end{cases}
$$

实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

# coding=utf-8

from numpy import *


# 读取数据
def loadDataSet(path):
    dataMat = []
    fr = open(path, 'rb')
    for line in fr.readlines():
        currLine = line.strip().split('\t')
        fltLine = map(float, currLine)
        dataMat.append(fltLine)
    return mat(dataMat)


# 以value为分界点切分数据集
def binSplitDataSet(dataSet, feature, value):
    mat0 = dataSet[nonzero(dataSet[:, feature] > value)[0], :]
    mat1 = dataSet[nonzero(dataSet[:, feature] <= value)[0], :]
    return mat0, mat1


# 计算叶子结点中数据的均值
def regLeaf(dataSet):
    return mean(dataSet[:, -1])


# 计算误差，即数据集的平方误差，这里使用方差乘以总个数计算
def regErr(dataSet):
    return var(dataSet[:, -1]) * shape(dataSet)[0]


# 二元切分选择分裂点
def chooseBestSplit(dataSet, leafType=regLeaf, errType=regErr, ops=(1, 4)):
    tolS = ops[0]
    tolN = ops[1]
    if len(set(dataSet[:, -1].T.tolist()[0])) == 1:
        return None, leafType(dataSet)
    m, n = shape(dataSet)
    S = errType(dataSet)
    bestS = inf
    bestIndex = 0
    bestValue = 0
    for featIndex in range(n - 1):
        for splitVal in set(dataSet[:, featIndex].T.tolist()[0]):
            # 切分数据集
            mat0, mat1 = binSplitDataSet(dataSet, featIndex, splitVal)
            # 判断切分后的数据集的条数是否满足要求
            if (shape(mat0)[0] < tolN) or (shape(mat1)[0] < tolN):
                continue
            # 计算两个数据集的误差
            newS = errType(mat0) + errType(mat1)
            # 选取误差最小的切分方式
            if newS < S:
                bestIndex = featIndex
                bestValue = splitVal
                bestS = newS
    # 如果误差已经小于要求的误差
    if S - bestS < tolS:
        return None, leafType(dataSet)
    mat0, mat1 = binSplitDataSet(dataSet, bestIndex, bestValue)
    # 如果不满足条数要求，不再继续分裂，返回结点的值
    if (shape(mat0)[0] < tolN) or (shape(mat1)[0] < tolN):
        return None, leafType(dataSet)
    return bestIndex, bestValue


def createTree(dataSet, leafType=regLeaf, errType=regErr, ops=(1, 4)):
    # 选择最佳分裂点
    feat, val = chooseBestSplit(dataSet, leafType, errType, ops)
    if feat == None:
        return val
    retTree = {}
    retTree['spInd'] = feat
    retTree['spVal'] = val
    lSet, rSet = binSplitDataSet(dataSet, feat, val)
    # 左子树
    retTree['left'] = createTree(lSet, leafType, errType, ops)
    # 柚子树
    retTree['right'] = createTree(rSet, leafType, errType, ops)
    return retTree


def main():
    data = loadDataSet('./data.txt')
    tree = createTree(data)
    print tree


if __name__ == '__main__':
    main()

代码详见CART

• 注：此处代码参考《机器学习实战》p161-p164，原书代码中存在错误，这里进行了修改。

分类树生成算法

  分类树使用基尼指数选择最优特征，同时决定最优二值切分点。

基尼指数

  基尼指数的定义如下，假设有$K$个分类，样本点属于第$k$类的概率为$p_k$，则概率分布的基尼指数为
$$
Gini(p) = \sum_{k=1}^{K}p_k(1-p_k)=1-\sum_{k=1}^{K}p^2_k
$$

  对于给定的样本集合$D$，其基尼指数为
$$
Gini(D)=1-\sum_{i=1}^{K}(\frac{|C_k|}{|D|})^2
$$
其中，$C_k$是$D$中属于第$k$类的样本子集，$K$是类的个数。

  如果样本集合$D$根据特征$A$是否取某一可能值$a$呗分割成$D_1$和$D_2$两部分，即
$$
D_1 = {(x,y) \in D|A(x)=a}\\
D_2 = D-D_1
$$
则在特征$A$的条件下，集合$D$的基尼指数定义为
$$
Gini(D,A)=\frac{|D_1|}{|D|}Gini(D_1)+\frac{|D_2|}{|D|}Gini(D_2)
$$
基尼指数$Gini(D)$表示集合$D$的不确定性，基尼指数$Gini(D,A)$表示经$A=a$分割后集合$D$的不确定性。基尼指数越大，样本集合的不确定性也就越大，这一点和熵一致。

• 注：CART分类树为什么要使用基尼指数作为选择最优特征的标准呢？因为基尼指数不仅拥有与熵类似的性质，而且计算简便，不用使用log函数。

生成算法

  输入：训练数据集$D$，停止计算的条件；

  输出：CART决策树

  根据训练数据集，从根节点开始，递归地对每个结点进行以下操作，构建二叉决策树

1. 设结点的训练数据集为$D$，计算现有特征对该数据集的基尼指数。此时对每个特征$A$，对其可能取值的每个值$a$，根据样本点对$A=a$的测试为“是”或”否“将$D$分割成$D_1$和$D_2$两部分，并计算$A=a$的基尼指数。
2. 在所有可能的特征$A$以及它们所有可能的切分点$a$中，选择基尼指数最小的特征及其对应的切分点作为最优特征和最优切分点。依最优特征与最优切分点，从现结点生成两个子结点，将训练数据依特征分配到两个子结点中。
3. 对两个子结点递归地调用1、2，直至满足停止条件
4. 生成CART决策树

  算法停止条件是结点中的样本个数小于预定阈值，或者样本集的基尼指数小于预定阈值（样本基本属于同一类），或者没有更多特征。

示例

  数据集合如下:

  分别以$A_1$、$A_2$、$A_3$、$A_4$表示年龄、有工作、有自己的房子和信贷4个特征。

  求特征$A_1$的基尼指数，以1，2，3表示年龄的值为青年、中年、老年

$$
Gini(D,A_1=1) = \frac{5}{15}[2 \times \frac{2}{5} \times (1-\frac{2}{5})] + \frac{10}{15}[2 \times \frac{7}{10} \times (1-\frac{7}{10})]=0.44\\
Gini(D,A_1=2) = \frac{5}{15}[2\times\frac{3}{5}\times(1-\frac{3}{5})] + \frac{10}{15}[2\times\frac{6}{10}\times(1-\frac{6}{10})]=0.48\\
Gini(D,A_1=3) = \frac{5}{15}[2\times\frac{4}{5}\times(1-\frac{4}{5})] + \frac{10}{15}[2\times\frac{5}{10}\times(1-\frac{5}{10})]=0.44
$$
其中$Gini(D,A_1=1)$和$Gini(D,A_1=3)$最小，所以$A_1=1$和$A_1=3$均可做为$A_1$的最优切分点

  求特征$A_2$的基尼指数，以1，2表示有工作和没有工作
$$
Gini(D,A_2=1) = \frac{5}{15}[2\times\frac{5}{5}\times(1-\frac{5}{5})] + \frac{10}{15}[2\times\frac{4}{10}\times(1-\frac{4}{10})]=0.32
$$
$A_2$只有一个切分点即$A_2=1$

  求特征$A_3$的基尼指数，以1，2表示有房子和没有房子
$$
Gini(D,A_3=1) = \frac{6}{15}[2\times\frac{6}{6}\times(1-\frac{6}{6})] + \frac{9}{15}[2\times\frac{3}{9}\times(1-\frac{3}{9})]=0.27
$$
$A_3$只有一个切分点即$A_3=1$

  求特征$A_4$的基尼指数，以1，2，3表示信贷情况为非常好、好和一般
$$
Gini(D,A_4=1) = \frac{4}{15}[2\times\frac{4}{4}\times(1-\frac{4}{4})] + \frac{11}{15}[2\times\frac{6}{11}\times(1-\frac{6}{11})]=0.36\\
Gini(D,A_4=2) = \frac{6}{15}[2\times\frac{4}{6}\times(1-\frac{4}{6})] + \frac{9}{15}[2\times\frac{5}{9}\times(1-\frac{5}{9})]=0.47\\
Gini(D,A_4=3) = \frac{5}{15}[2\times\frac{1}{5}\times(1-\frac{1}{5})] + \frac{10}{15}[2\times\frac{8}{10}\times(1-\frac{8}{10})]=0.32
$$

其中$Gini(D,A_4=3)$最小，所以$A_4=3$为$A_4$的最优切分点

  在4个特征中，$Gini(D,A_3=1)=0.27$最小，所以选择$A_3$为最优特征，$A_3=1$为其最优切分点。其中$A_3=1$的结点包含数据4、8、9、10、11，$A_3 \neq 1$结点包含数据1，2，3，5，6，7，13，14，15。$A_3=1$结点内数据类别相同，所以形成叶子结点。对$A_3 \neq 1$结点继续重复上述过程。

  数据如下

  求特征$A_1$的基尼指数
$$
Gini(D,A_1=1) = \frac{4}{10}[2\times\frac{1}{4}\times(1-\frac{1}{4})] + \frac{6}{10}[2\times\frac{3}{6}\times(1-\frac{3}{6})]=0.45\\
Gini(D,A_1=2) = \frac{2}{10}[2\times\frac{0}{2}\times(1-\frac{0}{2})] + \frac{8}{10}[2\times\frac{4}{8}\times(1-\frac{4}{8})]=0.4\\
Gini(D,A_1=3) = \frac{4}{10}[2\times\frac{3}{4}\times(1-\frac{3}{4})] + \frac{6}{10}[2\times\frac{1}{6}\times(1-\frac{1}{6})]=0.32
$$
其中$Gini(D,A_1=3)$最小，所以$A_1=3$是$A_1$的最优切分点

  求特征$A_2$的基尼指数
$$
Gini(D,A_2=1) = \frac{3}{10}[2\times\frac{3}{3}\times(1-\frac{3}{3})] + \frac{7}{10}[2\times\frac{1}{7}\times(1-\frac{1}{7})]=0.17
$$
$A_2$只有一个切分点即$A_2=1$

  求特征$A_4$的基尼指数
$$
Gini(D,A_4=1) = \frac{1}{10}[2\times\frac{1}{1}\times(1-\frac{1}{1})] + \frac{9}{10}[2\times\frac{3}{9}\times(1-\frac{3}{9})]=0.4\\
Gini(D,A_4=2) = \frac{5}{10}[2\times\frac{3}{5}\times(1-\frac{3}{5})] + \frac{5}{10}[2\times\frac{1}{5}\times(1-\frac{1}{5})]=0.4\\
Gini(D,A_4=3) = \frac{4}{10}[2\times\frac{0}{0}\times(1-\frac{0}{0})] + \frac{6}{10}[2\times\frac{4}{6}\times(1-\frac{4}{6})]=0.27
$$

其中$Gini(D,A_4=3)$最小，所以$A_4=3$为$A_4$的最优切分点

  在4个特征中，$Gini(D,A_2=1)=0.17$最小，所以选择$A_2$为最优特征，$A_2=1$为其最优切分点。其中$A_2=1$的结点包含数据3、13、14，$A_3 \neq 1$结点包含数据1、2、5、6、7、12、15。$A_2=1$结点内数据类别相同，所以形成叶子结点。对$A_2 \neq 1$结点内数据类别相同，所以形成叶子结点。

  到此，CART树生成完毕。可以发现与之前按照ID3算法所生成的决策树一样，如下图。

参考

1. 机器学习实战，人民邮电出版社
2. 统计学习方法，清华大学出版社

  决策树(decision tree)是一种多功能的机器学习算法，它可以实现分类和回归任务，同时也是随机森林的重要组成部分，本篇博客主要讨论用于分类的决策树。

  分类型决策树是一种描述对实例进行分类的树形结构。决策树由结点（node）和有向边（directed edge）组成。结点有两种类型：内部结点（internal node）和叶结点（leaf node）。内结点表示一个特征，叶子结点表示一个类别。

相关数学概念

信息熵

  在将实例分配到子结点的过程中需要选择一个度量条件来判断该实例属于哪个子结点，常见的方法是用信息熵作为度量条件。

  在信息论中，熵（entropy）表示随机变量不确定性的度量，设$X$是一个取有限个值的离散随机变量，其概率分布为
$$P(X=x_i)=p_i, i=1,2,…,n$$
则随机变量$X$的熵的定义为
$$
H(X)=-\sum_{i=1}^{n}p_ilog_{2}p_i
\tag{1}
$$
若$p_i=0$，则定义$0log0=0$。熵越大，说明随机变量的不确定性越大。

条件熵

  设有随机变量(X,Y)，其联合概率分布为
$$P(X=x_i,Y=y_i)=p_{ij}, i=1,2,…,n;j=1,2,…,m$$
条件熵H(Y|X)表示在已知随机变量$X$的条件下随机变量$Y$的不确定性。随机变量$X$给定的条件下随机变量$Y$的条件熵H(Y|X)，定义为$X$给定条件下$Y$的条件概率分布的熵对$X$的数学期望。
$$P(Y|X)=\sum_{i=1}^{n}p_iH(Y|X=x_i)$$
其中$p_i=P(X=x_i)$

  当熵和条件熵中的概率由数据估计（特别是极大似然估计）得到时，所对应的熵和条件熵分别称为经验熵(empirical entropy)和经验条件熵(empirical conditional entropy).

信息增益

  信息增益(information gain)表示得知特征$X$的信息而使得类别$Y$的信息不确定性减少的程度。

  特征A对训练数据集D的信息增益$g(D,A)$，定义为集合D的经验熵$H(D)$与特征A给定的条件下D的经验条件熵$H(D|A)$之差，即
$$g(D,A)=H(D)-H(D|A)$$

  一般地，熵$H(Y)$与条件熵$H(Y|X)$之差称为互信息，决策树中的信息增益等价于训练数据集中类别和特征的互信息。

信息增益比

  特征A对训练数据集D的信息增益比$g_R(D,A)$定义为其信息增益$g(D,A)$与训练数据集$D$的经验熵$H(D)$之比
$$g_R(D,A)=\frac{g(D,A)}{H(D)}$$

原理

  决策树在建模过程中使用信息增益准则选择特征，给定训练集$D$和特征$A$，经验熵$H(D)$表示对数据集$D$进行分类的不确定性。而经验条件熵$H(D|A)$表示在特征$A$给定的条件下对数据集$D$进行分类的不确定性，它们的差即信息增益，表示由于特征$A$而使得对数据集$D$的分类的不确定性减少的程度。显然，对于数据集$D$而言，信息增益依赖于特征，不同的特征具有不同的信息增益，信息增益大的特征具有更强的分类能力。

  根据信息增益准则的特征选择方法是：对训练数据集$D$，计算其每个特征的信息增益，并比较它们的大小，选择信息增益大的特征。

ID3算法

  设训练数据集为$D$，$|D|$表示其样本容量，即样本个数。设有$K$个类$C_k，k=1,2,…,K$，$|C_k|$为属于$C_k$的样本个数，$\sum_{k=1}^{K}|C_k|=|D|$。设特征$A$有$n$个不同的取值${a_1,a_2,…,a_n}$，根据特征$A$的取值将$D$划分为n个子集$D_1,D_2,…,D_n$，$|D_i|$为$D_i$的样本个数，$\sum_{k=1}^{n}|D_i|=|D|$。记子集$D_i$中属于类$C_k$的样本的集合为$D_{ik}$，记$D_{ik}=D_i \cap C_k$，$|D_{ik}|$为$D_{ik}$的样本个数。

  ID3算法如下

• 输入： 训练数据集$D$，特征集$A$，阈值$\epsilon$
• 输出： 决策树$T$

1. 若$D$中所有实例属于同一类别$C_k$，则$T$为单结点树，并将类$C_k$作为该结点的类标记，返回$T$；
2. 若$A=\emptyset$，则$T$为单点树，将$D$中实例数最多的类$C_k$作为该结点的类标记，返回$T$；

3. 否则，如下计算$A$中各特征对$D$的信息增益，选择信息增益最大的特征$A_g$

   1. 计算数据集$D$的经验熵
      $$H(D)=-\sum_{k=1}^{K}\frac{|C_k|}{|D|}log_2\frac{|C_k|}{|D|}$$
   2. 计算特征集$A$对数据集$D$的经验条件熵
      $$H(D|A)=\sum_{i=1}^{n}\frac{|D_i|}{|D|}H(D_i)=\sum_{i=1}^{n}[\frac{|D_i|}{|D|}\sum_{k=1}^{K}\frac{|D_{ik}|}{|D_{i}|}log_2\frac{|D_{ik}|}{|D_{i}|}]$$
   3. 计算信息增益
      $$g(D,A)=H(D) - H(D|A)$$

4. 如果$A_g$小于阈值$\epsilon$，则置$T$为单点树，将$D$中实例数最多的类$C_k$作为该结点的类标记，返回$T$；

5. 否则，对$A_g$的每一个值$a_i$，按照$A_g=a_i$将$D$分割为若干非空子集$D_i$，将$D_i$中实例最大的类作为标记，构建子结点，由结点和子结点构成树$T$，返回$T$
6. 对第$i$个子结点，以$D_i$为训练集，以$A_i-{A_g}$为特征集，递归地调用步骤1-5，得到子树$T_i$，返回$T_i$

C4.5算法

  C4.5算法和ID3算法类似，C4.5算法对ID3算法进行了改进，在生成模型的过程中采用信息增益比来选择特征。

示例

  以贷款申请数据为例讲解决策树的生成过程，采用ID3算法。

1. 计算经验熵$H(D)$
   $$H(D)=-\frac{6}{15}log_2\frac{6}{15}-\frac{9}{15}log_2\frac{9}{15}=0.971$$

2. 计算各特征对数据集$D$的信息增益

   • 年龄
     $$
     \begin{split}
     g(D,A_1)&=H(D) - [\frac{5}{15}H(D_1)+\frac{5}{15}H(D_2)+\frac{5}{15}H(D_3)]\\
     &= 0.971 - {\frac{5}{15}[-\frac{2}{5}log_2\frac{2}{5}-\frac{3}{5}log_2\frac{3}{5}]+
     \frac{5}{15}[-\frac{3}{5}log_2\frac{3}{5}-\frac{2}{5}log_2\frac{2}{5}]+
     \frac{5}{15}[-\frac{4}{5}log_2\frac{4}{5}-\frac{1}{5}log_2\frac{1}{5}]}\\
     &=0.971-0.888\\
     &=0.083
     \end{split}
     $$
   • 有工作
     $$
     \begin{split}
     g(D,A_2)&=H(D) - [\frac{5}{15}H(D_1)+\frac{10}{15}H(D_2)]\\
     &= 0.971 - {\frac{5}{15} \times 0+
     \frac{10}{15}[-\frac{6}{10}log_2\frac{6}{10}-\frac{4}{10}log_2\frac{4}{10}]}\\
     &=0.971-0.647\\
     &=0.324
     \end{split}
     $$
   • 有自己的房子
     $$
     \begin{split}
     g(D,A_3)&=H(D) - [\frac{6}{15}H(D_1)+\frac{9}{15}H(D_2)]\\
     &= 0.971 - {\frac{6}{15} \times 0+
     \frac{9}{15}[-\frac{3}{9}log_2\frac{3}{9}-\frac{6}{9}log_2\frac{6}{9}]}\\
     &=0.971-0.551\\
     &=0.420
     \end{split}
     $$
   • 信贷情况
     $$
     \begin{split}
     g(D,A_3)&=H(D) - [\frac{5}{15}H(D_1)+\frac{6}{15}H(D_2)+\frac{4}{15}H(D_2)]\\
     &= 0.971 - {\frac{5}{15}[-\frac{1}{5}log_2\frac{1}{5}-\frac{4}{5}log_2\frac{4}{5}] +
     \frac{6}{15}[-\frac{4}{6}log_2\frac{4}{6}-\frac{2}{6}log_2\frac{2}{6}] +
     \frac{4}{15} \times 0 }\\
     &=0.971-0.608\\
     &=0.363
     \end{split}
     $$
       综上，特征$A_3$（有自己的房子）的信息增益值最大，所以选择$A_3$作为根节点的特征。接下来将数据集$D$划分为两个子集$D_1$($A_3$=是)、$D_2$($A_3$=否)，分别进行下一轮的信息增益筛选。

3. 进入子集$D_1$($A_3$=是)，由于该子集中只有一个类别，所以形成叶子结点，结点的标记就是类别。

4. 进入子集$D_2$($A_3$=否)，数据如下

   • 计算经验熵
     $$H(D_2)=-\frac{6}{9}log_2\frac{6}{9}-\frac{3}{9}log_2\frac{3}{9}=0.918$$

   • 年龄特征信息增益
     $$
     \begin{split}
     g(D_2,A_1)&=H(D) - [\frac{4}{9}H(D_1)+\frac{2}{9}H(D_2)+\frac{3}{9}H(D_3)]\\
     &= 0.971 - {\frac{4}{9}[-\frac{3}{4}log_2\frac{3}{4}-\frac{1}{4}log_2\frac{1}{4}]+
     \frac{2}{9}\times0+
     \frac{3}{9}[-\frac{2}{3}log_2\frac{2}{3}-\frac{1}{3}log_2\frac{1}{3}]}\\
     &=0.971-0.667\\
     &=0.251
     \end{split}
     $$

   • 有工作特征信息增益

   $$
   \begin{split}
   g(D_2,A_2)&=H(D) - [\frac{6}{9}H(D_1)+\frac{3}{9}H(D_2)]\\
   &= 0.971 - [\frac{6}{9}\times0+\frac{3}{9} \times 0]\\
   &= 0.918-0\\
   &=0.918
   \end{split}
   $$

   • 信贷情况特征信息增益

     $$
     \begin{split}
     g(D_2,A_4)&=H(D) - [\frac{4}{9}H(D_1)+\frac{4}{9}H(D_2)+\frac{1}{9}H(D_3)]\\
     &= 0.971 - {\frac{4}{9}\times0+
     \frac{4}{9}[-\frac{2}{4}log_2\frac{2}{4}-\frac{2}{4}log_2\frac{2}{4}]+
     \frac{1}{9}\times0}\\
     &=0.918-0.444\\
     &=0.474
     \end{split}
     $$
       综上，特征$A_2$（有工作）的信息增益值最大，所以选择$A_2$作为该结点的特征。接下来将数据集$D_2$划分为两个子集$D_{21}$($A_2$=是)、$D_{22}$($A_2$=否)，分别进行下一轮的信息增益筛选。

5. 进入子集$D_{21}$($A_2$=是)，由于该子集中只有一个类别，所以形成叶子结点，结点的标记就是类别

6. 进入子集$D_{22}$($A_2$=否)，由于该子集中只有一个类别，所以形成叶子结点，结点的标记就是类别

  最终产生的决策树模型如下图所示

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

# coding=utf-8
from math import log
import pygraphviz as pgv


def calcShannonEnt(dataSet):
    # 数据集总条数
    numEntries = len(dataSet)
    labelCount = {}
    # 计算每个类别样本的数量
    for featVec in dataSet:
        currentLabel = featVec[-1]
        labelCount[currentLabel] = labelCount.get(currentLabel, 0) + 1
    shannonEnt = 0.0
    for key in labelCount:
        # 计算每个类别出现的概率
        prob = float(labelCount[key]) / numEntries
        # 计算信息熵
        shannonEnt -= prob * log(prob, 2)
    return shannonEnt


def splitDataSet(dataSet, axis, value):
    retDataSet = []
    for featVec in dataSet:
        # 挑选第axis个特征的值为value的数据
        if featVec[axis] == value:
            reducedFeature = featVec[:axis]
            reducedFeature.append(featVec[axis + 1:])
            # 注意，此处是extend，不是append
            retDataSet.extend(reducedFeature)
    return retDataSet


def chooseBestFeatureToSplit(dataSet):
    # 特征数量，由于最后一个是标签，所以减1
    numFeatures = len(dataSet[0]) - 1
    # 计算经验熵
    baseEntropy = calcShannonEnt(dataSet)
    bestInfoGain = 0.0
    bestFeature = -1
    for i in range(numFeatures):
        # 把第i个特征全部挑出
        featList = [example[i] for example in dataSet]
        # 第i个特征有几个不同的特征值
        uniqueVals = set(featList)
        # 条件经验熵
        newEntropy = 0.0
        for value in uniqueVals:
            # 把特征值等于value的数据挑出
            subDataSet = splitDataSet(dataSet, i, value)
            # 计算特征值等于value的数据占总样本的比例
            prob = len(subDataSet) / float(len(dataSet))
            # 计算经验条件熵
            newEntropy += prob * calcShannonEnt(subDataSet)
        # 计算信息增益
        infoGain = baseEntropy - newEntropy
        if infoGain > bestInfoGain:
            bestInfoGain = infoGain
            bestFeature = i
    return bestFeature


def majorityCnt(classList):
    classCount = {}
    for vote in classList:
        classCount[vote] = classCount.get(vote, 0) + 1
    sortedClassCount = sorted(classCount.iteritems(), key=lambda item: item[1], reverse=True)
    return sortedClassCount[0][0]


def createTree(dataSet, labels):
    classList = [example[-1] for example in dataSet]
    if classList.count(classList[0]) == len(classList):
        return classList[0]
    if len(dataSet[0]) == 1:
        return majorityCnt(classList)
    bestFeat = chooseBestFeatureToSplit(dataSet)
    bestFeatLabel = labels[bestFeat]
    myTree = {bestFeatLabel: {}}
    del (labels[bestFeat])
    featValues = [example[bestFeat] for example in dataSet]
    uniqueVals = set(featValues)
    for value in uniqueVals:
        subLabels = labels[:]
        myTree[bestFeatLabel][value] = createTree(splitDataSet(dataSet, bestFeat, value), subLabels)
    return myTree


def createDataSet():
    dataSet = [['1', '1', 'yes'],
               ['1', '1', 'yes'],
               ['1', '0', 'no'],
               ['0', '1', 'no'],
               ['0', '1', 'no']]
    labels = ['no surfacing', 'flippers']
    return dataSet, labels


def main():
    dataSet, labels = createDataSet()
    tree = createTree(dataSet, labels)
    print tree


if __name__ == '__main__':
    main()

参考

1. 机器学习实战，人民邮电出版社
2. 统计学习方法，清华大学出版社

  机器学习主要解决两类问题：回归和分类，上篇文章中的线性回归模型是回归模型的一种，但并不能用于分类问题，本文将介绍一种可以用于分类问题的回归算法–logistic回归。

  利用logistic回归进行分类的主要思想是：估算某个实例属于特定类别的概率，如果预估概率超过50%，则属于该类别；反之，属于其他类别，这样形成了一个简单的二元分类器。

原理

sigmoid函数

  我们希望能找到一种根据输入值输出0，1分类的函数，脑海中闪现的第一个函数应该就是单位阶跃函数，即
$$
\delta(x)=
\begin{cases}
0, & x \leq 0; \\
1, & x>0.
\end{cases}
$$
但是该函数在$x=0$处不可导，对于公式推导来说有很大问题，所以大牛们找到了sigmoid函数来代替阶跃函数。

  Sigmoid函数公式如下

$$
f(x)=\frac{1}{1+e^{-x}}
\tag{1}
$$

其函数图像如下所示

logistic回归模型

  将公式1推广到多维场景：给每个特征都乘以一个回归系数并相加，即$z=\theta_{0}+\theta_{1}x_{1}+\theta_{2}x_{2}+……+\theta_{n}x_{n}=\theta ^{T}x$（线性回归），然后将总和带入sigmoid函数中得到
$$
h_{\theta}\left(x \right)=\frac{1}{1+e^{-z}}=\frac{1}{1+e^{-\theta ^{T}x}}
\tag{2}
$$
  所以线性回归的的结果被映射到了sigmoid函数中。在sigmoid函数图像中可以看出，其函数值介于0～1，中间值为0.5，因此$h_{\theta}\left(x \right)$的输出可以看作分类的概率。当$h_{\theta}\left(x \right)>0.5$，说明$x$属于A类；反之，$h_{\theta}\left(x \right)<0.5$，则$x$属于B类。

参数估计

  对于给定的数据集$T={(x_1,y_1),(x_2,y_2),…(x_n,y_n)}$，其中$x_i \in R^n$，$y_i \in {0,1}$。分类结果为1的概率表示为
$$
P(y=1|x;\theta)=h_\theta(x)
$$
分类结果为0的概率表示为
$$
P(y=0|x;\theta)=1-h_\theta(x)
$$
将其合并则有
$$
p(y|x;\theta)=h_\theta(x_i) ^{y_i}(1-h_\theta(x))^{1-y_i}
$$
因为m个样本互相独立，所以其联合分布可以看作各个样本概率之积，则似然函数为
$$
L(\theta) = \prod_{i=1}^{n}p(y_i|x_i;\theta)=\prod_{i=1}^{n} [h_\theta(x_i)]^{y_i}[1-h_\theta(x_i)]^{1-y_i}
$$
取对数似然
$$
l(\theta) =\sum_{i=1}^{n} y_ilog[h_\theta(x_i)]+(1-y_i)log[1-h_\theta(x)]
\tag{3}
$$
当$l\left(\theta \right)$取最大值时，$\theta$的值为最终的解。求$l\left(\theta \right)$对$\theta$的导数

$$
\begin{split}
l^{‘}(\theta)&=\sum_{i=1}^{n} [\frac{y_i}{h_\theta(x_i)}-\frac{1-y_i}{1-h_\theta(x_i)}] \cdot h^{‘}_\theta(x_i)\\
&=\sum_{i=1}^{n} [\frac{y_i}{h_\theta(x_i)}-\frac{1-y_i}{1-h_\theta(x_i)}]\cdot[h_\theta(x_i)]^{‘}_{z} \cdot z^{‘}_\theta\\
&=\sum_{i=1}^{n}[\frac{y_i}{h_\theta(x_i)}-\frac{1-y_i}{1-h_\theta(x_i)}] \cdot h_\theta(x_i)[1-h_\theta(x_i)] \cdot z^{‘}_\theta \\
&=\sum_{i=1}^{n}[y_i-h_\theta(x_i)] \cdot x_i
\end{split}
\tag{4}
$$

* 注：sigmoid求导结果
$$
f’\left(x \right)=f\left(x \right)\left[1-f\left(x \right) \right]
$$

<\font>
对公式4进行迭代，可求得参数$\theta$。

思考

对偏置项b的处理

  在大部分资料中，会在输入$x$和权重$\theta$中增加一列常数列作为偏置项，这样做的好处是不需要对其进行单独处理，可随权重$\theta$一起求解，下面来讨论一下单独求解偏置项$b$的步骤。

  假设输入为$[x_1,x_2,…,x_n]$，权重为$[\theta_1,\theta_2,…,\theta_n]$，偏置项为$[b_1,b_2,…,b_n]$，则公式2变为
$$
h_{\theta}\left(x \right)=\frac{1}{1+e^{-z}}=\frac{1}{1+e^{-\theta ^{T}x-b}}
$$
接下来需让公示3对变量b求导
$$
\begin{split}
l^{‘}(b)&=\sum_{i=1}^{n} [\frac{y_i}{h_\theta(x_i)}-\frac{1-y_i}{1-h_\theta(x_i)}] \cdot h^{‘}_\theta(x_i)\\
&=\sum_{i=1}^{n} [\frac{y_i}{h_\theta(x_i)}-\frac{1-y_i}{1-h_\theta(x_i)}] \cdot [h_\theta(x_i)]^{‘}_{z} \cdot z^{‘}_b\\
&=\sum_{i=1}^{n}[\frac{y_i}{h_\theta(x_i)}-\frac{1-y_i}{1-h_\theta(x_i)}] \cdot h_\theta(x_i)[1-h_\theta(x_i)] \cdot z^{‘}_b \\
&=\sum_{i=1}^{n}[y_i-h_\theta(x_i)] \cdot x_i
\end{split}
\tag{5}
$$
对公式5进行迭代，可求得参数$b$。

代码实现

  为了加深对公示的理解，代码没有直接使用tensorflow内置的优化器，手动实现了逻辑回归梯度下降。示例如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

# coding=utf-8
import tensorflow as tf
from sklearn import datasets
import numpy as np
from sklearn.preprocessing import OneHotEncoder

iris = datasets.load_iris()
data_x = np.array(iris['data'])
data_y = np.array(iris['target']).reshape(-1, 1)
print data_x.shape
print data_y.shape
enc = OneHotEncoder()
enc.fit(data_y)
# 多分类问题，对label进行one-hot编码
targets = enc.transform(data_y).toarray()
# 输入行数不指定，列数为特征个数
X = tf.placeholder(dtype=tf.float32, shape=(None, data_x.shape[1]))
# 输出行数不指定，列数类别个数
y = tf.placeholder(dtype=tf.float32, shape=(None, 3))
# 学习率
alpha = 0.0001
# 迭代次数
epoch = 500
# 权重
theta = tf.Variable(tf.random_uniform([data_x.shape[1], 3], -1.0, 1.0), name='theta')
# 偏置
b = tf.Variable(tf.random_uniform([3], -1.0, 1.0), name='b')
# 预测值
predict_y = tf.nn.softmax(tf.matmul(X, theta) + b)
# 误差
error = tf.cast(y, tf.float32) - predict_y
# 代价函数，使用交叉熵函数，手动实现
cost = tf.reduce_mean(-tf.reduce_sum(y * tf.log(predict_y), reduction_indices=1))
# 权重的负梯度
theta_gradient = -tf.matmul(tf.matrix_transpose(X), error, name='theta_gradient')
# 偏置的负梯度
b_gradient = - tf.reduce_mean(tf.matmul(tf.transpose(X), error), reduction_indices=0, name='b_gradient')
# 权重迭代项
training_op1 = tf.assign(theta, theta - alpha * theta_gradient)
# 偏置迭代项
training_op2 = tf.assign(b, b - alpha * b_gradient)
# 全局初始化
init = tf.global_variables_initializer()
with tf.Session() as sess:
    init.run()
    for i in range(epoch):
        sess.run([training_op1, training_op2, cost], feed_dict={X: data_x, y: targets})
        print sess.run(theta, feed_dict={X: data_x, y: targets})
        print sess.run(b, feed_dict={X: data_x, y: targets})
        print sess.run(cost, feed_dict={X: data_x, y: targets})

代码详见逻辑回归代码

参考

1. Tensorflow实现梯度下降各种方法
2. 机器学习实战-基于Scikit-Learn和Tensorflow，机械工业出版社
3. 机器学习实战，人民邮电出版社

简介

  说到回归一般都指线性回归，回归的目的是预测数值型的目标值。线性回归模型的优点在于结果易于理解，计算简单，缺点在于对非线性的数据拟合不好。

原理

  线性回归就是对输入特征加权求和，在加上偏置项进行预测。公式如下所示：

$$
\widehat{h}=\theta_{0}+\theta_{1}x_{1}+\theta_{2}x_{2}+…+\theta_{n}x_{n}
\tag{1}
$$

其中，$\widehat{h}$代表预测值，$n$是特征的个数，${x_i}$代表第i个特征值，$\theta_{i}$ 代表第i个特征的权重，将上式向量化有

$$
\widehat{h}=h_{\theta}(X)=X\theta
\tag{2}
$$

其中$\theta$是模型的参数向量，X是实例的特征向量，$h_{\theta}$是使用模型参数$\theta$做参数的假设函数。

  为了衡量模型效果，一般使用MSE(均方误差，Mean Square Error)描述线性回归的预测值和真实值之间的偏差。假设输入的样本特征为$x_1,x_2,…,x_n$，对应的样本值为$y_1,y_2,y_3,…,y_n$，则MSE的定义为

$$
J(\theta)=\frac{1}{2}(\widehat{y}_i-y_i)^2=\frac{1}{2}\sum^n_{i=1}(h_{\theta}(x_i)-y_i)^2
\tag{3}
$$
写成向量形式有
$$
J(\theta)=\frac{1}{2}(X\theta-Y)^2
\tag{4}
$$
  以上就是线性回归的基本公式和损失函数。

求解

  我们的目标是找到合适的一组权重$\theta$让$J(\theta)$最小，这样就说明我们的模型越符合实际效果越好。一般有两种方式进行求解，直接求解和梯度下降。

解析求解

  我们已经拿到了损失函数的表达形式，可以从公式$(4)$中直接推到出$\theta$，过程如下。

$$
\begin{split}
J(\theta)=\frac{1}{2}(X\theta-Y)^2&=\frac{1}{2}(X\theta-Y)^T(X\theta-Y)\\
&=\frac{1}{2}(\theta^TX^T-Y^T)(X\theta-Y)\\
&=\frac{1}{2}(\theta^TX^TX\theta-\theta^TX^TY-Y^TX\theta-Y^TY)
\end{split}
\tag{5}
$$
  接下来对$J(\theta)$求偏导

$$
\begin{split}
\frac{\partial^{2}J(\theta)}{\partial\theta^{2}}&=\frac{1}{2}(2 \cdot X^TX\theta-2 \cdot X^TY)\\
&=X^TX\theta-X^TY
\end{split}
\tag{6}
$$

  令导数为0，有

$$\frac{\partial^{2}J(\theta)}{\partial\theta^{2}}=0$$

$$X^TX\theta-X^TY=0$$

$$X^TX\theta=X^TY$$

$$\theta=(X^TX)^{-1}X^TY \tag{7}$$

  最终推到出$(7)$，直接求解需要用到实例数据$X$及对应的预测值$Y$，在输入不大的情况下可以快速求得权重矩阵。但$(7)$中含有矩阵求逆运算，其算法复杂度通常在$O(n^2.4)$到$O(n^3)$之间，在高维大数据量的情况下直接求解的计算量很大；而且对于某些损失函数是凸函数的模型($X^TX$不可逆)来说，无法求得损失函数的偏导，该方法失效，只能使用梯度下降算法。

梯度下降

  梯度下降是一种通用的优化算法，梯度下降的中心思想就是迭代地调整参数从而使成本函数最小化。

  具体来说，首先使用一个随机的$\theta$值，即随机初始化，然后逐步改进，每次踏出一步，每一步都尝试降低成本函数，直到算法收敛到一个最小值。梯度下降中一个重要的参数是每一步的步长，该参数被称为学习率。如果学习率太低，算法需要经过大量迭代才能收敛，如果太高，有可能越过最小值甚至比上一步迭代的结果还大，导致模型无法收敛。

  对于线性回归，梯度下降的数学表达形式如下。假设第$i$次迭代的权重为$\theta_i$，学习率为$\eta$，第$i+1$次迭代的权重为
$$
\theta_{i+1}=\theta_i-\eta \cdot \frac{\partial^2J(\theta)}{\partial\theta^2}
\tag{8}
$$

实现

  下面使用tensorflow实现两种求解方式，以加州房价预测数据为例。

tensorflow实现解析求解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# coding=utf-8
import numpy as np
import tensorflow as tf
from sklearn.datasets import fetch_california_housing

n_epochs = 10000  # 把样本集数据学习1000次
learning_rate = 0.01  # 步长 学习率 不能太大 太大容易来回震荡 太小 耗时间，跳不出局部最优解
# 可以写learn_rate动态变化，随着迭代次数越来越大 ，学习率越来越小 learning_rate/n_epoches
housing = fetch_california_housing()
m, n = housing.data.shape
housing_data_plus_bias = np.c_[np.ones((m, 1)), housing.data]
housing_data_target = housing.target.reshape(-1, 1)
X = tf.constant(housing_data_plus_bias, tf.float32)
y = tf.constant(housing_data_target, tf.float32)
XT = tf.transpose(X)
theta = tf.matmul(tf.matmul(tf.matrix_inverse(tf.matmul(XT, X)), XT), y)
with tf.Session() as sess:
    theta_value = theta.eval()
    print theta_value

tensorflow实现梯度下降

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

# coding=utf-8
import tensorflow as tf
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.preprocessing import StandardScaler

n_epochs = 10000  # 把样本集数据学习1000次
learning_rate = 0.01  # 步长 学习率 不能太大 太大容易来回震荡 太小 耗时间，跳不出局部最优解
# 可以写learn_rate动态变化，随着迭代次数越来越大 ，学习率越来越小 learning_rate/n_epoches
housing = fetch_california_housing()
m, n = housing.data.shape
housing_data_plus_bias = np.c_[np.ones((m, 1)), housing.data]
# 可以使用TensorFlow或者Numpy或者sklearn的StandardScaler去进行归一化
# 归一化可以最快的找到最优解
# 常用的归一化方式：
# 最大最小值归一化 (x-min)/(max-min)
# 方差归一化 x/方差
# 均值归一化 x-均值 结果有正有负 可以使调整时的速度越来越快。
scaler = StandardScaler().fit(housing_data_plus_bias)  # 创建一个归一化对象
scaled_housing_data_plus_bias = scaler.transform(housing_data_plus_bias)  # 真正执行 因为来源于sklearn所以会直接执行，不会延迟。
housing_data_target = housing.target.reshape(-1, 1)

X = tf.placeholder(tf.float32, name='X')
y = tf.placeholder(tf.float32, name='y')

# random_uniform函数创建图里一个节点包含随机数值，给定它的形状和取值范围，就像numpy里面rand()函数
theta = tf.Variable(tf.random_uniform([n + 1, 1], -1.0, 1.0), name='theta')  #参数\theta，列向量，按照-1.0到1.0随机给
y_pred = tf.matmul(X, theta, name="predictions")  # 相乘 m行一列
error = y_pred - y  # 列向量和列向量相减 是一组数
mse = tf.reduce_mean(tf.square(error), name="mse")  # 误差平方加和，公式（3）的手动实现
gradients = 2.0 / m * tf.matmul(tf.transpose(X), error)  # 梯度公式
training_op = tf.assign(theta, theta - learning_rate * gradients)  # 即公式(8)，需对该式进行迭代

init = tf.global_variables_initializer()

with tf.Session() as sess:
    sess.run(init)  # 初始化

    for epoch in range(n_epochs):  # 迭代1000次
        if epoch % 100 == 0:
            print("Epoch", epoch, "MSE = ",
                  sess.run(mse, feed_dict={X: scaled_housing_data_plus_bias, y: housing_data_target}))  # 每运行100次的时候输出
        sess.run(training_op, feed_dict={X: scaled_housing_data_plus_bias, y: housing_data_target})

    best_theta = theta.eval()  # 最后的w参数值
    print(best_theta)

以上代码详见线性回归tensorflow实现

#参考

1. 【TensorFlow篇】–Tensorflow框架初始，实现机器学习中多元线性回归
2. 机器学习实战-基于Scikit-Learn和Tensorflow，机械工业出版社
3. 机器学习实战，人民邮电出版社

Topology创建过程

  Topology是storm的一个完整工作流，由spout、bolt等组件构成。下面我们来看一下Topology是如何被创建的。

入口函数

  我们一般会在storm的入口函数调用TopologyBuilder进行Topology的创建，如下所示

1
2
3
4
5
6
7

TopologyBuilder builder = new TopologyBuilder();
//设置spout
builder.setSpout(SENTENCE_SPOUT_ID, spout, 2);
//设置bolt
builder.setBolt(SPLIT_BOLT_ID, splitBolt, 2).setNumTasks(4).shuffleGrouping(SENTENCE_SPOUT_ID);
builder.setBolt(COUNT_BOLT_ID, countBolt, 4).fieldsGrouping(SPLIT_BOLT_ID, new Fields("word"));
builder.setBolt(REPORT_BOLT_ID, reportBolt).globalGrouping(COUNT_BOLT_ID);

TopolgyBuilder

  比较重要的实例变量

1
2
3
4
5
6

//所有提交的bolt放入_bolts中
private Map<String, IRichBolt> _bolts = new HashMap<>();
//所有提交的spout放入_spouts中
private Map<String, IRichSpout> _spouts = new HashMap<>();
//所有topology的spout和bolt放入_commons中
private Map<String, ComponentCommon> _commons = new HashMap<>();

setSpout

1
2
3
4
5
6
7
8
9

public SpoutDeclarer setSpout(String id, IRichSpout spout, Number parallelism_hint) throws IllegalArgumentException {
    //检测输入id是不是唯一的，主要是从实例变量的map里看有没有对应的key存在
    validateUnusedId(id);
    //构建ComponentCommon对象并进行初始化，最后放入_commons中
    initCommon(id, spout, parallelism_hint);
    //放入_spouts中
    _spouts.put(id, spout);
    return new SpoutGetter(id);
}

initCommon

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

   private void initCommon(String id, IComponent component, Number parallelism) throws IllegalArgumentException {
       ComponentCommon common = new ComponentCommon();
//设置消息流来源和分组方式
       common.set_inputs(new HashMap<GlobalStreamId, Grouping>());
//设置并行度
       if(parallelism!=null) {
           int dop = parallelism.intValue();
           if(dop < 1) {
               throw new IllegalArgumentException("Parallelism must be positive.");
           }
           common.set_parallelism_hint(dop);
       }
//设置组件的配置参数
       Map conf = component.getComponentConfiguration();
       if(conf!=null) common.set_json_conf(JSONValue.toJSONString(conf));
       _commons.put(id, common);
   }

setBolt

  与setSpout类似

1
2
3
4
5
6
7

public BoltDeclarer setBolt(String id, IRichBolt bolt, Number parallelism_hint) throws IllegalArgumentException {
    validateUnusedId(id);
    initCommon(id, bolt, parallelism_hint);
    _bolts.put(id, bolt);
    return new BoltGetter(id);
}

setSpout和setBolt区别

看上去二者完成的事情基本类似，但是返回值有区别

SpoutGetter的源码

1
2
3
4
5
6

protected class SpoutGetter extends ConfigGetter<SpoutDeclarer> implements SpoutDeclarer {
    public SpoutGetter(String id) {
        super(id);
    }        
}

而BoltGetter源码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

protected class BoltGetter extends ConfigGetter<BoltDeclarer> implements BoltDeclarer {
    private String _boltId;

    public BoltGetter(String boltId) {
        super(boltId);
        _boltId = boltId;
    }

    public BoltDeclarer fieldsGrouping(String componentId, Fields fields) {
        return fieldsGrouping(componentId, Utils.DEFAULT_STREAM_ID, fields);
    }

    public BoltDeclarer fieldsGrouping(String componentId, String streamId, Fields fields) {
        return grouping(componentId, streamId, Grouping.fields(fields.toList()));
    }

    public BoltDeclarer globalGrouping(String componentId) {
        return globalGrouping(componentId, Utils.DEFAULT_STREAM_ID);
    }

    public BoltDeclarer globalGrouping(String componentId, String streamId) {
        return grouping(componentId, streamId, Grouping.fields(new ArrayList<String>()));
    }

    public BoltDeclarer shuffleGrouping(String componentId) {
        return shuffleGrouping(componentId, Utils.DEFAULT_STREAM_ID);
    }

    public BoltDeclarer shuffleGrouping(String componentId, String streamId) {
        return grouping(componentId, streamId, Grouping.shuffle(new NullStruct()));
    }
	
    .....
}

可以发现BoltGetter还实现了不同的分组方式，如

1
2
3
4
5

private BoltDeclarer grouping(String componentId, String streamId, Grouping grouping) {
    _commons.get(_boltId).put_to_inputs(new GlobalStreamId(componentId, streamId), grouping);
    return this;
}

分组的本质是在_common中通过对应的boltId找到对应的ComponentCommon对象，对inputs属性进行设置。

createTopology()

  TopologyBuilder中还有一个比较重要的方法–createTopology()，其主要完成最后的封装工作。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

   public StormTopology createTopology() {
//bolt集合
       Map<String, Bolt> boltSpecs = new HashMap<>();
//spout集合
       Map<String, SpoutSpec> spoutSpecs = new HashMap<>();
       maybeAddCheckpointSpout();
       for(String boltId: _bolts.keySet()) {
    //通过boltId获取Bolt
           IRichBolt bolt = _bolts.get(boltId);
           bolt = maybeAddCheckpointTupleForwarder(bolt);
    //设置对应ComponentCommon对象的streams属性(输出的字段列表是否为直接流)
           ComponentCommon common = getComponentCommon(boltId, bolt);
           try{
               maybeAddCheckpointInputs(common);
	//把bolt和common一起放入bolt集合
               boltSpecs.put(boltId, new Bolt(ComponentObject.serialized_java(Utils.javaSerialize(bolt)), common));
           }catch(RuntimeException wrapperCause){
               if (wrapperCause.getCause() != null && NotSerializableException.class.equals(wrapperCause.getCause().getClass())){
                   throw new IllegalStateException(
                       "Bolt '" + boltId + "' contains a non-serializable field of type " + wrapperCause.getCause().getMessage() + ", " +
                       "which was instantiated prior to topology creation. " + wrapperCause.getCause().getMessage() + " " +
                       "should be instantiated within the prepare method of '" + boltId + " at the earliest.", wrapperCause);
               }
               throw wrapperCause;
           }
       }
//对spout的处理和bolt的处理基本一致
       for(String spoutId: _spouts.keySet()) {
           IRichSpout spout = _spouts.get(spoutId);
           ComponentCommon common = getComponentCommon(spoutId, spout);
           try{
               spoutSpecs.put(spoutId, new SpoutSpec(ComponentObject.serialized_java(Utils.javaSerialize(spout)), common));
           }catch(RuntimeException wrapperCause){
               if (wrapperCause.getCause() != null && NotSerializableException.class.equals(wrapperCause.getCause().getClass())){
                   throw new IllegalStateException(
                       "Spout '" + spoutId + "' contains a non-serializable field of type " + wrapperCause.getCause().getMessage() + ", " +
                       "which was instantiated prior to topology creation. " + wrapperCause.getCause().getMessage() + " " +
                       "should be instantiated within the prepare method of '" + spoutId + " at the earliest.", wrapperCause);
               }
               throw wrapperCause;
           }
       }

       StormTopology stormTopology = new StormTopology(spoutSpecs,
               boltSpecs,
               new HashMap<String, StateSpoutSpec>());

       stormTopology.set_worker_hooks(_workerHooks);

       return Utils.addVersions(stormTopology);

最终我们设置的bolt和spout都被封装到了StormTopology中。

总结

  总的来讲TopologyBuilder就是根据分组方式把spout和bolt节点连接起来形成一个拓扑结构。

并行度相关概念

  和spark driver、work、executor一样，storm也有一套自己的概念。Storm集群中的节点可以分为两类：主节点nimbus，从节点supervisor。nimbus主要负责分配计算资源，supervisor主要负责执行client提交的任务。supervisor节点在运行任务时，涉及到以下三个概念：worker、executor、task。

worker

  supervisor在执行任务时，会启动一个或者多个jvm进程，这些进程统称为worker。默认情况下一个supervisor最多可以启动4个worker，也可以通过修改storm.yaml参数进行修改。因为topology最终都是通过worker来执行的，所以nimbus并不关心有几个supervisor，只关心有多少个worker，至于worker分配在哪个节点上，nimbus是不关心的。需要说明的是，一个worker只能同时运行一个topology，如果在该worker运行时又有topology提交，nimbus会将其分配给其他空闲的worker。

executor

  在一个jvm进程中，有时候我们会开启多个线程执行任务，这里executor的概念就是线程的概念。

task

  executor中运行的topology的一个component，如spout、bolt，叫做task，task是storm中进行计算的最小运行单位，表示spout和bolt的运行实例。

资源设置的注意事项

worker

  worker个数的增加，会导致worker之间数据传输时间增加。如果程序瓶颈在于待处理的元组数据太多算力不足，那么通过可以通过增加worker个数提高计算效率。

executor

  executor是真正的并行度。executor初始数量=spout数量+bolt数量+acker数量，也就是task个数，默认一个executor对应一个task。其中spout、bolt、acker数量运行时是不会变化的，但是executor数量是可以变化的。

task

  task的存在是为了topology扩展的灵活性，与并行度无关。task在实际执行数据处理。如果单纯提高task个数，不增加executor个数，并不一定能提高性能。提高task任务数量，可以为后期进行弹性计算（rebalance）即后期动态调整某一组件的并行度。

并行度计算

计算并行度官网有个比较好的例子

上图中，有2个worker进程，

• 蓝色的BlueSpout有2个executor进程，每个executor有1个task，并行度为2；
• 绿色的GreenBolt有2个executor进程，每个executor有2个task，并行度为2；
• 黄色的YellowBolt有6个executor进程，每个executor有1个task，并行度为6；
• 上图总平行度是2+2+6=10，具体分配到每个worker上就是5个