Logistic回归代价函数的数学推导及实现

logistic回归的代价函数形式如下：

J(θ)=−1m[∑i=1my(i)loghθ(x(i))+(1−y(i))log(1−hθ(x(i)))]

可是这又是怎么来的呢？ 答：最大似然估计计算出来的

1.最大似然估计

我们先来简单的回顾一下最大似然估计(Maximum likelihood estimation),详细戳此处,见参数估计

所谓参数估计就是：对未知参数θ进行估计时，在参数可能的取值范围内选取，使“样本获得此观测值x1,x2...,xn”的概率最大的参数θ^作为θ的估计，这样选定的θ^有利于x1,x2...,xn”的出现。也就是说在已知数据集（结果）和模型（分布函数）的情况下，估计出最适合该模型的参数。

举个例子：

假设你有一枚硬币，随机抛10次；现在的结果是6次正面。我们都知道，抛一枚硬币，正面朝上和反面朝上的概率均是θ=0.5；但前提时，这是在大量的实验（抛硬币）情况下才有的结论。那在我们这个情况下，参数θ到底取何值时才能使得出现6次正面的肯能性最大呢？

我们知道，抛硬币是符合二项分布B(n,p)，也就是说我们现在已知样本结果以及函数分布，估计出使得该结果最大可能出现的参数θ^。则有：

L=P(X=6)=C610θ^6(1−θ^)4

而我们接下来要做的就是求当L取最大值时，θ^的值。我们很容易求得当θ^=0.6时L取得最大值0.25；而当θ^=0.5时，L=0.21

再假设你有一枚硬币，随机抛10次；现在的结果是7次正面。则此时使得该结果最大可能性出现参数θ^又是多少呢？按照上面的方法我们很容易求得当θ^=0.7时可能性最大。

再举个例子：

明显，在Logistic回归中，所有样本点也服从二项分布；设有x1,x2,x3三个样本点，其类标为1,1,0；同时设样本点为1的概率为P=hθ(x)，那么当P等于多少时，其结果才最可能出现1,1,0呢？于是问题就变成最大化：

P∗P(1−P)=hθ(x1)∗hθ(x2)∗(1−hθ(x3))

而这就是最大似然估计（求解参数估计的一个方法）

2.最大似然估计的数学定义

• 设总体X是离散型，其概率分布为P{X=x}=p(x;θ),θ为未知参数，X1,X2,...,Xn为X的一个样本，则X1,X2,...,Xn 取值为x1,...,xn的概率是：
  
  P{X1=x1,...,Xn=xn}=∏i=1nP{Xi=xi}=∏i=1np{xi;θ}

显然这个概率值是θ的函数，将其记为

L(θ)=L(x1,...,xn;θ)=∏i=1np{xi;θ}

称L(θ)为样本(x1,...,xn)的似然函数.
若θ^使得

L(x1,...,xn;θ^)=maxL(x1,...,xn;θ)

则称θ^=θ^(x1,...,xn)为未知参数θ的最大似然估计值.

求解步骤

a. 写出似然函数；

b. 方程两边同时取ln的对数；

c. 令∂lnL∂θi=0，求得参数

3.Logistic回归代价函的推导

设:

P(y=1|x;θ)=hθ(x)P(y=0|x;θ)=1−hθ(x)h(x)=g(z)=g(θTx)

将两者合并到一起为:

p(y|x;θ)=(hθ(x))y(1−hθ(x))1−y(3.1)

因此我们可以得到似然函数：

L(θ)=∏i=1mp(y(i)|x(i);θ)=∏i=1m(hθ(x(i)))y(i)(1−hθ(x))1−y(i)(3.2)

于是两边同时取自然对数得：

l(θ)=logL(θ)=∑i=1m[y(i)logh(x(i))+(1−y(i))log(1−h(x(i)))](3.3)

插播一个公式：

logabcd=logab+logcd=bloga+dlogc

且最大化(3.3)等价于最小化−1ml(θ)

下面我们梯度下降算法进行求解，所以就要对所有的参数求偏导，按照惯例，我们先考虑一个样本点(m=1)的情况，然后在矢量化：

∂l∂θj∂l∂b=−∂∂θjylog(hθ(x))+(1−y)log(1−hθ(x))=−∂∂θjylog(g(z))+(1−y)log(1−g(z))=−yg(z)∂g(z)∂θj−(1−y)(1−g(z))∂g(z)∂θj=−(yg(θTx)−(1−y)(1−g(θTx)))∂g(θTx)∂θTx∂θTx∂θj=−(y−g(z))xj=−(yg(θTx)−(1−y)(1−g(θTx)))∂g(θTx)∂θTx∂θTx∂b=−(y−g(z))∗1(3.4)(3.5)

这儿用的是g(z)=1/(1+exp(−z))，同样还可以是tanh(z),但是选取前者主要是因为其求导之后的结果为g′(z)=g(z)(1−g(z))便于化简。

于是根据梯度下降算法有：

θj=θj−1mα(h(x)−y)⋅xjb=b−1mα∑(h(x)−y)(3.6)(3.7)

注意，此时的x,y均为n个样本点的情况，与(3.4)中不同。xj是指，每个样本点的第j维所组成的列向量。

(3.6)矢量化后为:

θ=θ−1mαxT⋅(h(x)−y)

“⋅”表示点乘

4. python 实现

def gradDescent(X, y, W, b, alpha):    maxIteration = 500    for i in range(maxIteration):        z = X*W+b        error = sigmoid(z) - y        W = W - (1.0/m)*alpha*X.T*error        b = b - (1.0/m)*alpha*np.sum(error)    return W,b

源码