R 最適化計算 - optimize(), optim()
optimize(), optim() の使い方に関するメモ。
（１）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が、そのパラメータを代入したときの関数の値（最小値）。
（２）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が返る。
（３）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回とも同じ。
以上。
