从0手撸Gamma回归

模型公式

线性回归、逻辑回归、泊松回归。。。各种各样的回归可以抽象出一个通用线性回归

通用线性回归模型 [Generalized Linear Models (GLM)]：

$\hat{y}(w,X)=h(Xw)$

不同的回归模型所对应的函数$h$各不相同

标准的线性回归的$h$函数则为常数1：

${\hat{y}}_{linear}=Xw$

Gamma回归所对应的函数$h$则为$e$的${\hat{y}}_{linear}$次幂，完整的Gamma模型公式如下：

$\hat{y}=e^{Xw}$

损失函数

同样各种各样的回归都能抽象出来一个通用的损失函数，具体公式如下：

${\min }_w\frac{1}{2n_{\text{samples}}}{\sum }_{i}d(y_i,{\hat{y}}_i)+\frac{\alpha }{2}||w|{|}_2^2$

• 其中，$||w|{|}_2$为L2范数，$\alpha$则为L2的惩罚项

不同的回归模型所对应的函数$d$各不相同

标准的线性回归函数$d$则为，

$d(y_i,{\hat{y}}_i)=(y-\hat{y}{)}^2$

且$\alpha =0$

而Gamma回归模型所对应的函数$d$则为：

$d(y_i,{\hat{y}}_i)=2(\log \frac{\hat{y}}{y}+\frac{y}{\hat{y}}-1)$

代码实现Loss以及求导

当前章节则先用numpy实现求loss以及梯度 再借用Pytorch的自动求导机制实现求导 从而验证结果是否正确

import torch
import numpy as np

# 定义loss以及梯度的计算
def grad(coef, X, y, y_pred, alpha=1):
    # coef第0位存放bias
    offset = 1
    # 计算标准线性回归 WX + bais
    lin_pred = X @ coef[offset:] + coef[0]
    # 计算Gamma回归 基于标准线性回归 e^{WX + bais}
    y_pred = np.exp(lin_pred)

    # 计算梯度部分
    d1 = np.exp(lin_pred)
    coef_scaled = alpha * coef[offset:]
    temp = d1 * -2 * (y - y_pred) / np.power(y_pred, 2)
    devp = np.concatenate(([temp.sum()], temp @ X))
    grad = 0.5 * devp
    grad[offset:] += coef_scaled

    return grad

# 计算Loss部分
def loss(coef, y, y_pred, alpha=1):
    offset = 1
    # coef的第0位存放bias 详看后边完整代码
    coef_scaled = alpha * coef[offset:]
    # Unit Deviance d(y, y^hat)
    dev = np.sum(2 * (np.log(y_pred / y) + y / y_pred - 1))
    # 完整的loss函数 【The minimization problem】
    loss = 0.5 * dev + 0.5 * (coef[offset:] @ coef_scaled)

    return loss

X = np.array([[1, 2], [2, 3], [3, 4], [4, 3]])
y = np.array([19, 26, 33, 30])

"""numpy的实现"""
# 初始化模型
coef = np.zeros(X.shape[1] + 1)
bias = np.log(np.average(y))
coef[0] = bias
# 前传
linear_pred = X @ coef[1:] + coef[0]
gamma_pred = np.exp(linear_pred)
# 求loss
loss_np = loss(coef, y, gamma_pred, 1)
# 求梯度
grad_np = grad(coef, X, y, gamma_pred, 1)[1:]
print(f"Loss: {loss_np} Grad: {grad_np}")
# 由于用普通的SGD无法实现Gamma回归的迭代模型 故略去不写
# 详见后面完整代码

"""Pytorch实现进行验证"""
X = torch.from_numpy(X).double()
y = torch.from_numpy(y).double()
coef = torch.from_numpy(coef)
coef.requires_grad = True
# 前传
linear_pred = X @ coef[1:] + coef[0]
gamma_pred = torch.exp(linear_pred)
# 求loss
coef_scaled = 1 * coef[1:]
dev = (2 * (torch.log(gamma_pred / y) + y / gamma_pred - 1)).sum()
loss_pt = 0.5 * dev + 0.5 * (coef[1:] @ coef_scaled)
# 反向传播
loss_pt.backward()
# 查看loss以及梯度
grad_pt = coef.grad[1:].numpy()
print(f"Loss: {loss_pt} Grad: {grad_pt}")

assert loss_np == loss_pt
assert grad_np.all() == grad_pt.all()

完整代码实现

import numpy as np
from scipy.optimize import minimize


X = np.array([[1, 2], [2, 3], [3, 4], [4, 3]])
y = np.array([19, 26, 33, 30])


# 实现gamma回归
class GammaRegressor:

    def __init__(self, alpha=1):
        # L2范数惩罚系数
        self.alpha = alpha

    def init_model(self, X, y):
        _, n_features = X.shape

        # 生成 Weigh 与 Bias, 存储为coef 第0个为Bias
        coef = np.zeros(n_features + 1)
        coef[0] = np.log(np.average(y))

        return coef

    def show_coef(self, coef):
        print(f"[INFO] 当前模型: {coef}")

    def fit(self, X, y):
        """训练"""

        # 初始化函数
        coef = self.init_model(X, y)

        # 定义loss以及梯度的计算
        def cal_loss_and_grad(coef, obj, X, y):

            offset = 1
            # 计算标准线性回归 WX + bais
            lin_pred = X @ coef[offset:] + coef[0]
            # 计算Gamma回归 基于标准线性回归 e^{WX + bais}
            y_pred = np.exp(lin_pred)

            # 计算Loss部分
            coef_scaled = obj.alpha * coef[offset:]
            # Unit Deviance d(y, y^hat)
            dev = np.sum(2 * (np.log(y_pred / y) + y / y_pred - 1))
            # 完整的loss函数 【The minimization problem】
            loss = 0.5 * dev + 0.5 * (coef[offset:] @ coef_scaled)

            # 计算梯度部分
            d1 = np.exp(lin_pred)
            temp = d1 * -2 * (y - y_pred) / np.power(y_pred, 2)
            devp = np.concatenate(([temp.sum()], temp @ X))
            grad = 0.5 * devp
            grad[offset:] += coef_scaled
            print(f"[INFO] 当前梯度: {grad}")
            print(f"[INFO] 当前loss: {loss}")

            return loss, grad

        # 开始训练 优化方法采用 L-BFGS-B 采用SGD会梯度爆炸无法收敛
        args = (self, X, y)
        opt_res = minimize(
            cal_loss_and_grad,
            coef,
            method="L-BFGS-B",
            jac=True,
            options={
                "maxiter": 5,
                "disp": True
            },
            args=args,
            callback=self.show_coef
        )
        coef = opt_res.x

        return coef


gmr = GammaRegressor()
gmr.fit(X, y)