博士

やあ、ロボ子！今日はRコードの再現性に関する面白い記事を見つけたのじゃ。

ロボ子

博士、こんにちは。再現性ですか、科学の世界では非常に重要な概念ですね。

博士

そうじゃろう？特に多変量正規分布のサンプリングで問題が起きやすいらしいぞ。`set.seed()`を使っても、違うマシンだと違う乱数が生成されることがあるらしい。

ロボ子

それは困りますね。`set.seed()`は万能ではないということでしょうか？

博士

原因は浮動小数点演算にあるみたいじゃ。わずかな誤差が、OSやコンパイラの違いで変わってくるらしい。

ロボ子

浮動小数点演算は、厳密な算術演算とは異なるのですね。それで、`MASS::mvrnorm()`関数で特に顕著になると。

博士

`MASS::mvrnorm()`は固有値分解を使うからの。ちょっとした入力の変化で固有ベクトルの符号が反転したりして、結果が大きく変わることがあるんじゃ。

ロボ子

固有値分解は不安定なのですね。記事では、Cholesky分解を使う`mvtnorm::rmvnorm()`が推奨されているようですが。

博士

そうそう！Cholesky分解は正定値行列に対して一意な上三角行列を生成するから、固有値分解より安定しているんじゃ。

ロボ子

`MASS::mvrnorm()`と`mvtnorm::rmvnorm()`の違いは他にもあるみたいですね。乱数の入力方法が違うとか。

博士

`MASS::mvrnorm()`は乱数を列ごとに入力するから、サンプル数を変えると以前のサンプルも変わっちゃう。`mvtnorm::rmvnorm()`は行ごとだから、以前のサンプルは変わらないぞ。

ロボ子

なるほど、それは重要な違いですね。`mvtnorm::rmvnorm()`を使う際は、`method = "chol"`を設定すると良いのですね。

博士

その通り！あと、条件数は行列の反転における計算誤差の指標だけど、この問題には直接関係ないから、過信しないように、じゃ。

ロボ子

勉強になります。pivotオプションについても触れられていましたね。これは数値安定性のために行を入れ替える機能で、特異に近い行列の場合に有効なのですね。

博士

そうじゃ！まとめると、サンプリング直前に`set.seed()`を呼び出し、`MASS::rmvnorm()`の代わりに`mvtnorm::rmvnorm()`を使い、`method = "chol"`を設定する。これで再現性の問題はかなり軽減されるはずじゃ。

ロボ子

ありがとうございます、博士。これで安心して多変量正規分布のサンプリングができます。

博士

ところでロボ子、もしRコードが全然再現されなかったらどうする？

ロボ子

ええと…原因を徹底的に調査します。

博士

ざんねーん！　Re: 現象（再現）を諦める！…ってね！