この記事は

TeX ＆ LaTeX Advent Calendar 2016 - Adventar

の13日目の記事です．

はじめに

LaTeXでモンテカルロ積分を行うプログラムを作りました．Tikzで可視化もしました．

何を作ったのか

モンテカルロ積分によって円周率の近似値を求めるプログラムです．左にモンテカルロ積分を可視化した図，右にプロットした点の数や求まった円周率の近似解を表示するプログラムです．

モンテカルロ 積分について

モンテカルロ 積分

モンテカルロ積分とはモンテカルロ法によって定積分の近似値を求めるアルゴリズムです．モンテカルロ法は疑似乱数を用いてシミュレーションや数値計算を行う手法の総称です．つまりモンテカルロ積分は疑似乱数を用いて定積分の近似値を求めるアルゴリズムです．

モンテカルロ 積分で円周率を求めるには

積分値に ${\pi}$ が含まれるような関数の積分値の近似値を求めます．今回は以下の定積分を解いて円周率を解きます．

${ \displaystyle \int_0^1 \sqrt{1 - x^2} }$

この定積分の解は

${ \displaystyle \int_0^1 \sqrt{1 - x^2} = \frac{\pi}{4} }$

ですので定積分の近似解を4倍すれば円周率の近似値が得られそうです．

${ \displaystyle \int_0^1 4\sqrt{1 - x^2} = \pi }$

モンテカルロ 積分の手順

今回の関数は ${ \displaystyle \sqrt{1 - x^2} }$ の積分を0から1まで求めるので次のようなイメージになります． f:id:muscle_keisuke:20161205153744p:plain

この内部の面積を ${S}$ とします．

この関数に囲まれた部分の面積 ${S}$ を求めます．次にこの正方形内部にランダムに点を打ちます．次の図はプロット点数1000で描画した時のイメージです．青の点が内部の点，赤の点が外部の点です． f:id:muscle_keisuke:20161205155637p:plain

全体の点数と青の点数の比と正方形の面積と求めたい面積の比が近似することを利用して比の方程式で ${S}$ の近似値を求めます．
${全体の点数 : 青の点数 \approx 正方形の面積 : S}$
${1000 : 青の点数 \approx 1 : S}$
${S \approx \frac{青の点数}{1000}}$
${S}$ の本当の値は ${\frac{\pi}{4}}$ なので円周率の近似値は
${\pi \approx 4\frac{青の点数}{1000}}$ となります．

TeX+TikZでモンテカルロ 積分の実装

変数やマクロを定義する

積分範囲はTeXのマクロに定義し，赤青の点数などカウントが必要なものはLaTeXの変数に代入してしまいます．

\def\startpoint{0}
\def\endpoint{1}

\newcounter{redplot}
\newcounter{blueplot}
\setcounter{redplot}{0}
\setcounter{blueplot}{0}

プロット点数は標準入力で

プロット点数は入力できるようにしました．

\newcount\input

\message{プロット点数を入力してください}
\read-1 to \input

正方形や関数を描く

Tikzで正方形や関数を描画します．

% 枠を描画
\draw (\startpoint,\startpoint) rectangle (\endpoint,\endpoint);
% 関数を描画
\draw plot[domain=\startpoint:\endpoint,samples=100] (\x,{func(\x)});

これで準備は整いました．後は点をプロットして円周率の近似値を求めます．

fpパッケージを使いTeXで固定小数点演算を行う

0から1の積分なのでランダムにプロットする点は必ず小数になります．しかし，TeXは標準で四則演算くらいはできますが小数を扱うことができません．なので fpパッケージと呼ばれるTeXで固定小数演算を行うパッケージを使います．私の実行環境であるTeXLive2016では標準で入っていました． fpは小数演算の四則演算だけでなく指数計算や対数計算など様々な計算をすることが可能です．詳しくは以下の方の記事が分かりやすいです．

天地有情 [LaTeX] fp --- 固定小数点数演算

ランダムに点をプロットする

用意できたら点をランダムにプロットしていきます． fpパッケージには乱数を発生するマクロもあります．まずはプロットする点をランダムに決めていきます．

% 今日，0時から何分経過したかでシードを決定
\FPseed = \time
\begin{document}
.
.
.
% ランダムにx座標を決める
\FPrandom{\randomX}
% f(randomX)を求める
\FPpow{\randomFuncX}{\randomX}{2} % randomX^2
\FPsub{\randomFuncX}{1}{\randomFuncX} % 1 - randomX^2
\FProot{\randomFuncX}{\randomFuncX}{2} % √1 - randomX^2
% randomXをfp以外で評価できる形にする
\FPeval\plotX{\randomX}
% f(randomX)をfp以外で評価できる形にする
\FPeval\evalX{\randomFuncX}
% ランダムにy座標を決める
\FPrandom{\randomY}
% randomYをfp以外で評価できる形にする
\FPeval\plotY{\randomY}

本当は以下のようにも書けるのですがオーバーフローを引き起こすので上のような遠回りをしています．

\FPeval\plotX{random()}
\FPeval\evalX{root(1 - \plotX^2,2)}
\FPeval\plotY{random()}

点が関数の内側かどうかを判定する

ある点(X, Y)が関数の内側にあるかどうかは以下の不等式で判定可能です．
${f(X) > Y}$
これをTikZの条件式で判定します．

% 関数の求める範囲の外側の場合は赤，内側の場合は青で点をプロット
\pgfmathparse{\evalX < \plotY ? "red" : "blue"}
% blueかredを\colorに定義
\edef\color{\pgfmathresult}

次に定義した\colorを元に青の点と赤の点をカウントします．

% 赤の点と青の点を数える
\ifthenelse{\equal{\color}{red}}{\stepcounter{redplot}}{\stepcounter{blueplot}}

点をプロット

Tikzのコマンドで点をプロットします．

% 点をプロットする
\filldraw[fill=\color, ultra thin] (\plotX, \plotY) circle [radius=0.01];

入力分ループを回す

それではここまでの処理をTikZのforechで回していきます．

\foreach \i in {1,...,\input}
{
      % ランダムにx座標を決める
      \FPrandom{\randomX}
      \FPpow{\randomFuncX}{\randomX}{2}
      \FPsub{\randomFuncX}{1}{\randomFuncX}
      \FProot{\randomFuncX}{\randomFuncX}{2}
      \FPeval\plotX{\randomX}
      \FPeval\evalX{\randomFuncX}
      % ランダムにy座標を決める
      \FPrandom{\randomY}
      \FPeval\plotY{\randomY}
      % 関数の求める範囲の外側の場合は赤，内側の場合は青で点をプロット
      \pgfmathparse{\plotY > \evalX ? "blue" : "red"}
      \edef\color{\pgfmathresult}
      % 赤の点と青の点を数える
      \ifthenelse{\equal{\color}{red}}{\stepcounter{redplot}}{\stepcounter{blueplot}}
      % 点をプロットする
      \filldraw[fill=\color, ultra thin] (\plotX, \plotY) circle [radius=0.01];
}

コンパイル

それではコンパイルしていきます．私の環境ではplatex+dvipdfmxでコンパイルしています．コンパイルしたら「プロット点数を入力してください」と表示されるので 1000とか入力します．そうするとコンパイルされてdviが出てくるのでそれもpdfに変換して表示してみると以下のような図が出ると思います．

f:id:muscle_keisuke:20161205155637p:plain

近似値なども表示する

プロット点数や近似値も表示させます．全体の面積を求めます．

% 枠内全体の面積を変数に代入
% 四捨五入して整数にしている
\FPeval\overallS{round((\endpoint - \startpoint)^2,0)}

比の方程式を解いて，4倍して円周率の近似値を求めます．小数点以下が多いので第9位まで丸めます．

%\arabic{blueplot} ... 青の点数
%\overallS ... 正方形の面積
%\input ... 全体のプロット点数
\FPeval\result{round(4 * \arabic{blueplot} * \overallS / \input,9)}

さて，それでは全て表示させます．

\begin{tikzpicture}
    \draw (7,5) node {全体の点数 = \input};
    \draw (7,4) node {青い点数 = \arabic{blueplot}};
    \draw (7,3) node {赤い点数 = \arabic{redplot}};
    \draw (7,2) node {全体的な面積 = \overallS};
    \draw (7,1) node {$\displaystyle\int^1_0 4\sqrt{1-x^2}$ = \result};
\end{tikzpicture}