The One with ...

授業・研究用のメモなどを公開しています

オブジェクトやファイル名をループ内で割り当てる(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.m0appendしていきます.

モデルを変える場合は,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とします.