logistic回归与最大熵模型（二）

前言

上一篇文章对Logistic回归和最大熵模型有了初步的理解，但是其中还有些一知半解的地方，经过代码实现，自己对相关内容有了充分的理解，所以写出本总结（二）。

两个模型的共同之处

模型都是输入$x$，输出就是对应各类别的概率分布$p(y|x)$，不同的是两个模型的形式不同。
Logistic二分类回归模型是：
$$\begin{array}{r}P(y|x)=\{\begin{array}{ll}\frac{\exp (wx+b)}{1+\exp (wx+b)}& \text{ if }y=1\\ \frac{1}{1+\exp (wx+b)}& \text{ if }y=0\end{array}\end{array}$$

Logistic多分类回归模型是：
$$\begin{array}{r}P(y|x)=\{\begin{array}{ll}\frac{\exp (w_{1}x+b_1)}{1+{\sum }_{i=1}^n\exp (w_{i}x+b_i)}& \text{ if }y=1\\ \frac{\exp (w_{2}x+b_2)}{1+{\sum }_{i=1}^n\exp (w_{i}x+b_i)}& \text{ if }y=2\\ & ⋮\\ \frac{\exp (w_{n}x+b_n)}{1+{\sum }_{i=1}^n\exp (w_{i}x+b_i)}& \text{ if }y=n\\ \frac{1}{1+{\sum }_{i=1}^n\exp (w_{i}x+b_i)}& \text{ if }y=0\end{array}\end{array}$$

而最大熵二分类模型则是：
$$\begin{array}{r}P(y|x)=\{\begin{array}{ll}\frac{\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=1))}{{\sum }_{j=0}^1\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=j))}& \text{ if }y=1\\ \frac{\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=0))}{{\sum }_{j=0}^1\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=j))}& \text{ if }y=0\end{array}\end{array}$$

最大熵多分类模型就是：
$$\begin{array}{r}P(y|x)=\{\begin{array}{ll}\frac{\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=1))}{{\sum }_{j=0}^m\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=j))}& \text{ if }y=1\\ \frac{\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=2))}{{\sum }_{j=0}^m\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=j))}& \text{ if }y=2\\ ⋮& \\ \frac{\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=0))}{{\sum }_{j=0}^m\exp ({\sum }_{i=1}^n{w}_i{f}_i(x,y=j))}& \text{ if }y=0\end{array}\end{array}$$

一般来说其中的$f_i(x,y)$都是选定的函数，中间参数不用修改。比如我们可以选择$f_i(x,y)=(w_{i}x+b_i)y$。这样上面的最大熵二分类模型就是，在下面的方程中只考虑了一个f函数：
$$\begin{array}{rl}P(y|x)=& \{\begin{array}{ll}\frac{\exp (w_{\xi }(w_{1}x+b_1)\ast 1)}{{\sum }_{j=0}^1\exp (w_{\xi }(w_{1}x+b_1)\ast j)}& \text{ if }y=1\\ \frac{\exp (w_{\xi }(w_{1}x+b_1)\ast 0)}{{\sum }_{j=0}^1\exp (w_{\xi }(w_{1}x+b_1)\ast j)}& \text{ if }y=0\end{array}\\ =& \{\begin{array}{ll}\frac{\exp (w_{\xi }(w_{1}x+b_1))}{1+\exp (w_{\xi }(w_{1}x+b_1))}& \text{ if }y=1\\ \frac{1}{1+\exp (w_{\xi }(w_{1}x+b_1))}& \text{ if }y=0\end{array}\end{array}$$

上面这个方程跟logistic二分类回归模型很相似，不同的是，这里只有一个可变参数$w_{\xi }$，其他的两个参数$w_1,b_1$都是我们在设计方程的时候就已经确定下来的，所以我们接下来直接比较logistic回归跟最大熵二分类单方程回归之间的差异。
上述模型不好用，我该一下看看效果
$$\begin{array}{rl}P(y|x)=& \{\begin{array}{ll}\frac{\exp (w_{\xi }(w_{1}x+b_1)\ast 1)}{{\sum }_{j=0}^1\exp (w_{\xi }(w_{1}x+b_1)\ast j)}& \text{ if }y=1\\ \frac{\exp (w_{\xi }(w_{1}x+b_1)\ast 0)}{{\sum }_{j=0}^1\exp (w_{\xi }(w_{1}x+b_1)\ast j)}& \text{ if }y=0\end{array}\\ =& \{\begin{array}{ll}\frac{\exp (w_{\xi }(w_{1}x+b_1+w_{2}y))}{\exp (w_{\xi }(w_{1}x+b_1))+\exp (w_{\xi }(w_{1}x+b_1+w_{2}y))}& \text{ if }y=1\\ \frac{\exp (w_{\xi }(w_{1}x+b_1))}{\exp (w_{\xi }(w_{1}x+b_1))+\exp (w_{\xi }(w_{1}x+b_1+w_{2}y))}& \text{ if }y=0\end{array}\end{array}$$

示例代码

代码都以二分类为例，其他多分类也很容易拓展。

Logistic二分类回归代码：

#导入模块
import torch
from matplotlib import pyplot as plt
torch.set_default_tensor_type(torch.cuda.FloatTensor)
# 生成数据
torch.manual_seed(1)
w_0=torch.randn(3)
b_0=torch.randn(1)
x_0=torch.randn(15,3)
y_0=torch.where(torch.mv(x_0,w_0)+b_0>0,1,-1)
print(f"w_0,b_0,x_0,y_0:{w_0,b_0,x_0,y_0}")
# 设计模型
lr=torch.tensor(0.001)
temp1=torch.where(y_0<0,1,0)#负一的位置
temp2=torch.where(y_0>=0,1,0)#正一的位置
w1=torch.randn(3,requires_grad=True)
b1=torch.randn(1,requires_grad=True)
p=torch.exp(torch.mv(x_0,w1)+b1)/(1+torch.exp(torch.mv(x_0,w1)+b1))#模型对于分布的预测
p=(1-p)*temp1+temp2*p#正一的位置保留，负一的位置用1修正，这样的话p中数据就是对应位置出现对应类别的概率。
prod=-torch.sum(torch.log(p))
#print(torch.log(p),torch.sum(torch.log(p)))
flag=prod.data+1
while True:
    p=torch.exp(torch.mv(x_0,w1)+b1)/(1+torch.exp(torch.mv(x_0,w1)+b1))
    p=(1-p)*temp1+temp2*p
    prod=-torch.sum(torch.log(p))
    prod.backward()
    w1.data=w1.data-lr*w1.grad.data
    b1.data=b1.data-lr*b1.grad.data
    print(f"w1,b1,prod:{w1,b1,prod}")
    if(flag-prod.data<0.001):
        print(flag)
        print(prod.data)
        break
    flag=prod.data
    w1.grad.data.zero_()
    b1.grad.data.zero_()
print(f"w1,b1,prod:{w1,b1,prod}")
print(torch.exp(torch.mv(x_0,w1)+b1)/(1+torch.exp(torch.mv(x_0,w1)+b1)))
print(y_0)
#创建会有不同id，运算也会有不同id，除了赋值，其他的都会产生不同id。同时只有.data是在id不变的情况下运算的。
#创建时候的id变化、运算时候的id变化。

此代码逻辑上没有问题，但是执行中会有问题，最后一个分类错误，就很奇怪。想不明白为什么。
但是下面的代码就没有问题，仅仅是改动了生成数据的分布从torch.randn到torch.rand，可能是因为randn生成的数据恰好在局部最优点附近，所以无法找到全局最优？我好像知道为什么了，上面的代码没有梯度归零。添加梯度归零以后就好了。

#导入模块
import torch
from matplotlib import pyplot as plt
torch.set_default_tensor_type(torch.cuda.FloatTensor)
# 生成数据
torch.manual_seed(1)
w_0=torch.randn(3)
b_0=torch.randn(1)
x_0=torch.randn(15,3)
y_0=torch.where(torch.mv(x_0,w_0)+b_0>0,1,-1)
print(f"w_0,b_0,x_0,y_0:{w_0,b_0,x_0,y_0}")
lr=torch.tensor(0.001)
temp1=torch.where(y_0<0,1,0)#负一的位置
temp2=torch.where(y_0>=0,1,0)#正一的位置
#w1=torch.randn(3,requires_grad=True)
#b1=torch.randn(1,requires_grad=True)
w1=torch.rand(3)
w1.requires_grad_(True)
b1=torch.rand(1)
b1.requires_grad_(True)
print(w1,b1)
p=torch.exp(torch.mv(x_0,w1)+b1)/(1+torch.exp(torch.mv(x_0,w1)+b1))#模型对于分布的预测
p=(1-p)*temp1+temp2*p#正一的位置保留，负一的位置用1修正，这样的话p中数据就是对应位置出现对应类别的概率。
prod=-torch.sum(torch.log(p))
#print(torch.log(p),torch.sum(torch.log(p)))
flag=prod.data+1
while True:
    p=torch.exp(torch.mv(x_0,w1)+b1)/(1+torch.exp(torch.mv(x_0,w1)+b1))
    p=(1-p)*temp1+temp2*p
    prod=-torch.sum(torch.log(p))
    prod.backward()
    w1.data=w1.data-lr*w1.grad.data
    b1.data=b1.data-lr*b1.grad.data
    print(f"w1,b1,prod:{w1,b1,prod}")
    if(flag-prod.data<0.001):
        print(flag)
        print(prod.data)
        break
    flag=prod.data
    w1.grad.data.zero_()
    b1.grad.data.zero_()
print(f"w1,b1,prod:{w1,b1,prod}")
print(torch.exp(torch.mv(x_0,w1)+b1)/(1+torch.exp(torch.mv(x_0,w1)+b1)))
print(y_0)
#w1,b1,prod:(tensor([-2.2503,  5.0948, -1.0124], requires_grad=True), tensor([1.4083], requires_grad=True), tensor(2.3109, grad_fn=<NegBackward0>))

最大熵模型分类代码示例

又遇到问题了，真的离谱。

#导入模块
import torch
from matplotlib import pyplot as plt
torch.set_default_tensor_type(torch.cuda.FloatTensor)
# 生成数据
torch.manual_seed(1)
w_0=torch.randn(3)
b_0=torch.randn(1)
x_0=torch.randn(15,3)
y_0=torch.where(torch.mv(x_0,w_0)+b_0>0,1,0)
print(f"w_0,b_0,x_0,y_0:{w_0,b_0,x_0,y_0}")
# 设计模型
lr=torch.tensor(0.1)
temp1=torch.where(y_0<1,1,0)#负一的位置
temp2=torch.where(y_0>0,1,0)#正一的位置
print(f'lr,temp1,temp2:{lr,temp1,temp2}')
# tensor([0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0]))
# tensor([1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1])
# tensor([0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0])
w_1=torch.randn(3)
b_1=torch.randn(1)
w_xi=-torch.rand(1)
w_xi.requires_grad_(True)
print(f'w_1,b_1,w_xi:{w_1,b_1,w_xi}')
p=torch.exp((torch.mv(x_0,w_1)+b_1)*w_xi)/(1+torch.exp((torch.mv(x_0,w_1)+b_1)*w_xi))#模型对于分布的预测
print(p)
p=(1-p)*temp1+temp2*p#正一的位置保留，负一的位置用1修正，这样的话p中数据就是对应位置出现对应类别的概率。
print(p)
prod=-torch.sum(torch.log(p))
print(prod)
#print(torch.log(p),torch.sum(torch.log(p)))
flag=prod.data+1
print(flag)
for i in range(500):
#while True:
    p=torch.exp((torch.mv(x_0,w_1)+b_1)*w_xi)/(1+torch.exp((torch.mv(x_0,w_1)+b_1)*w_xi))
    print(p)
    p=(1-p)*temp1+temp2*p
    print(p)
    prod=-torch.sum(torch.log(p))
    print(prod)
    prod.backward()
    w_xi.data=w_xi.data-lr*w_xi.grad.data
    #if(flag-prod.data<0.000001):
    #    print(flag)
    #    print(prod.data)
    #    break
    flag=prod.data
    w_xi.grad.data.zero_()
    print(f'w_xi,prod:{w_xi,prod}')
#print(f"w1,b1,prod:{w1,b1,prod}")
print(torch.exp((torch.mv(x_0,w_1)+b_1)*w_xi)/(1+torch.exp((torch.mv(x_0,w_1)+b_1)*w_xi)))
print(y_0)
#创建会有不同id，运算也会有不同id，除了赋值，其他的都会产生不同id。同时只有.data是在id不变的情况下运算的。
#创建时候的id变化、运算时候的id变化。

试一试同一个$x$有不同的类别的情况。

关于最大熵模型

要想最大熵模型比较接近真实分布，最关键的就是要条件概率的函数尽可能的保留较多的信息，也就是：
$$E_{\widehat{P(x)}P(y|x)}(f(x,y))=E_{\widehat{P(x,y)}}(f(x,y))$$

这些条件要包含尽可能多的有关条件概率$P(y|x)$的信息，我们可以试想一下那些函数比较有效？需要多少个函数？

• 首先我们要确保函数应该是相互垂直的，也就是不相关的，这样能确保我们确定的函数数量是最优的，没有冗余的。

• 除此以外，我们针对单个样本x的条件概率分布去考虑它需要多少的不相关函数方程才能拥有足够的信息量。

这个方程的数量应该与样本的label标签数量是呈正相关的，也就是标签越多，对应的方程也应该越多。
比如说y对应的分类标签有：$1,2,3,4,5,6,\cdots ,m$，一共有m种标签，那么我们就需要有至少m个函数来确保拥有足够的信息量，最简单的就是：
$$\begin{array}{rl}P(y=1|x)& =p_1\\ ⋮& \\ P(y=m|x)& =p_m\end{array}$$

如果我们采用取均值的形式的函数，需要多少个函数？比如：
$$\begin{array}{rl}{f}_1(x,y=1)P(y=1|x)+\cdots +f_1(x,y=m)P(y=m|x)& =c_1\\ ⋮& \\ f_t(x,y=1)P(y=1|x)+\cdots +f_t(x,y=m)P(y=m|x)& =c_t\end{array}$$

这里的关键问题在于我们需要多大的t呢？同时每个$f_i(x,y)$应该是什么样的函数呢？
首先各函数应该是不相关的，对于单变量$x$来说，不相关的函数泰勒展开式告诉我们了有：$1,x,x^2,\cdots ,x^n,\cdots$
那么对于多变量来说呢？不相关的函数有哪些呢？多元泰勒展开式同样告诉了我们：$1,x_1,\cdots ,x_i,\cdots ,x_n,x_1^2,x_1{x}_2,\cdots ,x_i{x}_j,\cdots ,x_n^2,x_1^3,x_1^2{x}_2,\cdots ,x_i{x}_j{x}_k,\cdots ,x_n^3,\cdots$
足足有：$1+n+n^2+n^3+\cdots$个函数，那么我们需要多少个才可以充分表示上面的信息呢？
我们可以把上面的m个未知条件概率分布看作未知变量。所以我们只要建立一个$m\ast m$的方程应该就够了。并且他们需要是不相关的，也就是不同的方程说的是不同的事情，不能不同的方程说的同样的事情。
比如如果对于任何样本$x_1=x_2$，那么我们的两个方程：$f_1(x_1,y=1)P(y=1|x)+\cdots +f_1(x_1,y=m)P(y=m|x)=c_1$
跟另一个方程：$f_1(x_2,y=1)P(y=1|x)+\cdots +f_1(x_2,y=m)P(y=m|x)=c_1$其实就是一个方程，所以我们除了函数本身的不相关性以外，我们还需要考虑特征x的各个分量之间的相关性。所以为了解决这个问题，有两种方法：

• 提前对样本特征进行主成分分析，保留不相关的部分。这样我们可以选择前m个方程就可以解出条件概率分布
• 不对样本特征进行主成分分析，而是根据样本数据利用假设检验的方法检验各函数之间是否相关。

针对第二种方法，方法如下：
比如：假设$x_1,x_2$是相关的，然后我们利用统计量分析验证是否相关，例如$\text{COV}(x_1,x_2)$。或者$\text{COV}(x_i{x}_j{x}_k,x_p)$，不过这种方法需要不停寻找直到找到足够数量的不相关函数，时间消耗大。