【學習筆記】用numpy實現多項式邏輯回歸（PolynomialLogisticRegression）

多項式邏輯回歸就是在邏輯回歸的基礎上將高次項作為特征加進去，以實現高維特征的提取

一、模型構建

多項式邏輯回歸模型是由三個子模型組成：

（1）添加多項式特征

（2）標准化

（3）邏輯回歸

添加多項式特征

將各個特征之間相乘得到新的特征，比如原來的特征是\([x_0,x_1]\)

二次多項式特征是\([1,x_0,x_1,x_0^2,x_0x_1,x_1^2]\)

三次多項式特征是\([1,x_0,x_1,x_0^2,x_0x_1,x_1^2,x_0^3,x_0^2x_1,x_0x_1^2,x_1^3]\)

def poly_transform(X):
    XX=X.T
    ret=[np.repeat(1.0,XX.shape[1]),XX[0],XX[1]]
    for i in range(2,degree+1):
        for j in range(0,i+1):
            ret.append(XX[0]**(i-j)*XX[1]**j)
    return np.array(ret).T

超參數設置：

degree=3

標准化

將樣本縮放到均值為0，標准差為1，變換公式為\(\Large x=\frac{x-\mu}{\sigma}\)

坑點：

（1）標准差有可能為0，轉換的時候會出現分母為0的情況，需要將其轉化為非零的數

（2）只需要對訓練數據fit一遍，測試的時候按照訓練時的均值和標准差變換就行了，否則會出現偏差

def scaler_fit(X):
    global mean,scale
    mean=X.mean(axis=0)
    scale=X.std(axis=0)
    scale[scale<np.finfo(scale.dtype).eps]=1.0
def scaler_transform(X):
    return (X-mean)/scale

邏輯回歸

模型的核心部分，訓練過程由前向傳播、計算損失和反向傳播三部分構成：

前向傳播

由輸入數據計算輸出概率\(p(x)\)

\[p(x)=\sigma(w^Tx) \]

def sigmoid(x):
    return 1/(1+exp(-x))
def forward(x):
    return sigmoid(w@x.T)

從意義上來說\(p(x)\)代表\(x\)被判別為第1類的概率
按理說可以再加上一個偏置項\(b\)，對有些數據效果可能會更好一些

計算損失

假設上一步已經得到了訓練數據\(x\)的概率\(p=p(x)\)

對於單個訓練數據\([p,y]\)損失函數為：

\[J=-[y\log p+(1-y)\log(1-p)] \]

采用MBGD，一次取多個訓練數據，對於一組具有m個訓練數據的batch有：

\[J=-\frac{1}{m}\sum_{i=0}^{m-1}[y_i\log p_i+(1-y_i)\log(1-p_i)] \]

為防止\(w\)的權值過度膨脹（因為\(w\)各項權值的絕對值越大，\(p\)的值越接近0或1），加入L1正則項，因此最終的損失函數為：

\[J=-\frac{1}{m}\sum_{i=0}^{m-1}[y_i\log p_i+(1-y_i)\log(1-p_i)]+k||w||_1 \]

def compute_loss(p,y):
    return np.mean(-(y*log(p)+(1-y)*log(1-p)),axis=0)+k*np.linalg.norm(w,ord=1)

反向傳播+更新參數

計算損失函數\(J\)對\(w\)各個權值的偏導數

為方便思考，先只考慮單個數據\([p,y]\)時的情況：

\[\frac{\part J}{\part p}=-(\frac{y}{p}-\frac{1-y}{1-p})\\ \frac{\part p}{\part(w^Tx)}=p(1-p)\\ \frac{\part(w^Tx)}{\part w_i}=x_i \]

因此有

\[\frac{\part{J}}{\part{w_i}} =\frac{\part{J}}{\part{p}}\frac{\part{p}}{\part(w^Tx)}\frac{\part(w^Tx)}{\part w_i} =(p-y)x_i \]

對於m個訓練數據：

\[\frac{\part{J}}{\part{w_i}}=\frac{1}{m}\sum_{j=0}^{m-1}(p_j-y_j)x_{ji} \]

加入正則項之后：

\[\frac{\part{J}}{\part{w_i}}=\frac{1}{m}\sum_{j=0}^{m-1}(p_j-y_j)x_{ji}+k\text{sign}(w_i) \]

寫成向量的形式就是：

\[\frac{\nabla J}{\nabla w}=\frac{1}{m}\sum_{j=0}^{m-1}(p_j-y_j)x_{j}+k\text{sign}(w) \]

在反向傳播的同時更新參數，即：

\[w=w-\alpha\frac{\nabla J}{\nabla w} \]

\(\alpha\)代表學習率

def backward(p,x,y):
    global w
    w-=learning_rate*np.mean((p-y).reshape(-1,1)*x,axis=0)+k*sign(w)

訓練

采用MBGD進行訓練，設學習率為learning_rate，正則項權重為k，一個批次的訓練數據規模為batch_size，迭代輪次為epoch_num:

def logistic_fit(X,Y):
    global w
    XX=X.copy()
    YY=Y.copy()
    w=random.randn(XX[0].shape[0])*0.01
    I=list(range(len(XX)))
    LOSS=[]
    for epoch in range(epoch_num):
        random.shuffle(I)
        XX=XX[I]
        YY=YY[I]
        for i in range(0,len(XX),batch_size):
            x=XX[i:i+batch_size]
            y=YY[i:i+batch_size]
            p=forward(x)
            loss=compute_loss(p,y)
            LOSS.append(loss)
            backward(p,x,y)
    return LOSS

超參數設置：

learning_rate=0.9
k=0.005
batch_size=32
epoch_num=50

二、生成訓練數據

生成訓練數據集\(X\)和\(Y\)，其中\(X=\{x\}\)，\(Y=\{y\}\)，\(x=[x_0,F(x_0)]\)，\(x_0∈[l,r]\)，標簽為\(y\)，數據集大小為n

def generate_data(F,l,r,n,y):
    x=np.linspace(l,r,n)
    X=np.column_stack((x,F(x)))
    Y=np.repeat(y,n)
    return X,Y

三、模型訓練及預測

訓練

輸入訓練樣本\((X,Y)\)，更新scaler及logisticregression的參數並繪制出loss曲線，無返回值

def train(X,Y):
    XX=poly_transform(X)
    scaler_fit(XX)
    XX=scaler_transform(XX)
    LOSS=logistic_fit(XX,Y)
    plt.plot(list(range(len(LOSS))),LOSS,color='r')
    plt.show()

預測

輸入預測樣本\(X\)，返回預測結果\(\hat Y\)

def predict(X):
    XX=scaler_transform(poly_transform(X))
    return logistic_predict(XX)

繪制判別界面

輸入預測樣本\((X,Y)\)，散點圖、及判別界面

def plot_decision_boundary(X,Y):
    global predY1
    x0_min,x0_max=X[:,0].min()-1,X[:,0].max()+1
    x1_min,x1_max=X[:,1].min()-1,X[:,1].max()+1
    m=500
    x0,x1=np.meshgrid(
        np.linspace(x0_min,x0_max,m),
        np.linspace(x1_min,x1_max,m)
    )
    XX=np.c_[x0.ravel(),x1.ravel()]
    Y_pred=predict(XX)
    Z=Y_pred.reshape(x0.shape)
    plt.contourf(x0,x1,Z,cmap=plt.cm.Spectral)
    plt.scatter(X[:,0],X[:,1],c=Y,cmap=plt.cm.Spectral)
    plt.show()

四、主程序

random.seed(114514)
data_size=200
X1,Y1=generate_data(lambda x:x**2+2*x-2+random.randn(data_size)*0.8,-3,1,data_size,0)
X2,Y2=generate_data(lambda x:-x**2+2*x+2+random.randn(data_size)*0.8,-1,3,data_size,1)
X=np.vstack((X1,X2))
Y=np.hstack((Y1,Y2))
train(X,Y)
plot_decision_boundary(X,Y)

loss曲線：

判別界面：

五、完整代碼

導入依賴包

import numpy as np
from numpy import exp,log,sign,random
from matplotlib import pyplot as plt

模型構建

# LogisticRegression
learning_rate=0.9
k=0.005
batch_size=32
epoch_num=50
def sigmoid(x):
    return 1/(1+exp(-x))
def forward(x):
    return sigmoid(w@x.T)
def compute_loss(p,y):
    return np.mean(-(y*log(p)+(1-y)*log(1-p)),axis=0)+k*np.linalg.norm(w,ord=1)
def backward(p,x,y):
    global w
    w-=learning_rate*np.mean((p-y).reshape(-1,1)*x,axis=0)+k*sign(w)
def logistic_fit(X,Y):
    global w
    XX=X.copy()
    YY=Y.copy()
    w=random.randn(XX[0].shape[0])*0.01
    I=list(range(len(XX)))
    LOSS=[]
    for epoch in range(epoch_num):
        random.shuffle(I)
        XX=XX[I]
        YY=YY[I]
        for i in range(0,len(XX),batch_size):
            x=XX[i:i+batch_size]
            y=YY[i:i+batch_size]
            p=forward(x)
            loss=compute_loss(p,y)
            LOSS.append(loss)
            backward(p,x,y)
    return LOSS
def logistic_predict(X):
    return np.where(forward(X)>0.5,1,0)
    
# StandardScaler
def scaler_fit(X):
    global mean,scale
    mean=X.mean(axis=0)
    scale=X.std(axis=0)
    scale[scale<np.finfo(scale.dtype).eps]=1.0
def scaler_transform(X):
    return (X-mean)/scale

# PolynomialFeatures
degree=3
def poly_transform(X):
    XX=X.T
    ret=[np.repeat(1.0,XX.shape[1]),XX[0],XX[1]]
    for i in range(2,degree+1):
        for j in range(0,i+1):
            ret.append(XX[0]**(i-j)*XX[1]**j)
    return np.array(ret).T

模型訓練及預測

def train(X,Y):
    XX=poly_transform(X)
    scaler_fit(XX)
    XX=scaler_transform(XX)
    LOSS=logistic_fit(XX,Y)
    plt.plot(list(range(len(LOSS))),LOSS,color='r')
    plt.show()
def predict(X):
    XX=scaler_transform(poly_transform(X))
    return logistic_predict(XX)
def plot_decision_boundary(X,Y):
    global predY1
    x0_min,x0_max=X[:,0].min()-1,X[:,0].max()+1
    x1_min,x1_max=X[:,1].min()-1,X[:,1].max()+1
    m=500
    x0,x1=np.meshgrid(
        np.linspace(x0_min,x0_max,m),
        np.linspace(x1_min,x1_max,m)
    )
    XX=np.c_[x0.ravel(),x1.ravel()]
    Y_pred=predict(XX)
    Z=Y_pred.reshape(x0.shape)
    plt.contourf(x0,x1,Z,cmap=plt.cm.Spectral)
    plt.scatter(X[:,0],X[:,1],c=Y,cmap=plt.cm.Spectral)
    plt.show()

生成訓練數據

def generate_data(F,l,r,n,y):
    x=np.linspace(l,r,n)
    X=np.column_stack((x,F(x)))
    Y=np.repeat(y,n)
    return X,Y

主程序

random.seed(114514)
data_size=200
X1,Y1=generate_data(lambda x:x**2+2*x-2+random.randn(data_size)*0.8,-3,1,data_size,0)
X2,Y2=generate_data(lambda x:-x**2+2*x+2+random.randn(data_size)*0.8,-1,3,data_size,1)
X=np.vstack((X1,X2))
Y=np.hstack((Y1,Y2))
train(X,Y)
plot_decision_boundary(X,Y)

參考資料

• [1] 機器學習：邏輯回歸（使用多項式特征）
• [2] 深度學習PyTorch極簡入門教程