はじめに

この資料では、同一のデータに対して

  1. 線形回帰
  2. 2次関数による回帰
  3. 線形項と周期関数を組み合わせた回帰

を行い、回帰曲線の違いを可視化するとともに、
AIC（赤池情報量規準）を用いてモデルの良さを比較する。

データの読み込みと概要

library(this.path)
setwd(this.dir())
library(minpack.lm)
library(ggplot2)

# x軸とy軸の散布図用テストファイルの読み込み
tes = read.table("http://mit.sci.shizuoka.ac.jp/test.dat", header = F)
colnames(tes) = c("x", "y")

# データの簡易プロット
plot(tes$x, tes$y)

線形回帰

モデル： \[ y = a + b x \]

fi = nlsLM(
  y ~ a + b*x,
  data  = tes,
  start = list(a = mean(tes$y), b = 0)
)

a1 = coef(fi)["a"]
b1 = coef(fi)["b"]

tes$yhat1 = a1 + b1 * tes$x

2次関数による回帰

モデル： \[ y = a + c x^2 \]

fi2 = nlsLM(
  y ~ a + c*x^2,
  data  = tes,
  start = list(a = mean(tes$y), c = 0)
)

a2 = coef(fi2)["a"]
c2 = coef(fi2)["c"]

tes$yhat2 = a2 + c2 * tes$x^2

線形項と周期関数(周期6)を含む回帰

モデル： \[ y = a + b x + A\cos\!\left(\frac{2\pi x}{6}\right) + B\sin\!\left(\frac{2\pi x}{6}\right) \]

fi3 = nlsLM(
  y ~ a + b*x + A*cos(2*pi*x/6) + B*sin(2*pi*x/6),
  data  = tes,
  start = list(a = mean(tes$y), b = 0, A = 0, B = 0)
)

a3 = coef(fi3)["a"]
b3 = coef(fi3)["b"]
A3 = coef(fi3)["A"]
B3 = coef(fi3)["B"]

tes$yhat3 = a3 +
  b3 * tes$x +
  A3 * cos(2*pi*tes$x/6) +
  B3 * sin(2*pi*tes$x/6)

回帰結果の比較（可視化）

ggplot(tes, aes(x = x, y = y)) +
  geom_point(size = 2) +
  geom_line(aes(y = yhat1), linewidth = 1) +
  geom_line(aes(y = yhat2), linewidth = 1, colour = "red") +
  geom_line(aes(y = yhat3), linewidth = 1,
            colour = "blue", linetype = "dashed") +
  theme(
    panel.background = element_rect(
      fill = "white",
      colour = "black",
      linetype = "solid"
    ),
    panel.grid = element_blank()
  )

ggsave("comp.png", width = 4, height = 3, dpi = 400)

AICによるモデル評価

AIC（赤池情報量規準）に基づくモデル比較を行う。
AICとは、「当てはまりの良さ」と「モデルの複雑さ」の
バランスを評価するための指標であり、次式で定義される。 \[ \mathrm{AIC} = -2\log(L) + 2k \] ここで \(L\) は尤度（残差平方和の関数）、\(k\) は推定したパラメータ数である。
この式から分かるように、AICは

という性質を持つ。そのため、AICは「データへの当てはまり」と
「モデルの単純さ」のバランスを評価する指標であり、
値が小さいほど良いモデルと解釈される。

尤度だけでモデルを比較すると、パラメータ数の多いモデルが
常に有利となり、過剰適合が生じやすいが、AICではその点が補正される。

AIC(fi)
AIC(fi2)
AIC(fi3)

残差からAICを計算する方法（理論的背景）

最小二乗回帰との関係

最小二乗回帰では、観測値 \(y\) は

\[ y_i = f(x_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma^2) \]

という正規分布に従う誤差を持つと仮定する。ここで残差 \(r_i\) を

\[ r_i = y_i - f(x_i) \]

と定義する。このとき、各データ点に対する確率密度は

\[ p(y_i\mid x_i,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{r_i^2}{2\sigma^2}\right) \]

で与えられる。独立な \(n\) 個のデータが与えられたとき、尤度 \(L\) は

\[ L = \prod_{i=1}^{n} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{r_i^2}{2\sigma^2}\right) \]

と書ける。

計算を簡単にするため、通常は尤度の対数（対数尤度）を用いる。対数尤度は

\[ \log L = -\frac{n}{2}\log(2\pi) -\frac{n}{2}\log(\sigma^2) -\frac{1}{2\sigma^2}\sum_{i=1}^{n} r_i^2 \]

となる。ここで残差平方和を

\[ \mathrm{RSS} = \sum_{i=1}^{n} r_i^2 \]

と定義する。\(\sigma^2\) を最尤推定すると

\[ \hat{\sigma}^2 = \frac{\mathrm{RSS}}{n} \]

となり、これを代入すると

\[ -2\log L = n\left(\log(2\pi)+1\right) + n\log\!\left(\frac{\mathrm{RSS}}{n}\right) \]

が得られる。

最小二乗におけるAICの実用式

したがって、最小二乗回帰（正規誤差仮定）におけるAICは

\[ \mathrm{AIC} = n\left(\log(2\pi)+1+\log\!\left(\frac{\mathrm{RSS}}{n}\right)\right) + 2k \]

と書ける。ここで \(n\) はデータ数、\(\mathrm{RSS}\) は残差平方和、\(k\) は推定したパラメータ数である。

なお、R の AIC() 関数では、回帰係数に加えて誤差分散 \(\sigma^2\) も1つのパラメータとして数える。
したがって \(k\) は「回帰係数の数 + 1」という立場を取っている。

n = length(resid(fi))
rss = sum(resid(fi)^2)
k = length(coef(fi)) + 1
n * (log(2*pi) + 1 + log(rss/n)) + 2*k