《统计学习方法》1——逻辑斯蒂回归

1. 模型

二项逻辑斯蒂回归模型是由逻辑斯蒂分布的条件概率分布表示的分类模型。
逻辑斯蒂分布函数为

其中，u是位置参数，是中心点（对称点）的位置；r是形状函数。
其分布图形是一条S形曲线，也就是sigmoid曲线,S形曲线的范围是(0,1)。

二项逻辑斯蒂回归模型的输出只有0和1两个值，即一件事情发生或者不发生。定义基于x，事件发生的概率分布如下：

由于是二值的，当事件发生的概率为p，那么不发生的概率为1-p，因此

上面两个概率分布为二项逻辑斯蒂回归模型的分布。X是输入，Y是输出，只有0和1两个取值， w.x是w和x的内积，即对应分量相乘后相加。其中，P(Y=1|x)的分子和分母同除 ，就是sigmoid函数的原型，因此，P(Y=1|x)是关于w.x+b的sigmoid函数，相当于u=, r = -b。

为了方便，将w和x做扩充，将偏执统一进w和x中，扩充后，, 。得到

其中，w.x近似于一个线性回归。

2. 模型参数估计

对于训练数据集，其中，，y的取值为0或者1。令。
概率的对数似然函数为：

将对数似然函数作为目标函数，其梯度为：

对于批梯度下降，以上梯度可以直接于梯度下降，对于随机梯度下降，以上N可为1，且样本随机抽取。

根据以上理论,程序实现如下：

class logistic_regression(object):
    def __init__(self, init_w_b,
                 input_size,
                 input_data_path,
                 input_label_path,
                 output_path,
                 init_learningRate = 1e-3,
                 k = 8,
                 iter = 20000,
                 batch_size=8):
        self.w_b = init_w_b
        self.input_size = input_size
        self.init_learningRate = init_learningRate
        self.batch_size = batch_size
        self.one_vector = np.array([1] * batch_size)
        self.input_data_path = input_data_path
        self.input_label_path = input_label_path
        self.K = k   #每k份数据抽取一份验证数据。
        self.output_path = output_path
        self.iter = iter

    """***************************
                exp(w.x+b)
     p(Y=1|X) = ---------------
                 1+exp(w.x+b)
     w.x的sigmoid函数
    ***************************"""
    def sigmoid(self, x):
        w_x = np.dot(self.w_b, x.transpose())
        exp_w_x = np.exp(w_x)
        p = exp_w_x / (self.one_vector + exp_w_x)
        return p

    """*****************************************
    L = y*(w.x)-log(1+exp(w.x)
    *****************************************"""
    def Log_likelihood(self, p, y):
        one_y = np.array([1] * self.batch_size, dtype=float) - y
        one_p = np.array([1] * self.batch_size) - p
        log_p = np.log(p)
        log_one_p = np.log(one_p)
        L_array = np.multiply(y, log_p) + np.multiply(one_y, log_one_p)

        L = 0.0
        for l in L_array:
            L = L + l
        L = L / self.batch_size
        return L

    """****************************************
    dL                xj exp(w.x)
    --   =  y*xj - ---------------- = xj(y-p)
    dwj              1 + exp(w.x)
    其中，w.x表示w和x的内积。
    ***************************************"""
    def Jaccobian(self, x, y, p):
        delta_w = np.dot(x.transpose(), (y - p))
        delta_w /= self.batch_size
        """
        print p
        print y
        print delta_w
        """
        return delta_w


    def optimizer(self,delta_w, learning_rate):
        self.w_b -= learning_rate * delta_w
        #print self.w_b


    def accurancy(self, p, y):
        sum_correct = 0.0
        for i in range(self.batch_size):
            y_ = 0
            if p[i] > 0.6:
                y_ = 1
            if y_ == y[i]:
                sum_correct += 1.0
        sum_correct /= self.batch_size
        return sum_correct

    def split_data(self):
        train_x = []
        test_x = []
        train_y = []
        test_y = []
        x_lines = []
        y_lines = []
        with open(self.input_data_path, 'r') as f1:
            x_lines = f1.readlines()
        with open(self.input_label_path, 'r') as f1:
            y_lines = f1.readlines()
        counter = 0
        test_num = 0
        length = min(len(x_lines), len(y_lines))

        for i in range(length):
            Xis = x_lines[i].split(', ')
            yis = y_lines[i].strip()
            if counter % self.K == 0:
                test_num = random.randint(0, self.K)
                counter = 0
            if counter == test_num:
                test_x.append([int(Xis[0][1:]), float(Xis[1][1:]), float(Xis[2][:-2]), 1.0])
                test_y.append(int(yis))
            else:
                train_x.append([int(Xis[0][1:]),  float(Xis[1][1:]), float(Xis[2][:-2]), 1.0])
                train_y.append(int(yis))
            counter += 1

            """数据归一化"""
        train_x = np.array(train_x)
        test_x = np.array(test_x)
        for i in range(self.input_size - 1):
            divided_max_value = 1 / max(train_x[:,i])
            train_x[:,i] *= divided_max_value
            test_x[:, i] *= divided_max_value
        train_x = list(train_x)

        return train_x, train_y, test_x, test_y

    def train(self, train_x, train_y):
        length = len(train_x)
        learning_rate = self.init_learningRate
        acc_sum = 0
        for j in range(self.iter):
            rand_index_list = random.sample(range(length), self.batch_size)
            train_x_batch = []
            train_y_batch = []
            for i in rand_index_list:
                train_x_batch.append(train_x[i])
                train_y_batch.append(train_y[i])
            train_x_array = np.array(train_x_batch, dtype=float)
            train_y_array = np.array(train_y_batch, dtype=float)
            #train_x_array = np.array([train_x_batch.append(train_x[i]) for i in rand_index_list])
            #train_y_array = np.array([train_y_batch.append(train_y[i]) for i in rand_index_list])

            p = self.sigmoid(train_x_array)
            delta_w = self.Jaccobian(train_x_array, train_y_array, p)
            self.optimizer(delta_w, learning_rate)

            L = self.Log_likelihood(p, train_y_array)
            acc = self.accurancy(p, train_y_array)
            acc_sum += acc * acc
            if j % 1000 == 0:
                if learning_rate > 1e-5:
                    learning_rate -= (0.1 * (j/1000))
                    acc_mean = math.sqrt(acc_sum) / 1000
                    acc_sum = 0
                print self.w_b
                print "第%d次迭代，对数似然值为%.3f, 准确度为%.3f"%(j, L, acc)
            #plt.plot(L, i, '', marker='o', ms=20, c='r', label="对数似然估计的变化曲线")


    def offline_test(self, test_x, test_y):
        length = len(test_x)
        rand_index_list = random.sample(range(length), self.batch_size)
        test_x_batch = []
        test_y_batch = []
        for i in rand_index_list:
            test_x_batch.append(train_x[i])
            test_y_batch.append(train_y[i])
        test_x_array = np.array(test_x_batch, dtype=float)
        test_y_array = np.array(test_y_batch, dtype=float)
        p = self.sigmoid(test_x_array)
        print test_x_array
        print test_y_array
        print p
        sum_correct = 0
        for i in range(self.batch_size):
            y_ = 0
            if p[i] > 0.6:
                y_ = 1
            if y_ == test_y_array[i]:
                sum_correct += 1
        sum_correct /= self.batch_size
        print sum_correct
        return sum_correct

3. 多项式扩展

在实践过程中，发现用特征的线性函数(w.x)没有办法很好的做逻辑拟合。在没有更好的特征可以替换的时候，可以考虑把特征做多项式扩展。Sklearn库中的PolynomialFeatures类，可以做多项式扩展，可将如(a,b,c)这样的3维特征，做2阶扩展，得到这样的10维特征。

    poly = PolynomialFeatures(2)
    train_x_extend = poly.fit_transform(train_x)
    test_x_extend = poly.fit_transform(test_x)

4. 正则化

多特征进行多项式扩展后，得到更多维的衍生特征，这样，过拟合的风险升高，因此，通常需要引入正则化。

clf = LogisticRegression(penalty='l2', solver='sag', max_iter= 2000)
clf.fit(train_x, train_y)

以上程序中，LogisticRegression的第一个参数penalty为正则方法的选择，solver为优化方法的选择。