从伯努利分布、极大似然估计到逻辑回归

这篇文章也是想弄清楚一些自己没搞清楚的概念，也希望能帮助到同样有需求的朋友。

所有的统计学书籍都会从抛硬币试验开始讲起：随机变量 X=1 表示抛硬币正面朝上，设正面朝上的概率为 p ，那么随机变量 X 的 概率质量函数 （probability mass function, PMF）是：

P\left( x \right)=\left\{ \begin{align} & p,x=1 \\ & 1-p,x=0 \\ \end{align} \right.

这就是伯努利分布（Bernoulli distribution）概率质量函数的表达式了，需要指出的是，我们一般用大写字母 P 来表示概率质量函数，而用小写字母 p 或者 f 来表示概率密度函数（不是指这个式子中的 p ）。

这是很直观的有木有？在此基础上就可以引入二项分布的概念了：

对于 n 重伯努利试验，每次试验成功的概率为 p ，失败的概率为 1-p 。进行 n 次这样的试验，成功次数为 x ，失败次数为 n-x 的概率为：

P\left( x \right)=C_{n}^{x}{{p}^{x}}{{\left( 1-p \right)}^{n-x}}=\frac{n!}{x!\left( n-x \right)!}{{p}^{x}}{{\left( 1-p \right)}^{n-x}}

这就是二项分布（Binomial distribution）的概率质量函数。这个也很好理解有木有， n=1 时对应的就是伯努利分布。所以，我们可以将伯努利分布的概率质量函数写成：

P\left( x \right)={{p}^{x}}{{\left( 1-p \right)}^{1-x}}

其实和上面第一个式子是等价的有木有，但我们写成这种形式就可以保证和二项分布的统一。其实我们完全也可以把伯努利分布的概率质量函数写成其他的形式，例如：

P\left( x \right)=px+\left( 1-p \right)\left( 1-x \right)

但是这样并不能揭示其和二项分布的本质关联，所以我们通常将伯努利分布写成上面那种形式。接下来就可以推导对其参数 p 的极大似然估计了：

对于 极大似然估计 的本质，其实已经有网友讲得很清楚了：

极大似然估计是在极大似然原理的基础上衍生出来的一个统计方法。该方法用于通过观察试验数据来评估模型参数，即：“ 模型已定，参数未知 ”。

这里请大家参考： https://www. jianshu.com/p/48e2c0d4f dd8

从本质来讲，极大似然估计的过程是：通过若干次试验，观察结果并记录，利用试验的结果来推测某个参数，使得 样本出现的概率（这里体现为似然函数） 为最大，则称为极大似然估计。

其他网友也有一些更形象的解释：

https:// blog.csdn.net/zouxy09/a rticle/details/8537620

简单来讲，假设我们进行了 n 次扔硬币的实验，得到的样本集为：

D=\left\{ {{x}_{1}},{{x}_{2}},\ldots ,{{x}_{n}} \right\}

例如 n=10 ，而我们得到的结果是：(1,0,1,0,0,0,1,0,0,0)，由于每一次试验是独立的，也就是说这些样本是独立同分布的，那么这些样本同时出现的概率就是这些样本单独出现的概率的乘积，这就是似然函数（likehood function）:

L\left( p\left| {{x}_{1}},{{x}_{2}},\ldots ,{{x}_{n}} \right. \right)=\prod\limits_{i=1}^{n}{P\left( {{x}_{i}}\left| p \right. \right)}=\prod\limits_{i=1}^{n}{{{p}^{{{x}_{i}}}}{{\left( 1-p \right)}^{1-{{x}_{i}}}}}

这是很好理解的有木有。本质上来讲，就是在已知样本的情况下，计算这些样本同时出现的概率，这个概率是我们需要估计的参数的函数，如果这个函数能取得最大值，那么最大值情况下的参数则最有可能是真实的参数（或者最接近真实参数的情况，因为这些样本是已经真实发生的试验结果）。

由于上式是连乘运算，我们一般是取对数进行计算（取对数的一个原因是：容易超出计算机的浮点范围，造成下溢。而且取对数相当于把连乘转化成累加，更加方便于求偏导计算极值。）所以对数似然函数是：

L=\log \prod\limits_{i=1}^{n}{P\left( {{x}_{i}}\left| p \right. \right)}=\sum\limits_{i=1}^{n}{\log }P\left( {{x}_{i}}\left| p \right. \right)=\sum\limits_{i=1}^{n}{\left[ {{x}_{i}}\log p+\left( 1-{{x}_{i}} \right)\log \left( 1-p \right) \right]}

所以，要得到这个函数的极值，就应该在其对 p 导数为0的情况下得到：

\hat{p}=\arg \underset{p}{\mathop{\max }}\,L\left( p\left| X \right. \right)=\arg \underset{p}{\mathop{\max }}\,\sum\limits_{i=1}^{n}{\left[ {{x}_{i}}\log p+\left( 1-{{x}_{i}} \right)\log \left( 1-p \right) \right]}

\frac{\partial L}{\partial p}=\sum\limits_{i=1}^{n}{[\frac{{{x}_{i}}}{p}+\frac{1-{{x}_{i}}}{p-1}]=}\sum\limits_{i=1}^{n}{\frac{p-{{x}_{i}}}{p\left( p-1 \right)}}=0

那么，由于要求分子为0，我们可以得到：

\sum\limits_{i=1}^{n}{\left( p-{{x}_{i}} \right)}=0\Rightarrow p=\frac{1}{n}\sum\limits_{i=1}^{n}{{{x}_{i}}} ，具体到我们之前举的10次扔硬币的例子，可以很容易得到 p=0.3 ，虽然结论很显而易见，但是这里我们经过了严格的证明，并且在此基础上可以扩展到更为复杂的情况，例如，我们可以用同样的过程来计算逻辑回归的参数：

逻辑回归基于sigmoid函数（或者称为logistic函数的特殊形式）：

p\left( x \right)=\frac{1}{1+{{e}^{-x}}}

简单而言，逻辑回归就是计算：

P\left( Y=1\left| x \right. \right)=\frac{{{e}^{wx+b}}}{1+{{e}^{wx+b}}}

假设 P\left( Y=1\left| x \right. \right)\text{=}\pi \left( x \right) ，那么 P\left( Y=0\left| x \right. \right)=1-\pi \left( x \right) ，那么和上面的推导类似，逻辑回归的似然函数可以写成：

\prod\limits_{i=1}^{n}{{{\left[ \pi \left( {{x}_{i}} \right) \right]}^{{{y}_{i}}}}{{\left[ 1-\pi \left( {{x}_{i}} \right) \right]}^{1-{{y}_{i}}}}}

那么对数似然函数为：

\begin{align} & L\left( w \right)=\sum\limits_{i=1}^{n}{\left[ {{y}_{i}}\log \pi \left( {{x}_{i}} \right)+\left( 1-{{y}_{i}} \right)\log \left( 1-\pi \left( {{x}_{i}} \right) \right) \right]} \\ & =\sum\limits_{i=1}^{n}{\left[ {{y}_{i}}\log \frac{\pi \left( {{x}_{i}} \right)}{1-\pi \left( {{x}_{i}} \right)}+\log \left( 1-\pi \left( {{x}_{i}} \right) \right) \right]} \\ & =\sum\limits_{i=1}^{n}{\left[ {{y}_{i}}\left( w\cdot {{x}_{i}}+b \right)-\log \left( 1+{{e}^{w\cdot {{x}_{i}}+b}} \right) \right]} \\ \end{align}

在此基础上对 w 和 b 求偏导，再利用随机梯度下降，就可以训练出 w 和 b 了。

下面给出在PyTorch上的Logistic Regression实现，这里主要参考了廖星宇的《深度学习入门之PyTorch》上的数据和示例代码： L1aoXingyu/code-of-learn-deep-learning-with-pytorch

不过不知道为啥作者把原来书里的代码改掉了，我参考其他朋友的博客把原来书里的代码整理了一下，网上的代码不知道怎么回事总有些小问题：

import torch
from torch.autograd import Variable
from torch import nn
import numpy as np
import matplotlib.pyplot as plt
def get_parameter_number(net):#返回对网络参数的统计
    total_num = sum(p.numel() for p in net.parameters())
    trainable_num = sum(p.numel() for p in net.parameters() if p.requires_grad)
    return {'Total': total_num, 'Trainable': trainable_num}
with open('./data.txt', 'r') as f:
    data_list = [i.split('\n')[0].split(',') for i in f.readlines()]
    data = [(float(i[0]), float(i[1]), float(i[2])) for i in data_list]
# 标准化
x0_max = max([i[0] for i in data])
x1_max = max([i[1] for i in data])
data = [(i[0]/x0_max, i[1]/x1_max, i[2]) for i in data]
x0 = list(filter(lambda x: x[-1] == 0.0, data)) # 选择第一类的点
x1 = list(filter(lambda x: x[-1] == 1.0, data)) # 选择第二类的点
plot_x0 = [i[0] for i in x0]
plot_y0 = [i[1] for i in x0]
plot_x1 = [i[0] for i in x1]
plot_y1 = [i[1] for i in x1]
plt.plot(plot_x0, plot_y0, 'ro', label='x_0')
plt.plot(plot_x1, plot_y1, 'bo', label='x_1')
plt.legend(loc='best')
# plt.show()
np_data = np.array(data, dtype='float32') # 转换成 numpy array
x = torch.from_numpy(np_data[:, 0:2]) # 转换成 Tensor, 大小是 [100, 2]
y = torch.from_numpy(np_data[:, -1]).unsqueeze(1) # 转换成 Tensor，大小是 [100, 1]
class LogisticRegression(nn.Module):
    def __init__(self):
        super(LogisticRegression, self).__init__()
        self.lr = nn.Linear(2, 1)
        self.sm = nn.Sigmoid()
    def forward(self, x):
        x = self.lr(x)
        x = self.sm(x)
        return x
logistic_model = LogisticRegression()
if torch.cuda.is_available():
    logistic_model.cuda()
# 定义损失函数和优化器
criterion = nn.BCELoss()
optimizer = torch.optim.SGD(logistic_model.parameters(), lr=1e-3, momentum=0.9)
count=get_parameter_number(logistic_model)
print(count)
# 开始训练
for epoch in range(30000):
    if torch.cuda.is_available():
        x_data = Variable(x).cuda()
        y_data = Variable(y).cuda()
    else:
        x_data = Variable(x)
        y_data = Variable(y)
    out = logistic_model(x_data)
    loss = criterion(out, y_data)
    print_loss = loss.data.item()
    mask = out.ge(0.5).float()  # 以0.5为阈值进行分类
    correct = (mask == y_data).sum()  # 计算正确预测的样本个数
    acc = correct.item() / x_data.size(0)  # 计算精度
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
    # 每隔200轮打印一下当前的误差和精度
    if (epoch + 1) % 200 == 0:
        print('*'*10)
        print('epoch {}'.format(epoch+1)) # 训练轮数
        print('loss is {:.4f}'.format(print_loss))  # 误差
        print('acc is {:.4f}'.format(acc))  # 精度
w0, w1 = logistic_model.lr.weight[0]
b0 =logistic_model.lr.bias.item()
w0 = w0.item()
w1 = w1.item()
plot_x = np.arange(0.2, 1, 0.01)
plot_y = (-w0 * plot_x - b0) / w1