線形モデル（Linear Model）は統計的学習における最も基本的なモデルクラスであり、 一般化線形モデル（Generalized Linear Model, GLM）はその指数型分布族への体系的な拡張である。 両者は最尤推定・ベイズ推定・正則化・凸最適化と深く結びつき、 より複雑なモデルの理解における基礎を提供する。
入力（共変量）$\boldsymbol{x} \in \mathbb{R}^d$、 応答変数 $y \in \mathcal{Y}$、 パラメータ $\boldsymbol{\beta} \in \mathbb{R}^d$（切片を含む場合は $\boldsymbol{x}$ に 定数成分 $1$ を加えて $\boldsymbol{x} \in \mathbb{R}^{d+1}$）を考える。 i.i.d. 標本 $\mathcal{D}_n = \{(\boldsymbol{x}_i, y_i)\}_{i=1}^n$ が与えられているとする。 行列表記として $X \in \mathbb{R}^{n \times d}$（各行が $\boldsymbol{x}_i^\top$）、 $\boldsymbol{y} = (y_1, \ldots, y_n)^\top \in \mathbb{R}^n$ を用いる。
線形回帰モデルは
\[ y_i = \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i, \quad \varepsilon_i \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2), \quad i = 1, \ldots, n\]と定義される。行列形式では $\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}$、 $\boldsymbol{\varepsilon} \sim \mathcal{N}(\boldsymbol{0}, \sigma^2 I_n)$。 モデルの仮定を整理すると以下の通りである。
残差二乗和を最小化する最小二乗推定量（Ordinary Least Squares, OLS）は
\[ \hat{\boldsymbol{\beta}}^{\mathrm{OLS}} = \arg\min_{\boldsymbol{\beta} \in \mathbb{R}^d} \|\boldsymbol{y} - X\boldsymbol{\beta}\|_2^2\]と定義される。目的関数の勾配を零とおくと正規方程式（Normal Equations）
\[ X^\top X \hat{\boldsymbol{\beta}} = X^\top \boldsymbol{y}\]が得られ、$X^\top X$ が正則（$\mathrm{rank}(X) = d$）のとき一意解
\[ \hat{\boldsymbol{\beta}}^{\mathrm{OLS}} = (X^\top X)^{-1} X^\top \boldsymbol{y}\]が存在する。幾何学的には $\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}$ による予測値 $\hat{\boldsymbol{y}} = X\hat{\boldsymbol{\beta}} = H\boldsymbol{y}$ は $\boldsymbol{y}$ の列空間 $\mathrm{col}(X)$ への直交射影であり、 $H = X(X^\top X)^{-1}X^\top \in \mathbb{R}^{n \times n}$ を ハット行列（Hat Matrix）または射影行列と呼ぶ。 $H$ は対称・冪等（$H^2 = H$）であり、$\mathrm{rank}(H) = d$。
正規性の仮定を外し、$\mathbb{E}[\boldsymbol{\varepsilon}] = \boldsymbol{0}$、 $\mathrm{Cov}(\boldsymbol{\varepsilon}) = \sigma^2 I_n$ のみを仮定する。 Gauss–Markov 定理は以下を主張する： $\boldsymbol{c}^\top \boldsymbol{\beta}$（任意の線形結合）の 不偏推定量 $\boldsymbol{c}^\top \hat{\boldsymbol{\beta}}$ のうち、 OLS 推定量が最小分散を達成する（BLUE：Best Linear Unbiased Estimator）。
すなわち任意の線形不偏推定量 $\tilde{\boldsymbol{\beta}} = A\boldsymbol{y}$（$AX = I_d$）に対して
\[ \mathrm{Var}(\boldsymbol{c}^\top \tilde{\boldsymbol{\beta}}) \geq \mathrm{Var}(\boldsymbol{c}^\top \hat{\boldsymbol{\beta}}^{\mathrm{OLS}}) = \sigma^2 \boldsymbol{c}^\top (X^\top X)^{-1} \boldsymbol{c}\]が任意の $\boldsymbol{c} \in \mathbb{R}^d$ について成立する。 さらに正規性 $\boldsymbol{\varepsilon} \sim \mathcal{N}(\boldsymbol{0}, \sigma^2 I_n)$ を仮定すると、 OLS 推定量は最尤推定量（MLE）に一致し、 すべての不偏推定量のうち最小分散を達成する（UMVUE）。
正規性の仮定のもとで
\[ \hat{\boldsymbol{\beta}}^{\mathrm{OLS}} \sim \mathcal{N}\!\left(\boldsymbol{\beta},\, \sigma^2 (X^\top X)^{-1}\right)\]が成立する。$\sigma^2$ の不偏推定量として残差分散
\[ \hat{\sigma}^2 = \frac{\|\boldsymbol{y} - X\hat{\boldsymbol{\beta}}\|^2}{n - d} = \frac{\mathrm{RSS}}{n-d}\]を用いる（$\mathrm{RSS}$：残差二乗和、$n-d$ は自由度）。 $\hat{\sigma}^2$ は $\hat{\boldsymbol{\beta}}^{\mathrm{OLS}}$ と独立に $\sigma^2 \chi^2(n-d)/(n-d)$ に従う。 これらを用いて $\beta_j$ の $t$ 検定統計量
\[ t_j = \frac{\hat{\beta}_j}{\hat{\sigma}\sqrt{[(X^\top X)^{-1}]_{jj}}} \sim t(n-d)\]および $F$ 検定（モデル全体の有意性）が構成される。 また $\boldsymbol{\beta}$ の $95\%$ 信頼区間は $\hat{\beta}_j \pm t_{n-d,\,0.025} \cdot \hat{\sigma}\sqrt{[(X^\top X)^{-1}]_{jj}}$ として得られる。
モデルの当てはまりの指標として決定係数（Coefficient of Determination）
\[ R^2 = 1 - \frac{\mathrm{RSS}}{\mathrm{TSS}} = 1 - \frac{\|\boldsymbol{y} - \hat{\boldsymbol{y}}\|^2}{\|\boldsymbol{y} - \bar{y}\boldsymbol{1}\|^2} \in [0, 1]\]を用いる（$\mathrm{TSS}$：全二乗和、$\bar{y} = \frac{1}{n}\sum_i y_i$）。 $R^2$ は説明変数を追加するたびに単調増加するため、 モデルの複雑度（説明変数の数 $d$）を補正した自由度調整済み $R^2$
\[ \bar{R}^2 = 1 - \frac{\mathrm{RSS}/(n-d)}{\mathrm{TSS}/(n-1)} = 1 - (1-R^2)\frac{n-1}{n-d}\]がモデル比較に用いられる。 残差 $\hat{\varepsilon}_i = y_i - \boldsymbol{x}_i^\top\hat{\boldsymbol{\beta}}$ の 分散は均一でなく $\mathrm{Var}(\hat{\varepsilon}_i) = \sigma^2(1 - H_{ii})$ であるため、 スチューデント化残差 $r_i = \hat{\varepsilon}_i / (\hat{\sigma}\sqrt{1-H_{ii}})$ が診断に用いられる。
誤差の分散が均一でない場合（不均一分散）、 $\mathrm{Cov}(\boldsymbol{\varepsilon}) = \sigma^2 W^{-1}$（$W$ は既知の正定値対角行列）のとき、 加重最小二乗法（WLS）
\[ \hat{\boldsymbol{\beta}}^{\mathrm{WLS}} = \arg\min_{\boldsymbol{\beta}} (\boldsymbol{y} - X\boldsymbol{\beta})^\top W (\boldsymbol{y} - X\boldsymbol{\beta}) = (X^\top W X)^{-1} X^\top W \boldsymbol{y}\]が Gauss–Markov 定理の意味で BLUE となる。 さらに $\mathrm{Cov}(\boldsymbol{\varepsilon}) = \sigma^2 \Sigma$（$\Sigma$ は既知の正定値行列）の 一般的な相関誤差の設定では、一般化最小二乗法（GLS）
\[ \hat{\boldsymbol{\beta}}^{\mathrm{GLS}} = (X^\top \Sigma^{-1} X)^{-1} X^\top \Sigma^{-1} \boldsymbol{y}\]が BLUE である。GLS は変換 $\tilde{X} = \Sigma^{-1/2}X$、 $\tilde{\boldsymbol{y}} = \Sigma^{-1/2}\boldsymbol{y}$ による OLS として実行できる。
GLM の理論的基礎として、指数型分布族（Exponential Family）を導入する。 確率密度（または確率質量）関数が
\[ p(y ; \eta, \phi) = \exp\!\left(\frac{y\eta - b(\eta)}{a(\phi)} + c(y, \phi)\right)\]の形に書けるとき、この分布族を指数型分布族と呼ぶ。 $\eta$ を自然パラメータ（Natural Parameter）、 $\phi > 0$ を分散パラメータ（Dispersion Parameter）と呼ぶ。 $b(\eta)$ を対数分配関数（Log-Partition Function）、 $a(\phi)$ を分散関数、$c(y,\phi)$ を規格化定数と呼ぶ。
対数分配関数の微分は平均と分散を与えるという重要な性質がある：
\[ \mathbb{E}[Y;\eta] = b'(\eta), \qquad \mathrm{Var}(Y;\eta) = b''(\eta)\cdot a(\phi)\]この関係は対数尤度の微分（スコア関数）の期待値がゼロであることから導かれ、 GLM の理論全体の基礎となる。代表的な分布の対応を以下に示す。
| 分布 | $\mathcal{Y}$ | 自然パラメータ $\eta$ | $b(\eta)$ | $a(\phi)$ | 平均 $\mu = b'(\eta)$ |
|---|---|---|---|---|---|
| 正規分布 $\mathcal{N}(\mu,\sigma^2)$ | $\mathbb{R}$ | $\mu$ | $\eta^2/2$ | $\sigma^2$ | $\eta$ |
| ベルヌーイ分布 $\mathrm{Ber}(p)$ | $\{0,1\}$ | $\log\frac{p}{1-p}$ | $\log(1+e^\eta)$ | $1$ | $\sigma(\eta) = \frac{e^\eta}{1+e^\eta}$ |
| ポアソン分布 $\mathrm{Poi}(\lambda)$ | $\mathbb{Z}_{\geq 0}$ | $\log\lambda$ | $e^\eta$ | $1$ | $e^\eta$ |
| ガンマ分布 $\mathrm{Gamma}(\alpha,\beta)$ | $\mathbb{R}_{>0}$ | $-\beta$ | $-\log(-\eta)$ | $1/\alpha$ | $\alpha/\beta$ |
| 負の二項分布 | $\mathbb{Z}_{\geq 0}$ | $\log\frac{p}{1-p}$ | $-r\log(1-e^\eta)$ | $1$ | $r(1-p)/p$ |
一般化線形モデル（Nelder–Wedderburn, 1972）は以下の三要素で定義される。
自然パラメータ $\eta$ を線形予測子と直接等置する $g(\mu) = \eta$（自然パラメータが線形予測子に等しくなるリンク）を 正準リンク関数（Canonical Link Function）と呼ぶ。 正準リンクのもとでは MLE の正規方程式が特に単純な形になり、 数値的安定性も高い。
| モデル | 応答変数 | 分布 | 正準リンク $g(\mu)$ | 平均関数 $g^{-1}(\eta)$ |
|---|---|---|---|---|
| 線形回帰 | 連続値 | 正規分布 | $\mu$（恒等リンク） | $\eta$ |
| ロジスティック回帰 | 二値 $\{0,1\}$ | ベルヌーイ分布 | $\log\frac{\mu}{1-\mu}$（ロジット） | $\sigma(\eta) = \frac{1}{1+e^{-\eta}}$ |
| ポアソン回帰 | カウント $\mathbb{Z}_{\geq 0}$ | ポアソン分布 | $\log\mu$（対数リンク） | $e^\eta$ |
| ガンマ回帰 | 正の連続値 | ガンマ分布 | $1/\mu$（逆数リンク） | $1/\eta$ |
| プロビット回帰 | 二値 $\{0,1\}$ | ベルヌーイ分布 | $\Phi^{-1}(\mu)$（プロビット） | $\Phi(\eta)$ |
GLM の対数尤度は
\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \log p(y_i; \eta_i, \phi) = \sum_{i=1}^n \frac{y_i \eta_i - b(\eta_i)}{a(\phi)} + \text{const.}\]と書ける。$\boldsymbol{\beta}$ に関するスコア方程式は
\[ \frac{\partial \ell}{\partial \boldsymbol{\beta}} = \sum_{i=1}^n \frac{(y_i - \mu_i)}{a(\phi) b''(\eta_i)} \cdot \frac{d\eta_i}{d\mu_i} \cdot \boldsymbol{x}_i = \boldsymbol{0}\]と得られる。これは一般に閉形式解を持たないが、 反復重み付き最小二乗法（Iteratively Reweighted Least Squares, IRLS） によって効率的に解くことができる。
IRLS の各ステップでは現在の推定値 $\hat{\boldsymbol{\beta}}^{(t)}$ から 調整応答変数（Adjusted Dependent Variable）
\[ z_i^{(t)} = \hat{\eta}_i^{(t)} + (y_i - \hat{\mu}_i^{(t)}) \frac{d\eta_i}{d\mu_i}\bigg|_{\hat{\mu}_i^{(t)}}\]と作業重み（Working Weight）
\[ w_i^{(t)} = \frac{1}{b''(\hat{\eta}_i^{(t)}) \cdot \left(\dfrac{d\eta_i}{d\mu_i}\bigg|_{\hat{\mu}_i^{(t)}}\right)^2 \cdot a(\phi)}\]を計算し、加重最小二乗問題
\[ \hat{\boldsymbol{\beta}}^{(t+1)} = \arg\min_{\boldsymbol{\beta}} \sum_{i=1}^n w_i^{(t)} (z_i^{(t)} - \boldsymbol{x}_i^\top \boldsymbol{\beta})^2 = (X^\top W^{(t)} X)^{-1} X^\top W^{(t)} \boldsymbol{z}^{(t)}\]を解く（$W^{(t)} = \mathrm{diag}(w_1^{(t)},\ldots,w_n^{(t)})$）。 IRLS は Newton–Raphson 法と等価であり、 指数型分布族の対数尤度の凸性から大域収束が保証される （正準リンク使用時）。 正準リンクのもとでは $w_i = b''(\eta_i)/a(\phi) = \mathrm{Var}(Y_i)/a(\phi)^2$ が簡略化される。
正則条件のもとで GLM の MLE $\hat{\boldsymbol{\beta}}$ は 最尤推定の一般論（前節）と同様の漸近性質を持つ：
\[ \sqrt{n}(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}^*) \xrightarrow{d} \mathcal{N}(0,\, \mathcal{I}(\boldsymbol{\beta}^*)^{-1})\]フィッシャー情報行列は
\[ \mathcal{I}(\boldsymbol{\beta}) = \frac{1}{n} X^\top W X\]と表される（$W = \mathrm{diag}(w_1,\ldots,w_n)$ は IRLS の作業重みに対応）。 $\hat{\boldsymbol{\beta}}$ の漸近共分散行列の推定量は $(X^\top \hat{W} X)^{-1}$（$\hat{W}$ は収束時の作業重み行列）であり、 各係数の $t$ 検定・信頼区間の構成に用いられる。
ロジスティック回帰は GLM の中で最も広く用いられる二値分類モデルである。 $Y_i \in \{0,1\}$、$p_i = P(Y_i = 1 \mid \boldsymbol{x}_i)$ として、
\[ \log\frac{p_i}{1-p_i} = \boldsymbol{x}_i^\top \boldsymbol{\beta}, \qquad p_i = \sigma(\boldsymbol{x}_i^\top \boldsymbol{\beta}) = \frac{1}{1 + e^{-\boldsymbol{x}_i^\top \boldsymbol{\beta}}}\]ここで $\sigma(\cdot)$ はシグモイド関数（ロジスティック関数）である。
対数尤度は
\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \bigl[ y_i \log p_i + (1-y_i)\log(1-p_i) \bigr] = \sum_{i=1}^n \bigl[ y_i \boldsymbol{x}_i^\top \boldsymbol{\beta} - \log(1 + e^{\boldsymbol{x}_i^\top \boldsymbol{\beta}}) \bigr]\]であり、負の対数尤度（損失関数）はバイナリ交差エントロピー損失に一致する：
\[ -\frac{1}{n}\ell(\boldsymbol{\beta}) = \frac{1}{n}\sum_{i=1}^n \bigl[ -y_i \log p_i - (1-y_i)\log(1-p_i) \bigr]\]この損失関数は $\boldsymbol{\beta}$ について凸であるため（ヘッセ行列 $\nabla^2(-\ell) = \frac{1}{n}X^\top \mathrm{diag}(p_i(1-p_i)) X \succeq 0$）、 大域最適解が一意に存在する（$X$ が列フルランクのとき正定値となり狭義凸）。 ただし線形分離可能なデータでは $\hat{\boldsymbol{\beta}} \to \infty$ となり MLE が存在しない（分離問題）。この場合 $\ell_2$ 正則化が本質的に必要となる。
$K$ クラス分類（$Y_i \in \{1,\ldots,K\}$）への拡張として、 多項ロジスティック回帰（ソフトマックス回帰）は
\[ P(Y_i = k \mid \boldsymbol{x}_i) = \frac{\exp(\boldsymbol{x}_i^\top \boldsymbol{\beta}_k)} {\sum_{j=1}^K \exp(\boldsymbol{x}_i^\top \boldsymbol{\beta}_j)}, \quad k = 1, \ldots, K\]と定義される（$\boldsymbol{\beta}_K = \boldsymbol{0}$ と固定して識別可能性を保証）。 対応する損失関数はカテゴリカル交差エントロピーであり、 ニューラルネットワークの最終層と組み合わせて広く用いられる。 ソフトマックス関数 $\mathrm{softmax}(\boldsymbol{\eta})_k = e^{\eta_k}/\sum_j e^{\eta_j}$ は $\log\sum_j e^{\eta_j}$（Log-Sum-Exp）の勾配として得られ、 これは凸関数であるため多項ロジスティック回帰の負の対数尤度も凸である。
ポアソン回帰はカウントデータ（非負整数値の応答変数）のモデリングに用いられる。 $Y_i \mid \boldsymbol{x}_i \sim \mathrm{Poi}(\lambda_i)$、 $\log \lambda_i = \boldsymbol{x}_i^\top \boldsymbol{\beta}$（$\lambda_i > 0$）として、
\[ \lambda_i = \exp(\boldsymbol{x}_i^\top \boldsymbol{\beta})\]対数リンクにより $\lambda_i$ の非負性が自動的に保証される点が利点である。 対数尤度は
\[ \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \bigl[ y_i \boldsymbol{x}_i^\top \boldsymbol{\beta} - e^{\boldsymbol{x}_i^\top \boldsymbol{\beta}} - \log(y_i!) \bigr]\]であり、これも $\boldsymbol{\beta}$ について凸である。 ポアソン回帰の仮定（平均と分散が等しい $\mathbb{E}[Y_i] = \mathrm{Var}(Y_i) = \lambda_i$）が 成立しない過分散（Overdispersion）の場合は、 負の二項回帰（Negative Binomial Regression）が代替として用いられる。 また観測値に真のゼロが過剰に存在するゼロ過剰（Zero Inflation）の場合は、 ゼロ過剰ポアソンモデル（ZIP）が適用される。
GLM のモデル評価には、線形回帰の RSS に対応する逸脱度（Deviance）を用いる。 飽和モデル（各観測値に固有のパラメータを持つ完全適合モデル）の対数尤度を $\ell_{\mathrm{sat}}$、当てはめモデルの対数尤度を $\ell(\hat{\boldsymbol{\beta}})$ とするとき、
\[ D = 2\bigl[\ell_{\mathrm{sat}} - \ell(\hat{\boldsymbol{\beta}})\bigr] \geq 0\]を残差逸脱度（Residual Deviance）と定義する。 $D$ が小さいほどモデルの当てはまりが良い。 正則条件のもとで $D \xrightarrow{d} \chi^2(n-d)$（$n \to \infty$）が成立し、 モデルの適合度検定に利用できる。
二つのネストしたモデル（$\mathcal{M}_0 \subset \mathcal{M}_1$）の比較には 尤度比検定（前節）と等価な逸脱度の差
\[ \Delta D = D_0 - D_1 = 2[\ell(\hat{\boldsymbol{\beta}}_1) - \ell(\hat{\boldsymbol{\beta}}_0)] \xrightarrow{d} \chi^2(d_1 - d_0)\]を用いる。またモデル選択には AIC（前節）の GLM 版
\[ \mathrm{AIC} = -2\ell(\hat{\boldsymbol{\beta}}) + 2d\]が広く用いられる。ピアソン統計量 $X^2 = \sum_i (y_i - \hat{\mu}_i)^2 / \hat{V}_i$ （$\hat{V}_i = b''(\hat{\eta}_i) a(\phi)$ は分散の推定値）も 逸脱度と漸近的に等価な適合度統計量である。
高次元設定（$d \gg n$）や多重共線性の存在下では、 GLM の MLE が不安定または存在しないため正則化が必要となる。 正則化付き GLM の一般形は
\[ \hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \left[ -\frac{1}{n}\ell(\boldsymbol{\beta}) + \lambda\Omega(\boldsymbol{\beta}) \right]\]と書ける。$\Omega = \|\boldsymbol{\beta}\|_2^2$ の Ridge ロジスティック回帰は ロジスティック回帰の分離問題を解消しつつ係数を縮小し、 $\Omega = \|\boldsymbol{\beta}\|_1$ の Lasso ロジスティック回帰は スパースな特徴選択を行う。 目的関数は凸のまま保たれ（$-\ell$ が凸、$\Omega$ が凸の和）、 IRLS の各ステップに Ridge/Lasso の近接写像を組み合わせた 近接 Newton 法（glmnet アルゴリズム等）が効率的な解法として用いられる。
線形回帰は正規分布・恒等リンクによる GLM の特殊ケースであり、 OLS 推定量は Gauss–Markov 定理により BLUE、 正規性仮定のもとで MLE かつ UMVUE となる。 GLM は指数型分布族・リンク関数・線形予測子の三要素により、 正規・ベルヌーイ・ポアソン・ガンマ等の分布を統一的な枠組みで記述する。 MLE は IRLS（Newton–Raphson 法と等価）により計算され、 漸近正規性・効率性がフィッシャー情報行列を通じて保証される。 逸脱度は RSS の一般化として適合度評価・モデル比較・尤度比検定を提供し、 正則化付き GLM は高次元設定での推定安定性とスパース変数選択を可能にする。 GLM はニューラルネットワークの最終層・一般化加法モデル・ 混合モデルへと拡張され、現代の統計的機械学習の基礎的構成要素となっている。
Mathematics is the language with which God has written the universe.