正则化逻辑回归(2)

正则化逻辑回归

题目:

你是工厂主管,你有一些芯片在两次测试的结果,你需要决定它们是否合格,你手里有之前的结果,现在需要构建一个逻辑回归模型进行预测

代码:

1. 读取数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
path = 'ex2data2.txt'
df = pd.read_csv(path, header=None, names=['Test1', 'Test2', 'Accepted'])
df.head()
df.describe()

2. 绘制样本图像

positive = df[df['Accepted'].isin([1])]
negative = df[df['Accepted'].isin([0])]
fig, ax = plt.subplots(figsize=(12, 8))
ax.scatter(positive['Test1'], positive['Test2'], s=50, c='b', marker='o', label='Accepted')
ax.scatter(negative['Test1'], negative['Test2'], s=50, c='r', marker='x', label='Rejected')
ax.legend()
ax.set_xlabel('Test1 Score')
ax.set_ylabel('Test2 Score')
plt.show()

3. 运用特征映射 

特征映射是卷积神经网络CNN里的知识,我暂时跳过了

简单了解就是:

特征映射的目的是从原始输入数据中提取出有用的信息,并将其转化为一种适合后续处理(如分类、回归等)的形式

def feature_mapping(x, y, power, as_ndarray=False):
    data = {'f{0}{1}'.format(i-p, p): np.power(x, i-p) * np.power(y, p)
                for i in range(0, power+1)
                for p in range(0, i+1)
           }
    if as_ndarray:
        return pd.DataFrame(data).values
    else:
        return pd.DataFrame(data)
x1 = df.Test1.values
x2 = df.Test2.values
Y = df.Accepted
data = feature_mapping(x1, x2, power=6)
# data = data.sort_index(axis=1, ascending=True)
data.head()
data.describe()

 4. 使用代价函数

与上图即上次作业不同的是多了一个正则项

就是对于θ平方求和再乘了一个λ权重

 4.1 调用特征映射
theta = np.zeros(data.shape[1])
X = feature_mapping(x1, x2, power=6, as_ndarray=True)
X.shape, Y.shape, theta.shape
((118, 28), (118,), (28,))
 4.2 调用sigmoid函数,将模型的输出转换为概率值
def sigmoid(z):
    return 1 / (1 + np.exp(-z))
4.3 使用代价函数

l不是1,是字母l,代表的是λ

def cost(theta, X, Y):
    first = Y * np.log(sigmoid(X@theta.T))
    second = (1 - Y) * np.log(1 - sigmoid(X@theta.T))
    return -1 * np.mean(first + second)
def regularized_cost(theta, X, Y, l=1):
    theta_1n = theta[1:]
    regularized_term = l / (2 * len(X)) * np.power(theta_1n, 2).sum()
    return cost(theta, X, Y) + regularized_term

θ为偏置,一般不参加运算,所以上面 theta_1n = theta[1:],意思是从第一个开始,没有算第0个,也就是偏置

意思就是从第二个开始求和

而由于我们将θ设置为0,所以这个正则化代价函数与代价函数的值是相同的

 

?为什么要正则化代价函数

  • 正则化代价函数和传统的代价函数之间的主要区别在于正则化代价函数在原有的代价函数基础上增加了一个正则化项。这个正则化项的目的是为了限制模型的复杂度,防止模型过拟合,并帮助模型更好地泛化到未见过的数据。
    • L1正则化:L1正则化(也称为Lasso回归)通过在代价函数中添加权重的绝对值之和来惩罚模型复杂度。这有助于产生稀疏的权重矩阵,即许多权重为零。
    • L2正则化:L2正则化(也称为岭回归)通过在代价函数中添加权重的平方和来惩罚模型复杂度。这有助于减小权重的幅度,从而防止模型过于复杂。
    • 本题用到的是L2正则化,所以不用乘二分之一
    • 因为:L2正则化是在代价函数后面再加上一个正则化项λ,使得权重在更新的时候,乘以一个小于1的因子(1-a(λ/m)),这个可以防止W过大。正则化项里面有一个系数1/2,主要是为了后面求导的结果方便,因为后面那一项求导会产生一个2,与1/2相乘刚好凑整。所以,L2正则化本身并不包括二分之一这个系数。

5. 正则化梯度 

def gradient(theta, X, Y):
    return (1/len(X) * X.T @ (sigmoid(X @ theta.T) - Y))
def regularized_gradient(theta, X, Y, l=1):
    theta_1n = theta[1:]
    regularized_theta = l / len(X) * theta_1n
#     regularized_theta[0] = 0
#     np.concatenate是连接函数
    regularized_term = np.concatenate([np.array([0]), regularized_theta])
    
    return  gradient(theta, X, Y) + regularized_term
#     return  gradient(theta, X, Y) + regularized_theta

1. theta_1n = theta[1:]

这行代码将theta数组(通常是一个包含模型参数的向量)中的第一个元素去掉,并将其余部分赋值给theta_1n。假设theta是一个长度为n+1的向量,那么theta_1n就是一个长度为n的向量。

2. regularized_theta = l / len(X) * theta_1n

这行代码计算正则化项regularized_theta。这里,l是正则化参数(通常是一个很小的正数,用于控制正则化的强度),X是特征矩阵,len(X)是特征矩阵的行数(即样本数量)。theta_1n是除了第一个参数以外的所有参数。

正则化项的作用是对模型的复杂度进行惩罚,以防止过拟合。在这个例子中,正则化项与除了截距项(通常是theta的第一个元素)以外的所有参数成正比。

3. # regularized_theta[0] = 0

这行代码被注释掉了,但如果取消注释,它会将regularized_theta数组的第一个元素设置为0。这意味着截距项(通常是theta的第一个元素)不会被正则化。然而,由于上一行代码已经计算了regularized_theta,这行代码实际上是不必要的,因为regularized_theta已经不包括截距项了。

4. regularized_term = np.concatenate([np.array([0]), regularized_theta])

这行代码创建一个新的数组regularized_term,它通过将一个包含单个元素0的数组和regularized_theta数组连接起来而成。这样,regularized_term就变成了一个与原始theta数组具有相同长度的数组,其中第一个元素为0(对应于截距项),其余元素与regularized_theta相同。

总结起来,这段代码的目的是计算正则化项,并将其与原始参数向量结合,以创建一个包含正则化参数的新向量。

6. 拟合参数

此时通过使用Scipy.optimize.minimize拟合出了最优的θ(与上面0.6形成对比)

import scipy.optimize as opt
# func=cost,就是传的损失函数;x0初始点;args是传样本的输入和输出;jac传的梯度
res = opt.minimize(fun=regularized_cost, x0=theta, args=(X, Y), method='Newton-CG', jac=regularized_gradient)
res
 fun: 0.5290027297127388
     jac: array([ 7.26089191e-08,  4.22913232e-09,  8.15815876e-09,  6.15699190e-08,
        7.74567232e-09, -3.09360466e-08,  2.12821347e-08,  1.22156735e-08,
        1.96058084e-08, -3.19108791e-08, -4.39405717e-09, -2.76847096e-09,
       -2.77934021e-08,  1.23592858e-08, -7.14474161e-08,  8.98276579e-09,
        1.45962365e-08, -1.00120216e-08, -7.32796823e-09,  1.43317535e-08,
       -4.38679455e-08, -4.85023121e-09, -3.40732357e-10, -1.11668147e-08,
       -5.01047274e-09, -1.44326742e-08,  8.78794915e-09, -5.71951122e-08])
 message: 'Optimization terminated successfully.'
    nfev: 7
    nhev: 0
     nit: 6
    njev: 68
  status: 0
 success: True
       x: array([ 1.27273909,  0.62527214,  1.18108783, -2.01995993, -0.91742426,
       -1.43166279,  0.12400726, -0.36553444, -0.35723901, -0.17513021,
       -1.45815774, -0.05098947, -0.61555653, -0.27470644, -1.19281683,
       -0.24218793, -0.20600565, -0.04473137, -0.27778488, -0.2953778 ,
       -0.45635711, -1.04320321,  0.02777158, -0.29243198,  0.01556636,
       -0.32738013, -0.14388704, -0.92465213])

综上所述,这行代码的目的是使用牛顿共轭梯度方法从初始参数值 theta 开始,通过最小化正则化的代价函数 regularized_cost 来找到最优的参数值。这通常用于机器学习模型的训练,特别是当需要应用L2正则化来防止过拟合时。 

7. 预测分析 

def predict(theta, X):
    probability = sigmoid(X @ theta.T)
    return probability >= 0.5
    return [1 if x>=0.5 else 0 for x in probability]
#使用了一个方法
#support标签中出现的次数
#precision查准率,recall召回率,f1-score调和平均数
# final_theta = res.x
from sklearn.metrics import classification_report
Y_pred = predict(res.x, X)
print(classification_report(Y, Y_pred))
 precision    recall  f1-score   support
           0       0.90      0.75      0.82        60
           1       0.78      0.91      0.84        58
    accuracy                           0.83       118
   macro avg       0.84      0.83      0.83       118
weighted avg       0.84      0.83      0.83       118
真实类别
01
预测0TN(真反例)FN(假反例)
类别1FP(假正例)TP(真正例)

precision=tp / tp+fp  查准率 意思是:在所有正例中真正的正例占比,也就是准确率

recall=tp / tp+fn   召回率 意思是:你这个算法可以从真实的类别中回收回来多少正例

f1-score=2pr / p+r 调和平均数

8. 决策边界

#就是把上面一系列封装起来
# 得到最优的所有theta
#power就是用来控制高维特征的维度的
#l就是λ ,就是正则化项的那个权重
def find_theta(power, l):
    path = 'ex2data2.txt'
    df = pd.read_csv(path, header=None, names=['Test1', 'Test2', 'Accepted'])
    df.head()
    Y = df.Accepted
    x1 = df.Test1.values
    x2 = df.Test2.values
    X = feature_mapping(x1, x2, power, as_ndarray=True)
    theta = np.zeros(X.shape[1])
#     res = opt.minimize(fun=regularized_cost, x0=theta, args=(X, Y, l), method='Newton-CG', jac=regularized_gradient)
    res = opt.minimize(fun=regularized_cost, x0=theta, args=(X, Y, l), method='TNC', jac=regularized_gradient)
    return res.x
# 找决策边界,thetaX = 0, thetaX <= threshhold
def find_decision_boundary(density, power, theta, threshhold):
    #t1、t2是两个数组
    t1 = np.linspace(-1, 1.2, density)# 1000个样本
    t2 = np.linspace(-1, 1.2, density)
    #再将它们转化为坐标
    cordinates = [(x, y) for x in t1 for y in t2]
    #把坐标x和y独立出来,就是分隔开
    x_cord, y_cord = zip(*cordinates)
    #再特征映射(本来坐标是2维的,映射完就是高维的,也就是你本题theta的28维)
    mapped_cord = feature_mapping(x_cord, y_cord, power)
    
    pred = mapped_cord.values @ theta.T
    #abs取绝对值
    decision = mapped_cord[np.abs(pred) <= threshhold]
    
    return decision.f10, decision.f01
# 画决策边界
def draw_boundary(power, l):
    #1000个样本
    density = 1000
    # 阈值(就是0.002对应就是我们所说的ε)
    threshhold = 2 * 10**-3
    
    # 调用了上面的两个函数
    theta = find_theta(power, l)
    x, y = find_decision_boundary(density, power, theta, threshhold)
    
    # 画图
    positive = df[df['Accepted'].isin([1])]
    negative = df[df['Accepted'].isin([0])]
    fig, ax = plt.subplots(figsize=(12, 8))
    ax.scatter(positive['Test1'], positive['Test2'], s=50, c='b', marker='o', label='Accepted')
    ax.scatter(negative['Test1'], negative['Test2'], s=50, c='g', marker='x', label='Rejected')
    ax.scatter(x, y, s=50, c='r', marker='.', label='Decision Boundary')
    ax.legend()
    ax.set_xlabel('Test1 Score')
    ax.set_ylabel('Test2 Score')
    plt.show()

决策边界就是让令下面这个式子等于0

与之前不同的是,扩展到了高维(所以决策曲线就不光有一次项,即一条很显然的直线。它会更加复杂)

分析:

我们放宽条件(因为可能会有等于0.0001这样很接近的情况),所以,我们只需这个式子小于 ε 即可(等价于原来的等于0)👇

   具体实现这个的代码👇:

    #abs取绝对值

    decision = mapped_cord[np.abs(pred) <= threshhold]

np.linspace是NumPy库中的一个函数,用于生成等间隔的数值数组,也就是等差数列。等差数列是数学中的一个经典概念,指的是两个相邻的数之间的差是固定的。

使用np.linspace函数时,需要指定间隔的起始点、终止点,以及指定分隔值总数(包括起始点和终止点)。最终,函数会返回间隔类均匀分布的数值序列。

尽管arange函数也可以生成类似的结构,但与arange相比,许多人更愿意使用linspace,因为它更易于理解和使用。

 

     mapped_cord = feature_mapping(x_cord, y_cord, power)

    

    pred = mapped_cord.values @ theta.T

含义👇

 return decision.f10, decision.f01 

- 01就是原始的第一个特征

-10就是原始的第二个特征

返回的就是分界曲线上的x和y

 

 9. 通过调节正则项中 λ 权重的值,来调节拟合程度

拟合较好
draw_boundary(6, l=1)

 

 过拟合
draw_boundary(6, l=0)
欠拟合 
draw_boundary(6, l=100)