正規分布 vs t分布 - 外れ値に影響されやすい度グラフ化してみた

PRML 図2.16を再現してみた。

まずはグラフから
f:id:mathhun:20131129080215p:plain:w300
f:id:mathhun:20131129080218p:plain:w300

上は外れ値なし、下はあり。赤が正規分布で青はt分布。

外れ値なしでは2つのグラフはほぼ重なる。このグラフは乱数を何度か取り直してあえて重なりが少なくなるものを選んだが、たいていはほぼ完全に一致する。

Rコード

data <- rnorm(50, -1, 1.2)
noise <- rnorm(5, 10, 1) # 外れ値なしの場合は第一引数を0にする
data <- c(data, noise)

# norm
log.likelihood.norm <- function(x) {
	return(function(par) {
		mu <- par[1]
		sigma2 <- par[2]
		- length(x) / 2 * log(sigma2) - 1 / 2 * sum((x - mu)^2) / sigma2
	})
}

my.dt <- function(x, mu, lambda, v) {
	gamma(v/2 + 1/2) / gamma(v/2) * (lambda/pi/v)^(1/2) * (1 + lambda*(x - mu)^2 / v)^(-v/2 -1/2)
}
my.dt.log <- function(x, mu, lambda, v) {
	log(gamma(v/2 + 1/2)) - log(gamma(v/2)) + (1/2)*log(lambda/pi/v) + (-(v/2) - (1/2)) * log(1 + lambda*(x - mu)^2 / v)
}

# t
likelihood.t <- function(x) {
	return(function(par) {
		mu <- par[1]
		lam <- par[2]
		v <- par[3]
		
		sum(my.dt.log(x, mu, lam, v))
	})
}


# 解析的に求められる式で計算
mu <- mean(data)
sigma2 <- sum((data - mu)^2) / length(data)

# 数値的に求めた結果
# norm
opt <- optim(par=c(0,1), fn=log.likelihood.norm(data), control=list(fnscale=-1))
mu.guess <- opt$par[1]
sigma2.guess <- opt$par[2]

print(list(
  title="norm",
  mu=mu,
  sigma=sigma2,
  mu.guess=mu.guess,
  sigma2.guess=sigma2.guess
))

# t
opt <- optim(par=c(0,1,1), fn=likelihood.t(data), control=list(fnscale=-1))
mu.t.guess <- opt$par[1]
lambda.guess <- opt$par[2]
v.guess <- opt$par[3]

print(list(
  title="t",
  mu=mu.t.guess,
  lambda=lambda.guess,
  v=v.guess
))


# ---
# visualization
# ---

x.max <- 15
x.min <- -5
y.max <- 0.5

# raw data
hist(data, xlim=c(x.min, x.max), ylim=c(0, y.max), breaks=seq(x.min, x.max, 1), freq=F)
par(new=T)

# norm
f <- function(x) dnorm(x, mean=mu.guess, sd=sqrt(sigma2.guess))
curve(f, xlab='', ylab='', xlim=c(x.min, x.max), ylim=c(0, y.max), col="red")
par(new=T)

# t
f2 <- function(x) my.dt(x, mu.t.guess, lambda.guess, v.guess)
curve(f2, xlab='', ylab='', xlim=c(x.min, x.max), ylim=c(0, y.max), col="blue")

数式

t分布
本に載ってる定義
$St(x | \mu, \lambda, \nu)=\frac{\Gamma(\nu/2 + 1/2)}{\Gamma(\nu/2)} \biggl(\frac{\lambda}{\pi \nu} \biggr)^{1/2} \biggl[ 1 + \frac{\lambda(x - \mu)^2}{\nu} \biggr] ^{-\nu/2 - 1/2}$

はまったところ

この程度のグラフは簡単に描けると思ったがかなり苦戦した。t分布のパラメータ推定でいろいろ苦労したので書いておく。

t分布というと、「正規分布のパラメータ推定の時に脇役的に登場するあのヒトね」という感じだが、t分布自体のデータへの当てはめというのは知らなかった。ググっても情報少なめ。
- そもそも普通t分布のパラメータは1つ。正規分布のように平均などを指定できないの？と思った。
- よく見るとPRMLに載ってる式にはパラメータが3つあり、正規分布の2つのパラメータに対応するものがある。
推定値の計算
- 解析的に求める式が本に載っていないし自力で計算できるかも謎。Rのoptim関数というものがあることを知ったのでそれを使ってみることに。
- 最終推定でlogとらないバージョンで結果を求めると推定値が求まらない。(初期値として与えたのとほぼ同じ値が返ってくる)
- 調べてみると絶対値の小さな値を多量に掛け算しているので結果が0になっている(アンダーフロー)
- t分布の関数logバージョンで再度optimに掛けてそれらしき値が得られる