广义线性模型

广义线性模型（generalized linear model，GLM）用一个概率分布族描述响应的随机波动，再用连接函数把响应的条件均值与特征的线性组合联系起来。它既规定均值如何变化，也规定在均值给定后数据如何波动；“广义”并不表示可以任意选择一个非线性函数而忽略分布。

考虑三个合成观测：特征x为0、1、2，对应计数y为1、2、4。若采用泊松响应和对数连接，均值可按固定倍数随x增加，始终保持正数。这个例子用于展示GLM的组成、似然推导与实际迭代，数据规模本身不足以进行可靠的实际科学推断。

随机部分、线性预测与连接

给定固定设计向量 $x_i\in {ℝ}^p$，假设 $Y_1,\dots ,Y_n$ 条件于设计相互独立，密度或质量函数相对于共同参考测度写成 $$f(y_i;{\theta }_i,\varphi )=\exp \{\frac{y_i{\theta }_i-b({\theta }_i)}{a_i(\varphi )}+c(y_i,\varphi )\},a_i(\varphi )>0.$$ 这里θ为自然参数，φ为离散或尺度参数，b、a、c为已知函数。支撑不随θ变化，并假设在所讨论的参数内部可将归一化积分求导。对 $\int f\mathrm{d}y=1$ 求一次导数，有 $$0=𝔼\frac{Y_i-b^{\prime }({\theta }_i)}{a_i(\varphi )},\text{所以 }{\mu }_i=𝔼Y_i=b^{\prime }({\theta }_i).$$ 再求一次导数，得到 $$0=𝔼[\frac{(Y_i-b^{\prime }({\theta }_i){)}^2}{a_i(\varphi {)}^2}-\frac{b^{″}({\theta }_i)}{a_i(\varphi )}],\mathrm{Var}(Y_i)=a_i(\varphi )b^{″}({\theta }_i).$$ 这说明均值与方差来自同一个归一化分布，不能独立随意指定。

系统部分是 ${\eta }_i=x_i^{𝖳}\beta$，连接关系是 $g({\mu }_i)={\eta }_i$，g在有效均值区间上一一对应且足够可微。常数截距可作为x的第一列一。若连接使 ${\eta }_i={\theta }_i$，称为规范连接。正态均值的恒等连接、伯努利均值的logit连接、泊松均值的对数连接是不同的GLM组合；同一个分布也可考虑其他有效连接，选取它是一项建模决定。

泊松对数连接的得分方程

泊松质量为 $e^{-{\mu }_i}{\mu }_i^{y_i}/y_i!$，其中 $y_i\in {\mathbb{Z}}_{\ge 0}$。取 ${\theta }_i=\log {\mu }_i$、$b({\theta }_i)=e^{{\theta }_i}$、a为一，便有 $\mathrm{Var}(Y_i)={\mu }_i$。对数连接给出 ${\mu }_i=e^{x_i^{𝖳}\beta }$，所以独立样本的对数似然为 $$\ell (\beta )=\sum _i[y_i{x}_i^{𝖳}\beta -e^{x_i^{𝖳}\beta }-\log (y_i!)].$$ 令X的第i行为 $x_i^{𝖳}$，$W=\mathrm{diag}({\mu }_i)$，逐项求导得 $$U(\beta )=X^{𝖳}(y-\mu ),H(\beta )=-X^{𝖳}WX.$$ 设计满列秩且均值为正时，对任意非零v有 $v^{𝖳}{X}^{𝖳}WXv=\sum _i{\mu }_i(x_i^{𝖳}v{)}^2>0$，所以对数似然严格凹，有限驻点若存在就唯一且全局最优。

Newton步解 $X^{𝖳}WX\delta =U$，再令 ${\beta }_{\mathrm{n}\mathrm{e}\mathrm{w}}=\beta +\delta$。它也可写成加权最小二乘：令工作响应 $z=X\beta +(y-\mu )/\mu$（除法逐分量），则 $X^{𝖳}WX{\beta }_{\mathrm{n}\mathrm{e}\mathrm{w}}=X^{𝖳}Wz$。每一轮的权重和工作响应都依赖当前均值，这就是本例的迭代加权最小二乘。大步可能降低似然，实际实现可缩短步长。

三点计数的完整迭代

主例 $x_i^{𝖳}=(1,x_i)$，从 ${\beta }^{(0)}=(\log (7/3),0)$ 开始，此时三个均值均为7/3。得分为 $U=(0,3{)}^{𝖳}$，信息矩阵为 $$X^{𝖳}WX=\frac{7}{3}\left(\begin{array}{cc}3& 3\\ 3& 5\end{array}\right).$$ 解二元方程给出 $\delta =(-9/14,9/14)$，所以第一步参数约为 $(0.204440718,0.642857143)$，预测均值约为 $(1.226838724,2.333333333,4.437783337)$。

第二步约为 $(0.018902211,0.686213164)$，均值约为 $(1.019081989,2.024080200,4.020187482)$；再迭代趋向 $(0,\log 2)$，恰给出均值1、2、4。最后一点可独立验证：该参数使每个 $y_i-{\mu }_i=0$，所以得分为零；严格凹性证明它是唯一有限最大值。这里精确拟合是特意选取的教学数据，不是一般GLM都能逐点拟合。

点为合成计数1、2、4；最终均值 exp(x log2) 在本例恰经过它们。

Python 3

import math

def poisson_logline(x, counts, tol=1e-10, max_iter=100):
    x, y = tuple(x), tuple(counts)
    if len(x) != len(y) or len(x) < 2:
        raise ValueError("at least two paired observations required")
    if any(not math.isfinite(t) for t in x) or min(x) == max(x):
        raise ValueError("finite, nonconstant predictor required")
    if any(type(t) is not int or t <= 0 for t in y):
        raise ValueError("this teaching solver requires positive integer counts")
    if not math.isfinite(tol) or tol <= 0 or type(max_iter) is not int or max_iter < 1:
        raise ValueError("invalid convergence controls")
    try:
        xmean = math.fsum(x) / len(x)
        centered = [t - xmean for t in x]
        if not math.isfinite(xmean) or any(not math.isfinite(t) for t in centered):
            raise ArithmeticError("predictor centering exceeds floating range")
        xscale = max(abs(t) for t in centered)
        if xscale == 0:
            raise ArithmeticError("predictor scale is not representable")
        x = tuple(t / xscale for t in centered)
    except OverflowError as error:
        raise ArithmeticError("predictor preprocessing exceeds floating range") from error
    beta = [math.log(sum(y) / len(y)), 0.0]

    def original_coefficients(b):
        slope = b[1] / xscale
        intercept = b[0] - slope * xmean
        if any(not math.isfinite(v) for v in (slope, intercept)):
            raise ArithmeticError("original coefficients exceed floating range")
        if b[1] != 0 and slope == 0:
            raise ArithmeticError("original slope underflows")
        return intercept, slope

    def objective(b):
        try:
            eta = [b[0] + b[1] * t for t in x]
            value = sum(yi * z - math.exp(z) for yi, z in zip(y, eta))
            return value if math.isfinite(value) else -math.inf
        except OverflowError:
            return -math.inf

    for _ in range(max_iter):
        try:
            mu = [math.exp(beta[0] + beta[1] * t) for t in x]
        except OverflowError as error:
            raise ArithmeticError("mean exceeds floating range") from error
        if any(not math.isfinite(v) or v <= 0 for v in mu):
            raise ArithmeticError("mean is outside positive finite floating range")
        g0 = sum(yi - mi for yi, mi in zip(y, mu))
        g1 = sum(t * (yi - mi) for t, yi, mi in zip(x, y, mu))
        if any(not math.isfinite(g) for g in (g0, g1)):
            raise ArithmeticError("score exceeds floating range")
        if max(abs(g0), abs(g1)) <= tol * (1 + sum(y)):
            return original_coefficients(beta)
        h00 = sum(mu)
        h01 = sum(t * mi for t, mi in zip(x, mu))
        h11 = sum(t * t * mi for t, mi in zip(x, mu))
        det = h00 * h11 - h01 * h01
        if any(not math.isfinite(v) for v in (h00, h01, h11, det)):
            raise ArithmeticError("information matrix exceeds floating range")
        if det <= 1e-14 * h00 * h11:
            raise ArithmeticError("ill-conditioned design: rescale predictor")
        step = ((h11 * g0 - h01 * g1) / det,
                (h00 * g1 - h01 * g0) / det)
        old = objective(beta)
        if not math.isfinite(old) or any(not math.isfinite(v) for v in step):
            raise ArithmeticError("Newton calculation exceeds floating range")
        length = 1.0
        while length >= 2 ** -40:
            candidate = [b + length * d for b, d in zip(beta, step)]
            rounding = 1e-14 * (1 + abs(old))
            if objective(candidate) >= old + 1e-4 * length * (g0 * step[0] + g1 * step[1]) - rounding:
                beta = candidate
                break
            length *= 0.5
        else:
            raise ArithmeticError("line search failed")
    raise RuntimeError("iteration limit reached")

b0, b1 = poisson_logline([0, 1, 2], [1, 2, 4])
assert abs(b0) < 1e-7 and abs(b1 - math.log(2)) < 1e-7
assert all(abs(math.exp(b0 + b1 * x) - y) < 1e-6
           for x, y in zip([0, 1, 2], [1, 2, 4]))
print(b0, b1)

这份Python实现只拟合一个有限、非恒定特征加截距的泊松对数模型，要求正整数计数，以保持教学例子的有限最优解。它先将特征中心化并按最大绝对偏差缩放，在该无量纲坐标上检查得分收敛，再还原原单位的截距和斜率，避免仅因特征单位很小就提前停止。一般泊松模型允许零；若所有计数都为零，截距向负无穷时似然才达到上确界，根本没有有限MLE，不能用迭代停止值冒充解。程序显式报告病态设计、超出浮点范围、线搜索或迭代失败，不改写输入；二元方程逐步解析求解，每轮成本为 $O(n)$、存储为 $O(n)$，总成本乘实际迭代轮数。它没有包含任意多特征、偏移量或标准误计算。

纵轴省去共同的负log阶乘项，只比较同一数据内的迭代值。本例每次接受了完整Newton步。

观察时长与模型边界

若第i次记录持续 $e_i>0$ 个时间单位，应区分总计数均值与单位时间强度r：${\mu }_i=e_i{r}_i$，可写成 $\log {\mu }_i=\log e_i+x_i^{𝖳}\beta$，其中已知项 $\log e_i$ 的系数固定为一，称为偏移量。

变式。 主例计数仍为1、2、4，但每次都观察两个时间单位。按同一泊松率模型，单位时间拟合强度为0.5、1、2，参数改为 ${\beta }_0=-\log 2,{\beta }_1=\log 2$。加上已知偏移量log2后，总计数均值仍为1、2、4。若漏掉观察时长，数值上仍可拟合计数，却会把每时段数量误解释为单位时间强度。

对数连接作用于均值，不要求每个观测取对数；零计数在原质量函数中合法。泊松GLM还施加方差等于均值的条件，若数据来自共享随机强度、零膨胀或相关时段，常会违背这个条件。均值曲线写得合适并不使泊松标准误自动可靠，需检查分布与独立性，必要时采用其他响应族或另外论证的推断方法。GLM也不保证所有数据都有有限唯一的参数估计，满秩只能保证本例中严格凹的形状，不能独自保证最大值在有限位置取得。