対数尤度が正で、AICが負になるケース - StatsBeginner: 初学者の統計学習ノート

初歩的な話ですが、今日学生から訊かれたので例を考えてみました。
AICがマイナスの値になってもいいんだっけ？みたいな話です。

AICがマイナスに

統計モデルの最尤推定をする場合、尤度は1より小さい場合が多く、したがって対数尤度が負である場合が多いとしましょう*1。特に離散変数では必ずそうなります。その場合、AICは「-2×(対数尤度-パラメータ数)」なので、正になるでしょう。
でも、対数尤度が正でAICが負になることというのも、連続変数ではあり得ます。その場合AICは、負で絶対値が大きいほどよい、ということになる。

誤差の分布は正規分布だと考えているとして、たとえば、

モデルの当てはまりがめちゃめちゃいい
そもそも数値の単位が小さい

というような場合に、誤差の標準偏差（分散）の値が小さくなって、結果的にAICが負になりやすくなりますね。
　
　

正規分布の確率密度

統計モデルの誤差が平均ゼロ、標準偏差σの正規分布に従うとした場合、確率密度は誤差ゼロのときに最大になるわけですが、誤差ゼロのときの確率密度は、σが小さいと大きくなります。
Rでは、dnorm()という関数で正規分布上の確率密度が簡単に計算できますが、平均ゼロの正規分布上でゼロという値の確率密度を、標準偏差が1、0.1、0.01の3パターンで計算してみると、

> # sdは標準偏差
> dnorm(0, mean=0, sd=1)
[1] 0.3989423
> dnorm(0, mean=0, sd=0.1)
[1] 3.989423
> dnorm(0, mean=0, sd=0.01)
[1] 39.89423

というふうに、標準偏差に反比例して大きくなっていきます。正規分布の確率密度関数でいうと、 $x$ も $\mu$ （平均）も0で、右のところが $\exp(0)=1$ になるので、まぁ当たり前です。

$f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})$

平均値に一致する点（誤差ゼロの点）でなくても、たとえばプラス方向に1標準偏差分行ったところの値も、同じようになります。この場合はexpのところがexp(-0.5)に固定になりますからね。

> dnorm(1, mean=0, sd=1)
[1] 0.2419707
> dnorm(0.1, mean=0, sd=0.1)
[1] 2.419707
> dnorm(0.01, mean=0, sd=0.01)
[1] 24.19707

標準偏差（Sigma）が0.5の場合と0.1の場合について確率密度関数のグラフを書いてみると次のようになります。

x <- seq(from=-1, to=1, by=0.01)
Density <- dnorm(x, mean=0, sd=0.5)
plot(x, Density, type='l', col='blue', ylim=c(0,4))

Density <- dnorm(x, mean=0, sd=0.1)
par(new=T)
plot(x, Density, type='l', col='red', ylim=c(0,4))

text(x=0.1, y=3, labels='Sigma=0.1', col='red', pos=4)
text(x=0.2, y=1, labels='Sigma=0.5', col='blue', pos=4)

f:id:midnightseminar:20200123002137p:plain

標準偏差0.5の場合は、xがどんな値であっても確率密度が1より小さいので、それを掛け算してモデルの尤度が計算されるわけですから、1より小さい値でどんどん小さくなっていきます。当てはまりが良くないモデルだと、最尤法で計算する場合の母標準偏差の推定値も比較的大きくなって、こういう横に広がった分布が想定されることになります。
一方、標準偏差0.1の場合は、平均値まわりの割と広い範囲で確率密度が1を超えてますね。
　
　

当てはまりが良いとAICがマイナスになりやすくなる

切片ゼロ、傾き1（α=0、β=1）の線形モデルに、平均ゼロ、標準偏差sigmaの正規分布に従う乱数を加えてデータをつくって、単回帰分析を行い、そのモデルの決定係数、対数尤度、AICを計算してみます。まずsigma=1としてみます。

> # データを生成
> sigma <- 1
> x <- seq(from=0.1, to=10, by=0.1)
> y <- x + rnorm(n=100, mean=0, sd=sigma)
> 
> # 回帰モデルを推定して、対数尤度とAICを計算する
> lm1 <- lm(y~x)
> r2 <- summary(lm1)$r.squared
> lm1.loglik <- logLik(lm1)
> lm1.aic <- AIC(lm1)
> print(r2)
[1] 0.896859
> print(lm1.loglik)  # logLikの結果はlog Likという特殊なクラスになってる
'log Lik.' -136.302 (df=3)
> print(lm1.aic)
[1] 278.604
> # 対数尤度とAICの関係
> -2*(as.numeric(lm1.loglik)-3)
[1] 278.604

対数尤度（log Lik）は負で、AICが正ですね。df=3となっているのは、単回帰モデルなのでパラメータが「切片」「傾き」「誤差の標準偏差」の3つなのを表してるんでしょう。
以下のようにグラフにしてみます。

title <- paste('Sigma=',as.character(sigma),
               ', R2=',as.character(round(r2,2)),
               ', LogLik=', as.character(round(lm1.loglik,2)),
               ', AIC=', as.character(round(lm1.aic,2)),
               sep='')

plot(x, y, pch=20, main=title)  # 散布図
abline(lm1, col='red', lwd=2)  # 推定された回帰直線を引く

f:id:midnightseminar:20200123000600p:plain

同じやり方で、sigma=0.5、0.2、-0.1の場合も出してみると、次のようになります。

f:id:midnightseminar:20200123000645p:plain

誤差分散（標準偏差）が小さくて当てはまりがめちゃめちゃいい場合は、AICが負になることが分かりますね。
ちなみに、標準偏差を0.01から2まで少しずつ変えていって、その都度乱数でデータを生成して回帰モデルを当てはめ、AICを計算するという作業を繰り返してその結果をグラフにすると、次のように対数の形なります。

get_aic <- function(sigma){
  x <- seq(from=0.1, to=10, by=0.1)
  y <- x + rnorm(n=100, mean=0, sd=sigma)
  model.1 <- lm(y ~ x)
  return(AIC(model.1))
}

aics <- c()
sigmas <- seq(from=0.01, to=2, by=0.01)
for (i in sigmas){
  aic <- get_aic(i)
  aics <- c(aics, aic)
}

plot(x=sigmas, y=aics, type='l',xlab='Sigma', ylab='AIC')

f:id:midnightseminar:20200123001231p:plain
　
　

単純に、小さい数値を扱ってるとAICが負になりやすくなる

ところで、もう一つ、単にデータの値が小さい数字になるような単位だと、AICが負になりやすくなります。
これも同様に、誤差の分散が値としては小さくなるからですが、以下のように試してみます。
まず、誤差の標準偏差2のデータに単回帰モデルを当てはめて決定係数、対数尤度、AICを出してみて、元データを10で割ってから同じことをしてみます。

# 標準偏差2でやってみる
sigma <- 2
x1 <- seq(from=0.1, to=10, by=0.1)
y1 <- x1 + rnorm(n=100, mean=0, sd=sigma)
lm1 <- lm(y1~x1)
lm1.r2 <- summary(lm1)$r.squared
lm1.loglik <- logLik(lm1)
lm1.aic <- AIC(lm1)
title <- paste('Sigma=',as.character(sigma),
               ', R2=',as.character(round(lm1.r2,2)),
               ', LogLik=', as.character(round(lm1.loglik,2)),
               ', AIC=', as.character(round(lm1.aic,2)),
               sep='')
plot(x1, y1, pch=20, main=title)  # 散布図
abline(lm1, col='red', lwd=2)  # 推定された回帰直線を引く

f:id:midnightseminar:20200123011822p:plain

次に、単にデータを10で割り算してから同じことをやってみます。
Sigmaは、自分で指定してないですが、lmの結果の中に入っている推定値を取り出して見ておきます。（それもだいたい10分の1になります。）

# データを10分の1に
x2 <- x1/10
y2 <- y1/10
lm2 <- lm(y2~x2)
lm2.r2 <- summary(lm2)$r.squared
lm2.sigma.est <- summary(lm2)$sigma  # σは推定値を取り出す
lm2.loglik <- logLik(lm2)
lm2.aic <- AIC(lm2)
title <- paste('Sigma=',as.character(round(lm2.sigma.est,3)),
               ', R2=',as.character(round(lm2.r2,2)),
               ', LogLik=', as.character(round(lm2.loglik,2)),
               ', AIC=', as.character(round(lm2.aic,2)),
               sep='')
plot(x2, y2, pch=20, main=title)
abline(lm2, col='red', lwd=2)