optim()で階層構造を持った多変数の最適化をしてみる
最適化周りの処理について実装する必要が出たので，optim() を調べて使ってみましたよ，という話．
optim() って関数の形で表わせさえすれば，結構なんでも自由にできるっぽくて便利です．特に，階層的な構造を持ったデータの最適化もできるのは大きいです．
例
階層化ということで，ここで例に挙げるのは，比率の最適化です．イメージとしては，各クラスにおける学級委員長の選出を思い浮かべてもらえればと思います．
各クラスには複数の立候補者がいて，その支持率が既知のものとします*1．各クラスの生徒に対して，支持している候補者と，目指すクラスの方向性についてのアンケート聞いたとしてください．その得られたアンケート結果に対しなにがしかの係数をかけて，各立候補者の支持スコアを算出します．で，それを全立候補者のスコア合計で割ってあげれば，各立候補者の支持率が擬似的に出せます．
そうすると，クラスごとに算出した支持率と，実際の支持率の両者が得られます．アンケート結果にかける係数を調整することで，支持率の予測と実際の乖離を0に近づけましょう，ということになります．この場合，係数がかかるのは各クラスの生徒ですが，知りたいのはクラスごとの支持率という，クラス単位のパラメタだという階層構造になるわけです．
サンプルデータの作成
とりあえず，以下のようにクラスデータを作っておきます．わかりやすくするために，クラスは2つ，各クラスの候補者も2人としておきます．ここでポイントは，予め各候補者のカラムを作っておくことです．各生徒の支持対象はわかっているので，ここでは1または0の値で表しちゃいます*2．
library(dplyr) # 選択肢1, 2のデータを作成して結合 # 項目は以下の通り # クラス番号 # 候補者1支持ダミー # 候補者2支持ダミー # 価値観1 # 価値観2 location1 <- matrix(c(rep(1, 100), rep(1, 30), rep(0, 70), rep(0, 30), rep(1, 70), rep(1, 71), rep(0, 29), rep(0.3, 69), rep(0.7, 31)), 100, 5) location2 <- matrix(c(rep(2, 100), rep(1, 50), rep(0,50), rep(0, 50), rep(1,50), rep(0.1, 40), rep(0.8, 20), rep(0.5, 40), rep(0.2, 20), rep(1, 80)), 100, 5) d <- as.data.frame(rbind(location1, location2)) colnames(d) <- c("class_id", "person1", "person2", "value1", "value2") # クラス1, 2の支持率データ # 項目は以下の通り # クラス番号 # 候補者1支持率 # 候補者2支持率 s <- as.data.frame(t(matrix(c(1, 0.8, 0.2, 2, 0.6, 0.4), 3, 2))) colnames(s) <- c("class_id", "support1", "support2") # 結合してデータセットを作成 ds <- d %>% inner_join(s, by="class_id")
最適化対象の関数の作成
やりたいのは各候補者のスコアを出して，それを既知の支持率と比較することです．まずはわかりやすくするために，1クラス分のデータだけでやってみたいと思います．
# 1クラスぶんにデータを絞る ds1 <- ds %>% filter(class_id==1) # 最適化する関数 predict_rating <- function(b, d, print=FALSE) { d_tmp <- d d_tmp$predict1 <- b[1]*d_tmp$person1*b[3]*d_tmp$value1*b[4]*d_tmp$value2 d_tmp$predict2 <- b[2]*d_tmp$person2*b[3]*d_tmp$value1*b[4]*d_tmp$value2 # クラスidでグループ化して合計スコアを出した上で，比率に変換する d_tmp2 <- d_tmp %>% group_by(class_id, support1, support2) %>% summarize(., predict1=sum(predict1), predict2=sum(predict2)) total <- d_tmp2$predict1+d_tmp2$predict2 d_tmp2$predict1 <- d_tmp2$predict1/total d_tmp2$predict2 <- d_tmp2$predict2/total if (print) { print(d_tmp2) } # 最少化する関数は二乗誤差 sum((d_tmp2$support1-d_tmp2$predict1)^2+abs(d_tmp2$support2-d_tmp2$predict2)^2) }
optim() で実行
ここまで用意ができたら，実際にoptim() で実行します．この場合，最適化したい係数は非負で一定の範囲の実数に抑えたいとします．そのような指定ができる手法は，どうやら "L-BFGS-B*3" だけらしいので，これを選択します．各変数は0から10の範囲を取るという制約を与えます*4．各係数の初期値には，平均と分散が1の正規分布から乱数を生成し，負の値は正にするかたちで用いています．
res <- optim(abs(rnorm(4, 1, 1)), predict_rating, d=ds1, method="L-BFGS-B", lower=0, upper=10, control=list(maxit=10000)) res
これを実行すると，以下のように結果が得られ，実際に最適化がなされていることが見て取れます．
> res <- optim(abs(rnorm(4, 1, 1)), + predict_rating, + d=ds1, + method="L-BFGS-B", + lower=0, + upper=10, + control=list(maxit=10000)) > res $par [1] 3.2007261 0.5497433 0.5569501 0.7339547 $value [1] 1.603087e-14 $counts function gradient 10 10 $convergence [1] 0 $message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" > predict_rating(res$par, d=ds1, TRUE) Source: local data frame [1 x 5] Groups: class_id, support1 [?] class_id support1 support2 predict1 predict2 (dbl) (dbl) (dbl) (dbl) (dbl) 1 1 0.8 0.2 0.7999999 0.2000001 [1] 1.603087e-14
ここまで結果が得られたところで，最適化対象のクラス数を2にしてみましょう．ds1ではなくdsを使う，というだけですが．そうすると，見事に以下のように最適化された結果が得られています．もちろん2クラスに同じ係数を当てはめているため，クラス1については当てはまりが多少悪くなっています．それでもクラス1，2ともにそれなりによい当てはまりになっていることがわかるかと思います．
> res <- optim(abs(rnorm(4, 1, 1)), + predict_rating, + d=ds, + method="L-BFGS-B", + lower=0, + upper=10, + control=list(maxit=10000)) > res $par [1] 2.42237544 0.52544405 1.42063132 0.04326498 $value [1] 0.005155603 $counts function gradient 8 8 $convergence [1] 0 $message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" > predict_rating(res$par, d=ds, TRUE) Source: local data frame [2 x 5] Groups: class_id, support1 [?] class_id support1 support2 predict1 predict2 (dbl) (dbl) (dbl) (dbl) (dbl) 1 1 0.8 0.2 0.7600352 0.2399648 2 2 0.6 0.4 0.6313148 0.3686852 [1] 0.005155603
