= 3. 考察

== 課題(1)

=== 乱数の平均値が0.5になる理由、分散の値が1/12になる理由

  まず、一様乱数の確率密度関数を考える。一般に、一様乱数の確率密度関数は次のように表される。

  $ 
   f(x) = 
   cases(
    1 / (b-a) space ( a <= x <= b ),
    0 space ( "otherwise" )
   )
  $ <eq:uniform>

  これを用いて、平均値$E(X)$と分散$V(X)$を求めると、次のようになる。

  $
   E(X) &= integral_(-infinity)^infinity x f(x) d x = integral_(a)^b x f(x) d x \
    &= integral_(a)^b x / (b-a) d x \
    &= 1 / (b-a) [x^2 / 2]_(a)^b \
    &= (b + a) / 2
  $

  $
   E(X^2) &= integral_(-infinity)^infinity x^2 f(x) d x = integral_(a)^b x^2 f(x) d x \
    &= integral_(a)^b x^2 / (b-a) d x \
    &= 1 / (b-a) [x^3 / 3]_(a)^b \
    &= (b^2 + a b + a^2) / 3
  $

  $
   V(X) &= E(X^2) - E(X)^2 \
    &= (b^2 + a b + a^2) / 3 - (b + a)^2 / 4 \
    &= (b^2 - 2 a b + a^2) / 12 = (b - a)^2 / 12
  $

  よって、今回の場合、$a = 0$、$b = 1$ であるため、平均値$E(X)$、分散$V(X)$はそれぞれ$0.5$、$1/(12)$になることがわかる。


== 課題(2)

=== buffonの針で、何故$pi$の値が求められるのか

  まず、buffonの針問題とは、床の上に等間隔に引かれた平行線群があり、その距離を$2a$として，この床に長さ$2l( l < a)$の針をデタラメに落とすと，その針が平行線と交わる確率はいくらかという問題である。
  buffonの針で何故$pi$の値が求められるかについて以下の図を基に考える。

  #figure(
    image("../figs/buffon _fig.png", width: 60%),
    caption: [buffonの針問題],
  ) <fig:buffon_fig>

  針の中点を$y$、針と平行線のなす角を$theta$、針の長さを$2l$とする。対称性より、$y$の区間を$[0, a]$、$theta$の区間を$[0, pi/2]$とする。この時、針が平行線と交差する条件は、

  $ y <= l sin theta $

  と表される。ここで針をデタラメに落とすということは、$y$と$theta$は一様分布に従うと解釈できる。$(y, theta)$平面を考えると、針が平行線と交差する確率は、$y = l sin theta$を$theta$の領域で積分した面積を、全体の面積で割ったものになる。全体の面積は、$y$の区間が$[0, a]$、$theta$の区間が$[0, pi/2]$であるため、$(a pi)/2$になる。また、針が平行線と交差する領域は、

  $ integral_(0)^(pi/2) l sin theta d theta = l $

  となる。よって、針が平行線と交差する確率$P$は、

  $ P = l / ((a pi)/2) = (2l) / (a pi) $ <eq:prob>

  と求めることができる。
  今回の実験では、$y$の区間を$[a, 3a]$、$theta$の区間を$[0, pi]$とし、針の下端$y_1$、上端$y_2$を

  $ 
   cases(
    y_1 = y - l sin theta,
    y_2 = y + l sin theta
   )
  $ <eq:y>

  として、$y_1$と$y_2$の積が負になるとき、針が平行線と交差すると判定している。交わった回数を$n$、全試行回数を$N$とすると、@eq:prob より、

  $ pi = (2l) / (a P) = (2l N) / (a n) $

  と$pi$を求めることができる。[参考：@text2024e5]


=== フローチャート

  今回buffonの針問題を解くためのプログラムのフローチャートを以下に示す。ここで、$a = 1.0$、$l = 0.8$に設定した。そして、@eq:prob、@eq:y を用いて、$pi$の計算を行った。

  #figure(
    image("../figs/flowchart.png", width: 60%),
    caption: [buffonの針問題のフローチャート],
  ) <fig:buffon_flowchart>

== 課題(3)
=== 乱数から$pi$を求める⽅法とその精度

  $pi$を求めるのに用いられる方法として、入門的モンテカルロ法と当たり外れのモンテカルロ法がある。これら2つの方法について述べる前に、大数の法則について説明する。大数の法則とは、試行回数（一般的には標本すう）を十分大きくすると、確率変数の平均値が基の分布（母集団）の平均値に収束するという法則である。ここで、$[a, b]$で定義された関数$f(x)$と$[a, b]$に一様分布する確率変数$r$を考える。この時、$r$の確率密度関数$rho(r)$は、@eq:uniform と同様の式で表される。この時、$f(x)$の平均値$E(f(x))$は、

  $ E(f(x)) = integral_(a)^b f(x) rho(r) d r = I / (b - a) $ <eq:mean>

  となる。ここで、$I$は$f(x)$の区間$[a, b]$での積分値である。また、大数の法則より、十分大きい$N$に対して、$N$個の一様乱数$r_i$を発生させたとすると、$f(x)$の平均値は、

  $ E(f(r)) tilde.eq 1 / N sum_(i=1)^(N) f(r_i) $ <eq:mean2>

  となる。そして、@eq:mean と@eq:mean2 より、$I$は、

  $ I tilde.eq (b - a) / N sum_(i=1)^(N) f(r_i) $ <eq:mean3>

  と近似できる。これが入門的モンテカルロ法による求積である。

  次に、$0 <= f(x) <= c $として、$0 <= y <= c$、$a <= x <= b$の領域を考える。この領域に、$N$個の一様乱数$(x_i, y_i)$を発生させ、$f(x_i) >= y_i$の時の数を$n$とする。この時、領域と$f(x)$より下の領域の面積比は、

  $ 1 /(c (b - a)) integral_(a)^b f(x) d x tilde.eq n / N $ <eq:monte>

  となる。この式より、$I$は、

  $ I tilde.eq c (b - a) n / N $ <eq:monte2>

  となる。このように、$f(x_i) >= y_i$に当てはまるかどうかで積分値を求める方法を当たり外れのモンテカルロ法という。

  ここで、$f(x) = sqrt(1 - x^2)$、$a = 0$、$b = 1$とすることで、

  $ pi = 4 I $

  と求めることができる。

  次に、これら2つの方法の制度について考える。@fig:monte_hit_comp で見たように、どちらの手法も初めは回数が増えるにつれて精度が向上しているが、試行回数が$10^5$を超えてからは、精度の向上が$0$付近で緩やかになっている。大数の法則を考えると、初めは、試行回数が少ないため、平均値が母集団の平均値に収束していく様子がはっきりと見られるが、そこから試行回数が増えていくと、乱数の分布が理想的なものに近づくため誤差が$0$付近で一定になると考えられ、シミュレーション結果もそれを裏付けている。また、@fig:buffon_result、@fig:monte_result、@fig:hit_result を比べると、$30000$回の時に、$pi$の値の誤差が、入門的モンテカルロ法、当たり外れのモンテカルロ法の方がbuffonの針よりも小さく、精度が良いことがわかる。よって、小さい$N$($10^5$以下)なら、入門的モンテカルロ法を使い、大きい$N$($10^5$以上)なら、入門的モンテカルロ法、当たり外れのモンテカルロ法のどちらかを使うことで、精度良く$pi$の値を求めることができると考えられる。

== 課題(4)
=== ランダムウォークの理論値と実験値の比較

  @fig:randomwalk_comp を見ると、ランダムウォークの理論値と実験値は酷似していることがわかる。さらに、各$X$の値における理論値と実験値の相対誤差を計算した結果を@tab:randomwalk_error に示す。 

  #figure(
  caption: [ランダムウォークの理論値と実験値の比較],
  table(
  columns: 2,
  stroke: (none),
  table.header(
    [$X$], [相対誤差  "/%"]
  ),
  table.hline(),
  [-10], [3.02],
  [-8], [0.362],
  [-6], [0.363],
  [-4], [0.0192],
  [-2], [0.145],
  [0], [0.143],
  [2], [0.226],
  [4], [0.00469],
  [6], [0.367],
  [8], [1.01],
  [10], [2.41],
  table.hline(),
  )
  ) <tab:randomwalk_error>

  この表から、$X$の絶対値が大きくなるにつれ、相対誤差が大きくなるという傾向があることがわかる。これは、大数の法則を考えると、試行回数が多くなるにつれて平均値が理論値に収束していくため、平均値付近の値は理論値に近づくが、平均値から離れるにつれて、誤差が大きくなると考えられる。また、実際に平均値付近の$X$での誤差を確認してみると、$X$の絶対値が$4$の場所が誤差が最も小さくなっている。これは、発生させた乱数が、確率$1/2$で右、左に移動する様子を完全に再現できていないためであると考えられる。何故なら、確率が$1/2$からズレると、平均値も同時にずれてしまう。そのため、$X$の絶対値が$4$の場所付近が平均値となり、大数の法則により、誤差が最も小さくなっていると考えられる。

=== 二項分布に従う確率分布の理論値と実験値の比較

  @fig:binomial_comp だけを見ると、理論値と実験値は酷似しているように見える。より精度を確認するために、各患者数における理論値と実験値の相対誤差を計算した結果を@fig:binomial に重ねて表示させた。その結果を@fig:binomial_error に示す。この図から横軸の値が$0, 1, 2$の時は、相対誤差が$1$になっていることがわかる。しかし、この結果は妥当であると言える。なぜなら、患者数が$0, 1, 2$の時は、理論値が$0$に近く、オーダーが$10^(-9)$よりも小さいため、試行回数が$10^8$の時に実験値が$1$を超えることがないため、実験値が$0$として観測され、相対誤差が$1$になっていると考えられる。これを踏まえて、患者数が$0, 1, 2$の時の相対誤差を除いた結果を@fig:binomial_error2 に示す。この図から、期待値に近い患者数の時は、相対誤差が小さく、患者数が多くなるにつれて、相対誤差が大きくなるという傾向があることがわかる。しかし、患者数が$3$の時は、誤差が小さくなっている様子が見られる。これは、患者数が$0, 1, 2$の時と同様にオーダーの問題で、理論値が$10^(-8)$なため、試行回数が$10^8$の時に実験値が1の位になるため、ずれによる変動を受けにくいと考えられ、それに対し、患者数が$4,5$などの場所では、オーダーが$10^(-7)$なため、桁数が2桁になるため、1の位がずれの影響を受けやすいと考えられる。その結果、患者数が$3$の時の相対誤差が小さくなっていると考えられる。
  
  #figure(
  image("../figs/patients_error1.png", width: 90%),
  caption: [二項分布の理論値と相対誤差],
  ) <fig:binomial_error>


  #figure(
  image("../figs/patients_error2.png", width: 90%),
  caption: [二項分布の理論値と相対誤差(患者数が0, 1, 2の時を除く)],
  ) <fig:binomial_error2>

= 結論

  今回の実験では、乱数を用いて$pi$を求める方法、buffonの針問題、ランダムウォーク、二項分布に従う確率分布をシミュレーションするプログラムを実行し、その結果を比較した。その結果、$pi$を求める方法として、入門的モンテカルロ法、当たり外れのモンテカルロ法があることがわかった。また、これら2つの方法の精度について考察した結果、試行回数が$10^5$を超えると、精度の向上が緩やかになることがわかった。また、$30000$回の時に、$pi$の値の誤差が、入門的モンテカルロ法、当たり外れのモンテカルロ法の方がbuffonの針よりも小さく、精度が良いことがわかった。そして、ランダムウォーク、二項分布に従う確率分布の理論値と実験値の比較を行った結果、理論値と実験値は酷似していることがわかった。また、理論値と実験値の相対誤差を計算した結果、平均値から離れるにつれて、誤差が大きくなるという傾向があることがわかった。しかし、平均値付近の値は理論値に近づくため、誤差が小さくなることがわかった。