====== 【数値情報処理a】第5回 乱数とモンテカルロ法 ======

===== 前回の課題 =====

[[ans4]]

===== 乱数 =====

__**乱数**__はバラバラの数であり、理想的にはどの数が出るか予想ができず規則性も再現性も持たない。
例えばサイコロ1個を転がして出る目は 1 から 6 までの数字のどれがでるか予想できない乱数である。

コンピュータでは様々な場面で乱数が用いられる。
コンピュータで用いる乱数は、あるアルゴリズムで確定的に発生させるので__**規則性**__と__**再現性**__を持ち、理想的な乱数ではない。
これを__**疑似乱数**__という。
疑似乱数の発生アルゴリズムは色々と提案されている。

乱数は統計処理やシミュレーションでよく用いられるが、大量の理想的な乱数をすぐに得るのは難しい。
理想的な乱数ではないが、疑似乱数はコンピュータで短時間に大量に発生できるというメリットがある。

 
===== Rの疑似乱数 =====
R では様々な確率分布の疑似乱数を発生させる関数が用意されている。

^  確率分布                                  ^  乱数の関数           ^  引数                                 ^
| 一様分布（$0<x<1$）                    | ''runif(n)''          |                                       |
| 一様分布（min$<x<$max）                | ''runif(n,min,max)''  | ''min''：最小値, ''max''：最大値      |
| 二項分布 $\mathrm{B}(m,p)$             | ''rbinom(n,m,p)''     | ''m''：試行の回数, ''p''：成功の確率  |
| 標準正規分布 $\mathrm{N}(0,1)$         | ''rnorm(n)''          |                                       |
| 正規分布 $\mathrm{N}(m,s^{2})$         | ''rnorm(n,m,s)''      | ''m''：平均, ''s''：標準偏差          |
| 対数標準正規分布 $\Lambda(0,1)$        | ''rlnorm(n)''         |                                       |
| 対数正規分布 $\Lambda(m,s^{2})$        | ''rlnorm(n,m,s)''     |                                       |
| スチューデントのt分布 $\mathrm{t}(k)$  | ''rt(n,k)''           | ''k''：自由度                         |
| カイ二乗分布 $\chi^{2}(k)$             | ''rchisq(n,k)''       | ''k''：自由度                         |
| F分布 $\mathrm{F}(k_{1},k_{2})$        | ''rf(n,k1,k2)''       | ''k1'',''k2''：自由度                 |

これらの関数の引数に共通する ''n'' は出力する乱数の個数を指定するもので、一度に大量の乱数を発生させたい場合に便利である。

===== 一様乱数 =====

一様分布の乱数を__**一様乱数**__という。

$0<x<1$ の範囲の連続型の一様分布は以下の形をしている。
確率分布は面積が 1 になるように作られており、$0<x<1$ の範囲の実数値 $x$ のいずれかが等しい確率で得られる。

{{ dist_uniform.png?nolink |一様分布}} 

$0<x<1$ の一様乱数を 10 個発生させるには以下のようにする。
<code rsplus>
> runif(10)
 [1] 0.657080245 0.882882818 0.322379467 0.284648837 0.943793671 0.008148672
 [7] 0.762963097 0.677352903 0.925840073 0.281425054
</code>

次のプログラム uniform_dist.R では、一様乱数を 1000 個発生させて相対度数のヒストグラムを描く。

<file rsplus uniform_dist.R[enable_line_numbers="true"]>
# 一様乱数

# 乱数の個数
n <- 1000

# 一様乱数の発生
r <- runif(n)
print(r)

# ヒストグラム（相対度数）
hist(r, freq = FALSE, main = "一様乱数")
</file>

また、$1<x<7$ の一様乱数を 10 個発生させるには以下のようにする。
<code rsplus>
> runif(10, 1, 7)
 [1] 4.580631 1.360485 5.536313 3.697583 6.004030 6.948436 5.668975 5.357032
 [9] 6.656362 4.427038
</code>

整数の一様乱数が欲しければ、''as.integer()'' と組み合わせて小数点以下を切り捨てるとよい。

例えば、1 から 6 までの整数の一様乱数（つまり、サイコロの目）を 10 個発生させるには以下のようにする。
<code rsplus>
> as.integer(runif(10, 1, 7))
 [1] 4 1 5 3 6 6 5 5 6 4
</code>

^  runif(10,1,7) の範囲  ^  as.integer() の結果  ^
|  $1<x<2$               |  $1$                  |
|  $2\le x<3$            |  $2$                  |
|  $3\le x<4$            |  $3$                  |
|  $4\le x<5$            |  $4$                  |
|  $5\le x<6$            |  $5$                  |
|  $6\le x<7$            |  $6$                  |

===== 小テスト① =====
[[KMS>|Moodle Server（非公式）]]で第5回の小テスト①を受験しなさい。


===== 乱数の種 =====
疑似乱数には__**再現性**__がある。

例えば単純に ''runif(5)'' を繰り返し回実行すると、その度に異なる乱数列が得られる。
<code rsplus>
> runif(5)
[1] 0.2675082 0.2186453 0.5167968 0.2689506 0.1811683
> runif(5)
[1] 0.5185761 0.5627829 0.1291569 0.2563676 0.7179353
> runif(5)
[1] 0.9614099 0.1001408 0.7632227 0.9479664 0.8186347
</code>
__**乱数の種（seed）**__を ''set.seed(s)'' で設定すると、s の値に応じて毎回同じ乱数列が再現される。
<code rsplus>
> set.seed(10)
> runif(5)
[1] 0.50747820 0.30676851 0.42690767 0.69310208 0.08513597
> set.seed(10)
> runif(5)
[1] 0.50747820 0.30676851 0.42690767 0.69310208 0.08513597
</code>
<code rsplus>
> set.seed(999)
> runif(5)
[1] 0.38907138 0.58306072 0.09466569 0.85263123 0.78674676
> set.seed(999)
> runif(5)
[1] 0.38907138 0.58306072 0.09466569 0.85263123 0.78674676
</code>

理想的な乱数を用いた計算では、毎回異なる乱数なので計算を再現することはできないが、
疑似乱数を使った計算では、乱数の種を指定することで計算を再現できる。

===== 一様乱数の発生アルゴリズム =====
R では一様乱数の発生アルゴリズムとしてメルセンヌ・ツイスターが使用されている。

  * [[http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/mt.html|Mersenne Twister]]
  * [[http://www.math.sci.hiroshima-u.ac.jp/m-mat/TEACH/ichimura-sho-koen.pdf|あなたの使っている乱数、大丈夫？]]

メルセンヌ・ツイスターは松本眞氏と西村拓士氏によって考案されたアルゴリズムで、周期が $2^{19937}−1$ （$\simeq4.3\times 10^{6001}$）と非常に長く、さらに高速に疑似乱数を発生させられるというメリットがある。

一様乱数の発生アルゴリズムは''RNGkind()'' で変更できる。
どのようなアルゴリズムがあるかは ''%%help("RNGkind")%%'' で確認できる。

===== 正規乱数 =====

正規分布の乱数を__**正規乱数**__という。

標準正規分布 $\mathrm{N}(0,1)$ は以下の形をしている。
確率分布は面積が 1 になるように作られており、$x=0$ を中心に左右対称で $x=0$ から離れるにつれて確率が急激に小さくなる。

{{ dist_normal.png?nolink |正規分布}} 


$\mathrm{N}(0,1)$ に従う正規乱数を 10 個発生させるには以下のようにする。
<code rsplus>
> rnorm(10)
 [1] -0.46936885 -0.60760153 -0.62367666 -0.56246919  0.05175988 -1.21122009
 [7] -0.91658480  0.39809603  0.46113233 -0.32689456
</code>

次のプログラム normal_dist.R では、$\mathrm{N}(0,1)$ に従う正規乱数を 1000 個発生させて相対度数のヒストグラムを描く。
プログラムの最後では $\mathrm{N}(0,1)$ の確率密度関数 ''dnorm(x)'' を描いてヒストグラムと比較している。

<file rsplus normal_dist.R[enable_line_numbers="true"]>
# 正規乱数

# 乱数の個数
n <- 1000

# 正規乱数の発生
r <- rnorm(n)
print(r)

# ヒストグラム（相対度数）
hist(r, freq = FALSE, main = "正規乱数")

# N(0,1)の曲線との比較
curve(dnorm, add = TRUE, col = "red")
</file>


===== 中心極限定理による正規乱数 =====

正規乱数は ''rnorm()'' で簡単に発生させられるが、その仕組みはどうなっているのだろうか。
正規乱数を発生させる仕組みはいくつかあるが、一様乱数 ''runif()'' から__**中心極限定理**__を利用する方法がある。

__**中心極限定理**__とは以下のようなものである。

> 期待値 $\mu$、分散 $\sigma^{2}$ の、ある確率分布に従う $n$ 個の独立な確率変数 $X_{1}$, $X_{2}$, $\cdots$, $X_{n}$ があるとする。
> これらの平均 $\bar{X}$ の確率分布は $n$ が大きくなるにつれ正規分布 $\mathrm{N}(\mu, \sigma^{2}/n)$ に近づく。

例えば、連続型一様分布
\[
f(x)=1\quad (0<x<1)
\]
に従う $n$ 個の確率変数を $X_{i}$ ($i=1,2,\cdots,n$) とする。

$X_{i}$ の期待値$\mu$ と分散 $\sigma^{2}$ は以下のようになる。
\[
\mu=E(X_{i})=\int_{-\infty}^{\infty}xf(x)dx=\int_{0}^{1}xdx=\left[\frac{1}{2}x^{2}\right]_{0}^{1}=\frac{1}{2}
\]
\[
\sigma^{2}=E((X_{i}-\mu)^{2})=\int_{-\infty}^{\infty}(x-\mu)^{2}f(x)dx=\int_{0}^{1}(x-\mu)^{2}dx=\left[\frac{1}{3}(x-\mu)^{3}\right]_{0}^{1}=\frac{1}{12}
\]
中心極限定理によれば $X_{i}$ の平均
\[
\bar{X}=\frac{1}{n}\sum_{i=1}^{n}X_{i}
\]
の確率分布は $n$ が大きくなるにつれ $\mathrm{N}(1/2,1/12n)$ に近づくはずである。

この $\bar{X}$ を__**標準化**__によって次のように変数変換する。
\[
Z=\frac{\bar{X}-\mu}{\sqrt{\sigma^{2}/n}}=\frac{\bar{X}-1/2}{\sqrt{1/12n}}=\sqrt{12n}\left(\bar{X}-\frac{1}{2}\right)=\sqrt{\frac{12}{n}}\left(\sum_{i=1}^{n}X_{i}-\frac{n}{2}\right)
\]
標準化とは、期待値$\mu=0$ と分散 $\sigma^{2}=1$ の分布に変換することだから、$Z$ の確率分布は $n$ が大きくなるにつれ $\mathrm{N}(0,1)$ に近づく。


以上より、$Z$ を使うと $n$ 個の一様乱数 $X_{i}$ から $\mathrm{N}(0,1)$ に従う（近似的な）正規乱数を発生させることができる。
特に $n=12$ とすると、平方根を計算する必要がなくなるので都合がよい。
\[
Z=\sum_{i=1}^{12}X_{i}-6
\]

次のプログラム clt.R で一様乱数 ''runif()'' から正規乱数を発生できるのかを確認する。
<file rsplus clt.R[enable_line_numbers="true",highlight_lines_extra="7"]>
# 中心極限定理による正規乱数

# 乱数の個数
n <- 1000

# 正規乱数の発生
r <- replicate(n, sum(runif(12)) - 6)
print(r)

# ヒストグラム（相対度数）
hist(r, freq = FALSE, main = "正規乱数")

# N(0,1)の曲線との比較
curve(dnorm, add = TRUE, col = "red")
</file>

normal_dist.R とは 7行目が異なるだけである。

''sum(runif(12)) - 6'' では
\[
Z=\sum_{i=1}^{12}X_{i}-6
\]
を計算している。

''replicate(n, expr)'' は ''expr'' を ''n'' 回くり返し、各回の結果をまとめて ''n'' 個の数列のベクトルとして返す関数である。

===== 棄却法による乱数 =====
任意の確率密度関数 $f(x)$ に従う乱数を発生させる方法として__**棄却法**__がある。

$x$ の区間 $[a, b]$ で確率密度関数 $f(x)$ があるとする。
この区間で $f(x)\le c$ となる定数 $c$ が分かっているとする。

以下の処理を行うと $f(x)$ に従う乱数が得られる。

  - $a<x<b$ の一様乱数 $x$ を発生させる。
  - $0<y<c$ の一様乱数 $y$ を発生させる。
  - $y\le f(x)$ であれば $x$ を乱数として採択し、そうでなければ棄却する。
  - 1 に戻って繰り返す。

下図のように $y=f(x)$ のグラフで考えると、横が $a<x<b$、縦が $0$ から $0<y<c$ の長方形の領域内でランダムに打った点が $y=f(x)$ と $x$ 軸で囲まれた図形の中に入れば採択し、入らなければ棄却するということである。

$c$ はなるべく小さい方が棄却される点が少なくなり効率がよいので、$c$ は $f(x)$ の最大値にするのが望ましい。

{{ rejection.png?nolink |棄却法}} 

----
==== 例 ====
棄却法によって正規乱数を発生させる。

次のプログラム rejection.R は棄却法で正規乱数を発生させ、標準正規分布 $\mathrm{N}(0,1)$ と比較することで本当に正規乱数になっているのかを確認する。
確率密度関数 $f(x)$ は標準正規分布 ''dnorm()'' を用いている。

<file rsplus rejection.R[enable_line_numbers="true"]>
# 棄却法

# 乱数の個数
n <- 1000

# 確率密度関数 f(x)
f <- function(x) {
  return(dnorm(x))
}

# 棄却法
rejection <- function(a, b, c) {
  while (1) {
    x <- runif(1, a, b)
    y <- runif(1, 0, c)
    if (y <= f(x)) {return(x)}
  }
}

r <- replicate(n, rejection(-10, 10, 1))
print(r)

# ヒストグラム（相対度数）
hist(r, freq = FALSE, main = "棄却法による乱数")

# N(0,1)の曲線との比較
curve(dnorm, add = TRUE, col = "red")
</file>

関数 ''rejection(a,b,c)'' が棄却法で乱数を発生させているところである。

このプログラムでは $x$ の範囲を $[-10,10]$、$c=1$ としている。
この範囲は $\mathrm{N}(0,1)$ に比べて広すぎるので棄却される数が多くなり乱数の発生に時間がかかる。
そこで $c$ を $\mathrm{N}(0,1)$ の最大値 ''dnorm(0)'' にすると無駄を減らせる。

{{ rejection_normal.png?nolink |標準正規分布の範囲}}

rejection.R の確率密度関数 $f(x)$ を変えれば他の分布に従う乱数を発生させることができる。


===== モンテカルロ法 =====
乱数を用いて数値計算を行う方法を__**モンテカルロ法**__という。
[[https://www.monte-carlo.mc/en/|モンテカルロ]]というのはモナコ公国のカジノで有名な地名で、そこからこの方法の名前が付けられた。

{{ montecarlo.png?nolink |モンテカルロ法}}

ある図形の面積 $S$ を求めよう。

まず、この図形を囲む長方形を描く。
この長方形の各辺は $x$ 方向が $[x_{1}, x_{2}]$、$y$ 方向が $[y_{1}, y_{2}]$ とする。

そして以下の処理を行う。

  - $x_{1}<x<x_{2}$ の一様乱数 $x$ を発生させる。
  - $y_{1}<y<y_{2}$ の一様乱数 $y$ を発生させる。
  - このランダムな $(x,y)$ を $xy$ 平面上の座標として点をプロットする。
  - 点が図形の内部に入れば採択してカウントする。入らなければ棄却してカウントしない。
  - $n$ 個の点がプロットされるまで 1 に戻って繰り返す。

これは棄却法の処理と同じである。

図形の内部に点が入る確率は、図形の面積 $S$ と長方形の面積 $(x_{2}-x_{1})\times (y_{2}-y_{1})$ の比である。
図形の面積が大きいほど点は図形の内部に入りやすい。
図形の内部に入った点の数をカウントして $m$ とすると、図形の面積 $S$ は次のように近似できる。

\[
S\simeq \frac{m}{n}(x_{2}-x_{1})(y_{2}-y_{1})
\]

モンテカルロ法では、$n$ が大きいほど理論値に近い結果が出やすくなる。
理論値からのばらつき（標準偏差）は $1/\sqrt{n}$ に比例して小さくなっていく。

----
==== 例1 ====
半径 $r=1$ の円の面積 $S$ をモンテカルロ法で求める。
正確な値は円の面積の公式より $S=\pi r^{2}=\pi$ である。

{{ montecarlo_circle.png?nolink |半径 1 の円の図形}}

原点を中心とした半径 1 の円が入る $x_{1}=-1$, $x_{2}=1$, $y_{1}=-1$, $y_{2}=1$ の正方形を考える。

原点を中心とした半径 $r$ の円は次の式で表される。
\[
x^{2}+y^{2}=r^{2}
\]
この円の中の領域は不等式で表すと
\[
x^{2}+y^{2}<r^{2}
\]
である。
したがって、ランダムな点 $(x,y)$ が半径 $r=1$ の円の内部に入っているかは
\[
x^{2}+y^{2}<1
\]
で判断できる。

{{ montecarlo_circle_anim.gif?nolink |ランダムな点のプロット（赤が円内、青が円外）}}

次のプログラム montecarlo.R はモンテカルロ法で半径 1 の円の面積、つまり $\pi$ の近似値を求めるプログラムである。

<file rsplus montecarlo.R[enable_line_numbers="true"]>
# モンテカルロ法

# 判別関数
f <- function(x, y) {
  if (x ^ 2 + y ^ 2 < 1) return(1)
  return(0)
}

# 点の個数
n <- 1000

# 図形を囲む長方形の座標
x1 <- -1
x2 <- 1
y1 <- -1
y2 <- 1

# n個の一様乱数
x <- runif(n, x1, x2)
y <- runif(n, y1, y2)

m <- 0
for (i in 1 : n) {
  # 図形内部の点を判別関数でカウント
  m <- m + f(x[i], y[i])
}

# 面積
s <- m / n * (x2 - x1) * (y2 - y1)
print(s)
</file>

ランダムな点が図形の内部に入っているかどうかは判別関数で判断してカウントしている。
**<fc #ff0000>判別関数は、面積を計算したい図形の内部の点だと 1、外部の点だと 0 を返す関数として定義している。</fc>**

実行の度に結果が異なり、ばらつくことが分かる。

<hidden 短縮形>
このプログラムは次のようにかなり短くすることができる（ただし、分かりにくい）。
<code rsplus>
> n <- 1000
> sum(runif(n, -1, 1) ^ 2 + runif(n, -1, 1) ^ 2 < 1) / n * 4
</code>
</hidden>
\\

次のプログラム montecarlo_hist.R では、結果のばらつきがどのくらいかを相対度数のヒストグラムで確認できる。

<file rsplus montecarlo_hist.R[enable_line_numbers="true"]>
# モンテカルロ法のヒストグラム

# 試行回数
trial <- 100

# 点の個数
n <- 1000

# 面積をn回求める
r <- replicate(n, sum(runif(n, -1, 1) ^ 2 + runif(n, -1, 1) ^ 2 < 1) / n * 4)
print(r)

# 面積のヒストグラム（相対度数）
hist(r, freq = FALSE)
</file>

面積のばらつきを見るために、円の面積の計算を 100 回くり返している。

ここでは円で行ったが、円の対称性を考慮して四分円で行うこともできる。

----
==== 例2 ====

定積分
\[
S=\int_{-1}^{1}(-3x^{2}+3)dx
\]
をモンテカルロ法で求める。
正確な値は $S=4$ である。

{{ montecarlo_quad.png?nolink |定積分の図形}}

$y=-3x^{2}+3$ と $x$ 軸で囲まれた図形が入る $x_{1}=-1$, $x_{2}=1$, $y_{1}=0$, $y_{2}=3$ の長方形を考える。
この図形の内部は不等式では以下のように表される。
\[
y < -3x^{2}+3\ かつ\ y > 0
\]
これには以下のように montecarlo.R を書き換えればよい。

=== 判別関数 ===

4行目から7行目
<code>
f <- function(x, y) {
  if (y < -3 * x ^ 2 + 3 && y > 0) return(1)
  return(0)
}
</code>

=== 長方形の座標 ===
15,16行目
<code>
y1 <- 0
y2 <- 3
</code>

===== 小テスト② =====
[[KMS>|Moodle Server（非公式）]]で第5回の小テスト②を受験しなさい。


===== 課題 =====
[[KMS>|Moodle Server（非公式）]] で第5回の課題を行いなさい。

締め切り：2025年5月14日（水）20時