机器学习实战之逻辑回归logestic regression
机器学习实战这本书最大的缺点就是对代码注释的不够详细，尤其是很多运算过程直接省略。
最起码应该在开始代码前先将主要的数学结论说清楚。
逻辑回归logestic regression，
回归：指研究一组随机变量(Y1 ，Y2 ，…，Yi)和另一组(X1，X2，…，Xk)变量之间关系的统计分析方法
关系符合线性的叫线性回归，符合逻辑函数关系的叫逻辑回归。
这里用逻辑回归处理一些非线性的分类问题。
用到的函数sigmoid：近似一种阶跃函数，大于0.5归为1，小于0.5归为0；0.5作为分割线，
也可以指定别的值作为分割线。
这里涉及到矩阵乘法，矩阵相乘只有前列数和后行数一致才有意义，
得到一个前行后列的新矩阵。T表示矩阵的转置，行列互换。
w是回归系数向量，x是数据向量。
回归的目的就是为了求得最佳回归系数w，最佳回归系数可以用来测试新数据的分类。
最佳回归系数w符合其似然函数函数f(w)的最大值：
经过一系列的运算：f(w)对w的求导可以求得：
df(w)/dw=x*error
error是误差：error=y-sigmoid(z),就是error=y-sigmoid(w*x)
数据乘以回归系数后与实际分类值的误差。
为了求得f(w)的最大值，用到一种最优化算法：梯度上升算法，
还有梯度下降算法，可以求一个函数的最小值，比如最小二乘法的求解？不太懂，慢慢理解。
梯度上升算法的迭代公式：
w=w+step*df(w)/dw
step是步长值，可以自己设定，比如0.01；
在随机梯度算法中设定为动态的：step=0.01+4/（1.0+迭代次数n），随着迭代次数增多趋近于0.01；
其中：df(w)/dw称为梯度
根据上面已经算出的df(w)/dw=x*error
可得到：
w=w+step*x*error
# -*- coding: cp936 -*-
#logistics Regression,逻辑回归
#回归，指研究一组随机变量(Y1 ，Y2 ，…，Yi)和
#另一组(X1，X2，…，Xk)变量之间关系的统计分析方法
#最简单的情形是一元线性回归，由大体上有线性关系的一个自变量和一个因变量组成；
#模型是Y=a+bX+ε（X是自变量，Y是因变量，ε是随机误差）。
#这里logistics回归用到Sigmoid函数，Sigmoid函数计算的是某个向量属于1的概率。
#一般以0.5为分割线，大于0.5则认为这个向量属于1。可由此求出分割线方程。
#1.0/（1+exp(-z)）=0.5
#z=0
#w0x0+w1x1+w2x2+...+wnxn=0
from numpy import *
def sigmoid(inX):
return 1.0/(1+exp(-inX))
#sigmoid函数的输入记为z，即z=w0x0+w1x1+w2x2+...+wnxn，其中x0=1,w0表示常数项。
#如果用向量表示即为z=wTx，它表示将这两个数值向量对应元素相乘然后累加起来。
#其中，向量x是分类器的输入数据，
#w即为我们要拟合的最佳参数，从而使分类器预测更加准确。
def loadDataSet():
dataMat=[]
labelMat=[]
fr=open(r'H:\study\python\machine learning\machinelearninginaction\Ch05\testSet.txt')
for line in fr.readlines():
lineArr=line.strip().split()
dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
#testSet中，有三列数据，后一列是类标签。
#处理的时候，多加了一列数据在最前面，都为1.0
labelMat.append(int(lineArr[2]))
return dataMat,labelMat
def gradAscent(dataMatIn,classLabels):
#梯度上升，返回训练好的回归系数
#最小化损失函数，就用梯度下降，最大化似然函数，就用梯度上升。
dataMatrix=mat(dataMatIn)
labelMat=mat(classLabels).transpose()
#用到transpose变换函数，行列互换。
m,n=shape(dataMatrix)#n是列数，3列。
alpha=0.001#步长
maxCycles=500#迭代次数
weights=ones((n,1))#回归系数，3行1列
y1=[]
y2=[]
y3=[]#画图准备数据
for k in range(maxCycles):
h=sigmoid(dataMatrix*weights)
#此处是矩阵乘法，与数组乘法不同。
#大概理解是：数组乘法是行数据乘以行数据作为行数据
#矩阵乘法是第1行数据分别乘以第一列数据，然后相加作为第一行数据。
#矩阵乘法只有前的列数和后的行数相同才有意义，前后位置不可互换.
#转置 (AB)T=BTAT
error=(labelMat-h)
weights=weights+alpha*dataMatrix.transpose()*error
#回归系数为w，样本矩阵为x，误差为e，步长为alpha
#更新回归系数的公式就是：w=w+alpha*x*e
#矩阵要考虑行列的方向。
weights=array(weights)
y1.append(weights[0])
y2.append(weights[1])
y3.append(weights[2])
x=[k for k in range(1,maxCycles+1)]
return weights,y1,y2,y3,x
def plotBestFit(wei):
import matplotlib.pyplot as plt
weights=array(wei)
dataMat,labelMat=loadDataSet()
dataArr=array(dataMat)
n=shape(dataArr)[0]#得到行数，等同于n=len(dataArr)
xcord1=[]
ycord1=[]
xcord2=[]
ycord2=[]
for i in range(n):
if labelMat[i]==1:
xcord1.append(dataArr[i][1])#第1列是添加的数据1.0
ycord1.append(dataArr[i][2])
elif labelMat[i]==0:
xcord2.append(dataArr[i][1])
ycord2.append(dataArr[i][2])
plt.scatter(xcord1,ycord1,c='r',label='1')
plt.scatter(xcord2,ycord2,c='b',label='0')
xMin=dataArr.min(0)[1]#获取数据中第2列的最小值
xMax=dataArr.max(0)[1]
x=arange(xMin-1,xMax+1,0.1)
y=(-weights[0]-weights[1]*x)/weights[2]
#分割线方程：1.0/（1+exp(-z)）=0.5
#z=0
plt.plot(x,y,c='k',label='0.5')
plt.xlabel('X1')
plt.ylabel('X2')
plt.title('title')
plt.legend(loc='upper right')
plt.show()
def cal():
import matplotlib.pyplot as plt
'练习，用梯度上升求函数的最大值'
'f(x)=-x*x+3x-1，求导，df(x)=-2x+3,就是梯度'
'梯度上升迭代公式：w=w+step*df(w)'
#先划出函数f(x)
x=arange(-6,8,0.1)
y=-x*x+3*x-1
plt.plot(x,y,c='g')
step=0.01#设定步长
presision=0.00001#设定精度，应该指的是dw的精度
w_new=6#从6开始找寻
w_old=0
while abs(w_new-w_old)>presision:
w_old=w_new
w_new=w_old+step*(-2*w_old+3)
print(w_new)
plt.annotate('(%s,%s)' %(w_new,-w_new*w_new+3*w_new-1),xy=(w_new,-w_new*w_new+3*w_new-1))
plt.show()
def stocGradAscent0(dataMatrix,classLabels):
'随机梯度上升算法，单次'
dataMatrix=array(dataMatrix)
#对数据处理成数组，否则后面会报错。
m,n=shape(dataMatrix)
weights=ones(n)
alpha=0.01
for i in range(m):
h=sigmoid(sum(dataMatrix[i]*weights))
error=classLabels[i]-h
weights=weights+error*alpha*dataMatrix[i]
#在做数学运算的时候，一定要注意类型，是列表，数组还是矩阵，因为
#运算结果会很不一样。
return weights
def stocGradAscent1(dataMatrix,classLabels,num=150):
'多次迭代'
dataMatrix=array(dataMatrix)
m,n=shape(dataMatrix)
alpha=0.01
weights=ones(n)
y1=[]
y2=[]
y3=[]
for k in range(num):
for i in range(m):
h=sigmoid(sum(dataMatrix[i]*weights))
error=classLabels[i]-h
weights=weights+error*alpha*dataMatrix[i]
y1.append(weights[0])
y2.append(weights[1])
y3.append(weights[2])
x=[k for k in range(1,num*m+1)]
return weights,y1,y2,y3,x
def stocGradAscent2(dataMatrix,classLabels,num=150):
'随机选择,对算法效率改进明显'
dataMatrix=array(dataMatrix)
m,n=shape(dataMatrix)
alpha=0.01
weights=ones(n)
'''
y1=[]
y2=[]
y3=[]
'''
for k in range(num):
dataIndex=range(m)
for i in range(m):
alpha=0.01+4/(1.0+i+k)#每次迭代都调整步长，越来越小趋近于0.01
randIndex=int(random.uniform(0,len(dataIndex)))
#也可以写成：int(random.rand()*len(dataIndex))
h=sigmoid(sum(dataMatrix[randIndex]*weights))
error=classLabels[randIndex]-h
weights=weights+error*alpha*dataMatrix[randIndex]
del(dataIndex[randIndex])
'''
y1.append(weights[0])
y2.append(weights[1])
y3.append(weights[2])
'''
'x=[k for k in range(1,num*m+1)]'
return weights#',y1,y2,y3,x'
def classifyVector(inX,weights):
'逻辑回归分类函数'
prob=sigmoid(sum(inX*weights))
if prob>0.5:
return 1.0
else:
return 0.0
def colicTest():
frTrain=open(r'H:\study\python\machine learning\machinelearninginaction\Ch05\horseColicTraining.txt')
frTest=open(r'H:\study\python\machine learning\machinelearninginaction\Ch05\horseColicTest.txt')
trainingSet=[]
trainingLabels=[]
for line in frTrain.readlines():
currLine=line.strip().split('\t')#获得字符串列表
lineArr=[]
for i in range(21):
lineArr.append(float(currLine[i]))#将字符串转换成浮点数据
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[21]))
trainWeights=stocGradAscent2(trainingSet,trainingLabels,500)
#上行代码会出现runtime error，可能是运算量太大，可以忽略。
errorCount=0
numTestVec=0.0
for line in frTest.readlines():
numTestVec+=1.0
#直接len(frTest.readlines())
currLine=line.strip().split('\t')
lineArr=[]
for i in range(21):
lineArr.append(float(currLine[i]))
if int(classifyVector(array(lineArr),trainWeights))!= int(currLine[21]):
errorCount+=1
errorRate=float(errorCount)/numTestVec
print('错误率：%f' %errorRate)
return errorRate
#错误率：0.358209,。事实上，这个结果并不差，
#因为有3 0 %的数据缺失。当然，如果调整迭代次数和步长，]
#平均错误率可以降到20%左右。
if __name__=='__main__':
'''
dataArr,labelMat=loadDataSet()
weights,y1,y2,y3,x=gradAscent(dataArr,labelMat)
weights,error=gradAscent(dataArr,labelMat)
print(weights)
print('----------------')
print(error)
print('----------------')
h=sigmoid(mat(dataArr)*weights)
print(h)
import matplotlib.pyplot as plt
plt.subplot(3,1,1)
plt.plot(x,y1,'r')
plt.ylabel('w0')
plt.subplot(3,1,2)
plt.plot(x,y2,'y')
plt.ylabel('w1')
plt.subplot(3,1,3)
plt.plot(x,y3,'b')
plt.ylabel('w2')
plt.xlabel('numberofcycles')
plt.show()
'''
'''
import matplotlib.pyplot as plt
dataArr,labelMat=loadDataSet()
weights,y1,y2,y3,x=gradAscent(dataArr,labelMat)
weights,yy1,yy2,yy3,xx=stocGradAscent1(dataArr,labelMat,500)
weights,yyy1,yyy2,yyy3,xxx=stocGradAscent2(dataArr,labelMat)
plt.subplot(3,3,1,)
plt.plot(x,y1,'r')
plt.ylabel('w0')
plt.subplot(3,3,4)
plt.plot(x,y2,'y')
plt.ylabel('w1')
plt.subplot(3,3,7)
plt.plot(x,y3,'b')
plt.ylabel('w2')
plt.xlabel('numberofcycles')
plt.subplot(3,3,2,)
plt.plot(xx,yy1,'r')
plt.ylabel('w0')
plt.subplot(3,3,5)
plt.plot(xx,yy2,'y')
plt.ylabel('w1')
plt.subplot(3,3,8)
plt.plot(xx,yy3,'b')
plt.ylabel('w2')
plt.xlabel('numberofcycles')
plt.subplot(3,3,3,)
plt.plot(xxx,yyy1,'r')
plt.ylabel('w0')
plt.subplot(3,3,6)
plt.plot(xxx,yyy2,'y')
plt.ylabel('w1')
plt.subplot(3,3,9)
plt.plot(xxx,yyy3,'b')
plt.ylabel('w2')
plt.xlabel('numberofcycles')
plt.show()
#比较可知，随机梯度上升算法的效率最高，用较少的迭代次数获得参数的稳定值。
#前两个本质效果是一样的。
'''
errorSum=0.0
numTests=10
for i in range(numTests):
errorSum+=colicTest()
arrageError=errorSum/float(numTests)
print('测试%d 次的平均错误率是：%f' %(numTests,arrageError))