読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

R 最適化計算 - optimize(), optim()

optimize(), optim() の使い方に関するメモ。

 

(1)optimize()による1パラメータ推定

ベクトルXの平均値を(ややこしく)求めてみる。

X <- c(1, 2, 3, 5, 8)
mean(X)     # 3.8

optimize()を適用し、パラメータの最適値を求める対象の関数をまず定義する。

find_mean <- function(x){
    return(function(m){
      return( sum( (x-m)^2 ) )
    })
}

xを外側関数、mを内側関数が受け取る2重関数。

find_mean(x)(m)と実行すると計算結果が返ってくる。

ベクトルxを渡し、その平均mを最適化計算の関数(optimize、optim)で返す処理を行う場合に定義する。

sum( (x-m)^2 )は、xをベクトルとしたとき、x-mの二乗和を最小にするmを求めることは、xの平均値を求めることに等しくなる。

1-10の間でfind_mean(X)関数を最小とするmの解(Xの平均)を探す。

optimize(find_mean(X), interval=c(1,10), tol=0.0001, maximum=FALSE)

$minimum
[1] 3.8

$objective
[1] 30.8

 minimunは推定したパラメータ値。

objectiveが、そのパラメータを代入したときの関数の値(最小値)。

 

 

(2)optim()による2パラメータ推定

最尤法により、正規分布の平均と分散を

Y <- rnorm(50)
mean(Y)     # -0.02269269
sd(Y)          # 1.032003

 

ベクトルx(正規分布を仮定)の平均mu、標準偏差sigの2つのパラメータを最尤推定する処理を行う場合に定義する。

find_norm_para <- function(x){
    return(function(para){
    mu <- para[1]
    sig <- para[2]

    # 対数尤度を返す
    return( - length(x)/2 * log(sig) - 1/2 * sum( (x - mu)^2 ) / sig )
    })
}

 

para[1]=mu=-1、para[2]=sig=3を初期値(適当に設定)として、それらの真の解(mean(Y)、sd(Y))を探す。

求めるパラメータが2つ以上の場合、optim(初期値, 関数. ...)を使用。

optim(par=c(-1,3), fn=find_norm_para(Y), control=list(fnscale=-1))
parが適当に設定する初期値。ベクトルで与える。
control=list(fnscale = -1) で最大化を指示できる。デフォルトは最小化。

$par
[1] -0.02262334  1.04377987

$value
[1] -26.07003

$counts
function gradient 
      67       NA 

$convergence
[1] 0

$message
NULL

parがパラメータの推定結果。

valueが推定したパラメータの場合の式の値(最大対数尤度)。

countsのfunctionが計算繰り返し回数。

convergenceには、収束したなら0が返る。

 

 

(3)optim()を使用したその他の例

data1

   score1 score2 score3 ans
1       0      0      1   0
2       1      1      1   1
3       2      1      0   0
4       1      0      0   0
5       0      0      0   0
6       0      1      2   1
7       1      0      0   0
8       2      0      0   0
9       0      1      2   0
10      0      0      2   1
11      0      0      0   0

このようなdata1があったとする。

各スコアは0,1,2を取る、よってscore1+2+3の最大は6となる場合がある。

scoreの重み付き線形和とansの差の二乗和が最小となる重みを求めたいとする。
重みの範囲は[0,1]とする。

要するに、p1、p2、p3が推定するパラメータだったとき(max_score=6)、solとdata1のans列の差の2乗和を最小にする問題。

sol = (p1*score1 + p2*score2 + p3*score3) / max_score
sum( (sol-ans)^2 ) ⇒ 最小化

 

これまで見てきたように、optim_fn(データ)(パラメータ)形式の関数と定義する。
max_score <- 6
optim_fn <- function(data){
    score1 <- data[,1]
    score2 <- data[,2]
    score3 <- data[,3]
    ans <- data[,4]

    return(function(para){
      p1 <- para[1]
      p2 <- para[2]
      p3 <- para[3]

      sol <- (p1*score1 + p2*score2 + p3*score3) / max_score
      return( sum((sol-ans)^2) ) # want to minimize this
    })
}

 

optim(par=para_list, fn=optim_fn(data1), method="L-BFGS-B", lower = 0, upper = 1)
methodに"L-BFGS-B"を指定すると、探索値の範囲指定(lower=0、upper=1)が行える。

$par
[1] 0.1818182 1.0000000 1.0000000

$value
[1] 1.434343

$counts
function gradient 
       6        6 

$convergence
[1] 0

$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

 

最適解がいくつあるか分からない問題(これはシンプルだが...)の場合、ランダムに複数回初期値を選んで計算させてみる。

10回実施。

n_try <- 10

n_para <- 3

res_list <- list()

count <- 1
for(i in 1:n_try){
    para_list <- runif(n_para)
    res <- optim(par=para_list, fn=optim_fn(data1), method="L-BFGS-B", lower = 0, upper = 1)

    res_listcount <- c(para_list, res$par, res$value)
    count <- count+1
}

結果をデータフレームにまとめる。
mtx1 <- matrix( unlist(res_list), ncol=7, byrow=TRUE )
res_data <- data.frame(st_para1=mtx1[,1], st_para2=mtx1[,2], st_para2=mtx1[,3],
end_para1=mtx1[,4], end_para2=mtx1[,5], end_para3=mtx1[,6], value=mtx1[,7])

res_data

     st_para1   st_para2 st_para2.1 end_para1 end_para2 end_para3    value
1  0.15277128 0.14019332  0.3100025 0.1818182         1         1 1.434343
2  0.96345322 0.29611349  0.6449935 0.1818182         1         1 1.434343
3  0.57237208 0.52925936  0.3130487 0.1818182         1         1 1.434343
4  0.52376898 0.34242409  0.1968113 0.1818182         1         1 1.434343
5  0.13490861 0.70259612  0.9709726 0.1818182         1         1 1.434343
6  0.42852674 0.74335816  0.6959915 0.1818182         1         1 1.434343
7  0.55659350 0.42611310  0.5567565 0.1818182         1         1 1.434343
8  0.63883635 0.99873597  0.5545664 0.1818182         1         1 1.434343
9  0.11357835 0.05827815  0.1648482 0.1818182         1         1 1.434343
10 0.08554884 0.15794897  0.7886483 0.1818182         1         1 1.434343

end_para1、end_para2、end_para3が推定した結果となる。今回は10回とも同じ。

 

以上。