waterMelon
Author Ron
本文主要讲述对数几率回归的原理以及代码实现,数据集用的西瓜数据集3.0
α
\alpha
α
1.算法原理
考虑单调可微函数
g
(
)
g()
g(),令
y
=
g
−
1
(
ω
T
x
+
b
)
(1)
y=g^{-1}(\omega^Tx+b)\tag{1}
y=g−1(ωTx+b)(1)这样得到的模型称为广义线性模型
显然单位阶跃函数不连续,因此我们采用"Sigmoid函数"
y
=
1
1
+
e
−
(
ω
T
x
+
b
)
(2)
y=\frac{1}{1+e^{-(\omega^Tx+b)}}\tag{2}
y=1+e−(ωTx+b)1(2)
作为替代函数,(2)式可变化为
ln
y
1
−
y
=
ω
T
x
+
b
(3)
\ln\frac{y}{1-y}=\omega^Tx+b\tag{3}
ln1−yy=ωTx+b(3)
将(3)中的y视为类后概率估计
p
(
y
=
1
∣
x
)
p(y=1|x)
p(y=1∣x),则有:
ln
p
(
y
=
1
∣
x
)
p
(
y
=
0
∣
x
)
=
ω
T
x
+
b
(4)
\ln\frac{p(y=1|x)}{p(y=0|x)}=\omega^Tx+b\tag{4}
lnp(y=0∣x)p(y=1∣x)=ωTx+b(4)
显然有
p
(
y
=
1
∣
x
)
=
e
ω
T
x
+
b
1
+
e
ω
T
x
+
b
(5)
p(y=1|x)=\frac{e^{\omega^Tx+b}}{1+e^{\omega^Tx+b}}\tag{5}
p(y=1∣x)=1+eωTx+beωTx+b(5)
p
(
y
=
0
∣
x
)
=
1
1
+
e
ω
T
x
+
b
(6)
p(y=0|x)=\frac{1}{1+e^{\omega^Tx+b}}\tag{6}
p(y=0∣x)=1+eωTx+b1(6)
于是我们可通过极大似然法估计
ω
\omega
ω和
β
\beta
β,即令每个样本属于其真实标记的概率越大越好。令
β
=
(
ω
;
b
)
\beta=(\omega;b)
β=(ω;b),则有:
l
(
β
)
=
∑
i
=
1
m
ln
p
(
y
i
∣
x
i
;
β
)
(7)
l(\beta)=\sum_{i=1}^m\ln{p(y_i|x_i;\beta)}\tag{7}
l(β)=i=1∑mlnp(yi∣xi;β)(7)
p
(
y
i
∣
x
i
;
β
)
=
y
i
p
1
(
x
;
β
)
+
(
1
−
y
i
)
p
0
(
x
;
β
)
(8)
p(y_i|x_i;\beta)=y_ip_1(x;\beta)+(1-yi)p_0(x;\beta)\tag{8}
p(yi∣xi;β)=yip1(x;β)+(1−yi)p0(x;β)(8)
将(8)代入(7),结合(5)(6)可知最大化式(7)等价于最小化式(9):
l
(
β
)
=
∑
i
=
1
m
(
−
y
i
β
T
x
i
+
ln
(
1
+
e
β
T
x
i
)
)
(9)
l(\beta)=\sum_{i=1}^m(-y_i\beta^Tx_i+\ln{(1+e^{\beta^Tx_i})})\tag{9}
l(β)=i=1∑m(−yiβTxi+ln(1+eβTxi))(9)
对(9)求导:
∂
l
(
β
)
∂
β
=
∑
i
=
1
m
x
i
(
−
y
i
+
p
(
y
=
0
∣
x
i
;
β
)
)
(10)
\frac{\partial{l(\beta)}}{\partial{\beta}}=\sum_{i=1}^mx_i(-y_i+p(y=0|xi;\beta))\tag{10}
∂β∂l(β)=i=1∑mxi(−yi+p(y=0∣xi;β))(10)
=
∑
i
=
1
m
x
i
(
−
y
i
+
S
i
g
m
o
i
d
(
x
;
β
)
)
(11)
=\sum_{i=1}^mx_i(-y_i+Sigmoid(x;\beta))\tag{11}
=i=1∑mxi(−yi+Sigmoid(x;β))(11)
再对式(11)进行梯度下降Gradient Descent求最优解
2.代码
import numpy as np
import csv
import pandas as pd
from sklearn.cross_validation import train_test_split
def getData(fileName):
data = pd.read_csv(fileName)
#方便矩阵运算,加入一列1,用来和bias相乘
data['one'] = 1.0
x = np.mat(data[['one','density','sugar']])
y = np.mat(data[['label']])
x_train,x_test,y_train,y_test = train_test_split(x,y,test_size=0.25,random_state=33)
return x_train,x_test,y_train,y_test
def sigmoid(z):
return 1/(1 + np.exp(-z))
def gradDescent(x,y,iter):
lr = 0.1# 学习率
# initialization,通常随机数效果会好一点
weights = np.random.rand(3,1)
for i in range(iter): # 迭代次数
#具体推导看课本笔记
p = sigmoid(x.dot(weights))
error = p -y
d = x.T.dot(error)
weights -= lr * d
return weights
def testWeights(x,y,weights):
pre=np.zeros([x.shape[0],1])
for idx,value in enumerate(np.dot(x,weights)):
if sigmoid(value)>=0.5:#sigmoid(value)>=0.5视为正例
pre[idx] = 1
else:
pre[idx] = 0
print('The predictions are:')
print(pre)
print('The real labels are:')
print(y)
fileName = 'waterMalon.csv'
x_train,x_test,y_train,y_test = getData(fileName)
weights = gradDescent(x_train,y_train,1000)
testWeights(x_test,y_test,weights)