回归就是对已知公式的未知参数进行估计。比如已知公式是y=a∗x+by=a∗x+b，未知参数是a和b，利用多真实的(x,y)训练数据对a和b的取值去自动估计。估计的方法是在给定训练样本点和已知的公式后，对于一个或多个未知参数，机器会自动枚举参数的所有可能取值，直到找到那个最符合样本点分布的参数（或参数组合）。
logistic分布
设X是连续随机变量，X服从logistic分布是指X具有下列分布函数和密度函数：
F(x)=P(x≤x)=11+e−(x−μ)/γ f(x)=F′(x)=e−(x−μ)/γγ(1+e−(x−μ)/γ)2F(x)=P(x≤x)=11+e−(x−μ)/γ f(x)=F′(x)=e−(x−μ)/γγ(1+e−(x−μ)/γ)2
其中，μμ为位置参数，γγ为形状参数。
f(x)f(x)与F(x)F(x)图像如下，其中分布函数是以(μ,12)(μ,12)为中心对阵，γγ越小曲线变化越快。
logistic回归模型
二项logistic回归模型如下：
P(Y=1|x)=exp(w⋅x+b)1+exp(w⋅x+b) P(Y=0|x)=11+exp(w⋅x+b)P(Y=1|x)=exp(w⋅x+b)1+exp(w⋅x+b) P(Y=0|x)=11+exp(w⋅x+b)
其中，x∈Rnx∈Rn是输入，Y∈{0,1}Y∈{0,1}是输出，w称为权值向量，b称为偏置，w⋅xw⋅x为w和x的内积。
参数估计
假设：
P(Y=1|x)=π(x),P(Y=0|x)=1−π(x)P(Y=1|x)=π(x),P(Y=0|x)=1−π(x)
则似然函数为:
∏_i=1N[π(x_i)]y_i[1−π(x_i)]1−y_i∏_i=1N[π(x_i)]y_i[1−π(x_i)]1−y_i
求对数似然函数：
L(w)=∑_i=1N[y_ilogπ(x_i)+(1−y_i)log(1−π(x_i))] =∑_i=1N[y_ilogπ(x_i)1−π(x_i)+log(1−π(x_i))]L(w)=∑_i=1N[y_ilogπ(x_i)+(1−y_i)log(1−π(x_i))] =∑_i=1N[y_ilogπ(x_i)1−π(x_i)+log(1−π(x_i))]
从而对L(w)L(w)求极大值，得到ww的估计值。
求极值的方法可以是梯度下降法，梯度上升法等。
梯度上升确定回归系数
logistic回归的sigmoid函数：
σ(z)=11+e−zσ(z)=11+e−z
假设logstic的函数为:
y=w_0+w_1x_1+w_2x_2+...+w_nx_ny=w_0+w_1x_1+w_2x_2+...+w_nx_n
可简写为:
y=w_0+wTxy=w_0+wTx
梯度上升算法是按照上升最快的方向不断移动，每次都增加α∇_wf(w)α∇_wf(w)，
w=w+α∇_wf(w)w=w+α∇_wf(w)
其中，∇_wf(w)∇_wf(w)为函数导数，αα为增长的步长。
训练算法
本算法的主要思想根据logistic回归的sigmoid函数来将函数值映射到有限的空间内，sigmoid函数的取值范围是0~1，从而可以把数据按照0和1分为两类。在算法中，首先要初始化所有的w权值为1，每次计算一次误差并根据误差调整w权值的大小。
- 每个回归系数都初始化为1
- 重复N次
- 计算整个数据集合的梯度
- 使用α⋅∇f(x)α⋅∇f(x)来更新w向量
- 返回回归系数
#!/usr/bin/env python
# encoding:utf-8
import math
import numpy
import time
import matplotlib.pyplot as plt
def sigmoid(x):
return 1.0 / (1 + numpy.exp(-x))
def loadData():
dataMat = []
laberMat = []
with open("test.txt", 'r') as f:
for line in f.readlines():
arry = line.strip().split()
dataMat.append([1.0, float(arry[0]), float(arry[1])])
laberMat.append(float(arry[2]))
return numpy.mat(dataMat), numpy.mat(laberMat).transpose()
def gradAscent(dataMat, laberMat, alpha=0.001, maxCycle=500):
"""general gradscent"""
start_time = time.time()
m, n = numpy.shape(dataMat)
weights = numpy.ones((n, 1))
for i in range(maxCycle):
h = sigmoid(dataMat * weights)
error = laberMat - h
weights += alpha * dataMat.transpose() * error
duration = time.time() - start_time
print "duration of time:", duration
return weights
def stocGradAscent(dataMat, laberMat, alpha=0.01):
start_time = time.time()
m, n = numpy.shape(dataMat)
weights = numpy.ones((n, 1))
for i in range(m):
h = sigmoid(dataMat[i] * weights)
error = laberMat[i] - h
weights += alpha * dataMat[i].transpose() * error
duration = time.time() - start_time
print "duration of time:", duration
return weights
def betterStocGradAscent(dataMat, laberMat, alpha=0.01, numIter=150):
"""better one, use a dynamic alpha"""
start_time = time.time()
m, n = numpy.shape(dataMat)
weights = numpy.ones((n, 1))
for j in range(numIter):
for i in range(m):
alpha = 4 / (1 + j + i) + 0.01
h = sigmoid(dataMat[i] * weights)
error = laberMat[i] - h
weights += alpha * dataMat[i].transpose() * error
duration = time.time() - start_time
print "duration of time:", duration
return weights
start_time = time.time()
def show(dataMat, laberMat, weights):
m, n = numpy.shape(dataMat)
min_x = min(dataMat[:, 1])[0, 0]
max_x = max(dataMat[:, 1])[0, 0]
xcoord1 = []; ycoord1 = []
xcoord2 = []; ycoord2 = []
for i in range(m):
if int(laberMat[i, 0]) == 0:
xcoord1.append(dataMat[i, 1]); ycoord1.append(dataMat[i, 2])
elif int(laberMat[i, 0]) == 1:
xcoord2.append(dataMat[i, 1]); ycoord2.append(dataMat[i, 2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcoord1, ycoord1, s=30, c="red", marker="s")
ax.scatter(xcoord2, ycoord2, s=30, c="green")
x = numpy.arange(min_x, max_x, 0.1)
y = (-weights[0] - weights[1]*x) / weights[2]
ax.plot(x, y)
plt.xlabel("x1"); plt.ylabel("x2")
plt.show()
if __name__ == "__main__":
dataMat, laberMat = loadData()
#weights = gradAscent(dataMat, laberMat, maxCycle=500)
#weights = stocGradAscent(dataMat, laberMat)
weights = betterStocGradAscent(dataMat, laberMat, numIter=80)
show(dataMat, laberMat, weights)未优化的程序结果如下，
随机梯度上升算法（降低了迭代的次数，算法较快，但结果不够准确）结果如下，
对αα进行优化，动态调整步长（同样降低了迭代次数，但是由于代码采用动态调整alpha，提高了结果的准确性），结果如下
本文 由 cococo点点 创作，采用 知识共享 署名-非商业性使用-相同方式共享 3.0 中国大陆 许可协议进行许可。欢迎转载，请注明出处：
转载自：cococo点点 http://www.cnblogs.com/coder2012