オブジェクトやファイル名をループ内で割り当てる(R)
ファイル名を直書きするのがめんどうくさい
50個のデータセットを使い複数のモデルで推定して自由エネルギーを計算する必要があったのでコードを書いてみました.
多分忘れるだろうから,ここに貼り付けておきます.
fevalues.m0 <- c() for (i in 1:50) { stanfit <-stan(file='model0.stan',iter=10500,warmup=500,data=f2(i),chains=4) saveRDS(stanfit,c(paste('fit/fe', i,'_0.rds', sep = ''))) fevalues.m0 <- append(fevalues.m0,-logml(bridge_sampler(stanfit,method="warp3")) ) }
ここで関数f2( )
はstan用に各データセットをlist型に整形する関数です.
stan( )
の引数data
に関数の出力を渡しても問題ないようです.
途中でRがコケてもstanfitが残るようにsaveRDS
で逐次保存します.このときpaste( )
をfor
ループで回してデータセットに対応したファイル名を自動作成します(10個くらいなら直打ちのほうが早いでしょう.50×3個のファイル名はさすがに書く気がしない).
saveRDS
内でpaste( )
を使っても問題ないようです.
stanfit
は上書きされて1つしか残らないため,メモリ不足でRがクラッシュすることを防ぎます(多分)
計算した統計量(ここではブリッジサンプリングから計算した自由エネルギー)は空のベクトルfevalues.m0
にappend
していきます.
モデルを変える場合は,stan( )
のモデル部分だけ変え,上記のfor文をコピペして使い回します.
1回のループにまとめることもできますが,かえって可読性が悪くなると判断しました.
このくらいで分割した方が,数ヶ月後の自分には理解しやすいでしょう.
呼び出したstanfitからの計算
readRDS
で保存したstanfitを使ってbridge_sampler( )
を適用すると,「modelがないぞ」とrstan
に怒られます.
「コンパイルしたモデルは保存する設定にしとるやろがい!」と悪態をつきたくなる気持ちをぐっとこらえ,stanfitに対応したモデルをchains=0
で再コンパイルします(これは時間はかかりません).
mod<-stan(file='model0.stan',iter=10500,warmup=500,data=f2(i),chains=0) -logml(bridge_sampler(stanfit,mod,method="warp3"))
これなら計算してくれます.
他に解決策があるのかもしれませんが,今回はこれでOKとします.