モンテカルロシミュレーションによる円周率πの計算

滋賀大学教育学部紀要
モンテカルロシミュレーションによる円周率πの計算
119
No.65, pp. 119-128, 2015
モンテカルロシミュレーションによる円周率πの計算
水　上　善　博
On Calculations of π by Monte Carlo Simulations
Yoshihiro MIZUKAMI
Abstract
Monte Carlo simulations are performed on personal computers to estimate π . RND function of
VBA(Visual Basic for Application) in Microsoft Excel® is used for pseudorandom number generation.
Mersenne Twister is used for pseudorandom number generation, too. Central limit theorem is
adapted to estimate accuracy of Monte Carlo simulations. The limit of frequency of random number
generation by RND function is discussed. Monte Carlo simulations are performed on various types of
personal computers to check machine time. Personal computers, which are commonly used in class
room, will be available for Monte Carlo simulation of π in reasonable machine time in lessons.
キーワード：円周率、モンテカルロシミュレーション、中心極限定理、計算時間
１．はじめに
ることから、児童・生徒から興味をもってもらえ
高等学校の教科「情報」には、「モデル化と
る対象であると思われる。大学の演習授業におい
シミュレーション」という単元 があり、シミュ
てエクセルを用いてπのモンテカルロシミュレー
レーションの例として、乱数を用いたモンテカ
ションを実践した試みでは 1)、学生たちがシミュ
ルロ法がしばしば取り上げられている。乱数を
レーションに興味を抱いたとの報告がなされてい
利用したシミュレーションは、パーソナルコン
る。
ピュータ ( パソコン ) に標準的にインストール
本研究では、高等学校の教科「情報」や中学
されているアプリケーション、例えば、エクセ
校の技術科において、円周率πの値をシミュ
ル（Microsoft Excel®）のような表計算ソフト
レーションで計算する授業を実施する際に必
でも実行できるので、比較的手軽に授業で実施
要となる、精度や計算時間などの基礎的な事
できる。しかし、パソコンの性能によってシミュ
項を検証することを目的とする。具体的には、
レーションに必要な時間が異なり、古いパソコ
エクセルの乱数を利用して、円周率πの値を
ンでは、計算時間が長くなって授業で行うには
求めるモンテカルロシミュレーションのプロ
適さないこともある。
グ ラ ム を エ ク セ ル の VBA（Visual Basic for
また、シミュレーションを実行する対象の選
Applications）で作成し、シミュレーションの
択は重要である。円周率（3.14）は高等学校の生
回数とπの計算結果の精度について中心極限定
徒のみならず、小学校の高学年の児童や中学校の
理などの統計数学の観点から検証し、適切なシ
生徒にとっても比較的身近な定数であり、シミュ
ミュレーションの回数を探った。さらに、様々
レーションの結果得られる円周率πの値を正確
な性能のパソコンでシミュレーションを実施
なπの値（3.14159…）と比較することによって、
し、結果が得られるまでの計算時間を比較した。
どのくらいの規模のシミュレーションでどのくら
これらを基に、現在、教育現場で一般的に使用
いの精度の値が得られるかを実感することができ
されているパソコンの性能でどのくらいのシ
計算する授業を実施する際に必要となる、精度や計算時間などの基礎的な
ことを目的とする。具体的には、エクセルの乱数を利用して、円周率 π の
120
水　上　善　博
テカルロシミュレーションのプログラムをエクセルの
VBA
（Visual Basic for
作成し、シミュレーションの回数と π の計算結果の精度について中心極限
数学の観点から検証し、適切なシミュレーションの回数を探った。さらに、
ミュレーションが可能であるかを考察した。
とπの近似値が求まる。　ソコンでシミュレーションを実施し、結果が得られるまでの計算時間を比
を基に、現在、教育現場で一般的に使用されているパソコンの性能でどの
２．方法
VBA のプログラミング
レーションが可能であるかを考察した。
モンテカルロ法による円周率πの計算
具体的なシミュレーションはエクセルの
円周率πの計算の歴史については、ベックマ
VBA の乱数を用いたコンピュータシミュレー
2)
ンの著作「πの歴史」 に詳しい。その中の第 15
ションで行う。エクセルのマクロで VBA のプ
章「モンテ・カルロ法」では、シミュレーショ
ログラムを作成したが、エクセルのマクロを実
による円周率 π の計算
ンによってπを求める方法が紹介されている。い
行するには、メニューに「開発タブ」を表示さ
算の歴史については、ベックマンの著作「π の歴史」2)に詳しい。その中の
わゆる「ビュフォンの針」というもので、
「長さ
せる必要がある。Excel 2010 以降のバージョン
テ・カルロ法」では、シミュレーションによって
π
を求める方法が紹介され
L の針を、L より大きい間隔 d で平行な多数の直
では、「オプション　→　リボンのユーザー設
る「ビュフォンの針」というもので、
「長さ L の針を、L より大きい間隔
d
線が引かれている床の上に、でたらめに落下させ
定　→　開発　にチェックを入れる
→　OK」
直線が引かれている床の上に、でたらめに落下させる。このとき、針が直
る。このとき、針が直線のどれかと交わる確率を
で、開発タブを表示できる。マクロで VBA の
プログラム（例えば、pai という名前の関数）
とすると、円周率
わる確率を P P
とすると、
円周率   2 L / Pd と求まる」とい
と求まる」というものである。
うものである。本で紹介されている実験結果では
を作成するには、「開発　→　Visual Basic　→　いる実験結果では
12000 回のシミュレーションで π として 3.14 程度の値が
12000 回のシミュレーションでπとして 3.14 程度
挿入　→　標準モジュール」として、図 2 のよ
詳しくは「π の歴史」2)を参照のこと。
の値が得られている。詳しくは「πの歴史」2) を
うなプログラムを入力する。本来ならば、変数
参照のこと。
のデータ型の指定を厳密に行うべきであるが、
中学生や高校生にわかりやすいように、巨大な
整数（例えば、百万：106 ）を代入する num と
x
いう変数のみ Long（長整数型）を指定している。
また、このプログラムでは、繰り返し処理を実
行する、For ～ Next ステートメントや、条件
分岐を実行する、If 条件式 Then ステートメン
トなどの基本的な文法が用いられており、高校
生のみならず、中学校 3 年生程度がプログラミ
y
ングの基本を理解する教材としても適している
と思われる。
プログラムの入力が終わったら、「ファイル　→　終了して Microsoft Excel へ戻る」でマク
ロを閉じて、エクセルのシートに戻って計算を
図 1　円周率πのシミュレーションに用いられたモン
実行する。シミュレーションの実行は、例えば、
テカルロ法
図 1 円周率 π のシミュレーションに用いられたモンテカルロ法
100 回ランダムな点を発生してπを計算する場
本研究では、コンピュータシミュレーション
合、セルに　= pai(100)　と入力する。あるい
でより実行しやすい方法を用いる。図 1 に示す
は、セル A1 に 100 と入力して、セル B1 に　=
ように、一辺が長さ 1 の正方形の内側に内接す
pai(A1)　と入力してもよい。すなわち、エクセ
る四分の一の円を描く。正方形内でランダムな
ルで標準に用意されている関数と同様に、関数
本研究では、コンピュータシミュレーションでより実行しやすい方法を用いる。図 1 に
点を N 個発生させて、そのうち四分の一の円
pai が利用できるということである。もし、う
示すように、一辺が長さ
1
の正方形の内側に内接する四分の一の円を描く。正方形内でラ
の中に入る点の数 S を数える。ランダムな点を
まく実行できない場合は、プログラムの入力に
多数発生させれば、N
は正方形の面積 1 に対応
ミスがないかどうかを確認するとともに、マク
ンダムな点を N 個発生させて、そのうち四分の一の円の中に入る点の数
S を数える。ラン
ロが有効になっているかどうかをチェックする
し、S
は半径
1
の単位円の四分の一の面積
π
/
4
ダムな点を多数発生させれば、N は正方形の面積 1 に対応し、S は半径 1 の単位円の四分の
とよい。
に対応する。よって、 S / N ≅ (π / 4) / 1 = π / 4
一の面積  / 4 に対応する。よって、 S / N  ( / 4) / 1   / 4 となり、
となり、
注意が必要なのは作業が終わってファイルを
保存するときに、通常のエクセルファイルとし
（１）
てではなく、「ファイルの種類　→　Excel マ
  4S / N
（１）
モンテカルロシミュレーションによる円周率πの計算
121
クロ有効ブック　→　保存」のようにマクロ有
題があることが指摘されている 3)。まず、単精
効ブックとして保存することである。
度なので 10 進数表現で 7 ケタまでしか精度が
なお、本研究では、計算開始と計算終了の時刻
保障されない。また、乱数を生成する周期が
を表示するマクロを作成してシミュレーション
1677 万（周期は 10 進法で 7 桁程度）なので、
にかかる計算時間を算出した。
107 回～ 108 回以上の乱数を発生させると同じ
数が出現して、乱数としての意味がなくなると
乱数について
思われる。すなわち 106 回程度の乱数を発生さ
エクセルの VBA では乱数（疑似乱数）を生
せるシミュレーションが限界であると予想され
成 す る の に Visual Basic の RND 関 数 を 用 い
る。また、65536 種類の系列の乱数しか得られ
るが、RND 関数が生成する乱数には様々な問
ないので、かなりの確率でまったく同じ系列の
図 2　VBA で作成した、モンテカルロ法で円周率πを計算する関数のプログラムの一例
図 2 VBA で作成した、モンテカルロ法で円周率 π を計算する関数のプログラムの一例
表 1　乱数の発生回数が異なるモンテカルロシミュレーションによる π の値。上は
表１　乱数の発生回数が異なるモンテカルロシミュレーションによる
π の値。上は RND 関数で乱数
RND 関数で乱数を生成したシミュレーションの結果。下はメルセンヌ・ツイスタ
を生成したシミュレーションの結果。下はメルセンヌ・ツイスタで乱数を生成したシミュレーション
で乱数を生成したシミュレーションの結果。平均値、標準偏差、最大値、最小値
の結果。平均値、標準偏差、最大値、最小値は100回実行して求めた統計量。
は 100 回実行して求めた統計量。
誤差 εは、ε
μ－π ) /はエクセルの関数
π ただし、π pi()
はエクセルの関数
pi() で得られる円
誤差 εは、ε
= (μ－π) /=π ( ただし、π
で得られる円周率の値
周率の値 3.1415926535897。
3.1415926535897。
RND
乱数の発生回数（N ） 平均値（μ ） 標準偏差（σ ） 最大値
最小値
2
3.1316 0.148947776
3.48
2.72
10
3.12788 0.052139674
3.24
2.98
103
3.140552 0.014351003
3.1704
3.1048
104
3.1415944 0.005419461
3.15452
3.12572
105
6
3.1415902
0.001621922
3.14514
3.13792
10
7
3.141582696
0.000190778
3.1419424
3.1412132
10
3.141587 0.000014490 3.14162292 3.14155796
108
誤差（ε ） ％
-0.318076043
-0.436487320
-0.033125033
0.000055590
-0.000078100
-0.000316960
-0.000179959
メルセンヌ・ツイスタ
乱数の発生回数（N ） 平均値（μ ）
3.1636
102
3.14108
103
3.141376
104
3.1409432
105
3.14144408
106
7
3.14162714
10
8
3.141596535
10
誤差（ε ） ％
0.700515593
-0.016318271
-0.006896298
-0.020672750
-0.004729244
0.001097737
0.000123555
標準偏差（σ ） 最大値
最小値
0.148711264
3.44
2.84
0.051851650
3.288
2.984
0.014010318
3.1688
3.0976
0.004876898
3.154
3.1284
0.001561359
3.145032
3.137992
0.000512842 3.1427024
3.140452
0.000157429 3.14203112 3.14122144
122
水　上　善　博
乱数が生成される危険性がある。そこで、比較
値に近づいて 105 から 108 では有効数字が小数
のため、現在、最も信頼性が高い疑似乱数を生
点以下 4 桁まで真値に一致する。一方、メルセ
成すると言われているメルセンヌ・ツイスタ法
ンヌ・ツイスタでは発生回数が 102 のシミュレー
（周期は 10 進法で 6 千桁以上）によって生成し
ションではあまり良いπの値は得られないが、
た乱数によるモンテカルロシミュレーションも
103 からは真値に近づいて 104 から 106 では有効
あわせて実施した。なお、エクセルでのメルセ
数字が小数点以下 2 桁、107 では 3 桁、108 では
ンヌ・ツイスタ法による乱数生成は、和田維作
4 桁まで真値に一致する。
4)
氏によって作成された、DLL をエクセルにア
授業中に児童・生徒が 1 回だけシミュレー
ドインして実行した。
ションをしたときに得られるπの値の目安とし
て最大値、最小値を見ると、表 1 より、RND
３．結果と考察
の場合は、乱数発生が 106 のシミュレーション
シミュレーションの結果
の 場 合 3.14 ～ 3.15、107 で 3.141 ～ 3.142、108
2
表 1 に 10 から 10 ま で 10 倍 ず つ 乱 数 の 発
で 3.1416 程度となっている。一方、メルセンヌ・
生回数（N）を変えてモンテカルロシミュレー
ツイスタでは、106 で 3.14 ～ 3.15 程度、107 で 3.14
ション実行してπの値を求めた結果を示す。表
程度、108 でも 3.14 程度となっている。
の上部分が RND 関数によって乱数を生成した
文献 5) では、本研究と同様のシミュレーショ
シミュレーションの結果、表の下部分がメルセ
ンによるπの計算結果が示されており（乱数は
ンヌ・ツイスタによって乱数を生成したシミュ
C++ の RAND 関数で生成）、11111 回の乱数発
レーションの結果である。πの平均値、標準偏
生で誤差 0.44％、111111 回の乱数発生で誤差
差、最大値および最小値は 100 回のシミュレー
0.01% であったと報告されている。これを本研
ションから求めた統計量である。log10N を横軸
究の結果と比較する。表１で 103 の乱数を発生
に、πの平均値を縦軸に描いたグラフを図３に
させたシミュレーションを 100 回実行して得ら
示す。エクセルの pi 関数で求めた小数点以下
れたπの平均値はトータルで 105 回乱数を発生
13 桁までのπ（3.1415926535897：これを真値
したシミュレーション 1 回分に相当するので、
と呼ぶ）を合わせて示す。RND の乱数では発
表１の 103 の時の誤差（RND で 0.44%、メルセ
2
8
3
生回数が 10 と 10 のシミュレーションではあ
ンヌ・ツイスタで 0.02%）が、文献 5) の 11111
4
まり良いπの値は得られないが、10 からは真
回の乱数発生の誤差 0.44％と比較できる。また、
3.17
3.165
π の平均値
3.16
3.155
RND
3.15
3.1415926535897
3.145
メルセンヌ・
ツイスタ
3.14
3.135
3.13
3.125
0
1
2
3
4
5
6
7
8
9
log10N
図 3　N 回の乱数発生のモンテカルロシミュレーションを 100 回実行して得られたπの平均値 .　黒丸は RND 関数の乱数を用いたシミュレーション、×はメルセンヌ・ツイスタによる乱
数を用いたシミュレーション .
図 3 N 回の乱数発生のモンテカルロシミュレーションを 100 回実行して得られた π の平均値.
黒丸は RND 関数の乱数を用いたシミュレーション、×はメルセンヌ・ツイスタによる
乱数を用いたシミュレーション.
log10 σ log10 σ
中心極限定理は「母集団の分布がどのようなものであっても（ただし、期待
規分布で近似できる」というものである
。
中心極限定理は「母集団の分布がどのようなものであっても
規分布で近似できる」というものである
中心極限定理は「母集団の分布がどのようなものであっても（た
る。また、表１の
104 の時の誤差（RND で 0.03%、
6)
定義できない分布を除く）
、確率変数の和の確率分布の形は標本の大
規分布で近似できる」というものである
規分布で近似できる」とい
。中心
定義できない分布を除く）
、確率変数の和の確率分布の形は標本の大きさが大
定義できない分布を除く）
、確率変数の和の確率分布の形は標本
定義できない分布を除く）
、確率変数の和の確率分布の形は標本の大
11/ /1 m
ると標準偏差が
m倍になることがわか
倍になることがわか
ると標準偏差が
数を発生したシミュレーション
表
6) 回分に相当するので、
献 5)の 111111 回の乱数発生の誤差
0.01％と比較でき
6)
6)
規分布で近似できる」というものである
。中心極限定理によると
6)
規分布で近似できる」というものである
。中心極限定理によると、標本の数
1 / mると標準偏差が
倍になることがわかって
1 / m 倍に
ると標準偏差が
表１のシミュレーションは乱数の発生回
規分布で近似できる」というものである
。中心極限定理による
表１のシミュレーションは乱数の発生回
規分布で近似できる」というものである
。中心極限定理によると
メルセンヌ・ツイスタで
0.02%）が、文献 5)の 11111 回の
を用
モンテカルロシミュレーションによる円周率πの計算
123
11//これらの結果を総合的に見ると、一見、RND
倍になることがわかっている。
ると標準偏差が
表１のシミュレーション
N
1 / m 倍になることがわかっている。
ると標準偏差が
10
ので、標本の数は
10倍ずつ増えている。よ
倍ずつ増えている。よ
1m
/ 表１のシミュレーションは乱数の発生回数
mので、標本の数は
倍になることがわかっている。
ると標準偏差が
4
m
倍になることがわかっている。
ると標準偏差が
る。また、表１の
10
の時の誤差（RND
で
0.03%、メルセ
8
ンヌ・ツイスタのそれよりも精度がよいか、少なく
2
8 22から
2 1010
8
表１のシミュレーションは乱数の発生回数
N
を
10
10
8まで
ので、標本の数は
10
倍ずつ増えている。よって
ので、標本の数は
倍ずつ
表１のシミュレーションは乱数の発生回数
N
を
10
から
10
まで
倍ずつ
表１のシミュレーションは乱数の発生回数
N 10
を 10
から
ま
表１のシミュレーションは乱数の発生回数
Nを
から
10 10
まで
2 0.01％と比較できる。
献 5)の
111111の標準偏差
回の乱数発生の誤差
る。そこで、より詳しい計算精度のチェックを行う
3 3はは10102 の標準偏差
の標準偏差
の標準偏差2 2のの 11/ / 1
ので、標本の数は
10 倍ずつ増えている。よって、中心極限定理より
ので、標本の数は
10
倍ずつ増えている。よって、中心極限定理より、乱数の発
ので、標本の数は
10
倍ずつ増えている。よって、中心極限定理よ
ので、標本の数は
10
倍ずつ増えている。よって、中心極限定理より
これらの結果を総合的に見ると、一見、RND
の標準偏差  3 は 102 の標準偏差
の標準偏差
は 10
1を用いた
/ 2 の標準
10 に
表１の 104 の時の誤差（RND で 0.03%、メルセ
中心極限定理より、乱数の発生回数が
103 の23の
析を行った。
2
2  は 102 2の標準偏差
4 同様に 1044の
2
の標準偏差

の
1
/
10
になる。
ンヌ・ツイスタのそれよりも精度がよいか、少なくとも
ンヌ・ツイスタで 0.01%）が、文献 5) の 111111
の標準偏差の標準偏差
 3標準偏差
は 10 の標準偏差
 2の標準偏差
の
1 / 10
になる。
同様に
10 の標準偏差

の
になる。これを、式で
準偏差
の標準偏差
は の標準偏差
10

110
/ になる。
10
になる。
同様に
3 23の
1/1// 10
になる。これを、式
準偏差
32 の
33は
は310
10
の標準偏差
の
10
になる。
同様に
10 1の
21の
4
3
回の乱数発生の誤差 0.01％と比較できる。
同様に 10
の標準偏差
る。そこで、より詳しい計算精度のチェックを行うため
の 10
1 / の標準偏差
10 準偏差
になる。これを、式で表す
 33の
の 1 / 10 にな
準偏差  34は
中心極限定理
3 になる。これを、式で表すと、
の
1の
になる。これを、式で表すと、
準偏差
これらの結果を総合的に見ると、一見、RND
/ 10
になる。これを、式で表すと、
準偏差  の
n n11(1(1/ / 10
析を行った。
/中心極限定理は「母集団の分布がどのようなもの
10
になる。これを、式で表すと、
3 1準偏差
10))n n ただし、n
ただし、n==2,2,3,・
3,
1// 110
10
になる。これを、式で表すと、
3 の3（RND
数を発生したシミュレーション 1 回分に相当するので、表１の3 10準偏差
の時の誤差
で
0.44%、
を用いたシミュレーションのほうがメルセン
 n 1  (1 / 10 ) n、確率変数の和の確率分布
ただし、n
 n 1  (1 / = 10
2, 3,・・・
) n た
定義できない分布を除く）
メルセンヌ・ツイスタで 0.02%）が、文献 5)の 11111 回の乱数発生の誤差 0.44％と比較でき
を作用して、整理する
となる。両辺に
log
を作用して、整理する
となる。両辺に
log
1010
6)
ただし、n
=
2,
3,・・・,
8
ヌ・ツイスタのそれよりも精度がよいか、少な


(
1
/
10
)

ただし、n
=
2,
3,・・・,
8
中心極限定理


(
1
/
10
)

ただし、n
=
2,
3,・・・,
8
n

1
n
。中心極
/ )10
ただし、n
2, 3,・・・,
 n 1 nn1(10.01%）が、文
/ (1規分布で近似できる」というものである
10
 n) nただし、n
= 2,=3,・・・,
8 8
n 1
る。また、表１の 104 の時の誤差（RND で 0.03%、メルセンヌ・ツイスタで
となる。両辺に
log10 を作用して、整理すると、
となる。両辺に log10 を作用
くとも同程度の精度があるように思われる。そ
中心極限定理は「母集団の分布がどのようなものであ
1 / m 倍になることがわかっている
ると標準偏差が
を作用して、整理すると、
となる。両辺に
log10log
献 5)の 111111 回の乱数発生の誤差 0.01％と比較できる。
となる。両辺に
log10 を作用して、整理すると、
log
log1010n n11log
log1010n n00.5.5
を作用して、整理すると、
となる。両辺に
を作用して、整理すると、
となる。両辺に
log10log
10
こで、より詳しい計算精度のチェックを行うた
となる。両辺に
10 を作用して、整理すると、
定義できない分布を除く）
、確率変数の和の確率分布の形
表１のシミュレーションは乱数の発生回数
Nを
これらの結果を総合的に見ると、一見、RND を用いたシミュレーションのほうがメルセ
log


log

 100.56)n 1  log10 （
10
n 1
10
n log
n
めに、統計数学的な手法を用いた解析を行った。
規分布で近似できる」というものである
。中心極限定
ので、標本の数は
10 倍ずつ増えている。よって、中
log


log



0
.
5
（２）
ンヌ・ツイスタのそれよりも精度がよいか、少なくとも同程度の精度があるように思われ
log 10  n 1　log
log



0
.
5
（２）
10
n

1
10
n
n

log
N
なので、乱数の発生回
また、
log


log



0
.
5
（２）
n

log
N
なので、乱数の発生回
また、
10n 1n n1log10 10n n  0.5 （２）
（２）
1010
10 10
ると標準偏差が 1 / m 倍になることがわかっている。
る。そこで、より詳しい計算精度のチェックを行うために、統計数学的な手法を用いた解
 log
N なので、乱数の発生回数
N
なので
N
また、
また、
を描くと、理論的には傾きが－0.5
になると
の標準偏差
n3 は
10210の標準偏差
 2nの log
1 / 1010
にな
中心極限定理
を描くと、理論的には傾きが－0.5
になる
2
表１のシミュレーションは乱数の発生回数
n N を 10 か
n
n の log
n

log
N
なので、乱数の発生回数
N（=10
また、
n）と標準偏差
析を行った。
n

log
N
なので、乱数の発生回数
N（=10
）と標準偏差
10
また、
を描くと、理論的には傾きが－0.5
を描くと、理論的には傾き
になると予想
プロットを示す。
図図44ににNNととn nののlog-log
 log
N なので、乱数の発生回数
なので、乱数の発生回数
N（=10
）と標準
また、
log-log
プロットを示す
中心極限定理は「母集団の分布がどのような
n nlog
N（=10
）と標準偏差
また、
10
n
また、
N
10 なので、乱数の発生回数
10 N
ので、標本の数は
10 倍ずつ増えている。よって、中心極
n
を描くと、理論的には傾きが－0.5
になると予想される。

の
log-log
プロットを示す。

の
log-log
上図
図
4
に
N
と
図
4
に
N
と
を描くと、理論的には傾きが－0.5
になると予想される。
（=10
）と標準偏差
の
log-log
プロットを描く
ものであっても（ただし、期待値や分散が定義
N1
nnnになると予想される。
)が
22から
を描くと、理論的には傾きが－0.5
になると予想される。
の
になる。これを、式で表すと
準偏差 log
n
N(/ (==10
)が
から66までのデータは、直
までのデータは、直
を描くと、理論的には傾きが－0.5
1010
3nlog
2

の
log-log
プロットを示す。
上図は
RND
関数で乱数を
図
4
に
N
と
中心極限定理
と、理論的には傾きが－
の
上図は
RND
関数で乱数を生成した
nプロットを示す。
図 4 に N と図
は
の標準偏差
の
10関数で乱数を
になる。
(プロットを示す。
=10
n0.5
)が
2近い値となった。しかし、n
から
6
までのデータは、直線で
( =1RND
n/ )が
2 から
6がが
まで
log
になると予想され
できない分布を除く）、確率変数の和の確率分
 log
の
log-log
プロットを示す。
上図は
RND
関数で乱
図
にの標準偏差
Nと
10N
10
値の－0.5
近い値となった。しかし、n
7同
上図は
に4log-log
N
と
n4
3 値の－0.5
2N
n のn log-log
N
(
=
n
)が
2
から
6
までのデータは、直線で近似でき、その傾き
log
中心極限定理は「母集団の分布がどのようなものであっても（ただし、期待値や分散が
10
布の形は標本の大きさが大きいときは正規分布
が 7 と88
)がる。
2Nから
までのデータは、直線で近似でき、その傾きは－0.491
log10N ( = n log
図はメルセンヌ・ツイスタで乱数を生成した
n )が
2 値の－0.5
から
log
6n までのデータは、直線で近似でき、その傾き
までのデータは、直線で近似でき、その傾
(1近い値となった。しかし、n
/ 10 ) n 値の－0.5
ただし、n近い値となった
= 2, 3,・・・,
図はメルセンヌ・ツイスタで乱数を生成し
(N= (n6=)が
2 から
10
10
16
6)
値の－0.5
近い値となった。しかし、n
が
7
と
8
では直線から外れ下
定義できない分布を除く）
、確率変数の和の確率分布の形は標本の大きさが大きいときは正
図
4
に
N
と
の
log-log
プロットを示す。上
で近似できる」というものである
。中心極限

の
1
/
10
になる。これを、式で表すと、
準偏差
図はメルセンヌ・ツイスタで乱数を生成した結果
図はメルセンヌ・ツイスタ
値の－0.5 近い値となった。しかし、n
が
7
と
8
では直線から外れ下落してい
n
で近似でき、その傾きは－0.496
となり、理
値の－0.5
近い値となった。しかし、n
8 では直線から外
で近似でき、その傾きは－0.496
となり、
値の－0.5
近い値となった。しかし、n
が 7がと7 8とでは直線から外れ下
3
となる。両辺に log10 を作用して、整理すると、
88
図はメルセンヌ・ツイスタで乱数を生成した結果である。
N8が
2が
から
７
図は
RND
関数で乱数を生成した結果である。
定理によると、標本の数が m 6)倍になると標準
規分布で近似できる」というものである
。中心極限定理によると、標本の数が
m
倍にな
７
で近似でき、その傾きは－0.496
で近似でき、その傾きは－
となり、理論か
図はメルセンヌ・ツイスタで乱数を生成した結果である。
N
が
2
から
までの
RND
の乱数発生回数が
10
および
の
図はメルセンヌ・ツイスタで乱数を生成した結果である。
N10
2
RND
の乱数発生回数が
10
および
図はメルセンヌ・ツイスタで乱数を生成した結果である。N が
210から
8
７
で近似でき、その傾きは－0.496
となり、理論から予想される－0.5
log
N
(
=
n
)
が
2
から
6
までのデータは、直線
偏差が
倍になることがわかっている。
10
ると標準偏差が 1 / m 倍になることがわかっている。
 n 1 RND
(1ぎるためと考えられる。π
/ぎるためと考えられる。π
10
) nとなり、理論から予想される－0.5
ただし、n
= 2,
3,・・・,
の乱数発生回数が
RND
10
の乱数発生回数が
および
10 8のデー
で近似でき、その傾きは－0.496
となり、理論から予想される－0.5
に近い値と
の有効数字は乱数
で近似でき、その傾きは－0.496
となり、理論から予想される－
の有効数字は乱
で近似でき、その傾きは－0.496
log 10 10n７81および
７0.491
log1010
88nのデータの下落は標準偏差
8 8 0.5
（２）
8
７
で近似でき、その傾きは－
となり、理論
表１のシミュレーションは乱数の発生回数
N 10
RND
の乱数発生回数が
７
表１のシミュレーションは乱数の発生回数
N を 102 から
まで
10
倍ずつ増やしている
７
７
ぎるためと考えられる。π
ぎるためと考えられる。π
の有効数字は乱数発生
RND
の乱数発生回数が
10
および
10
のデータの下落は標準偏差が理論値
発生回数が
101010
、10
のシミュレーションの
RND
の乱数発生回数が
および
10
のデータの下落は標準
発生回数が
、10
のシミュレーションの
RND
の乱数発生回数が
10 10
および
のデータの下落は標準偏差
を作用して、整理すると、
となる。両辺に log
2
8
10
3
値の－ 0.5 近い値となった。しかし、n
が 7 と10７、108 のシミ
を 10 から
10 まで 10 倍ずつ増やしているの
ぎるためと考えられる。π
の有効数字は乱数発生の回数とほぼ一致す
ので、標本の数は
10 倍ずつ増えている。よって、中心極限定理より、乱数の発生回数が
10
発生回数が
10７、108 のシミュレーションの有効
発生回数が
ぎるためと考えられる。π
の有効数字は乱数発生の回数とほぼ一致するので、
ぎるためと考えられる。π
の有効数字は乱数発生の回数とほぼ一
ぎるためと考えられる。π
の有効数字は乱数発生の回数とほぼ一致
n

log
N
なので、乱数の発生回数
N（=
また、
8
７
10
8 10 ７、10
8７、10
では直線から外れ下落している。一方、下図
で、標本の数は 10 倍ずつ増えている。よって、
8
発生回数が
7
７ 8のシミュレーションの有効数字は、それぞれ
発生回数が
10
のシミュレーションの有効数字は、それぞれ
7 桁、8 桁と
log
、10
10
4N発生回数が
3 のシミュレーションの有効数字は、それぞれ
発生回数が
10 10
、10
のシミュレーションの有効数字は、それぞれ
7
の標準偏差  3 は 102 の標準偏差  2 の 1 / 10 になる。同様に 10
の標準偏差
は
10
の標
log


log



0
.
5
（２）
4 10
を描くと、理論的には傾きが－0.5
になると予想され
n 1
10
n
0
log10 N
図 4 に N と  n の log-log プロットを示す。上図は R
0
1
2
3
4
5
6
7
8
9
準偏差  3 の 1 / 10 になる。これを、式で表すと、
N（=10n）
また、
n )が
から 6 までのデータは、直線で近似
log10nN ( =log
10 N2なので、乱数の発生回数
0
-1
0
1
2
3
4
5
6
7
8 近い値となった。しかし、n
9
を描くと、理論的には傾きが－0.5
になると予想される。
値の－0.5
が 7 と 8 では
RND
 n 1  (1 / 10 ) n ただし、n-1= 2, 3,・・・, 8
図 4図はメルセンヌ・ツイスタで乱数を生成した結果で
に N と  n の log-log プロットを示す。上図は RND 関
-2
( = n )が 2 から 6 までのデータは、直線で近似でき
log10N
で近似でき、その傾きは－0.496
となり、理論から予
RND
となる。両辺に log10 を作用して、整理すると、
y = -0.491x + 0.1619
７
-2
値の－0.5
近い値となった。しかし、n
が 7 10
と8 8のデータの
では直線
RND
の乱数発生回数が 10 および
-3
2
R = 0.9986
y = -0.491x
+ 0.1619
図はメルセンヌ・ツイスタで乱数を生成した結果である。
ぎるためと考えられる。π の有効数字は乱数発生の
log 10  n 1  log10  n   0.5 -3
（２）
で近似でき、その傾きは－0.496
となり、理論から予想さ
-4
2
発生回数が 10７、108 のシミュレーションの有効数字
R = 0.9986
RND の乱数発生回数が 10７および 108 のデータの下落
N（=10n）と標準偏差  n の log-log プロット
また、 n  log10 N なので、乱数の発生回数
-4
-5
ぎるためと考えられる。π の有効数字は乱数発生の回数と
を描くと、理論的には傾きが－0.5 になると予想される。
発生回数が 10７、108 のシミュレーションの有効数字は、
-5
上図は RND 関数で乱数を生成した結果である。
図 4 に N と  n の log-log プロットを示す。
log10 N
log10N ( = n )が 2 から 6 までのデータは、直線で近似でき、その傾きは－0.491 となり、理論
0
log10 N
値の－0.5 近い値となった。しかし、n が 7 と 8 では直線から外れ下落している。一方、下
0
1
2
3
4
5
6
7
8
9
図はメルセンヌ・ツイスタで乱数を生成した結果である。
N が 2 から 8 までのデータは直線
0
-1
0
1
2
3
4 に近い値となった。
5
6
7
8
9
で近似でき、その傾きは－0.496 となり、理論から予想される－0.5
メルセンヌ・ツイスタ
log10 σ log10 σ
8
RND の乱数発生回数が 10７および-1
-210 のデータの下落は標準偏差が理論値よりも小さす
メルセンヌ・ツイスタ
ぎるためと考えられる。π の有効数字は乱数発生の回数とほぼ一致するので、例えば乱数の
y = -0.496x + 0.1689
-2
8
-3
発生回数が 10７、10 のシミュレーションの有効数字は、それぞれ
7 桁、8 桁となるが、RND
R² = 0.9996
y = -0.496x
+ 0.1689
-3
-4
R² = 0.9996
図 4　乱数の生成回数（N）と標準偏差（σ）
の log-log プロット . 上は RND 関数の結果 . 下はメルセンヌ・
ツイスタの結果 .-4
図4
乱数の生成回数（N）と標準偏差（σ）の log-log プロット. 上は RND 関数
124
水　上　善　博
はメルセンヌ・ツイスタで乱数を生成した結果
れていると考えられる。
である。N が 2 から 8 までのデータは直線で近
エクセルの RND 関数は手軽にモンテカルロ
似でき、その傾きは－ 0.496 となり、理論から
シミュレーションができて便利であるが、乱数
予想される－ 0.5 に近い値となった。
の発生回数としては 106 程度までが適当であり、
RND の乱数発生回数が 107 および 108 のデー
107 回や 108 回の乱数発生のシミュレーション
タの下落は標準偏差が理論値よりも小さすぎる
で良い結果が得られる場合は同一の乱数が発生
ためと考えられる。πの有効数字は乱数発生の
しているための見かけ上のものであると考えら
回数とほぼ一致するので、例えば乱数の発生回
れる。すなわち、エクセルの RND 関数で生成
数が 107、108 のシミュレーションの有効数字は、
する乱数を用いたシミュレーションは最大でも
それぞれ 7 桁、8 桁となるが、RND 関数の精度
106 程度までの乱数発生にとどめるべきである
が単精度 10 進数でいえば 7 桁程度までとなっ
といえる。
ているため、10 回以上乱数を発生させると、
図５に乱数の発生回数が 106 でπを求めるモ
同一の乱数が生じて、標準偏差が小さくなると
ンテカルロシミュレーションを 100 回実行した
7
考えられる。すなわち、RND 関数を用いて 10
結果のヒストグラムを示す。上図が RND、下
回以上乱数を発生させるシミュレーションは信
図がメルセンヌ・ツイスタの結果である。RND
頼性に乏しいと言える。
はあまりきれいな正規分布には見えないが、メ
一方で、メルセンヌ・ツイスタは少なくとも
ルセンヌ・ツイスタでは比較的きれいな正規分
今回研究を行った乱数の発生回数が 108 までの
布を示している。発生回数が 106 の場合でも、
シミュレーションでは、中心極限定理で予想さ
メルセンヌ・ツイスタの乱数を用いたシミュ
れるゆらぎをもつバランスの良い乱数が生成さ
レーションの方が信頼性が高いと考えられる。
頻度
頻度
7
30
30
25
25
20
20
15
15
10
10
5
5
0
0
3.138 3.139 3.14 3.141 3.142 3.143 3.144 3.145 3.146
3.138 3.139 3.14 3.141 3.142 3.143 3.144 3.145 3.146
頻度
頻度
π
π
6
35
35
30
30
25
25
20
20
15
15
10
10
5
5
0
0
3.138 3.139 3.14 3.141 3.142 3.143 3.144 3.145 3.146
3.138 3.139 3.14 3.141 3.142 3.143 3.144 3.145 3.146
π
π
回の乱数発生のモンテカルロシミュレーションを 100 回実行したときのπの値のヒストグ
図 5　10
図 5 106 回の乱数発生のモンテカルロシミュレーションを 100 回実行したときの π の値の
6 回の乱数発生のモンテカルロシミュレーションを
ラム。上図は
VB
の RND 関数で乱数を生成した場合
. 下図はメルセンヌ・ツイスタで乱数を
図 5 10
100 回実行したときの
π の値の
ヒストグラム。上図は
VB の RND 関数で乱数を生成した場合.
下図はメルセンヌ・
生成した場合 ツイスタで乱数を生成した場合.
.
ヒストグラム。上図は
VB の RND 関数で乱数を生成した場合. 下図はメルセンヌ・
ツイスタで乱数を生成した場合.
モンテカルロシミュレーションによる円周率πの計算
125
6
表 2　様々な性能のパソコンにおいてエクセルの RND 関数で 10 回乱数を発生させて円周率πを計
算したモンテカルロシミュレーションの結果。
（#1
と #2 は 10 回、#3 ～ #14 は 100 回のシミュ
6
回乱数を発生させて円周率πを計算したモンテカルロ
表2　様々な性能のパソコンにおいてエクセルのRND関数で10
レーションを実行した）
シミュレーションの結果。（#1と#2は10回、#3～#14は100回のシミュレーションを実行した）
Microsoft
WindowsのOSの
バージョン
Microsoft
Excel
のバージョン
パソコン
タイプ
#1
Windows 95
Excel 97
デスクトップ
Intel®Pentium®
64MB
133MHz
49.2
3.1415308 0.00180338
#2
Windows 95
Excel 97
ノート
Intel®Pentium®
32MB
120MHz
25.0
3.1414592 0.00180756
#3
Windows XP
Excel 2007
デスクトップ
Intel®Pentium®4
1GB
2.60GHz
1.66
3.1415766 0.00167359
#4
Windows XP
Excel 2007
デスクトップ
Intel®Pentium®4
2GB
3.00GHz
1.36
3.1415766 0.00167359
#5
Windows XP
Excel 2007
デスクトップ
Intel®Pentium®4
3GB
3.40GHz
1.21
3.1415766 0.00167359
#6
Windows Vista
Excel 2007
ノート
TM
Intel®Core 2
2.53GHz
4GB
1.11
3.1415679 0.00168244
#7
Windows Vista
Excel 2010
ノート
TM
Intel®Core 2
Duo
2.00GHz
2GB
0.88
3.1415766 0.00167359
TM
Intel®Core 2
デスクトップ
Duo
2.66GHz
4GB
0.76
3.1415766 0.00167359
4GB
0.70
3.1415766 0.00167359
24GB
0.57
3.1415766 0.00167359
12GB
0.56
3.1415766 0.00167359
32GB
0.47
3.1415902 0.00162192
8GB
1.00
3.1415679 0.00168244
4GB
0.43
3.1415672 0.00169303
#8
Windows 7
Excel 2007
CPU
RAM
平均計算
πの平均値
時間（秒）
標準偏差
TM
#9
Windows 7
Excel 2013
ノート
Intel®Core
i5-3337
1.80GHz
#10
Windows 7
Excel 2013
デスクトップ
Intel®Core
i7-960
3.20GHz
#11
Windows 7
Excel 2013
デスクトップ
Intel®Core
i7-960
3.20GHz
#12
Windows 7
Excel 2013
デスクトップ
Intel®Core
i7-3970X
3.50GHz
ノート
Intel®Core
i3-3217U
1.80GHz
デスクトップ
Intel®Core
i5-3570
3.40GHz
TM
TM
TM
TM
#13
Windows 8.1
Excel 2013
TM
#14
Windows 8.1
Excel 2013
126
水　上　善　博
計算時間
6 月に閣議決定された「世界最先端 IT 国家創
様々な性能のパソコン上においてエクセル
造宣言」8) では、「初等・中等教育段階でのプロ
の RND 関数を用いて 106 回の乱数を発生させ
グラミング、情報セキュリティ等の IT 教育の
たモンテカルロシミュレーションでπの計算を
充実」および「高等教育段階での産業界と教育
した結果を表２に示す。Windows 95 のパソコ
現場との連携の強化の推進」が述べられている。
ンでは数十秒かかるが、Windows XP 以降のパ
ICT（情報通信技術）は重要な教育課題となっ
ソコンでは 2 秒以内で計算が終了することがわ
ており、小学校や中学校でもプログラミングの
かる。また、RAM の大きさよりも CPU の性
基礎について学ぶことが推奨されつつある。そ
能の方が計算時間の短縮に重要であることがわ
こで、教育現場が現有している設備を利用して
かった。さらに、シミュレーションを 100 回実
手軽にプログラミングの基礎を学べるような教
行して平均をとったπの値を見ると、#3, 4, 5, 7,
材開発が望まれる。本研究で示したようなエク
8, 9, 10, 11 で同じ値となっている。これは、初
セルの VBA のプログラムを用いたシミュレー
期化せずに乱数を発生させたため同一の系列の
ションは比較的簡単に実施できるので、小、中
乱数が生成された結果である。#6 と #13 も同
学校におけるプログラミングの教材としても適
様である。なお、Windows XP 以降のパソコン
切であると思われる。
のπの平均値は小数点以下 5 桁目を四捨五入す
我々は人間同士のコミュニケーションの道具
る と、3.1416 に な り、 真 値（3.1415926535897）
として日本語や英語などの自然言語を学ぶが、
の小数点以下 5 桁目を四捨五入した値と一致す
コンピュータとコミュニケーションをとるため
る。
の人工言語でプログラミングを学ぶ意義のひと
教 育 情 報 化 振 興 協 会 が 2013 年 に 実 施 し た
7)
つは、プログラミング言語の文法やアルゴリ
調査報告によると 、小学校 632、中学校 351
ズム ( 問題を解くための手順、ものの考え方 )
の合計 983 校で児童・生徒が使っているコン
を学ぶことにある。プログラミング言語には、
ピュータの OS の割合は、Windows 98 が 0.3%、
Basic、C 言 語、Fortran、Java な ど 数 多 く 存
Windows XP が 31.4%、Windows Vista が
在するが、現在の小、中学校におけるパソコン
18.5％、Windows 7 が 50.1%、Windows 8 が 3.2%
環境で手軽に実施できるプログラミング教材と
となっている。小、中学校の授業用のパソコン
しては、一つは、HTML 言語によるホームペー
の OS は 99％以上が Windows XP 以降のバー
ジ作成がある。HTML 言語は文法がシンプル
6
ジョンとなっており、本研究の結果より 10 回
で扱いやすく、プログラムもパソコンに標準装
の乱数発生のモンテカルロシミュレーションに
備されている「メモ帳」で書ける。また、作成
かかる計算時間は 2 秒以内であるので、授業で
したプログラム ( すなわち作成したホームペー
十分に実行可能なものであると思われる。
ジ ) は、やはりパソコンに標準装備されている
インターネットのホームページを閲覧するブラ
学校の授業での活用
ウザを用いて確認できる。
文献 7) によると、小、中学校におけるコン
もう一つのプログラミング教材としては、中
ピュータ活用の目的としては、「インターネッ
学校の技術科で以前から取り上げられてきた
トを活用して「調べ学習」などをさせる」が
Basic がある。Basic は数年前までは大学入試
93.2% と最も多く、次いで、
「文字を入力するな
センター試験の数学の分野でも出題されていた
どの基本的な操作を行わせる」65.5%、「プレゼ
が、新課程の数学の内容からコンピュータが削
ンテーション能力を高める」46.0%、「情報化社
減されたことを受けて、2015 年実施の大学入
会の特徴（利点、危険、セキュリティ、個人情
試センター試験からは出題されなくなった。し
報保護など）を学習させる」44.2% となっており、
かし、プログラミングを通じてアルゴリズムを
シミュレーションやプログラミングにつながる
学ぶことは重要なので、今後は、Basic が発展
「計算能力など児童・生徒の基礎学力を向上さ
してビジュアル環境に適応した、Visual Basic
せる」は 12.7% にとどまっている。平成 27 年
でプログラミングを学ぶことが有意義であると
モンテカルロシミュレーションによる円周率πの計算
127
考えられる。Visual Basic は表計算ソフトのエ
4）　和田維作氏によって作成された、エクセ
クセルにおいてマクロ機能を利用したプログ
ルのためのダイナミック・リンク・ライ
ラミング言語 (VBA) としても用いられており、
ブ ラ リ (DLL) 一 式（libZMT.zip）　http://
本研究で示した例のように、手軽に興味深いシ
www001.upp.so-net.ne.jp/isaku/rand2.html
ミュレーションを実行できる。
(2015 年 9 月 11 日閲覧 )
5）　“Estimation of the value of
π using
大学の演習授業での実践
Monte-Carlo Method and Related Study
大学のコンピュータ実習におけるモンテカ
of Errors, Udayan Singh, Mathematics in
ルロシミュレーションの実践として、滋賀大
School, 21-23 (2013)
学教育学部で開講している教科「情報」の教
員免許の取得に必要な科目の一つである「モデ
6）　「統計学入門」，東京大学教養学部統計学教
室編，東京大学出版会 (1991)
ル化とシミュレーション」の授業において、エ
7）　「第９回教育用コンピュータ等に関するア
クセルによるπのモンテカルロシミュレーショ
ンケート調査報告書」一般社団法人教育情
ンを実施した。乱数の発生回数は 1 から 10 倍
報化振興協会（2014）
7
8
ず つ 増 や し、10、100、・・・、10 、10 ま で
として、それそれ、1 回ずつシミュレーション
8
を実行した。乱数発生の回数が 10 回のシミュ
レーションでは、結果が出るまでに 44 秒程度
http://www2.japet.or.jp/info/japet/report/
ICTReport9.pdf　(2015 年 9 月 11 日閲覧 )
8）　「世界最先端 IT 国家創造宣言」の変更につ
いて（平成 27 年 6 月 30 日　閣議決定）
かかった。結果が出るまでコンピュータが比較
https://www.kantei.go.jp/jp/singi/it2/
的長い時間計算する姿はあまり目にしたことが
kettei/pdf/20150630/siryou1.pdf　(2015 年
ないようで、学生たちは興味を示していた。な
9 月 11 日閲覧 )
お、誤って乱数の発生回数の入力を 1 桁間違え
て 109 として、シミュレーションを実行してし
まうと、計算結果が得られるまで、数十分程度
コンピュータがかたまってしまうので、乱数の
発生回数の入力には注意を促す必要がある。ま
た、RND 関数で乱数を発生するシミュレーショ
ンの場合、乱数の発生回数が 107 回以上の結果
は信頼性が低いことにも言及する必要がある。
謝辞
有益な助言をいただいた丸目諒氏と高田優介
氏に感謝いたします。
文献
1）
「Excel を用いたシミュレーション演習授
業 」， 松 永 豊，Bulletin of Aichi Univ. of
Education, 60 (Education Science), 187-190
(2011)
2）　「πの歴史」，ペートル・ベックマン著，田
尾陽一，清水韶光訳，ちくま学芸文庫（2006）
3）　「良い乱数・悪い乱数」
http://www001.upp.so-net.ne.jp/isaku/
rand.html (2015 年 9 月 11 日閲覧 )