-
Notifications
You must be signed in to change notification settings - Fork 3
/
logistic.py
155 lines (127 loc) · 4.47 KB
/
logistic.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
# ----------------------------------------------------
# Copyright (c) 2017, Wray Zheng. All Rights Reserved.
# Distributed under the BSD License.
# ----------------------------------------------------
# -*- coding: utf-8 -*-
import numpy as np
DEBUG = False
#############################################
# 调试信息输出
#############################################
def debug(*args, **kwords):
if DEBUG:
print(*args, **kwords)
#############################################
# sigmoid 函数（Y = 1 的概率)
#############################################
def sigmoid(theta, X):
result = 1.0 / (1 + np.exp(-theta * X.T))
return result
#############################################
# 计算梯度
# X 为特征矩阵，Y 为标签，theta 为参数
#############################################
def gradient(X, Y, theta):
return (sigmoid(theta, X) - Y) * X / Y.size
#############################################
# 计算 Hessian 矩阵
# X 为特征矩阵，Y 为标签，theta 为参数
#############################################
def hessianMatrix(X, Y, theta):
h = sigmoid(theta, X)
t = np.multiply(h, (1-h))
M = np.mat(X) / Y.size
for i in range(theta.size):
M[:,i] = np.multiply(M[:,i], t.T)
return X.T * M
#############################################
# 牛顿法估计参数
# X 为特征矩阵，Y 为标签，iterNum 为迭代次数
#############################################
def newtonMethod(X, Y, iterNum=10):
# 在矩阵 X 的第一列插入 1
X = np.insert(X, 0, 1, 1)
# m 是训练样本数，n-1 是样本的特征数
m, n = X.shape
# 初始化 theta 值
theta = np.mat(np.zeros(n))
# 迭代求解theta
for iterIndex in range(iterNum):
g = gradient(X, Y, theta)
H = hessianMatrix(X, Y, theta)
theta -= (H.I * g.T).T
debug("theta({}):\n{}\n".format(iterIndex + 1, theta))
return theta
#############################################
# 正则化牛顿法估计参数
# X 为特征矩阵，Y 为标签，iterNum 为迭代次数
#############################################
def regularizedNewtonMethod(X, Y, iterNum=10):
# 在矩阵 X 的第一列插入 1
X = np.insert(X, 0, 1, 1)
# m 是训练样本数，n-1 是样本的特征数
m, n = X.shape
# 初始化 theta 值
theta = np.mat(np.zeros(n))
# 惩罚系数
lambda_ = 0.01 / Y.size
# hessian 矩阵正则项
A = np.eye(theta.size)
A[0,0] = 0
# 迭代求解theta
for iterIndex in range(iterNum):
g = gradient(X, Y, theta) + lambda_ * theta
H = hessianMatrix(X, Y, theta) + lambda_ * A
theta -= (H.I * g.T).T
debug("theta({}):\n{}\n".format(iterIndex + 1, theta))
return theta
#############################################
# 梯度下降法估计参数
# X 是样本特征矩阵，每一行表示一个样本
# Y 是样本 X 中对应的标签，数组类型
# alpha 表示步长
#############################################
def gradientDescent(X, Y, alpha=0.001):
# 在矩阵 X 的第一列插入 1
X = np.insert(X, 0, 1, 1)
# m 是训练样本数，n-1 是样本的特征数
m, n = X.shape
# 初始化 theta 值
theta = np.mat(np.zeros(n))
# 迭代求解theta
for i in range(200000):
theta -= alpha * gradient(X, Y, theta)
debug("theta:", theta)
return theta
#############################################
# 正则化的梯度下降法估计参数
# X 是样本特征矩阵，每一行表示一个样本
# Y 是样本 X 中对应的标签，数组类型
# alpha 表示步长
#############################################
def regularizedGradientDescent(X, Y, alpha=0.001):
# 在矩阵 X 的第一列插入 1
X = np.insert(X, 0, 1, 1)
# m 是训练样本数，n-1 是样本的特征数
m, n = X.shape
# 初始化 theta 值
theta = np.mat(np.zeros(n))
# 惩罚系数
lambda_ = 0.01 / Y.size
# 迭代求解theta
for i in range(200000):
theta -= alpha * gradient(X, Y, theta)
theta[0,1:] -= lambda_ * theta[0,1:]
debug("theta:", theta)
return theta
#############################################
# 预测样本对应 Y=1 的概率
# X 是测试样本矩阵
# theta 是牛顿法估计出的参数
#############################################
def predict(X, theta):
# 在矩阵 X 的第一列插入 1
X = np.insert(X, 0, 1, 1)
# 计算 X 样本对应的 Y = 1 的概率
predictedY = np.array(sigmoid(theta, X)).flatten()
return predictedY