IRIS: R 質問チャンネルで再現性を担保するための標本抽出器

概要

R の初心者の質問などで役立つように再現性を担保できるテストデータを生成する関数を思いついた. なお, R Advent Calendar とは全く関係ない.

イントロダクション

r-wakalang では日々 R に関する質問と回答が行われている.

以下のスライドでも言及されているように, 質問に対して正確に問題の原因を把握し有効な回答を提示するには実行時の状況を再現できることが重要である.

再現性に関する解説

冒頭に示したスライドでも紹介されているように, データフレームからコードを逆生成する datapasta, 環境情報を出力する seessionInfo()関数は必要である*1. また, より厳格に再現性を担保したいなら packrat パッケージも有効だろう.

しかし, 中には再現に注意を払わねばならない話がある. 例えば, 疑似乱数は実行するたびに違う結果を返すため, 乱数生成を伴うプログラムは set.seed() 関数で乱数を固定する必要がある. しかし実は乱数シード値が保証する再現性は同一環境での乱数生成であり, 異なるシステム上では異なる結果を返すことがある*2.

R でも, 実は一部の乱数が環境依存である. 具体的には sample() がそれで, 公式のニュースリリースにあるように R 3.6.0 以前と以降では, 離散乱数生成器はシード値を固定してもバージョンによって出力される並びがかわる*3*4. packrat はパッケージのバージョン管理であり, R 本体のバージョン管理はできない.

set.seed(42)
sample(1:10, 10)

##  [1]  1  5 10  8  2  4  6  9  7  3

「初心者の質問」特有の問題

r-wakalang では必然的に初心者の質問が多い. そのため, プログラムの再現性とはなんなのかという認識がまだできていない*5. 現在は再現可能性のための様々なRパッケージが開発されているが, それらを初心者にインストールしておけと一喝するのは酷であるし, OS 依存の問題は解決が難しい. また質問によっては外部に公開できないデータセットが絡むこともある.

r-wakalang に寄せられる R の問題の多くは, データフレームを意図した形式に加工するにはどうするかというものである. R の特徴を考えれば自然なことである.

さて, このデータフレームに関する問題に, OS やバージョンによって変化せず, 再現性をもたせる方法とは何か. その1つは, iris である.

data(iris)
head(iris)

Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
5.1	3.5	1.4	0.2	setosa
4.9	3.0	1.4	0.2	setosa
4.7	3.2	1.3	0.2	setosa
4.6	3.1	1.5	0.2	setosa
5.0	3.6	1.4	0.2	setosa
5.4	3.9	1.7	0.4	setosa

iris は現代統計学の基礎を作った1人であるフィッシャーの利用した有名なデータセットで, R ならば追加パッケージもなく利用できる. しかし, 150件しかないため, 大きなデータを扱う問題で使用しづらい. よって, ある程度ランダムでありながら, 環境に依存せず再現性を担保し, かつデータの大きさに柔軟に対応できる方法が求められている.

補足:「再現性」と「再現可能性」

再現性と再現可能性はどちらも reproducibility の訳語であるが, この単語は文脈によって意味合いが変わってくるので注意が必要だ. プログラミングで言えば再現性*6とは単に同じプログラムを実行して同じ結果を得られることを意味することが多いが, R の活躍する分野である実証研究ではより広い意味を持つ*7. 後者の文脈では同じ結果が出るとは言っても, 実際には同じデータでなくても同じ結果がでるということであり, 偶然性のゆらぎを超えた法則の一般性, いわば外的妥当性のことも含めることが多い(この点はコンセンサスが取れてるか怪しいので, 言い切るのは早計かもしれない). 今回はプログラミングの質問の話なので, 前者である.

再現可能性のすゝめ (Wonderful R 3)

不変的で再現可能な入門用標本抽出器, IRIS の提案

そこで今回, 初心者の質問に役に立つ「不変的で再現可能な入門用標本抽出器」を考案した. この名称はiris データセットにちなんで IRIS: Invariant Reproducible Introductory Sampler とした. 疑似乱数生成器と同様, シード値の固定で同じデータセットを作成できる上に, 異なるOSでも再現性が保証される. 実装には組み込み関数のみを使用しているので, 将来のアップデート後も再現できる可能性が期待できる.

rand_iris <- function(n = 150, n_numeric = 4, n_character = 1,
                      seed = NULL, M = 2^31-1, a = 16807){
  if(is.null(seed)) seed <- as.numeric(Sys.time())
  else seed <- as.integer(seed)
  
  rng_iris <- function(n, seed, M, a){
    x <- rep(NA, n + 1)
    x[1] <- seed
    for(i in 2:(n+1)){
      x[i] <- (a * x[i-1]) %% M
    }
    return(x[-1] %% 150)
  }
  d <- list()
  if(n_numeric > 0){
    iris_cols <- rng_iris(n_numeric, seed, M, a) %% 4 + 1
    for(i in 1:n_numeric){
      d[[paste0("Num", i)]] <- iris[rng_iris(n, seed, M, a), iris_cols[i]]
    }
  }
  if(n_character > 0){
    iris$Species <- as.character(iris$Species)
    for(i in 1:n_character){
      d[[paste0("Chr", i)]] <- iris$Species[rng_iris(n, seed, M, a)]
    }
  }
  as.data.frame(d)
}

rand_iris(n = 10)

Num1	Num2	Num3	Num4	Chr1
3.3	2.5	3.3	6.0	virginica
3.0	1.5	3.0	4.5	versicolor
2.8	2.1	2.8	5.6	virginica
2.9	1.3	2.9	4.3	versicolor
4.4	0.4	4.4	1.5	setosa
2.6	2.3	2.6	6.9	virginica
3.2	1.8	3.2	4.8	versicolor
3.6	0.2	3.6	1.0	setosa
2.7	1.4	2.7	3.9	versicolor
3.0	1.4	3.0	4.6	versicolor

解説

n=150 は出力するデータフレームの行数である. n_numeric=5 は, 出力するデータセットの数値変数の列数, n_character は文字変数の列數である. そして, seed=NULL は, この関数独自のシード値である. rand_iris() は通常の set.seed() には依存しないし, rnorm() などの従来の乱数生成関数にも影響しない. 環境依存しないように乗算合同法 (multiplicative congruential method) という古い疑似乱数生成アルゴリズムを新たに実装した. このアルゴリズムは現在よく使われるメルセンヌ・ツイスタ法よりも周期が短いため, 厳密なシミュレーションには不適だが, この通り簡単に実装できる. さらに馴染みのあるデータセット iris の値を参照することで, 数値にもカテゴリカル変数にも対応できる.

乗算合同法は線形合同法に比べるとマイナーかもしれない. 乗算合同法に関する説明は以下が参考になる.

意味がわからない人へのネタバレ

iris のデータセットをバージョンやOSをまたいで再現できる乱数に基づいておもちゃにできるような関数を書いた. すでに述べたように, R は古いバージョンだと乱数の結果が変わることがあるので*8, より再現性を徹底するなら, みんなよく知っていてどのバージョンでも内容の変らない iris を使うべきでは, と考えこういうネタを思いついた. 乗算合同法は実装は簡単だが現在は役に立つ場面がほとんどない. しかし初心者に長いコードをコピペしろ, XXというパッケージをインストールしてロードしろ, というのは酷なので, 組み込みの関数のみでかつ短くした. しかし作っているうちに「これ別にあんまり面白くないな」と気づいたので日和見して, バージョン依存のケースも実際ありうるので set.seed() だけじゃ再現性の担保が不十分じゃないの, というちょっと役に立つ話も載せてみたり「初心者向け」と称してクオリティのなさをごまかしりした.

その他

この記事は rmakrdown/knitr の練習も兼ねて書いてみた. 私はブログを現状 LyX で書いており, Word との機能比較を書いていたりする.

その環境がだいぶ整備されているのでながらく rmarkdown を使ってこなかったが, knitr もバックグラウンドは LaTeX + pandoc で回っており, RStudio上では数式も図表も表示できるため, 細かい体裁の指定にこだわらなければ pdf 生成にそれなりに便利であることがわかった. LyX より明確に優れた点としては, knitr::kable() を使えば表の作成が非常に簡単なこと, プログラムのシンタックスハイライトが見やすい (listings のハイライトはかなり変だ.) ことが挙げられる. 一方で, 私は tex としてもブログに投稿する html としても1つのソースから生成できるような機能がほしかったのだが, table, figure 環境したりに付番したり, 相互参照したりするのは難しそうだ. さらに具体的な課題を書くと, 私が投稿する可能性のある媒体ははてなブログとqiitaで, どちらもかなり独自の仕様にしている*9ので, それらの対応も大変そうだ.

---

*1:このデータをこう変形したい, 程度の質問ならデータフレームだけでもよい. 以下2つに書かれているように, 情報が多すぎるとかえって煩わしくなることもある. 問題と関係のありそうな箇所を自分で試行錯誤しながら切り詰めるのも, 自己解決能力の鍛錬になる.

*2:C++ でも, メルセンヌ・ツイスタ乱数生成器の実装が環境依存であることが指摘されている. C++のmt19937は環境非依存だがdistributionは環境依存 - Qiita

*3:ダウングレードして確認するのが億劫なら, 現時点では paiza.io のRのバージョンが 3.4 なので確認してみると良い

*4:私は以前気づいたのだが原因までは特定していなかった. R-wakalang でこのことをなんとなく投稿したら tobcap 氏に指摘された.

*5:インターネット上で初心者の質問を受け付けるコミュニティでは昔から「エスパー」という言葉が使われた. 初心者は問題の原因をさがすとっかかりすら全くわからないので, 原因特定に役立つ情報を提示しないことが多い. こういうなかで原因が特定できる情報が出揃う前に原因を言い当てる者をエスパーと呼んだり, 情報を提示しない質問者に対して「俺達はエスパーじゃないんだ」返信が付いたりした.

*6:reproducibility というよりは reproducing とかそういう表現のほうが多いかもしれない.

*7:例えば高橋他 (2018) 『再現可能性のすゝめ』の意図するところは前者に近い.

*8:実際, R-wakalang 質問者の中には古いバージョンを使っていたのかこちらと結果が変わる場合があった. その質問では乱数の結果の違いが問題の原因ではなかったが.

*9:前者なら「見出し」が実際には h3 タグだとか, 後者なら id 属性によるアンカーが無効になるとか