廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

左切断右打ち切りデータからの最尤推定

切断と打ち切り

観測対称がとる値がある範囲を超えたとき、まったく観測されない場合を切断(truncated)と呼ぶ。

それに対して、観測されない対象の個数はわかっている場合を打ち切り(censored)と呼ぶ。

観測範囲がある値以上に限定されている場合、左切断という。

観測値がある値を超えたことはわかっているが、値そのものはわからない場合を右打ち切りと呼ぶ。

尤度関数

d_i を打ち切りの有無を表すインジケータ(x_i が右打ち切りされた観測値の場合に 1、そうでなければ 0 の値をとる)とすると、右打ち切りデータの尤度関数は、

\displaystyle \prod_{i=1}^{n} \{g(x_i)\}^{d_i} \{1-G(x_i)\}^{1-d_i}

である。

ここで、g は左切断分布の密度関数、

 \displaystyle g(x)=f(x|a < X) = \frac{f(x)}{1-F(a)}

G は左切断分布の分布関数、

 \displaystyle G(x)=\int^{x}_{a} g(u) \, du = \frac{F(x)-F(a)}{1-F(a)}

である。

R によるシミュレーション

0 で左切断、3 で右打ち切りされた正規分布のデータから、最尤推定を行うには以下のようにする。

dtnorm <- function(x,mu,sigma){
  dnorm(x,mu,sigma)/pnorm(0,mu,sigma,lower.tail = FALSE)
}
ptnorm <- function(x,mu,sigma){
  (pnorm(x,mu,sigma)-pnorm(0,mu,sigma))/pnorm(0,mu,sigma,lower.tail = FALSE)
}
ll <- function(par){
  sum(d*log(dtnorm(x,par[1],par[2])))+sum((1-d)*log(1-ptnorm(x,par[1],par[2])))
}
x <-rnorm(100,2,1)
x <- x[x>0] #左切断
d <- x<3
x[!d] <-3 #右打ち切り
optim(c(1,1),ll,control = list(fnscale=-1))