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 -