Stataでモンテカルロシミュレーション
（円周率のシミュレーション）

2.1　事前準備（点数の設定、観測数の設定、乱数の種の設定）

はじめにclearコマンドでデータセットをクリアします。次にlocalコマンドでマクロ変数myobsに100を代入し、setコマンドでobsにmyobsを設定します。この操作によりデータセットの観測数（行数）が100行に設定されます。さらに、seedとして1を設定します。この後の無作為な数字を生成する際に使用されますが、この値を設定することで再現性を確保できます。以上で無作為に点を打つ準備が整いました。

2.2　一様分布から乱数を生成（“無作為な点”のX座標、Y座標として利用）

genコマンドで変数x, y, rを生成します。変数x, yには、runiform関数により0以上1以下の値が無作為に抽出され、代入されます。この値をXY座標上の点とみなし、正方形の範囲内に無作為に点を打ったと想定します。変数rにはxとyの二乗和の平方根を代入します。この値が1以下である点(x, y)は四分円の内側にある点とみなすことができます。

2.3　四分円の内側の点のみ描画（四分円の外側の点は描画しない）

twowayコマンドでグラフを作成します。プロットタイプとしてscatterを指定し散布図を作成します。Y座標を変数y、X座標を変数xとし、条件「rが1以下」により四分円の内側にある点のみを描画します。オプションysize(5) xsize(5)で5×5インチのグラフを出力します。

続いてsuコマンドで条件「rが1以下」を満たすrの統計量を出力します。全体で100の点を打ちましたが、四分円の内側にある点は85であることが読み取れます。

2.4　四分円の内側の点数の割合の4倍を出力

suコマンドを実行すると内部変数が作成されr(N) に観測数が格納されるので、マクロ変数aに代入します。四分円の内側の点の割合を4倍した値を計算し、displayコマンドで出力します。計算結果として3.4と出力されますが、この値は半径が1の円の面積すなわち円周率になることが期待されます。

上記の結果は観測数を100としているため一桁目のみ合っているという結果でしたが、観測数を指定しているmyobsの値を1,000や10,000などのより大きな数字に変更して再度実行することで円周率に近づいていくことを確認してみましょう。

4　まとめ

• Stataを用いて円周率をシミュレーションする方法を紹介しました。モンテカルロシミュレーションはこのように、試行回数を増やしていくと真の分布に近づきます。
• 観測数や繰り返しの回数を変更して再度実行し、変数pisの平均値が円周率に近い値となるか、ヒストグラムの形状がどのように変化するかなど、doファイルを編集してぜひお試しください。