逻辑回归实现（python）

逻辑回归

g(z)=11+ez

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
z=np.arange(-5,5,0.1)
g=[]
for i in z:
g.append(1/(1+np.exp(-i)))
plt.xlabel('Input')
plt.ylabel('Output')
plt.title('Sigmoid')
plt.xlim(-5,5)
plt.ylim(0,1)
plt.grid(True)
plt.plot(z,g)

hθ(x)=g(θTx)

path='C:/Users/logit_data/data1.txt'
data=np.loadtxt(path,delimiter=',')
#查看数据维数
data.shape
(100, 3)

#查看数据内容
data[:5,:]
array([[34.62365962, 78.02469282,  0.        ],
[30.28671077, 43.89499752,  0.        ],
[35.84740877, 72.90219803,  0.        ],
[60.18259939, 86.3085521 ,  1.        ],
[79.03273605, 75.34437644,  1.        ]])

#定义一个可视化函数
def plotData(data,label_x,label_y,label_pos,label_neg,axes=None):
#获取不同标签
neg=data[:,2]==0
pos=data[:,2]==1
if axes==None:
axes=plt.gca()
axes.scatter(data[pos][:,0],data[pos][:,1],color='g',marker='+',s=60,linewidth=2,label=label_pos)
axes.scatter(data[neg][:,0],data[neg][:,1],color='r',s=60,linewidth=2,label=label_neg)
axes.set_xlabel(label_x)
axes.set_ylabel(label_y)
axes.legend()
plotData(data,'Exam 1 score','Exam 2 score','Pass','Fail')

l(θ)=1mi=1m[y(i)log(hθ(x(i)))+(1y(i))log(1hθ(x(i)))]

l(θ)=1m((log(g(Xθ))Ty+(log(1g(Xθ))T(1y))

#首先定义sigmoid函数
def sigmoid(x):
return (1/(1+np.exp(-x)))
#定义似然函数
def likelihoodFunction(theta,x,y):
m=y.size #m为样本数
h=sigmoid(x.dot(theta))
l=-1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y)) #目标似然函数是要最大化，我们仍以最小化求解，前面加负号

if np.isnan(l[0]):
return (np.inf)
return l[0]

δl(θ)δθj=1mi=1m(y(i)hθ(x(i)))x(i)j

δJ(θ)δθj=1mXT(yg(Xθ))

#梯度函数
m=y.size
h=sigmoid(x.dot(theta.reshape(-1,1)))
return (grad.flatten()) #降到一维
#将样本和目标值取出
x=np.c_[np.ones((data.shape[0],1)),data[:,0:2]] #前面添加一列1，实际为2+1维矩阵
y=np.c_[data[:,2]]
initial_theta=np.zeros(x.shape[1])
likelihood=likelihoodFunction(initial_theta,x,y)
print('likelihood:',likelihood)
print('grad:',grad)
likelihood: 0.6931471805599452


from scipy.optimize import minimize
res
      fun: 0.20349770158950983
hess_inv: array([[ 2.85339493e+03, -2.32908823e+01, -2.27416470e+01],
[-2.32908823e+01,  2.04489131e-01,  1.72969525e-01],
[-2.27416470e+01,  1.72969525e-01,  1.96170322e-01]])
jac: array([-2.68557631e-09,  4.36433478e-07, -1.39671758e-06])
message: 'Optimization terminated successfully.'
nfev: 34
nit: 25
njev: 30
status: 0
success: True
x: array([-25.16131634,   0.2062316 ,   0.20147143])

#构建一个预测函数
def predict(theta,x,threshold=0.5):
p=sigmoid(x.dot(theta.T))>=threshold
return (p.astype('int'))
#预测一位同学考试1得了45，考试2得了85的通过情况
predict(res.x,np.array([1,45,85]))
1


plt.scatter(45,85,s=60,c='r',marker='v',label='(45,85)')
plotData(data,'Exam 1','Exam 2','Pass','Fail')
x1_min,x1_max=x[:,1].min(),x[:,1].max()
x2_min,x2_max=x[:,2].min(),x[:,2].max()
X1,X2=np.meshgrid(np.linspace(x1_min,x2_max),np.linspace(x2_min,x2_max))
h=sigmoid(np.c_[np.ones((X1.ravel().shape[0],1)),X1.ravel(),X2.ravel()].dot(res.x))
h=h.reshape(X1.shape)
plt.contour(X1,X2,h,[0.5],linewidths=1,colors='b')

#导入学生成绩数据
path2='C:/Users/logit_data/data2.txt'
data2[:5,:]
array([[ 0.051267,  0.69956 ,  1.      ],
[-0.092742,  0.68494 ,  1.      ],
[-0.21371 ,  0.69225 ,  1.      ],
[-0.375   ,  0.50219 ,  1.      ],
[-0.51325 ,  0.46564 ,  1.      ]])

plotData(data2,'Test 1','Test 2','y=1','y=0')

l(θ)=1mi=1m[y(i)log(hθ(x(i)))+(1y(i))log(1hθ(x(i)))]+λ2mj=1nθ2j

l(θ)=1m((log(g(Xθ))Ty+(log(1g(Xθ))T(1y))+λ2mj=1nθ2j

#做一个6阶多项式
from sklearn.preprocessing import PolynomialFeatures
poly=PolynomialFeatures(6)
xx=poly.fit_transform(data2[:,0:2])
xx.shape
(118, 28)

#定义似然函数
def likelihoodFunctionReg(theta,reg,*args):
m=y.size
h=sigmoid(xx.dot(theta))
l=-1.0*(1.0/m)*(np.log(h).T.dot(y)+np.log(1-h).T.dot(1-y))+(reg/(2*m))*np.sum(np.square(theta[1:]))
if np.isnan(l[0]):
return (np.inf)
return (l[0])

δl(θ)δθj=1mi=1m(y(i)hθ(x(i)))x(i)j+λmθj

δl(θ)δθj=1mXT(yg(Xθ))+λmθj

#梯度函数
m=y.size
h=sigmoid(xx.dot(theta.reshape(-1,1)))
return (grad.flatten())
#标签
y=np.c_[data2[:,2]]
initial_theta=np.zeros(xx.shape[1])
likelihoodFunctionReg(initial_theta,1,xx,y)
0.6931471805599453

fig,axes=plt.subplots(1,3,figsize=(17,5))
'''

lambda=0没有正则化
lambda=1：正常
lambda=100：过大
'''
for i,C in enumerate([0.0,1.0,100.0]):
#准确率
accuracy=100.0*sum(predict(res2.x,xx)==y.ravel())/y.size

plotData(data2,'Test1','Test2','y=1','y=0',axes.flatten()[i])
#画出决策边界
x1_min,x1_max=xx[:,1].min(),xx[:,1].max()
x2_min,x2_max=xx[:,2].min(),xx[:,2].max()
XX1,XX2=np.meshgrid(np.linspace(x1_min,x1_max),np.linspace(x2_min,x2_max))
h=sigmoid(poly.fit_transform(np.c_[XX1.ravel(),XX2.ravel()]).dot(res2.x))
h=h.reshape(XX1.shape)
axes.flatten()[i].contour(XX1,XX2,h,[0.5],linewidths=1,colors='g')
axes.flatten()[i].set_title('Train accuracy{}% with lambda={}'.format(np.round(accuracy,decimals=2),C))

• 广告
• 抄袭
• 版权
• 政治
• 色情
• 无意义
• 其他

120