はじめに：アクチュアリーの「飛行シミュレーター」

CS1のカリキュラムの中でも、最も実用的な章へようこそ！Paper Bでは、単に確率を計算するだけでなく、モデルを構築することになります。確率変数のシミュレーションは、アクチュアリーにとっての「飛行シミュレーター」だと考えてください。保険会社が新商品を発売する前に、何千もの支払備金請求のパターンをシミュレーションして、会社が支払余力を維持できるか（ソルベンシー）を「テスト飛行」して確認します。この章では、Rの組み込み関数や逆関数法（Inverse Transform Method）を使って、これらの「テスト飛行」を行う方法を学びます。

1. 「r」ファミリーの関数

Rでは、ほとんどの確率分布に乱数を生成するための専用関数が用意されています。これらの関数は常に、"random"（ランダム）を意味する文字 r で始まります。

基本的な構造は通常、r[分布名](n, [パラメータ]) となります。ここで \(n\) は生成したい観測値の数です。

代表的な離散分布

  • 二項分布: rbinom(n, size, prob) — \(X \sim Bin(size, prob)\) から \(n\) 個のサンプルを生成します。
  • ポアソン分布: rpois(n, lambda) — \(X \sim Pois(\lambda)\) から \(n\) 個のサンプルを生成します。
  • 幾何分布: rgeom(n, prob) — 最初の成功が起こるまでの「失敗」の回数を生成します。
  • 負の二項分布: rnbinom(n, size, prob) — 指定された回数の成功が起こるまでの「失敗」の回数を生成します。

代表的な連続分布

  • 正規分布: rnorm(n, mean, sd) — Rでは分散 \(\sigma^2\) ではなく、標準偏差 \(\sigma\) を使用することに注意してください。
  • 指数分布: rexp(n, rate) — 強度パラメータ \(\lambda\) を使用します。
  • ガンマ分布: rgamma(n, shape, rate) — scale (\(1/rate\)) を使用することもできます。
  • 一様分布: runif(n, min, max) — 下限と上限の間の値を生成します。
  • 対数正規分布: rlnorm(n, meanlog, sdlog)。
  • ベータ分布: rbeta(n, shape1, shape2)。

クイック復習: 平均が5のポアソン分布から100件の請求をシミュレートする必要がある場合は、rpois(100, 5) を使用します。

2. 再現性：set.seed() 関数

コンピュータは、実は本当の意味での「ランダム」な数字を生成しているわけではありません。アルゴリズムを使用して「擬似乱数」を作成しています。そのため、rnorm(5, 0, 1) を2回実行すると、異なる結果が得られます。これは、あなたの答案を採点しなければならない試験官にとっては問題となります！

解決策： set.seed() を使いましょう。シミュレーションコードの前に set.seed(123)（または任意の整数）と入力することで、コードを実行するたびにRが「全く同じ」乱数を生成するように設定できます。

試験での重要なヒント： 試験問題で「シード値（seed）」が指定されているかどうか、必ず確認してください。「シード値42を使用せよ」と指示がある場合、コードの最初の行は必ず set.seed(42) でなければなりません。

3. 逆関数法 (ITM)

特定の「r」関数を使わずに、手動で確率変数を生成するように求められることがあります。逆関数法は、離散変数と連続変数の両方でこれを実現するための標準的な手法です。

逆関数法のロジック

すべての累積分布関数 (CDF) \(F(x)\) は、0から1の間の値を出力します。もし 一様分布 Uniform(0,1) から乱数 \(U\) を取得すれば、その確率に対応する \(x\) の値を「逆算」して見つけることができます。

ステップ・バイ・ステップの手順：

  1. runif(1) を使って、\(U \sim Uniform(0,1)\) から乱数 \(u\) を生成する。
  2. \(F(x) = u\) と置く。
  3. 逆関数を計算して \(x\) について解く： \(x = F^{-1}(u)\)。

指数分布の例：
累積分布関数は \(F(x) = 1 - exp(-\lambda x)\) です。
1. \(u = 1 - exp(-\lambda x)\) と置く。
2. \(1 - u = exp(-\lambda x)\)
3. \(\ln(1 - u) = -\lambda x\)
4. \(x = -\frac{1}{\lambda} \ln(1 - u)\)

Rで、一様乱数のベクトル u がある場合、これをコードにすると次のようになります： x <- - (1/lambda) * log(1 - u)

知っていましたか？ \(U\) と \(1-U\) はどちらも \(Uniform(0,1)\) 分布に従うため、アクチュアリーは式を簡略化して \(x = -\frac{1}{\lambda} \ln(u)\) とすることもあります。

4. 既存データからのサンプリング

シラバスでは、理論的な分布からだけでなく、既存のデータセットからサンプルを生成することについても触れています。これには sample() 関数を使用します。

sample() の主な引数：

  • x: サンプリング元のデータまたはベクトル。
  • size: 取り出す個数。
  • replace: 取り出した後に元に戻すかどうか（TRUE = 復元抽出、FALSE = 非復元抽出）。

実社会での応用：

ブートストラップ法： 推定量（バイアスや分散など）の性質を推定するために、「ブートストラップ法」を用いることがあります。これは、手元のデータから復元抽出（replace = TRUE）を行う手法です。これにより、母集団から新しいサンプルを取り直すことを模倣します。

置換検定（パーミュテーション・テスト）： ノンパラメトリックな仮説検定を行うために、しばしば非復元抽出（replace = FALSE）を行ってデータを「シャッフル」し、得られたパターンが偶然によるものかどうかを確認します。

5. シミュレーションと理論の比較

シミュレーションデータが得られたら、それを既知の分布と比較することがシラバスで求められています。特に中心極限定理 (CLT) に関連して、正規分布と比較することが多いです。

任意の分布から1,000個の標本平均をシミュレートした場合、中心極限定理によれば、それらの平均の分布は正規分布の曲線に似てくるはずです。Rでは、以下の手順でこれを確認できます：

  1. シミュレーションデータのヒストグラムを描く： hist(sim_data)
  2. 密度曲線を重ねる： lines(density(sim_data))
  3. Q-Qプロットを使用して、点が直線上に並ぶか確認する（詳細は「探索的プロット」の章で扱います）。

重要なポイント： シミュレーションを使うことで、中心極限定理を視覚的に検証できます。元のデータが（指数分布のように）大きく歪んでいても、サンプルサイズが大きくなるにつれて、それらのシミュレーションの「平均」は「正規分布」に近づいていくのです。

まとめ ＆ クイック復習

回避すべきよくある落とし穴：

  • 標準偏差 vs 分散： rnorm() には \(\sigma\) が必要ですが、多くの試験問題では \(\sigma^2\) が与えられます。必要に応じて sqrt() を使いましょう！
  • set.seed() の忘れ： シード値を設定しないと、数値の答えが採点基準（マークスキーム）と一致しなくなります。
  • 逆関数法 (ITM) のロジック： 逆関数法は常に runif() から始まることを覚えておきましょう。「確率」を「値」に変換しているのです。

注：これらのシミュレーションサンプルの特定の確率や要約統計量の計算については、「Rにおける確率、分位点、および要約統計量」の章を参照してください。