Python进行逻辑回归并给出OR值置信区间 逻辑回归正则化python

通过芯片在两种测试中的测试结果，建立正则化逻辑回归算法模型，来决定芯片是否通过测试（拥有过去芯片测试的数据集）

一、导入库和数据集

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['font.sans-serif']=['SimHei']#用黑体显示中文
matplotlib.rcParams['axes.unicode_minus']=False # 可以显示负号
df=pd.read_csv('C:/Users/Administrator/AppData/Local/Temp/Temp1_machine-learning-ex2.zip/machine-learning-ex2/ex2/ex2data2.txt',
              names=['测试一','测试二','接受'])

二、可视化

positive=df[df['接受'].isin([1])]
negetive=df[df['接受'].isin([0])]
fig,ax=plt.subplots(figsize=(8,5))
ax.scatter(positive['测试一'],positive['测试二'],label='接受')
ax.scatter(negetive['测试一'],negetive['测试二'],marker='x',label='未接受')
box=ax.get_position()
ax.set_position([box.x0,box.y0,box.width,box.height*0.8])
ax.legend(loc='center left',bbox_to_anchor=(0.2,1.12),ncol=3)
ax.set_xlabel('测试一')
ax.set_ylabel('测试二')
plt.show()

结果为

*由图可以发现，该数据集没有线性的决策边界，所以不能直接用logistic回归（theta.T X=0得到的是线性的决策边界）

三、特征映射

为了能很好的拟合，所以我们要进行特征映射，从而得到非线性的决策边
1.定义特征映射函数

def feature_mapping(x1,x2,power):
    data={}
    for i in np.arange(power+1):
        for p in np.arange(i+1):
            data["f{}{}".format(i-p,p)]=np.power(x1,i-p)*np.power(x2,p)
    return pd.DataFrame(data)

与下边的等价

def feature_mapping(x1,x2,power):
    data={"f{}{}".format(i-p,p):np.power(x1,i-p)*np.power(x2,p)
          for i in np.arange(power+1)
          for p in np.arange(i+1)}
    return pd.DataFrame(data)

2.将特征变量进行特征映射

x1=df['测试一']
x1=x1.values
x2=df['测试二']
x2=x2.values
df2=feature_mapping(x1,x2,6)

$x_1,x_2$经映射变为：$1,\quad x_1,x_2,\quad x_1^2,x_2^2,x_1x_2,\quad x_1^3,x_2^3,x_1^2x_2,x_1x_2^2,\quad x_1^4,x_2^4,x_1^2x_2^2,x_1^3x_2,x_1x_2^3,\quad x_1^5,x_2^5,x_1^4x_2,x_1x_2^4,x_1^3x_2^2,x_1^2x_2^3,\quad x_1^6,x_2^6,x_1^5x_2,x_1x_2^5,x_1^4x_2^2,x_1^2x_2^4,x_1^3x_2^3$

3.重新选取变量

X=df2.values  # 将特征映射后得到的变量定义为X
y=df['接受'].values
theta=np.zeros(X.shape[1])

注意：
特征映射将低维特征向量（本例中为二维）转化为高维特征向量（本例中为28维），因为特征向量的个数变多，所以需要用正则化，否则容易出现过拟合。

四、正则化logistic回归

正则化后的代价函数：

1.正则化后的代价函数

#定义sigmoid函数
def sigmoid(z):
    return 1/(1+np.exp(-z))
#定义代价函数
def cost(theta,X,y):
    first=[e-mail](sigmoid(X@theta))
    second=(1-y)@np.log(1-sigmoid(X@theta))
    return (first+second)/(-len(X))
#定义正则化代价函数
def costReg(theta,X,y,lam=1):
    _theta=theta[1:]
    reg=(lam/(2*len(X)))*(_theta@_theta)
    return cost(theta,X,y)+reg
 
costReg(theta,X,y)    # 0.6931471805599454

2.定义正则化后的梯度函数（偏导数）

def gradient(theta,X,y):
    return (X.T@(sigmoid(X@theta)-y))/(len(X))
def gradientReg(theta,X,y,lam=1):
    reg=(lam/len(X))*theta
    reg[0]=0    # 第一项没有惩罚因子
    return gradient(theta,X,y)+reg
(gradientReg(theta,X,y)).shape   #(28,)

3.使用高级优化算法

import scipy.optimize as opt
result=opt.fmin_tnc(func=costReg,x0=theta,fprime=gradientReg,args=(X,y,2))

与下边是等价的

opt.minimize(fun=costReg,x0=theta,args=(X,y,2),method='TNC',jac=gradientReg)

4.定义预测函数,将求得的参数finally_theta和X（特征映射后的）带入预测函数中，并计算预测精度

def predict(theta,X):
    probability=sigmoid(X@theta)
    return [1 if x>=0.5 else 0 for x in probability ]
finally_theta=result[0]
predictions=predict(finally_theta,X)
correct=[1 if a==b else 0 for (a,b) in zip(predictions,y)]
accuracy=sum(correct)/len(X)
accuracy     # 0.8305084745762712

或者用sklearn中的方法进行评估

from sklearn.metrics import classification_report 
print(classification_report(predictions,y))

结果如下：

5.决策边界

x = np.linspace(-1, 1.5, 250)
xx, yy = np.meshgrid(x, x)   #生成网格点坐标矩阵
# 得到的xx,yy是array形式，维度为(250,250)
z = feature_mapping(xx.ravel(), yy.ravel(), 6).values  # xx.ravel()将xx一维化(62500,)
z = z @ finally_theta  #z.shape (62500,)
z = z.reshape(xx.shape)

fig,ax=plt.subplots(figsize=(8,6))
ax.scatter(positive['测试一'],positive['测试二'],label='接受')
ax.scatter(negetive['测试一'],negetive['测试二'],marker='x',label='未接受')
# 设置图例显示在图的上方
box = ax.get_position()
ax.set_position([box.x0,box.y0,box.width,box.height*0.8])
ax.legend(loc='center left',bbox_to_anchor=(0.2,1.12),ncol=3)
ax.set_xlabel('测试一')
ax.set_ylabel('测试二')
ax.set_title('决策边界lambda为2')
plt.contour(xx,yy,z,0) # 生成等高线，0是指图像一分为二
plt.show()