公式解释与python代码一、二分类
1、逻辑回归（logistic回归）
首先给出线性回归模型：
写成向量形式为：
同时“广义线性回归”模型为：
可以从线性回归的回归模型引出logistic回归的分类模型，logistic回归是处理二分类问题的，所以输出的标记y={0,1}，并且线性回归模型产生的预测值z=wx+b是一个实值，所以我们将实值z转化成0/1值便可，这样有一个可选函数便是“单位阶跃函数”：
这种如果预测值大于0便判断为正例，小于0则判断为反例，等于0则可任意判断！
但是单位阶跃函数是非连续的函数，我们需要一个连续的函数，“Sigmoid函数”便可以很好的取代单位阶跃函数：
这样我们在原来的线性回归模型外套上sigmoid函数便形成了logistic回归模型的预测函数，可以用于二分类问题：
对上式的预测函数做一个变换为：
接下来求w：
1、为求解参数w，我们需要定义一个准则函数 J(w)，利用准则函数求解参数w
2、我们通过最大似然估计法定义准则函数J(w)
3、接下来通过最大化准则函数J(w)便可求出参数w的迭代表达式
4、为了更好地使用数据求出参数w，我们将第三步得到的w的迭代时向量化。自此便完成了对于参数w的推导过程，接下来便可以进行实例应用了！
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
from numpy import *
import matplotlib.pyplot as plt
#从文件中加载数据：特征X，标签label
def loadDataSet():
dataMatrix=[]
dataLabel=[]
#这里给出了python 中读取文件的简便方式
f=open('testSet.txt')
for line in f.readlines():
#print(line)
lineList=line.strip().split()
dataMatrix.append([1,float(lineList[0]),float(lineList[1])])
dataLabel.append(int(lineList[2]))
#for i in range(len(dataMatrix)):
# print(dataMatrix[i])
#print(dataLabel)
#print(mat(dataLabel).transpose())
matLabel=mat(dataLabel).transpose()
return dataMatrix,matLabel
#logistic回归使用了sigmoid函数
def sigmoid(inX):
return 1/(1+exp(-inX))
#函数中涉及如何将list转化成矩阵的操作：mat()
#同时还含有矩阵的转置操作：transpose()
#还有list和array的shape函数
#在处理矩阵乘法时，要注意的便是维数是否对应
#graAscent函数实现了梯度上升法，隐含了复杂的数学推理
#梯度上升算法，每次参数迭代时都需要遍历整个数据集
def graAscent(dataMatrix,matLabel):
m,n=shape(dataMatrix)
matMatrix=mat(dataMatrix)
w=ones((n,1))
alpha=0.001
num=500
for i in range(num):
error=sigmoid(matMatrix*w)-matLabel
w=w-alpha*matMatrix.transpose()*error
return w
#随机梯度上升算法的实现，对于数据量较多的情况下计算量小，但分类效果差
#每次参数迭代时通过一个数据进行运算
def stocGraAscent(dataMatrix,matLabel):
m,n=shape(dataMatrix)
matMatrix=mat(dataMatrix)
w=ones((n,1))
alpha=0.001
num=20 #这里的这个迭代次数对于分类效果影响很大，很小时分类效果很差
for i in range(num):
for j in range(m):
error=sigmoid(matMatrix[j]*w)-matLabel[j]
w=w-alpha*matMatrix[j].transpose()*error
return w
#改进后的随机梯度上升算法
#从两个方面对随机梯度上升算法进行了改进,正确率确实提高了很多
#改进一：对于学习率alpha采用非线性下降的方式使得每次都不一样
#改进二：每次使用一个数据，但是每次随机的选取数据，选过的不在进行选择
def stocGraAscent1(dataMatrix,matLabel):
m,n=shape(dataMatrix)
matMatrix=mat(dataMatrix)
w=ones((n,1))
num=200 #这里的这个迭代次数对于分类效果影响很大，很小时分类效果很差
setIndex=set([])
for i in range(num):
for j in range(m):
alpha=4/(1+i+j)+0.01
dataIndex=random.randint(0,100)
while dataIndex in setIndex:
setIndex.add(dataIndex)
dataIndex=random.randint(0,100)
error=sigmoid(matMatrix[dataIndex]*w)-matLabel[dataIndex]
w=w-alpha*matMatrix[dataIndex].transpose()*error
return w
#绘制图像
def draw(weight):
x0List=[];y0List=[];
x1List=[];y1List=[];
f=open('testSet.txt','r')
for line in f.readlines():
lineList=line.strip().split()
if lineList[2]=='0':
x0List.append(float(lineList[0]))
y0List.append(float(lineList[1]))
else:
x1List.append(float(lineList[0]))
y1List.append(float(lineList[1]))
fig=plt.figure()
ax=fig.add_subplot(111)
ax.scatter(x0List,y0List,s=10,c='red')
ax.scatter(x1List,y1List,s=10,c='green')
xList=[];yList=[]
x=arange(-3,3,0.1)
for i in arange(len(x)):
xList.append(x[i])
y=(-weight[0]-weight[1]*x)/weight[2]
for j in arange(y.shape[1]):
yList.append(y[0,j])
ax.plot(xList,yList)
plt.xlabel('x1');plt.ylabel('x2')
plt.show()
if __name__ == '__main__':
dataMatrix,matLabel=loadDataSet()
#weight=graAscent(dataMatrix,matLabel)
weight=stocGraAscent1(dataMatrix,matLabel)
print(weight)
draw(weight)二、Fisher判断分析
fisher原理
费歇（FISHER）判别思想是投影，使多维问题简化为一维问题来处理。选择一个适当的投影轴,使所有的样品点都投影到这个轴上得到一个投影值。对这个投影轴的方向的要求是：使每一类内的投影值所形成的类内离差尽可能小，而不同类间的投影值所形成的类间离差尽可能大。
两类问题Fisher线性分类器计算步骤如下：
三、SVM支持向量模型
SVM（Support Vector Machines）——支持向量机是在所有知名的数据挖掘算法中最健壮，最准确的方法之一，它属于二分类算法，可以支持线性和非线性的分类。一个能使两类之间的空间大小最大的一个超平面。
这个超平面在二维平面上看到的就是一条直线，在三维空间中就是一个平面...，因此，我们把这个划分数据的决策边界统称为超平面。离这个超平面最近的点就叫做支持向量，点到超平面的距离叫间隔。支持向量机就是要使超平面和支持向量之间的间隔尽可能的大，这样超平面才可以将两类样本准确的分开，而保证间隔尽可能的大就是保证我们的分类器误差尽可能的小，尽可能的健壮。
如何确定这个超平面呢？从直观上而言，这个超平面应该是最适合分开两类数据的直线。而判定“最适合”的标准就是这条直线离直线两边的数据的间隔最大。所以，得寻找有着最大间隔的超平面。
如下图所示，中间的实线便是寻找到的最优超平面（Optimal Hyper Plane），其到两条虚线边界的距离相等，这个距离便是几何间隔
#svm算法的实现
from numpy import*
import random
from time import*
def loadDataSet(fileName):#输出dataArr(m*n),labelArr(1*m)其中m为数据集的个数
dataMat=[];labelMat=[]
fr=open(fileName)
for line in fr.readlines():
lineArr=line.strip().split('\t')#去除制表符，将数据分开
dataMat.append([float(lineArr[0]),float(lineArr[1])])#数组矩阵
labelMat.append(float(lineArr[2]))#标签
return dataMat,labelMat
def selectJrand(i,m):#随机找一个和i不同的j
j=i
while(j==i):
j=int(random.uniform(0,m))
return j
def clipAlpha(aj,H,L):#调整大于H或小于L的alpha的值
if aj>H:
aj=H
if aj<L:
aj=L
return aj
def smoSimple(dataMatIn,classLabels,C,toler,maxIter):
dataMatrix=mat(dataMatIn);labelMat=mat(classLabels).transpose()#转置
b=0;m,n=shape(dataMatrix)#m为输入数据的个数，n为输入向量的维数
alpha=mat(zeros((m,1)))#初始化参数，确定m个alpha
iter=0#用于计算迭代次数
while (iter<maxIter):#当迭代次数小于最大迭代次数时（外循环）
alphaPairsChanged=0#初始化alpha的改变量为0
for i in range(m):#内循环
fXi=float(multiply(alpha,labelMat).T*\
(dataMatrix*dataMatrix[i,:].T))+b#计算f(xi)
Ei=fXi-float(labelMat[i])#计算f(xi)与标签之间的误差
if ((labelMat[i]*Ei<-toler)and(alpha[i]<C))or\
((labelMat[i]*Ei>toler)and(alpha[i]>0)):#如果可以进行优化
j=selectJrand(i,m)#随机选择一个j与i配对
fXj=float(multiply(alpha,labelMat).T*\
(dataMatrix*dataMatrix[j,:].T))+b#计算f(xj)
Ej=fXj-float(labelMat[j])#计算j的误差
alphaIold=alpha[i].copy()#保存原来的alpha(i)
alphaJold=alpha[j].copy()
if(labelMat[i]!=labelMat[j]):#保证alpha在0到c之间
L=max(0,alpha[j]-alpha[i])
H=min(C,C+alpha[j]-alpha[i])
else:
L=max(0,alpha[j]+alpha[i]-C)
H=min(C,alpha[j]+alpha[i])
if L==H:print('L=H');continue
eta=2*dataMatrix[i,:]*dataMatrix[j,:].T-\
dataMatrix[i,:]*dataMatrix[i,:].T-\
dataMatrix[j,:]*dataMatrix[j,:].T
if eta>=0:print('eta=0');continue
alpha[j]-=labelMat[j]*(Ei-Ej)/eta
alpha[j]=clipAlpha(alpha[j],H,L)#调整大于H或小于L的alpha
if (abs(alpha[j]-alphaJold)<0.0001):
print('j not move enough');continue
alpha[i]+=labelMat[j]*labelMat[i]*(alphaJold-alpha[j])
b1=b-Ei-labelMat[i]*(alpha[i]-alphaIold)*\
dataMatrix[i,:]*dataMatrix[i,:].T-\
labelMat[j]*(alpha[j]-alphaJold)*\
dataMatrix[i,:]*dataMatrix[j,:].T#设置b
b2=b-Ej-labelMat[i]*(alpha[i]-alphaIold)*\
dataMatrix[i,:]*dataMatrix[i,:].T-\
labelMat[j]*(alpha[j]-alphaJold)*\
dataMatrix[j,:]*dataMatrix[j,:].T
if (0<alpha[i])and(C>alpha[j]):b=b1
elif(0<alpha[j])and(C>alpha[j]):b=b2
else:b=(b1+b2)/2
alphaPairsChanged+=1
print('iter:%d i:%d,pairs changed%d'%(iter,i,alphaPairsChanged))
if (alphaPairsChanged==0):iter+=1
else:iter=0
print('iteraction number:%d'%iter)
return b,alpha
#定义径向基函数
def kernelTrans(X, A, kTup):#定义核转换函数（径向基函数）
m,n = shape(X)
K = mat(zeros((m,1)))
if kTup[0]=='lin': K = X * A.T #线性核K为m*1的矩阵
elif kTup[0]=='rbf':
for j in range(m):
deltaRow = X[j,:] - A
K[j] = deltaRow*deltaRow.T
K = exp(K/(-1*kTup[1]**2)) #divide in NumPy is element-wise not matrix like Matlab
else: raise NameError('Houston We Have a Problem -- \
That Kernel is not recognized')
return K
class optStruct:
def __init__(self,dataMatIn, classLabels, C, toler, kTup): # Initialize the structure with the parameters
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataMatIn)[0]
self.alphas = mat(zeros((self.m,1)))
self.b = 0
self.eCache = mat(zeros((self.m,2))) #first column is valid flag
self.K = mat(zeros((self.m,self.m)))
for i in range(self.m):
self.K[:,i] = kernelTrans(self.X, self.X[i,:], kTup)
def calcEk(oS, k):
fXk = float(multiply(oS.alphas,oS.labelMat).T*oS.K[:,k] + oS.b)
Ek = fXk - float(oS.labelMat[k])
return Ek
def selectJ(i, oS, Ei):
maxK = -1; maxDeltaE = 0; Ej = 0
oS.eCache[i] = [1,Ei]
validEcacheList = nonzero(oS.eCache[:,0].A)[0]
if (len(validEcacheList)) > 1:
for k in validEcacheList: #loop through valid Ecache values and find the one that maximizes delta E
if k == i: continue #don't calc for i, waste of time
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if (deltaE > maxDeltaE):
maxK = k; maxDeltaE = deltaE; Ej = Ek
return maxK, Ej
else: #in this case (first time around) we don't have any valid eCache values
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
def updateEk(oS, k):#after any alpha has changed update the new value in the cache
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
def innerL(i, oS):
Ei = calcEk(oS, i)
if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
j,Ej = selectJ(i, oS, Ei) #this has been changed from selectJrand
alphaIold = oS.alphas[i].copy(); alphaJold = oS.alphas[j].copy()
if (oS.labelMat[i] != oS.labelMat[j]):
L = max(0, oS.alphas[j] - oS.alphas[i])
H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
else:
L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
H = min(oS.C, oS.alphas[j] + oS.alphas[i])
if L==H: print("L==H"); return 0
eta = 2.0 * oS.K[i,j] - oS.K[i,i] - oS.K[j,j] #changed for kernel
if eta >= 0: print("eta>=0"); return 0
oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta
oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
updateEk(oS, j) #added this for the Ecache
if (abs(oS.alphas[j] - alphaJold) < 0.00001): print("j not moving enough"); return 0
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])#update i by the same amount as j
updateEk(oS, i) #added this for the Ecache #the update is in the oppostie direction
b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,i] - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[i,j]
b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.K[i,j]- oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.K[j,j]
if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
else: oS.b = (b1 + b2)/2.0
return 1
else: return 0
#smoP函数用于计算超平的alpha,b
def smoP(dataMatIn, classLabels, C, toler, maxIter,kTup=('lin', 0)): #完整的Platter SMO
oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler, kTup)
iter = 0#计算循环的次数
entireSet = True; alphaPairsChanged = 0
while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
alphaPairsChanged = 0
if entireSet: #go over all
for i in range(oS.m):
alphaPairsChanged += innerL(i,oS)
print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
iter += 1
else:#go over non-bound (railed) alphas
nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
for i in nonBoundIs:
alphaPairsChanged += innerL(i,oS)
print("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
iter += 1
if entireSet: entireSet = False #toggle entire set loop
elif (alphaPairsChanged == 0): entireSet = True
print("iteration number: %d" % iter)
return oS.b,oS.alphas
#calcWs用于计算权重值w
def calcWs(alphas,dataArr,classLabels):#计算权重W
X = mat(dataArr); labelMat = mat(classLabels).transpose()
m,n = shape(X)
w = zeros((n,1))
for i in range(m):
w += multiply(alphas[i]*labelMat[i],X[i,:].T)
return w
#值得注意的是测试准确与k1和C的取值有关。
def testRbf(k1=1.3):#给定输入参数K1
#测试训练集上的准确率
dataArr,labelArr = loadDataSet('testSetRBF.txt')#导入数据作为训练集
b,alphas = smoP(dataArr, labelArr, 200, 0.0001, 10000, ('rbf', k1)) #C=200 important
datMat=mat(dataArr); labelMat = mat(labelArr).transpose()
svInd=nonzero(alphas.A>0)[0]#找出alphas中大于0的元素的位置
#此处需要说明一下alphas.A的含义
sVs=datMat[svInd] #获取支持向量的矩阵，因为只要alpha中不等于0的元素都是支持向量
labelSV = labelMat[svInd]#支持向量的标签
print("there are %d Support Vectors" % shape(sVs)[0])#输出有多少个支持向量
m,n = shape(datMat)#数据组的矩阵形状表示为有m个数据，数据维数为n
errorCount = 0#计算错误的个数
for i in range(m):#开始分类，是函数的核心
kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))#计算原数据集中各元素的核值
predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b#计算预测结果y的值
if sign(predict)!=sign(labelArr[i]): errorCount += 1#利用符号判断类别
### sign（a）为符号函数：若a>0则输出1，若a<0则输出-1.###
print("the training error rate is: %f" % (float(errorCount)/m))
#2、测试测试集上的准确率
dataArr,labelArr = loadDataSet('testSetRBF2.txt')
errorCount = 0
datMat=mat(dataArr)#labelMat = mat(labelArr).transpose()此处可以不用
m,n = shape(datMat)
for i in range(m):
kernelEval = kernelTrans(sVs,datMat[i,:],('rbf', k1))
predict=kernelEval.T * multiply(labelSV,alphas[svInd]) + b
if sign(predict)!=sign(labelArr[i]): errorCount += 1
print("the test error rate is: %f" % (float(errorCount)/m))
def main():
t1=time()
dataArr,labelArr=loadDataSet('testSet.txt')
b,alphas=smoP(dataArr,labelArr,0.6,0.01,40)
ws=calcWs(alphas,dataArr,labelArr)
testRbf()
t2=time()
print("程序所用时间为%ss"%(t2-t1))
if __name__=='__main__':
main()四、多分类（多分类Logistic回归模型）
Logistic回归分析用于研究X对Y的影响，并且对X的数据类型没有要求，X可以为定类数据，也可以为定量数据，但要求Y必须为定类数据，并且根据Y的选项数，使用相应的数据分析方法。
只要是logistic回归，都是研究X对于Y的影响，区别在于因变量Y上，logistic回归时，因变量Y是看成定类数据的，如果为二元（即选项只有2个），那么就是二元logistic回归； 如果Y是多个类别且类别之间无法进行对比程度或者大小，则为多分类logistic回归；如果Y是多个类别且类别之间可以对比程度大小（也称为定量数据，或者有序定类数据），此时则使用有序logistic回归。
多分类logistic回归的难点在于：因变量为类别数据，研究X对Y的影响时，如果为类别数据，那么不能说越如何越如何，比如不能说越满意越愿意购买；而只能说相对小米手机来说，对于手机外观越满意越愿意购买苹果手机。这就是类别数据的特点，一定是相对某某而言。这就导致了多分类logistic回归分析时，文字分析的难度加大，最好是使用SPSSAU的智能文字分析对应查看。
单独进行多分类logistic回归时，通常需要有以下步骤，分别是数据处理，模型似然比检验，参数估计分析和模型预测效果分析共4个步骤。
多类别逻辑回归是一个将逻辑回归一般化成多类别问题得到的分类方法。用更加专业的话来说，它就是一个用来预测一个具有类别分布的因变量不同可能结果的概率的模型。
# coding=utf-8
import random
import matplotlib.pyplot as plt
import numpy as np
x, y = [], []
x_test1, x_test2, x_test3 = [], [], []
# 随机生成3种不同分类的点，分别打上标签存在y中
for i in range(0, 20):
x1 = random.random()
x2 = random.random()
if x1 + x2 < 1:
x.append([x1, x2, 1])
x_test1.append([x1, x2])
y.append([1, 0, 0])
x.append([x1 * 2, x2 + 1, 1])
x_test2.append([x1 * 2, x2 + 1])
y.append([0, 1, 0])
if x1 > x2:
x.append([x1 + 1, x2, 1])
x_test3.append([x1 + 1, x2])
y.append([0, 0, 1])
# 将list转换为numpy array
x = np.array(x)
y = np.array(y)
# 学习率
lr = 0.5
w = np.ones((3, 3))
m = len(x)
for i in range(0, 1000):
w_gard = 0
for j in range(0, m):
# sigmoid函数分类，这里用到了矩阵乘法
w_gard += -x[j] * np.atleast_2d(1 / (1 + np.exp(-np.matmul(w, x[j]))) - y[j]).T
w += lr * w_gard / m
# 画点
plt.plot(np.array(x_test1)[:, 0], np.array(x_test1)[:, 1], 'ro')
plt.plot(np.array(x_test2)[:, 0], np.array(x_test2)[:, 1], 'gx')
plt.plot(np.array(x_test3)[:, 0], np.array(x_test3)[:, 1], 'b.')
# 画分类直线
plt.plot([0, 2], [-w[0][2] / w[0][1], -(w[0][2] + 2 * w[0][0]) / w[0][1]], 'r')
plt.plot([0, 2], [-w[1][2] / w[1][1], -(w[1][2] + 2 * w[1][0]) / w[1][1]], 'g')
plt.plot([0, 2], [-w[2][2] / w[2][1], -(w[2][2] + 2 * w[2][0]) / w[2][1]], 'b')
plt.show()
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
