R言語で対数尤度関数を作り最適化を行う

Rで最尤法を行おうと思ったのですが、ネットで調べてもいまいち分からなかったため自分で作ってみました。自然な設定として線形回帰モデルを考え、最尤推定量を求めます。

モデル：

$\displaystyle y = X\mu + \epsilon$

ここで $y \in \mathbb{R}^n, X \in \mathbb{R}^n \times \mathbb{R}^q, \mu\in \mathbb{R}^q, \epsilon \in \mathbb{R}^n$ . 誤差項の分布として $\epsilon \sim$ $N(0, \sigma^{2})$ を仮定。

設定:

N <- 1000 # サンプル数
mu <- c(5,2,3) # 回帰係数
sig <- 3 # 分散
X <- matrix(runif(3*N,min=1,max=5), nr=N, nc=3) # 説明変数行列
par.init <- c(2,2,2,1) # 最適化の初期値

対数尤度関数:

$\displaystyle \ell_n(\theta):= \sum_{j=1}^{n} \log \frac{1}{\sigma}f \left(\frac{y_j - X_j \mu}{\sigma}\right)$

ここで $f(x)$ は $N(0,1)$ の確率密度関数、 $y_j: y$ の $j$ 番目の成分、 $X_j: X$ の $j$ 行。

log.likelihood <- function(y){
  return(function(par){
    mu1 <- par[1]
    mu2 <- par[2]
    mu3 <- par[3]
    mu <- c(mu1, mu2, mu3)
    sig <- par[2]
    return(sum(log(dnorm((y-X%*%mu)/sig, 0,1)/sig)))
  })
}

推定:

~
t <- proc.time()
y0 <- X%*%mu + rnorm(N,0,sig)  
mle.sample <- optim(par.init, log.likelihood(y0), method = "L-BFGS-B", lower = c(-Inf,-Inf,-Inf,0), upper = c(Inf, Inf,Inf,Inf), control = list(fnscale = -1))$par
proc.time()-t
mle.sample
~

計算時間と推定結果

> proc.time()-t
   user  system elapsed 
  0.006   0.001   0.006 
> mle.sample
[1] 4.670545 2.559163 2.805950 1.000000

とげまる日記

九州の大学で数学の博士課程に通ってます。自分の日記として使っています。

R言語で対数尤度関数を作り最適化を行う