...

セル・オートマトン セル・オートマトン(cellular automata)の一例 として

by user

on
Category: Documents
4

views

Report

Comments

Transcript

セル・オートマトン セル・オートマトン(cellular automata)の一例 として
環境シミュレーション
セル・オートマトン
セル・オートマトン(cellular automata)の一例
として,ここではイギリスの科学者スティーブン・
ウォルフラムが開発したセル・オートマトンの一次
元モデルを紹介する(ウォルフラムは数学ソフト
"Mathematica"の作者としても名高い).
1)4 行以下のシートには実際には「1」または「0」
の数値が入っている.これを[セルの書式設定][ユー
ザ定義]を用いて「[青][=1]"●";[=0]"";G/標準」と指
定したものである.また全ての列幅を 4 mm に統一
している.
ウォルフラムのセル・オートマトンは,一つひと
つのセルの白黒が,そのセルの上段の三つのセル
(真上,右上,左上)の状態(白黒)のみによって
規定されるという単純なものである.しかし,規則
の設定次第では,フラクタル図形を描く不思議なオ
ートマトン(自動機械)である.
2)プログラムとしては先ず A1 セルであるが,こ
こは 10 進法で 0 から 255 を入力する場所になって
いる.これは,セルが白黒(このプログラムでは青
●か空欄)のいずれかになるかのルールの組合せが,
256 通りあることに対応している.256 通りとなる
理 由 は , 上 段 三 つ の セ ル の 白黒の組み合わせが
2x2x2 = 8 通りであり,それぞれの組合せにおいて
下のセルが白黒のいずれかとなる 2 通りがあるので,
ルールとしては全部で 28 = 256 通りとなるのであ
る.
スピンボタン
ちなみに上図は[表示][ズーム]倍率指定 10%とし
たもので,これを通常の 100%で見ると下図のよう
になっている.
3)B1:I1 には下表のような式が入っている.この
B1:I1 は A1 の値を 2 進法で表したものである.
下セル○か●か
右上
真上
左上
- 11 -
環境シミュレーション
B1
C1
D1
E1
F1
G1
H1
I1
=INT(A1/2^7)
=INT((A1-2^7*B1)/2^6)
=INT((A1-2^7*B1-2^6*C1)/2^5)
=INT((A1-2^7*B1-2^6*C1-2^5*D1)/2^4)
=INT((A1-2^7*B1-2^6*C1-2^5*D1-2^4*E1)/2^3)
=INT((A1-2^7*B1-2^6*C1-2^5*D1-2^4*E1-2^3*F1)/2^2)
=INT((A1-2^7*B1-2^6*C1-2^5*D1-2^4*E1-2^3*F1-2^2*G1)/2)
=A1-2^7*B1-2^6*C1-2^5*D1-2^4*E1-2^3*F1-2^2*G1-2*H1
「=INT()」は Excel 関数で,()内の
数の小数点以下を切り捨て整数を
返す機能をもつ
4)つづいて 4 行(A4:IO4)(ここでの IO4 はア
イ・オー・4)には「1」または「0」の値が 50%ず
つの確率でランダムに入っている.設定方法として
は 先 ず A4:IO4 の 全 て の セ ル に
「=INT(RAND()+0.5)」を入力,これで「1」また
は「0」のランダムな数値を発生する.しかしこの
ままでは,シートの再計算のたびにこの範囲の値が
変化するため,それを防ぐために,A4:IO4 をすべ
てコピーし,同じ場所に[形式を選択して貼り付
け][値]でペーストする.
こでは,スピンボタンを設定する.スピンボタンを
利用すれば,マウスクリックで A1 の値を変化させ
ることが可能となる.
5)つづいて B5 に
を選択.「+」字
ツールバーからスピンボタン
カーソルが画面に現れるので,これをドラッグして
スピンボタンの場所と大きさを確定する.じゃまに
ならないよう,画面右下角がいいだろう.次に,ス
ピンボタンを右クリック,現れるメニューから[コン
トロールの書式設定]を選択.
スピンボタンの挿入には,[表示][ツールバー]から
[フォーム]を選ぶ.すると画面に下のようなフォー
ムのツールバーが現れる(縦横長はドラッグで自由
に変更できる).
=IF(A4=1,IF(B4=1,IF(C4=1,$B$1,$C$1),IF(C4=1,
$D$1,$E$1)),IF(B4=1,IF(C4=1,$F$1,$G$1),IF(C4
=1,$H$1,$I$1)))
の式を入力する.この式は上段三つのセルの白黒の
組み合わせ 8 通りのそれぞれにおいて下のセルが
「1」となるか「0」となるかのルールを B1:I1 の
1/0 で 指 定 し た も の . あ と は , こ の 式 を 全 範 囲
(A5:IO274)にコピーしてやればよい.注意:少
しでもファイルサイズを小さくするために,倍率
10%ズームで画面上で見えるセル範囲(A5:IO274)
だけに限定している.
ただし,A 列と最後の IO 列については注意が必
要である.なぜなら A 列なら左上のセルが,IO 列
なら右上のセルに値がないからである.このため例
えば A5 には
=IF(IO4=1,IF(A4=1,IF(B4=1,$B$1,$C$1),IF(B4=1
,$D$1,$E$1)),IF(A4=1,IF(B4=1,$F$1,$G$1),IF(B4
=1,$H$1,$I$1)))
IO5 には
=IF(IN4=1,IF(IO4=1,IF(A4=1,$B$1,$C$1),IF(A4=
1,$D$1,$E$1)),IF(IO4=1,IF(A4=1,$F$1,$G$1),IF(
A4=1,$H$1,$I$1)))
[コントロール]タブで最小値 0,最大値 255,変化
の増分 1,リンクするセルとして「A1」を指定して
完成である.
の式を入れ,これらの式をそれぞれ A 列と IO 列の
最後の行までコピーしてやる.上記の 2 式が意味す
るところは,A 列なら最後の IO 列 1 行前のセルを
左上のセルとして,IO 列なら最初の A 列 1 行前の
セルを右上のセルと見なすということである.した
がって,画面上は平面に見えるが,A 列と IO 列が
つながった円筒になっていると考えればよい.
ちなみにこのスピンボタンは,A1 に手入力した
値から開始することも可能である.
6)基本的には5)までのステップでセル・オート
マトンは完成である.A1 に 0-255 の値を代入すれ
ば,描画の変化を見ることができる.ただし,いち
いち A1 に値を入力するのは手間である.そこでこ
- 12 -
環境シミュレーション
回数を n としたならば,その 2 つの回数の比率はど
うなるだろう.もし N が十分大きければ,この N 対
n の比は,正方形の面積 1 と 4 分の 1 円の面積 π / 4
との比に,ほぼ等しくなると予想される.
3.モンテカルロで乾 π(かんぱい)
ここでは,モンテカルロ法と呼ばれているシミュ
レーション手法によって円周率 π を求めてみる.
先ずは,題名「モンテカルロで乾 π」についてであ
るが,もともと「乾」という文字は「高い空(天)」を意
味するものであった(通常,われわれが「乾」の字から
思い浮かべる「乾かす」の意は,高く空に掲げて干す
から転じたもの).これに地を意味する「坤」をつなげ
れば「乾坤(けんこん)」となる.よく聞く「乾坤一擲
(けんこんいってき)」とは,天地を賭けてサイコロを
ふること(要するに,一か八かの大勝負にでること).
モンテカルロ法とは,サイコロをふって答えを見つ
けるような方法だから「乾 π」とシャレてみた.(庄
野真代の唄った「モンテカルロで乾杯」を知らない
人は,お父さんか,お母さんに聞いてみよう.)
近代シミュレーションは,J.von Neumann(ゲー
ム理論でも出てきたノイマン)と S.M.Ulam の両名
が水爆の研究開発の中で,物質中における中性子の
拡散現象をモンテカルロ法によってコンピュータ
で模擬したことに発すると言われている.高速にラ
ンダムな方向へ飛び回る中性子の動きをコンピュ
ータ内で乱数を使って数値的に再現しようとした
のである(大気汚染の拡散シミュレーションに用い
られているランダムウォークモデルもランダムな
方向へ飛び回っている汚染物質の運動をモンテカ
ルロ法で再現しようとしたもの).
N: n ≅1:
π
4
上式を変形すれば, π ≅ 4n N となり, π は 4n N
で近似できることがわかる.直感的にもわかるよう
に,この近似は N が大きければ大きいほど,その精
度が増す.
シミュレーション・プログラムの作り方
ここではモンテカルロ法による π 計算を 3 種類
の方法で試みる.これら 3 種類に,今後の授業で用
いる基本的な Excel VBA のテクニックがすべて含
まれているので十分に習熟してほしい.
(1) Excel 関数だけを用いた方法
先ず下図が Excel 関数のみを用いて,モンテカル
ロ法(100 回の試行)によって円周率 π を求めて
みたワークシート(1)である.
このようにコンピュータ内で乱数を発生させる
ことによって確率現象を模擬して,もとの現象を解
明しようとするシミュレーション手法を,乱数を用
いることから,カジノで有名なモナコの都市モンテ
カルロから名をとって「モンテカルロ・シミュレー
ション」と呼ぶ(もともと,このノイマンらのプロジ
ェクトが最高軍事機密であったことから,軍内部で
の隠語としてそう呼ばれていたようである).以前は,
シミュレーションといえばモンテカルロ・シミュレ
ーションを指す時代もあった.
A1:E5 を拡大すると下のようになる.
では,モンテカルロ法による π の求めかたを見
てみよう.先ずは,下図のように,半径 1 の 4 分の
1 円と,その 4 分の 1 円を囲む一辺 1 の正方形を考
えてみる.
それぞれのセルに入っている式は下表のとおり.
B3
C3
D3
B5
C5
D5
E5
もしこの図の上にランダムな点 ( x, y ) を打ち,そ
の点が正方形の中に入る回数を N ,さらにその内で
4 分の 1 円の中に入る,すなわち x 2 + y 2 ≤ 1 となる
- 13 -
=COUNT(B4:B104)
=COUNT(C4:C104)
=4*C3/B3
=RAND()
=IF(B5^2+E5^2<=1,E5,"")
=IF(B5^2+E5^2>1,E5,"")
=RAND()
環境シミュレーション
また,A5:A104 には 1 から 100 までの通し番号
が,B6:E104 には,B5:E5 をコピー・ペーストした
式が入っている.
ここで B5 と E5 は RAND 関数を用いて 2 つの乱
数(0-1 の一様乱数)を発生させているセルであ
る.それぞれが x と y に相当する.つぎに C5 では,
x 2 + y 2 (B5^2+E5^2)を計算して,もし(IF)
x 2 + y 2 ≤ 1(B5^2+E5^2<=1)ならば y の乱数(E5)
を,そうでなければ空の文字列「""」を,また D5
では逆に x 2 + y 2 > 1 (B5^2+E5^2>1)ならば y の
乱数(E5)を,そうでなければ空の文字列「""」を
打ち出す.
B3 と C3 は Excel の COUNT 関数を用いて,括
弧内に指定した B4:B104 と C4:C104 の範囲にある
数値の個数をそれぞれカウントしている.例示され
ている B3 の 100 と C3 の 82 は,N =100 組の ( x, y )
を生成させたところ x 2 + y 2 ≤ 1 となったものが n =
下図のような系列設定になっていることを確認
する.もしそうなっていない場合は,系列の削除と
追加を使って,2 つの系列を下図のように指定する.
82 組あったことを表す.
最後に D3 の「=4*C3/B3」は,π の近似値とし
て 4n N を求めたものである.例では 3.28 となって
いる.
次に乱数で生成した ( x, y ) の点をグラフ上に描い
てみる.
1) B4:D104 を指定して,グラフウィザード
をク
リック,ウィザード画面で「散布図」を指定する.
出 来 上 が っ た グ ラ フ で は x2 + y2 ≤1 の 点 と
x 2 + y 2 > 1 の点が色分けされて表示されるはずであ
る.
(2) Sub プロシージャを用いた方法
つぎに 2 つ目の方法として下図が Sub プロシージ
ャという VBA プログラムを用いて,モンテカルロ
法で円周率 π を求めてみたワークシート(2)である.
コマンドボタン
[次へ]をクリックして,[系列]タブへ.
- 14 -
環境シミュレーション
シート上の「スタート」と書かれたコマンドボタ
ンをクリックすると,VBA プログラムが起動しワ
ークシート上に 3000 組の乱数 ( x, y ) を発生させ,
グラフにプロットするとともに,ワークシート(1)
と同様に D3 セルで π を計算するようになっている.
作成方法としては,先ず前掲のワークシート(1)
のコピーを作成する.作成方法はシート(1)を開き
[編集][シートの移動またはコピー]を選び,[コピー
を作成する]のチェックボックスを して[OK]ボタ
ンを押す.
続いて Sub プロシージャ(VBA プログラム)を
書き込むための,モジュールシートを挿入する.
モジュールの挿入
Sub プロシージャや後述する Function プロシー
ジャ(両者をあわせて VBA プログラムと呼ぶ)を
作成するためには,モジュールと呼ばれるシートが
必要である.モジュールを挿入するには,先ず[ツー
ル]から[マクロ],[Visual Basic Editor]を選ぶ.(注
意:Excel 2003 以降では VBA プログラムを使用す
るために,モジュールの挿入の前にセキュリティレ
ベルの設定変更が必要となる.そのためには 1)新し
いブックを開く,2)[ツール][マクロ][セキュリティ]
を選び,セキュリティレベルを「中」に設定,3)そ
の後,Excel を一旦終了して,再起動する,4)マク
ロの有効無効を聞いてくるので,[マクロを有効にす
る]をクリックする.これで VBA が使用可能とな
る.)
実は先ほどのオプション設定の意味は,この文頭
の一行をモジュールシートに自動的に挿入するた
めのものであった.この「Option Explicit」の意味は,
もし,そのシートに書かれるプログラム中に,デー
タ型を宣言していない変数が使われた場合には,必
ずそのことをエラーとして警告するように,という
ものである(データ型宣言の意味については,また後
述する).ここでは簡単に,変数のデータ型を一つ一
つ明示しておいたほうが,プログラムのエラーが見
つけやすくなる,とだけ言っておこう.
モジュールシートにプログラムを書いてゆく前
に,もう一つだけ設定をしておこう.[表示][ツール
バー]から[デバッグ]を する.これによって下図の
ようなツールバーが画面上に現れる.
[Visual Basic Editor]を選ぶと下の画面のような
[Visual Basic]編集画面が現れる.
このツールバーは,VBA プログラムを作成する
際に便利な機能をツールボタンとしてまとめたも
ここでモジュールシートを挿入する前に,若干の
オマジナイをしておく.それは,[ツール][オプショ
ン]の[編集]で[変数の宣言を強制する]のチェックボ
ックスを しておくことだ.この設定の意味は後述
する.編集画面でさらに[挿入]から[標準モジュール]
を選ぶと,文頭に「Option Explicit」と書かれた白紙
のシート(モジュールシート)が現れる.(またこ
こで[挿入]から[ユーザーフォーム]を選ぶと後述す
るフォーム画面が現れる.)
(F8 キー)を使え
のである.特にステップイン
ば,プログラムを 1 行 1 行実行させることができる.
Excel 2003 以降では,途中,変数にマウスカーソル
を合わせるだけでその変数の保持している値がウ
ォッチできるので便利である.また,このツールバ
ーは,ウィンドウ内のどこにでも移動させることが
できる.
さてツールバーの用意までできたところで,モジ
ュールシートに以下の文字列を一字一句間違えず
に打ち込んでいこう.ちなみに VBA では,大文字
小文字の区別は原則としてない.ただし,モジュー
ルの編集機能として,最初に使用した用法に合わせ
て,大文字小文字を自動的に修正する機能がある.
- 15 -
環境シミュレーション
Sub pi_dot()
Dim i As Integer
Dim x, y As Double
Range("A5:E3004").Clear
For i = 1 To 3000
Cells(4 + i, 1) = i
x = Rnd: y = Rnd
'[0,1]の一様乱数を返す
Cells(4 + i, 2) = x
If x ^ 2 + y ^ 2 <= 1 Then
Cells(4 + i, 3) = y
Else
Cells(4 + i, 4) = y
End If
ActiveSheet.Calculate
Next i
最後にボタンの名称をクリックして「スタート」
と書き直せば完成である.ちなみに一度作成したコ
マンドボタンを編集するためには,先ず<Ctrl キー>
を押しながらそのボタンをクリックして選択して
やる.また登録マクロの変更にはボタンを右クリッ
クして[マクロの登録]を選択してやる.
End Sub
以上の一連の文が pi_dot と名づけた Sub プロシ
ージャである.Sub プロシージャの基本的構文は以
下のようになっている.Sub 名称...End Sub.
なお Sub プロシージャの起動自体はワークシー
トの[ツール][マクロ][マクロ]で実行可能なマクロ
(プロシージャ)一覧がでるので,そこから pi_dot
を選択,[実行]ボタンをクリックしてやればよい.
Sub pi_dot()
...
End Sub
[Visual Basic Editor]に戻って,Sub プロシージ
ャ pi_dot を1行ずつ見ていこう.
ここでは「pi_dot」がこの Sub プロシージャの名
称である.Sub プロシージャは,引数を名称の後ろ
の()の中にとることもできるが,ここでは取ってい
ない.Sub プロシージャと,後述する Function プ
ロシージャとの違いは,後者が値を返すのに対して
前者は値を返さない点にある.また,Sub プロシー
ジャだけは,Excel のシートからマクロとして呼び
出せる特徴を持っている(ただし,Sub プロシージャ
も引数をとると,マクロとして呼び出せなくなる).シ
ート上の「スタート」コマンドボタンをクリックす
ると,この「pi_dot」がマクロとして起動することに
なる.
Dim i As Integer
Dim x, y As Double
最初にプログラムの冒頭に現れる「Dim」につい
て説明する.これは変数のデータ型を宣言し,型に
応じたメモリ領域を変数に割り当てるための宣言
ステートメント(命令)である. ここでは i を整数
(Integer),x と y を倍精度浮動小数点数型(Double)
と宣言している.これらのデータ型の違いの意味は
後述する.TIPS:モジュールの中では,調べたい文字に
カーソルを置き,[F1]キーを押せば,ヘルプ画面が現れ,
それらの説明文を見ることができる.
コマンドボタンへのマクロの登録
本来なら Sub プロシージャを完成させてから,コ
マンドボタンへマクロ登録するのであるが,ここで
は最初に登録方法を説明しておこう.コマンドボタ
ンを描画するためには先ず,ワークシートに戻り,
[表示][ツールバー]から[フォーム]を選択,ツールバ
ーからボタン
を選択する.すると「+」字カーソ
ルが画面に現れるので,これをドラッグしてコマン
ドボタンの場所と大きさを確定する.ここで[マクロ
の登録]ダイアログボックスが現れ,起動できるマク
ロの一覧がでてくるので,「pi_dot」という Sub プロ
シージャを選択する.これによってボタンに
「pi_dot」がリンクされたことになる.
Range("A5:E3004").Clear
続く上記の1行はワークシートの A5:E3004 の範
囲(Range)をクリアするための命令文である.
For i = 1 To 3000 Step 1
...
Next i
本文中の上記の部分は,VBA によるフロー制御
ステートメント(For...Next ステートメント)である.
- 16 -
環境シミュレーション
ここでは,変数 i が 1 から 3000 までの間(For),そ
の下のステートメント群を処理し,Next 行に到達
したら,次(Next)の i で処理を繰り返す.For...Next
ステートメントで囲まれる部分は,プログラムを見
やすくするために,インデントすることが推奨され
ている.
上ステートメントは,ランダムな点 ( x, y ) が半径
1 の 4 分の 1 円内にあるかどうかを判断している.
もし(if)(x ^ 2 + y ^ 2)が 1 より小さか,あるいは
等しければ(<= 1)(4 分の 1 円内ならば),その場合
(Then)は,4+i 行 3 列のセルに先ほど生成した乱数
y を打ち出す.もしそうでなければ(Else)(4 分
の 1 円外ならば)4+i 行 4 列のセルに y を打ち出す.
ちなみに VBA におけるべき乗演算には「^」(ハット
マーク)をつかう.
「Step」は繰り返す際,変数 i の増加の幅を指定す
る.ただし,1 の場合(Step 1)は省略可能である.
If ステートメントの終了を宣言するのが「End
if」である.
Cells(4 + i, 1) = i
続く Cells は VBA の関数で,現在開いているワ
ークシート上のセルに値を打ち出すためのコマン
ドである.この例では,4+i 行 1 列のセルに i を打
ち出すことになる.
ActiveSheet.Calculate
ワークシート上の計算をやり直し,グラフを描き
直すための VBA コマンドである.このコマンドが
ないと,Sub プロシージャがすべて終了してからで
ないとグラフは描き直されない.
x = Rnd: y = Rnd
'[0,1]の一様乱数を返す
End Sub
次に「x = Rnd: y = Rnd」についてである.先ず行
中の「:」は,本来は 2 行に書かれるべきプログラムが
1 行に詰め込まれていることを示している.TIPS:
Sub プロシージャの終了を示す.
VBA は固定フォーマット言語で,ひとつのステートメン
トは 1 行で完結しなければならない.もし複数行に渡る場
合には「_」(アンダースコア)を用い,複数ステートメント
を 1 行に納めるためには「:」(コロン)を用いる.
ちなみに,以上の Sub プロシージャは,ワークシ
ート(1)において手作業で行ったことを再現してい
るにすぎない.2 つの乱数 ( x, y ) を発生させて, x
その下の 1 行はコメント(注釈)行である.「'」は,
それ以下の行末までの文がコメント文であること
を示す.コメント文はプログラムの実行時に無視さ
れる.
ここで Rnd は Rnd 関数と呼ばれ,0 以上 1 未満
の範囲で一様乱数値を返す VBA 関数である.先に
出てきた RAND 関数(RAND())は,機能はまったく
同じであるが Excel 関数であり,Excel 関数は,モ
ジュールシート内で,そのままでは使用できないの
で注意が必要である.
の乱数をワークシート B 列に打ち出し, x 2 + y 2 ≤ 1
ならば y の乱数を C 列に,x 2 + y 2 > 1 ならば D 列に
打ち出している.
最後にシート(1)の 100 回の試行をシート(2)で
3000 回に増やしたことによる修正を施す.
1) B3 と C3 の COUNT 関数の指定範囲をそれぞ
れ「B4:B3004」と「C4:C3004」に変更する.
2) グラフを右クリック,[元のデータ]で[データ範
囲 ] を 「 ='(2)'!$B$4:$D$3004 」 , [ 系 列 ]
「 x^2+y^2<=1 」 の 「 x の 値 」 を
「 ='(2)'!$B$5:$B$3004 」 「 y の 値 」 を
「='(2)'!$C$5:$C$3004」に,[系列]「>1」の「x
の値」を「='(2)'!$B$5:$B$3004」「y の値」を
「='(2)'!$D$5:$D$3004」に変更する.
すなわち,ここでは Rnd 関数によって ( x, y ) を生
成している.
Cells(4 + i, 2) = x
ただし例に示されるように 3000 回の試行でも π
の近似値は 3.121333 とまだ理論値 3.141593...との
間には開きが見られる.より多くの試行回数が必要
のようであるが,ワークシート上に値を打ち出すこ
の方法には限界がある.一般のワークシートの行数
は 65536 が限界であり,またこれだけの数値を発生
させると Excel ファイルが大きくなりすぎてしまう.
そこで次には Function プロシージャを用いた方法
を試してみる.
続く行は,Cells によって 4+i 行 2 列のセルに先
ほど生成した乱数 x を打ち出している.
次に If...Then...Else ステートメントについて述
べる.
If x ^ 2 + y ^ 2 <= 1 Then
Cells(4 + i, 3) = y
Else
Cells(4 + i, 4) = y
End If
(3) Function プロシージャを用いた方法
- 17 -
環境シミュレーション
下図が Function プロシージャ(ユーザ定義関
数)という VBA プログラムを用いて,モンテカル
ロ法で円周率 π を求めてみたワークシート(3)であ
る.このシートでは,異なる回数 N に対して,
Function プロシージャによって π の近似値を求め,
N と近似値の関係をグラフ化している.ちなみに図
の中央付近にある(E10)の値は,PI 関数(Excel 関数
「=PI()」)で求めた π の値である. N の値が大き
くなるに連れて,4n N の値が π の理論値に近づい
ている様子がわかる.
試行回数 N )を渡し,π の近似値( 4n N )を演算結果
として返させようとしているのである.ただし,こ
こでの pi_cal という関数は,Excel に初めから定義
(用意)されていた関数ではない.これはユーザ定義
関数(Function プロシージャ)と呼ばれ,このシ
ミュレーションのために定義・作成した関数である.
この関数(Function プロシージャ)は,Sub プ
ロシージャと同じく,モジュールシート上で定義さ
れる.1 つのモジュールシート上に複数の Sub プロ
シージャや Function プロシージャを定義すること
ができる.
モジュールシートに以下の文字列を一字一句間
違えずに打ち込んでいこう.
ではこのシートの作成方法を順に説明していこう.
Function pi_cal(n As Long)
Dim ni, i As Long
Dim x, y As Double
ni = 0
Randomize
'乱数ジェネレータを初期化(乱数系列を再設定)する数値
'演算ステートメント
For i = 1 To n
x = Rnd: y = Rnd '[0,1)の一様乱数を返す
If x ^ 2 + y ^ 2 <= 1 Then ni = ni + 1
Next i
pi_cal = ni / n * 4
End Function
以上の一連の文が VBA によって書かれたユーザ
定義関数 pi_cal である.
オブジェクトの挿入
シートの左上に挿入された図は,ペイント
(Windows にアクセサリとして付いてくる簡単なお
絵描きソフト)で描いたものである.Excel のシート
上でペイントの絵を新規に作成するためには,[挿
入][オブジェクト][ビットマップイメージ]を選ぶ.
挿入された図は,ダブルクリックすることによって
いつでも編集可能であり,図のサイズも,図の 4 辺
および 4 角にある小さな黒いボックスをドラッグす
ることで変更可能である(ペイントによる作画の方
法は省略する).
ユーザ定義関数のプログラミングの基本は次の
ように,Function...(...)で始まり,End Function で
終わるサンドウィッチ構文である.
Function sample(x,y)
...
sample=...
...
End Function
上の例は,名前が sample の Function プロシー
ジャで,引数として x と y を取り,返り値として
「sample=」以下の値を返す関数となっている.いっ
たんこの形で関数をモジュール内で定義すれば,あ
と は 通 常 の シ ー ト の セ ル に 数 式
「=sample(A1,1.22)」のように入力するだけである
(モジュール内の他のプロシージャからこの関数を呼び出
すことも可能である).
罫線
表中の縦横 2 本の罫線は,[セルの書式設定]の[罫
線]タブで設定したもの.
ユーザ定義関数
さて π の近似値が示されている最も肝心な
E4:E8 について説明しよう.サンプルプログラムの
これらのセルには,例えば E4 になら,
では 1 行 1 行見ていく.
=pi_cal(D4)
Function pi_cal(n As Long)
という数式が入っている.要するに,pi_cal という
関数に引数として,その数式左側のセルの値(D4:
ユーザ定義関数 pi_cal の 1 行目は pi_cal(n As
Long)から始まる.これは pi_cal が引数として変数
- 18 -
環境シミュレーション
n を取り,そのデータ型が長整数型(Long)である
と宣言したものである.ここで長整数型にした理由
は,pi_cal の引数(ここでは,ランダム試行回数 N )
が最大で 1,000,000 を取るのに対して,整数型
(Integer)の記憶領域サイズは 2 バイトで,値と
して-32,768~32,767 までしか使用できないからで
ある.一方,長整数型なら 4 バイトまで割り当てら
れており,-2,147,483,648~2,147,483,647 の整数を
表わすことが可能である.
( x, y ) の原点からの距離の 2 乗(x ^ 2 + y ^ 2)がも
し(If)1 より小さか,あるいは等しければ(<= 1)(半径
1 の 4 分の 1 円の中ならば),その場合(Then)は,変
数 ni に 1 を加えて(ni + 1),それを変数 ni に再び代
入する(つまり ni を 1 つ増やす),というものである.
これによって 4 分の 1 円の中に入る ( x, y ) の数をカ
ウントしている.
pi_cal = ni / n * 4
Dim ni, i As Long
最終的には pi_cal = ni / n * 4 として,ni(4 分円の
中に入る回数を n )を n(正方形の中に入る回数 N )
で割ったものに 4 を掛け(すなわち 4n N ),その結
Dim x, y As Double
プログラムの 2-3 行は変数型宣言文である.ここ
では ni と i を長整数型(Long)に,x と y を倍精
度浮動小数点数型(Double)と宣言している.ちな
みにここで変数 ni は,4 分の 1 円の中に入る回数 n
を表す.
果をこの関数の返り値(すなわち π の近似値)として
返す.
このようにユーザ定義関数では必ずその返り値
を「関数名=...」の形で定義しておかなければならな
い.
ni = 0
グラフの作成
ワークシート(3)のような表がいったん完成す
れば,後は,結果をグラフ化するだけである.作
成の手順は先ず,1)グラフ化するデータの範囲
4 行目の「ni = 0」は,変数 ni を 0 で初期化したも
のである.
(D3:E8)を指定し,[グラフウィザード]ボタン
をクリックする.2)グラフの種類から[標準][折れ
線]を選ぶ.3)後はウィザードの指示に従ってその
他のオプションを設定するだけである.ちなみに
シート(3)のグラフは x 軸を対数目盛とし,x
軸が y 軸と 3.141593 で交差するように指定してい
る.
Randomize
次 に Randomize と だ け 書 か れ た 一 行 は ,
Randomize ステートメントと呼ばれ,その下で使
われている乱数ジェネレータ(Rnd)を初期化(乱数
系列を再設定)するための数値演算ステートメント
である.
For i = 1 To n
...
Next i
続く上記の部分は,For...Next ステートメントで
ある.ここでは,For と Next に囲まれたステート
メント群を引数の n 回繰り返す.
x = Rnd: y = Rnd
次の「x = Rnd: y = Rnd」は先の例と全く同じであ
る.ここで Rnd 関数で一辺 1 の正方形内にランダ
ムな点 ( x, y ) を生成している.
If x ^ 2 + y ^ 2 <= 1 Then ni = ni + 1
続く If ステートメントも,やや簡略形ではあるが,
やっていることは先の例とそれほど変わらない.点
- 19 -
環境シミュレーション
RAND()関数の応用例1
RAND 関数自体は 0-1 の一様乱数を生成するが,
これを利用して正規分布データを生成することも
できる.下記がワークシート上で生成した例である.
セル上に
=SQRT(-2*LN(RAND()))*COS(2*PI()* RAND())
偏差 1 の正規分布データであるが,この関数を σ 倍
して, x を加えてやれば,平均 x ,標準偏差 σ の
正規分布データとなる.
ページ下のヒストグラムは[ツール][分析ツール]
から[ヒストグラム]を選んで描いた.なお[分析ツー
ル]で[乱数発生]を選べば,同様に正規分布の乱数や
一様乱数を発生することができる.
ただし分析ツールを使用するためには,Excel を
フルインストールしていること,アドインによる分
析ツール使用の設定が必要である.
と入力すればよい.
ちなみにこの関数で生成されるのは平均 0,標準
RAND()を使った正規分布データの発生
"=SQRT(-2*LN(RAND()))*COS(2*PI()*RAND())"
28
-0.14491
29
1.150287
30
-0.79687
31
1.601563
32
-0.60391
33
0.498597
1
-0.99881
2
0.257725
3
-0.62075
-3.5
0
34
-0.77858
4
1.393733
-2.5
2
35
1.093049
5
1.322116
-1.5
6
36
0.752229
6
-0.81586
-0.5
25
7
-0.93269
0.5
38
98
-2.09073
8
0.63843
1.5
21
99
-0.74471
9
0.349568
2.5
8
100
0.966562
10
0.261899
3.5
0
11
1.288906 次の級
12
-0.49166
13
0.318864
14
-0.17804
15
0.076804
16
0.146498
17
0.42808
18
-0.61577
19
1.642217
20
-1.00306
21
0.45726
22
-0.25134
23
-0.67187
24
1.191922
25
0.872832
26
-3.04242
27
-0.2946
データ区間
頻度
0
- 20 -
環境シミュレーション
RAND()関数の応用例2(乱数を使った確率計算)
理論解:
男女が同数の学校において,任意の8人の男女を
選んで,横一列に並ばせた.この時,ある男子生徒
の両脇が女子生徒となる確率は? 男を 0,女を1
とおく.
A
7
8
1
B
C
D
E
F
G
H
I
1
2
3
4
5
6
7
8
0
0
0
1
1
0
1
1
1
1
⋅ 6 ⋅ = 0.09375
8
8
=INT(RAND()+0.5)
9
2
1
1
0
1
1
0
1
1
10
3
1
0
1
1
1
1
0
1
105
98
0
0
1
0
0
1
0
0
106
99
1
0
0
1
1
1
0
1
107
100
0
1
1
1
1
0
0
1
Q
K
L
M
N
O
P
7
2
3
4
5
6
7
8
0
0
0
0
1
0
1
=SUM(K8:P8)
=IF(C8=0,IF(B8+D8=2,1,0),0)
9
0
1
0
0
1
0
2
10
1
0
0
0
0
1
2
105
0
0
0
0
0
0
0
106
0
0
0
0
0
1
1
107
0
0
0
0
0
0
0
84=SUM(Q8:Q107)
0.105=Q108/800
参考図書
「Excel で学ぶ理工系シミュレーション入門」CQ 出版社(2003)
- 21 -
環境シミュレーション
今回登場した Excel 関数
・COUNT 関数
・PI 関数
今回登場した VBA 関数
・Rnd 関数
課題
次のいずれかの課題を実施して次週までに提出
せよ.定期試験には,このうち 1 問を出題する.
1) Excel 関数を使って,サイコロを投げて出る目
をシミュレートせよ(ヒント:「=RAND()*6」
で 0 から 6 の一様乱数を生成できる).次に 3
つのサイコロを投げて出た目を a,b,c とする
とき i)積 abc が偶数である確率を求めよ(答え
7/8),ii)(a-b)(b-c)(c-a)=0 である確率を求めよ
(答え 4/9)(類 岐阜経大)
2) Excel 関数を使って,3/5 の確率で 1,2/5 の確
率で 0 となる事象をシミュレートせよ.次にセ
ルの書式設定で 1 なら赤い「●」,0 なら「○」
と表示するようにせよ.以上を踏まえて,次の
問題に答えよ.袋の中に赤球 6 個,白球 4 個が
入っている.この袋から任意に 2 球を取り出す
とき,その 2 球が同じ色である確率を求めよ.
(答え 7/15)
3) センター試験において英語と歴史の平均値と
標準偏差がそれぞれ 60±40 点,65±35 点のと
き,両科目の合計点が 180 点を超える人の割合
を求めよ.ただし両科目とも 100 点満点とす
る.(約 9%)
4) 以下の課題のプログラムを完成させて,ワーク
シート(グラフを含む)と Sub プロシージャ
をプリントアウトしたものを提出せよ.問題
y = sin x の 曲 線 と x 軸 で 囲 ま れ る 面 積
( 0 ≤ x ≤ π )を Sub プロシージャを使ったモ
ンテカルロ法でもとめよ.(答え 2)
5) 以下の課題のプログラムを完成させて,ワーク
シート(グラフを含む)と Sub プロシージャ
をプリントアウトしたものを提出せよ.問題
y = cos x の 曲 線 と x 軸 で 囲 ま れ る 面 積
( 0 ≤ x ≤ π 2 )を Sub プロシージャを使った
モンテカルロ法でもとめよ.(答え 2)
- 23 -
Fly UP