数値演算言語OCTAVEによるプログラム入門

競馬@Wikiへ戻る

OCTAVEとは? †

OCTAVEとは理工系の技術者/研究者に愛用されている数値演算用インタプリタ言語MATLABのクローンソフトです。しかもフリーウェアです。
元々UNIX上で開発されたものですが、現在はWINDOWS上でも動作して、世界中の大学生/研究者に愛用されています。
残念ながら、日本ではそれ程メジャーではありませんが、非常に直感的に扱える為、電卓代わりとして、また最初のコンピュータ･プログラミングを覚える為にもとても有用なソフトウェアだと思います。

何故OCTAVE? †

現在はコンピュータ言語としてはVISUAL BASICを初めとして、C言語、Delphi(Pascal言語の一種)が有名です。JRA-VANでも開発サポートしている言語は大体この3種類ですし、世の中のソフトウェアも大体その3つで作られています*1。当然MATLABやOCTAVEで作られた商用ソフトなんてありません。
では何故OCTAVEなのでしょうか?
それは言語としての開発思想が他の汎用プログラミング言語とまるっきり違うからです。
OCTAVEは数値演算言語と言う、かなり特殊な用途に用いられる為に設計されています。よって、数学的なライブラリが豊富で、またその手の特殊な計算に対して非常に高速に計算結果を返してくれるように設計されています。また、数学的な一連の手続きを思いついた時に、簡単に記述できてその計算結果を瞬時に確認する事が出来ます。反面、汎用プログラムのように、例えばビデオゲームを作りたいと言ったような要求には丸っきり応じられません。
一方、VISUAL BASICやC言語、Delphi等は何でも出来る反面、例えば高度で複雑な数学的計算のようなかなり特殊なプログラムを組もうと思っても、ライブラリが豊富ではなかったり、また、場合によってはライブラリそのものを自作しなければならない可能性もなきにしもあらず、です。また、かなり様式も厳密に決められているので、思いついたままにコードを記述して即実行と言うラフ･スタイルによるプログラミングにも向いていません。キチンと各言語仕様における構文に則った設計/記述じゃないと、すぐバグ報告を返してきます。反面、ゲームを作ったりするのには向いてるんですね。そう言う画像系のライブラリは豊富でしょうから。
要するに餅は餅屋なんです。例えばモンテカルロ法なんかを簡易に行うなら、OCTAVEの方が得意。反面、競馬予想ソフトを作りたいのなら汎用言語がイイ、と。それだけなんですね。
ただし、VB/VBA等のインタプリタ言語と、OCTAVEは基本的なプログラムの作り方に関する発想は似ています。そして、仮に競馬予想ソフトで数学的ロジックを取り入れたい場合でも、最初にOCTAVEで計算過程を見てみるのもいいでしょう。それからVBに移植しても遅くないと思います。特にOCTAVEはオープンソースなんで、内臓のライブラリにたくさんの統計関数等が仕込まれているので、それを参考にして自分で予想ソフト用の統計解析実行エンジンを組み立てたりもできると思います。

OCTAVEはどこでダウンロードできるの? †

コレがインストーラです。←コレは山梨大学 工学部 電気電子システム工学科のサーバーに置いてあるOCTAVEの一種、OCTAVE Workshopをダウンロード/インストールする為の日本に於けるミラーサイトです。

OCTAVE Workshopって何? †

OCTAVE Workshopとは、最新ヴァージョンのOCTAVEにGUI*2環境を組み込んで、また、色々なOCTAVEの拡張ソフト*3をも組み込んだヴァージョンです。
何でこんなOCTAVE別ヴァージョンが存在するのでしょうか?
と言うのも、先に説明したとおり、元々OCTAVEと言うのはUNIX上で開発されたソフトなんですね。つまりオリジナルは文字だらけのプログラム言語なんです。
仮に、WINDOWSで起動すると、下の画面のようなDOS プロンプトで起動します。

こんな画面は怖すぎます。黒字に白なんてのは見ていられません(泣)。

こんなホラーな画面を見てて平気なのは、

NEC PC 6001かSHARP MZ-80の経験者だけです(怒)。

普通の感覚だったら、

天地鳴動して雷鳴轟き、嫁は逃げ出し、娘はギャルサー加入、一家離散で宇多田ヒカルは太って豚田ヒカルと化し*4ます。


こんな悪魔的かつダミアン的でオーメン的な悲惨な状況を迎えるくらいでしたら、素直にWindows世代としてはOCTAVE Workshopを手に入れるべきだと思います。
ちなみに画面写真は下のようになります。

こっちの方が

Windows世代らしい

ですね。一家離散の危機も回避できそうですし、これだったら

キャバ嬢に疎まれるウンチクオジサン化の危機も脱出できそうです。

やはり、人間、若者の話にもついて行けるような

若い感性を持っていたいものですね。

いや、ホントにそのテの話はどーでもイイんですが(笑)。*5

OCTAVEを電卓として使う †

ではキャバ嬢にモテモテオジサンを目指す為に簡単なOCTAVE Workshopの使い方を見てみましょう。
OCTAVE Workshopをインストールした後、デスクトップにアイコンが出来る筈なんですが、それをダブルクリックして起動すると、次のようなGUI環境が立ち上がります。
まずはその画面構成を覚えましょう。

OCTAVE Workshopは基本的には次の3つの画面で構成されています。

1. 右側全面:プロンプト(コマンド入力画面)。ここに様々な計算命令を入力します。
2. 左上側:ディレクトリ表示。ここにOCTAVE Workshop内の様々なフォルダ(ディレクトリ)/ファイルが表示されています。元々UNIX用として設計されているので、アイコン表示ではなく、こんな味も素っ気もないような文字表示になっています。が、しかし、文字表示でもキチンとマウスのクリックには反応します。
3. 左下側:何だかいまだよく分かりません(笑)。

つまり、右側のプロンプト(＞)以降に入力さえ行えば、最低限のOCTAVE Workshopの使用法はマスターした事になります。
では、お約束ですが、基本的な計算方法を記述していきましょう。

OCTAVEの基本的な使用方法 †

早速バカバカしい例ですが、

を次のようにOCTAVE Workshopのプロンプト以降に入力してみましょう。

そしてリターンキーを押します。

このように電卓として充分機能します*6。 一つの計算が終わればまた新しいプロンプト(＞)がどんどん下に生成されていくので、色々な計算式を入力して試してみてください。

例えば括弧()を利用して次のような計算も一気に行えます。

関数電卓としてのOCTAVE †

もうちょっと複雑な計算、例えば

なんかはsqrt*7と言ったコマンドで簡単に計算できます。

また、

のような複雑な計算も次のようにして簡単に計算出来ます。

OCTAVE Workshopには次のような代表的な数学関数が用意されています。

代入/付値 †

ではいよいよプログラミング言語としてのOCTAVE Workshopの特徴である初歩的な計算法に付いて解説します。
単なる電卓とOCTAVE Workshopの大きな違いは、文字に特定の数値を代入できると言った部分です。
これは専門的には付値、またはプログラミング言語としては定義と呼ばれる部分です。これが出来るか出来ないかが、単なる電卓とプログラミング言語の違いなのです。
実例を見てみましょう。次のように入力してみましょう。

そしてリターンキーを押しますと次のようになる筈です。

と入力したのと全く同じ式が返されてきました。しかしこれは今後の作業を司る大事な一歩なんです。
今行ったのは、OCTAVE Workshopに、xと言う文字を1として認識させたと言った作業に該当するんです。
これは数学的にはまさしくxと言う文字に1を代入したと言った意味ですし、また、コンピュータ･プログラミングで言うとxと言う変数に1を格納したまたはxと言う文字を1と定義したと言う意味なんです。
ちょっと分かり辛いかもしれないでしょうから、今度はOCTAVE Workshopに次のように打ち込んでみましょう。

今度はyを2として定義しています。
さて、xを1として定義、yを2として定義すれば、ではxとyの間の計算を色々行えばどんな結果になるんでしょうか?
これは算数で考えれば凄くバカバカしい愚問にも思えますが、一応OCTAVE Workshopで計算結果を見てみましょう。

計算の内訳はバカバカしいので解説しませんが、

全部x＝1、y＝2として計算されているのがお分かりになると思います。つまり確かにx、yにはそれぞれの数字が代入されているんです。
ところで、ここの代入/付値の項目を読んで、

だから何なの?そんな計算結果は算数と同じだし、代入なんてメンドくさいだけ。
これなら電卓と大して変わんないよ。

って思う人も多い筈です。そしてその指摘はその通りなんです。
はっきり言うと、この文字(変数)への数字の代入って作業だけではあまり意味はありません。
しかし、この作業が単なる電卓とプログラミング言語を分かつ最大の理由は、大まかに言うと、

｢コンピュータ･プログラム｣とは、こちらで設定した条件に従って変数に自動で
数値を代入する過程を記述する事を言う。
また、｢プログラムを実行する｣とはプログラムによる計算結果を表示する事を指す。

プログラミング作業とは原理的には上の事を指すんです。
単純に言うと、例えマリオでさえ、あるプログラムによって定義された変数の集合体であり、ゲームの流れも連続した計算結果の表示なのです。

その他の関数 †

他にもOCTAVE Workshopには、次のような便利な関数が定義されています。

上の4つの命令は特に競馬関係のプログラムには関係ありませんが複素数計算でよく使う命令です。
OCTAVE Workshopでは虚数単位をiかjで記述します。
例えば、

と定義して、次のような命令を実行してみましょう。

それぞれの計算結果が返されてきますね。
また、最後のdisp関数は、表示に関わる関数です。
例えば、良くプログラミング言語入門書の最初で、

Hello, World!

と表示させる例がありますが、それをOCTAVE Workshopでやってみましょう。

そしてリターンキーを押すと、

となり、確かにHello, World!と表示されます。
このように文字列を表示させたい場合は、"文字列"と挟みます。 また、自分で定義した変数(文字)や、次項のようなOCTAVE Workshopの内部で定義された特殊な定数の内訳を見る場合もdisp関数を用います。

特殊な定数 †

OCTAVE Workshopでは円周率πはpi、ネピア数eはeとして、内部で定義されています。
disp関数を用いてその内訳を覗いてみましょう。

行列の計算 †

元々OCTAVEとその元になったMATLABは、行列を用いた線形代数の数値計算を素早く行える為に設計された数値演算言語です*10。その為、行列の計算が大の得意で、かつ、その入力も他のソフトに見られるような煩わしさがありません。大変直感的に扱えるようになっています。
ここでは基本的な行列の計算法を紹介したいと思います。

行列の入力 †

例えば、以下のようなn×n行列

があった場合、OCTAVE Workshopではスペースとセミコロン(;)を利用して次のように入力します。

ちょっと実例をあげてみましょう。
例えば、次のような行列aを入力してみたいとします。

これをOCTAVE Workshopでは次のように入力します。

そしてリターンキーを押すと、

となって、確かに行列として認識されています。
では列ベクトル

はどのようにして入力すればいいのでしょうか?

このように、セミコロン(;)を上手く利用して、3行1列の行列として捉えれば入力は成功します。
同様に、セミコロン(;)を使わずにスペースを使って入力すれば以下のような

行ベクトルとして認識されます。

逆行列の計算 †

では行列

の逆行列mはどのようにして計算すればいいのでしょうか?
OCTAVE Workshopの逆行列を計算するコマンドはinv()*11です。
では、先程定義した行列aの逆行列を実際求めてみましょう。

従って、aの逆行列mは、

である事が分かります。

投資金配分の問題 †

詳しくはクラメールの公式のページに譲りますが、オッズa、b、c、の3点の馬券があって、投資予算がS円、どの馬券が当たっても払戻金を同じにしたいとすると、それぞれの馬券への投資金は行列を利用して

と定義され、これらの行列の関係は

で書き表されると言うものでした。
また、この式は必ず解を持ち、それは逆行列を用いて

となる、と言うものです。
さて、そこで次の問題を考えてみましょう。

(問)今、手元に資金1万円がある。この元手で3.5倍、6.4倍、8.3倍の馬券を3点買って、
どの馬券が当たっても払戻額を全て同じにしたい。
それぞれの馬券を一体いくらづつ買えば良いのか?
ただし、馬券の最低購入額は100円なので、答えは100円未満四捨五入で近似値でも良し
とする。

そうすると、問題設定としては次の3つの行列が問題となるわけです。

では実際OCTAVE Workshopを用いて各馬券に対する投資金配分をしてみましょう。
まずは行列Aと列ベクトルBを、OCTAVE Workshopに入力します。

次にAの逆行列Mを計算します。

最後に行列Mと列ベクトルBの積を求めます。

これが列ベクトルXとなり、解になります。
すなわち、

1. オッズ3.5の馬券には5,100円投資する
2. オッズ6.4の馬券には2,800円投資する
3. オッズ8.3の馬券には2,100円投資する

このようにOCTAVE Workshopを用いれば、至極簡単に投資金配分をする事が出来ます。

OCTAVEで初級プログラミング講座 †

ではいよいよOCTAVE Workshopでプログラムを組んでみましょう。
通常この手の数値演算言語でのプログラムは関数と呼ばれ、これは例えばMicrosoft Excelなんかで言うと、VBAを用いたマクロにあたるかもしれません。
要するに OCTAVE Workshop本体から関数名で呼び出し、様々な処理を行うようになっているのです。
いずれにせよ、OCTAVE WorkshopはVISUAL BASICと同じインタプリタ言語であり、基本はほぼ同じです。それどころかBASICよりは遥かに簡単です。よって、汎用プログラミング言語を覚えたいお方でも、最初にOCTAVE Workshopによるインタプリタ言語の大まかな記述法になれておけば、他の言語(例えばBASICとかDelphiとか)に移ってもそんなに戸惑う事はないとは思います。

OCTAVE Workshopのテキストエディタ †

では、実際にOCTAVE Workshopでプログラムを書き始めようと思うのですが、その前にテキストエディタに付いてちょっとお話しましょう。
テキストエディタとは、Windowsで言うと、OSの付属品であるNotePad?(メモ帳)やWordPad?の類のモノです。

｢ああ、あのワープロのショボいヤツか｣

そうです(笑)。そのワープロのショボいヤツです(笑)。
しかしながら、これはある意味強力なワープロのショボいヤツで(笑)*12、恐ろしい事に殆どのファイルやプログラムを開いてしまうのも事実です。Microsoft WORDや一太郎なら開けないようなファイルでもお構いなしだったりします。一体何故なんでしょう?
実はテキストエディタと言うのは｢開く為にある｣のではなくって｢記述する為にある｣からこそ強力なんです。何の記述かって?何でもです。例えばコンピュータ･プログラムやHTML言語を。それらをフォーマット通りに書き込む為に存在しているのですね。
例えばホームページ(ブログではない)を持ってる人ならご存知でしょうが、どんなHPをデザインするにせよ、基本は全てテキスト(文章)です。これら、例えばメモ帳で記述したHTML形式の文書をサーバー上にアップロードすればキレイなHPで自分で書いた文書が普通にIEで読めるようになります。また、CGI等のカスタマイズもこのテキストエディタで修正してアップロードします。僕がなかやさんのトコにしている投稿も基本的には全てこのテキストエディタで製作しています。ワープロのショボいヤツは意外に使えるのです。
同様に、適切なインタプリタ、またはコンパイラ*13さえあれば、テキストエディタは強力なプログラム言語記述用ツールに変身します。適切なファイル方式でセーヴして、インタプリタ/コンパイラで呼び出せばプログラムは走っちゃうのです。これでマリオカートも思いのまま作れちゃうんですね(笑)。
つまり、VBでのプログラムだろうと、はたまたC言語のプログラムだろうと、OCTAVE Workshop上のプログラムであろうと、基本的にはWindows付属のNotePad?(メモ帳)やWordPad?で記述できちゃうんです。素晴らしいですね。
しかしながら、何で各商用パッケージに専用のエディタがあるのかと言うと･･････やはり、使いやすさと言うか｢××専用｣に特化した方がユーザが使うのにラクだろうから、と言うのがあるんですね。色々な標準的じゃない機能を搭載できますし、しかも作ったプログラムをセーヴするのに一々エクステンション*14を弄らなくっていい。それが理由なんです。
さて、オリジナルのOCTAVEはUNIXで生まれた為かこの専用のエディタと言うのが存在しません。と言うのもUNIXは元々文字だらけの世界なんで、エディタはそこらに普通にあるのが当たり前なのです。従って、各プログラム言語専用のエディタが別々に存在する事は殆どなく、代わりに汎用エディタが幅を利かせています。
一方それがWindowsの世界の住人には至極不安なわけです。確かにNotePad?(メモ帳)やWordPad?でプログラムを記述できるのは分かった、でもイヤなんですよね(笑)。大体そもそも一々エクステンション変更ってのがUNIXと違い気に食わないわけです。
そこで。実はオリジナルのOCTAVEと違い、OCTAVE Workshopには専用のテキストエディタが付属で付いてきます。エライですね(笑)。
では、OCTAVE Workshop専用テキストエディタを開いてみましょう。
それはココ↓にあります。

OCTAVE Workshopの左側上方にFileと言うプルダウンメニューがあります。それを開けばNew、Openと言う二つの項目があります。Newは新しいファイルを作る為、Openは以前作った、または既存のファイルを開く為存在しています。
取りあえずここではNewを選択してみましょうか。

な〜んの特徴もない真っ白なEditorが現れたと思います。これがOCTAVE Workshop専用テキストエディタです。
これからここにプログラムを記述していくのです。
そして用語としては、OCTAVE Workshopの為に記述されたプログラムが書かれたテキストファイルを特にm.fileと呼びます。これはOCTAVE Workshop用のプログラムファイルのエクステンションが.mだから、という事に由来します。
つまり、Editorでm.fileを作れば、それはOCTAVE Workshop でプログラムを書いた、と言うのと同義なわけです。*15

最初のプログラムと関数定義 †

ではいよいよOCTAVE Workshopを使って最初のプログラムを組んでみましょう。その最初のプログラムは思い切って前出の計算法を利用して6点買いの資金配分のプログラムを作ってみたいと思います。ちょっと難しく思うかもしれませんが、これなら既に計算ロジックは組んでいるので、比較的簡単に作れると思うからです。
さて、先程の単なる計算と今から組むプログラムの違いとは何でしょうか?
それは

単なる計算の場合は計算過程をいちいち書かなければならないが、プログラムの場合は
一旦計算過程さえ組んでしまえば、あとはオッズさえ入力すればユーザーが計算の中身
を見る必要性はない

と言った部分です。つまり、重要なのは適切なオッズさえあれば、中のルーチンに従って、資金配分を勝手に計算してくれれば成功、と言う事です。
そこで最初に、先程書いた資金配分の計算過程をOCTAVE Workshopのエディタ上にそのまま書き写してみましょう。

たった4行で済むんですね(笑)。
さて、これを改造して行こうと思うんですが、問題点が2つあります。

1. 元々の計算は3点買いの為のものなんで、これを6点買いに改造しなければならない。
2. オッズが3.5、6.4、8.3と言う固定された数値だったのを、オッズ入力によって別の答えを返す汎用プログラムにしたい。

意識するのはこの2点です。
1番目の問題は、これはプログラムの問題と言うより数学の問題です。詳しくはクラメールの公式に譲りますが、自分で連立方程式を組むなり、もしくは勘のいい人だったら3点買いの行列設定を見ただけで帰納的に6点買いの行列構成を予測できると思います。
一応答えを書いておきますが、6点買いの行列A、Bの設定は、

となります*16
さて、上の行列A、Bは数学の部分なので、a、b、c、d、f、g、Sと言った7つの文字によって一般的な6点買いの資金配分の為の表現になっているのですが、要するにこの7つの文字をプログラム内部でも具体的な数字を与えず、外部から入力を促すような設計にしたら当初の目的を達する事が出来るわけです。
このような、プログラムに入力してもらうa、b、c、d、f、g、Sのような変数を引数*17と言います。
そして、OCTAVE Workshopでのプログラムフォーマットは関数と言って次の形で定義されるのです。

function 返り値 = 関数の名前(引数)
関数定義本体

これだけです。functionと言うのはこれはプログラム(関数)ですよと言った宣言です。返り値は、最終的に答えとしたい部分(この例で言うと資金配分の内訳ですね)の名前で、yでもzでも何でも構いません。関数の名前も自由に設定して構いません。ただ、OCTAVE Workshopで呼び出す時、この名前で呼び出すので、分かり易い名前にした方がいいと思います。そして引数を設定すればいいんですね。
ちょっとまだピンと来ないかもしれないので、ここで返り値をy、関数名をshikinhaibunとして設定してみましょう。必要とする引き数は7つあったので次のようになります。

function y = shikinhaibun(a,b,c,d,f,g,S)

後は、関数定義を今までの議論を省みれば、次のようにEditorに書き込めばいいのが分かるでしょうか?

これだけなんです(笑)。簡単でしょ?これが貴方がOCTAVE Workshopで書いたはじめてのプログラムです。
なお、最後にy＝M*Bとしているのは、ここのM*Bの値(要するに資金配分の解ですよね)が最終的に欲しい/知りたい答えなので、この部分が返り値なんです。そこで最初の設定とつじつまを合わせる為、y＝と付け足してるわけですね。

プログラムのセーブとディレクトリ †

では前項で書いたプログラム(m.file)を保存しましょう。
OCTAVE Workshopでは、決められた場所(ディレクトリ。Windowsではフォルダと呼ぶ。)にm.fileをsaveしておけば、好きな時に呼び出す事が出来ます。その一連の流れがプログラムを実行すると言う事なんですね。
ところで、OCTAVE Workshop専用Editor上で書いたプログラムはそのまま関数定義名で自動に名付けられて保存されますし、またエクステンションも自動でm.fileとして保存されますので、その辺りは心配しなくて構いません。問題はm.fileをどこへ保存するのか?だけです。これはOCTAVE Workshopの仕様で殆ど厳密に決められているのですね。
さて、OCTAVE Workshopの左上のスペースをちょっと見て下さい。

これがWindowsで言うところのOCTAVE Workshopフォルダ内のフォルダ構成なんです。個人ホームページ(ブログではありません)をお持ちのお方ならご存知でしょうが、これはUNIX系の表現で、いわゆる相対パスの表現になっているんですね。
この味も素っ気もない文字はクリックするとキチンと開く事が出来ます。そしてフォルダの階層構造が把握する事が出来ると思います。
それで本題に戻りますと、m.fileをsaveするべき場所は、OCTAVE Workshop内のフォルダ構成から始まって、

msys/⇒octave2.9.4/⇒share/⇒octave/⇒2.9.4/⇒m/⇒プログラム保管場所

となります。6階層先になるわけですが、この場所を良く覚えておいて下さい。
実際問題、m.fileをsaveするのは以下のようなEditorのFileプルダウンメニューでSave asを利用して行います。

これでsaveを行う際、どのフォルダにsaveするのか尋ねるポップアップが出てきますので、上で説明したとおり、6階層先のmフォルダ内に作成したプログラム(m.file)を保存して下さい。
これで準備は完了です。

プログラムの呼び出し †

では作成/保存したプログラムshikinhaibunを実際使ってみましょう。
仮に、予算一万円で、オッズ7.4倍、8.2倍、11.3倍、14.2倍、15.2倍、17.1倍の馬券の計6点買いをするとしましょう。そこでOCTAVE Workshopの本体のプロンプトに次のように入力します。

shikinhaibun(7.4, 8.2, 11.3, 14.2, 15.2, 17.1, 10000)

これが作ったプログラムを実行する命令です。
そして実行すると、次のような計算過程がズラズラと出てくるはずです。

確かにプログラムは上手く動いたようです。
そして資金配分は四捨五入すると、

7.4倍の馬券･･････2,500円
8.2倍の馬券･･････2,300円
11.3倍の馬券･･････1,600円
14.2倍の馬券･･････1,300円
15.2倍の馬券･･････1,200円
17,1倍の馬券･･････1,000円

となる事が分かります。

プログラムの改良 †

もうちょっとプログラムの手直しをしてみましょうか?
確かにプログラムshikinhaibunは上手く動いてますが、毎回毎回計算過程が表示されるのも、初めはカッコイイんですが、目がチカチカする気もします。
そこで、計算過程をユーザー側から見えないようにしてみましょう。
OCTAVE WorkshopのFileプルダウンメニューからOpenを選んで、先程saveしたshikinhaibunプログラムを6階層先のm/フォルダから開きます。そして、次のように修正してみます。

｢え?どこを変えたの?｣

って思うかもしれません。が良く見て下さい。行列A、列ベクトルB、逆行列M、解M*Bの内訳のうしろ全てにセミコロン(;)が付け加えられているのに注目して下さい。
OCTAVE Workshopではこのように、数式のうしろにセミコロン(;)を付け加える事によって、プログラム実行時に計算過程を表示させない、と言ったような約束事があるのです。
取りあえず上のようにプログラムshikinhaibunを書き換えてSaveして下さい。そして、今度はオッズ9.2倍、9.8倍、12.6倍、13.5倍、17.7倍、18.2倍の6点の馬券を1万円で購入する資金配分をOCTAVE Workshopで計算させてみましょう。

ご覧のように、今度は計算過程は表示されずに単に計算結果だけ返してくれます。
この場合の資金配分は四捨五入すれば

9.2倍の馬券･･････2,300円
9.8倍の馬券･･････2,100円
12.6倍の馬券･･････1,700円
13.5倍の馬券･･････1,600円
17.7倍の馬券･･････1,200円
18.2倍の馬券･･････1,200円

である事が分かります。
今度はもうちょっと改良してどれか馬券が当たった場合の期待回収値も表示できるようにしてみましょうか?
詳しくはクラメールの公式のページに譲りますが、例えばオッズa、b、c、d、f、gの6点買いをした場合、合成オッズと言われる公式で、期待回収率αは次のように計算されます。

つまりこの公式で計算された値を表示させれば済むわけです。
そこでshikinhaibunプログラムを次のように書き換えてみましょう。

最後の返り値の指定y=M*Bの前に一行kaishuurituの式を加えただけです。しかもこれはセミコロン(;)を行末に加えていないので、これだけは計算過程として計算表示されるわけです*18。
そこで、先程行った、オッズ9.2倍、9.8倍、12.6倍、13.5倍、17.7倍、18.2倍の6点の馬券を1万円で購入する資金配分をもう一度OCTAVE Workshopで計算させてみましょう。

この1行を加えただけで、この6点買いの期待回収値も一目で分かるわけです。この例ですと、2.1倍の馬券を1点買いする回収率と変わらない結果だという事が分かりますね。
OCTAVE Workshopではこんな感じで資金配分の為のプログラムも簡単に組む事が出来るのです。10点買い専用プログラムとか、予算といつも買う点数に合わせてプログラムを色々と組んでいただけたらな、と思います。

プログラムを組む為のテクニック †

前項で資金配分の為のプログラムを組みましたが、OCTAVE Workshopでは結構簡単に組めるので驚いた方々も多いと思います。特に汎用言語を使った経験がある方々なら余計に、だと思います。
最大の違いは、VB/VBA、もしくはC言語等と言う汎用言語では変数(前項の例では引数ですが)の型を厳密に定義しなければいけない事、また、前項では行列を利用したプログラムだったのですが、この行列の為のライブラリが存在しなければ、たったの4行前後ではとてもじゃないけどプログラムが組めない、と言う部分でです。これも数値演算用の特殊な言語であるOCTAVE Workshopならではの簡易さ、だと思います。
しかしながら、確かに前項の例は簡単過ぎますし、他のプログラムと同様、もうちょっと複雑な計算過程を扱うとなると、それなりのプログラムの考え方に慣れないといけません。
つまり、制御命令と言われるテクニックに付いて、です。

条件分岐:if, else †

ある条件で場合分けをして処理を行いたい場合はif文、else文を使います。
一番簡単なif命令の形は

if (条件文)
処理
endif

と言うようなものです。基本的には処理をif〜endifで挟むと言うのが構文になります。
この命令を実行すると、｢条件文｣が真である時だけ、｢処理｣が実行されます。
｢条件文｣が偽である時実行する処理を指定するには次のようにします。

if (条件文1)
処理1
else
処理2
endif

例として、入力xが負なら正にして出力、xが負でなければそのまま出力、つまり絶対値を計算する関数myabsを作ってみましょう。

function y = myabs(x)
     if (x < 0) 
         y = -x;
     else
         y = x;
     endif
     y;

では実際プログラムmyabsをOCTAVE Workshopで使ってみましょう。

確かに絶対値が計算されているのが分かると思います。
条件文が2つある(処理自体は3つ)場合は、次のように記述します。

if (条件文1)
　　処理1
elseif (条件文2)
　　処理2
else
　　処理3
endif

同様に、条件文が3つある(処理自体は4つ)ある場合は

if (条件文1)
　　処理1
elseif (条件文2)
　　処理2
elseif (条件文3)
　　処理3
else
　　処理4
endif

となり、このようにelseifを使えばいくつでも条件を増やす事が出来ます。
なお、条件文設定の為の演算子には次のようなものがあります。

また、複数の条件文を組み合わせる複合条件文の為の論理演算子には以下のようなものがあります。

繰り返し文:for †

ある処理を繰り返し行いたい場合には、同じ文を何回も書く代わりにfor文を用います。いわゆるループと言われる処理ですね。
for命令は次のような構文です。

for ループ変数 = リスト
    繰り返し処理
endfor

このように繰り返しさせたい処理をfor〜endforで挟みます。
例として、0が入ってる変数aに5回だけ1を足すプログラムmyloopを作ってみましょう。

function y = myloop()
 a=0;
  for i = 1:5
   a = a + 1;
  endfor
 y = a;       

これは引数を特に定義しないので、関数名のうしろの中身をブランクにしています。
また、1から5までと言うのはコロン(:)を使って1:5と記述します。これをリストと呼びます。
プログラムmyloopを実行すると、次のような結果になります。

計算結果は5となります。
これは当たり前と言えば当たり前で、0から初めて1づつ5回足すと5になるんですね。
上のプログラムは次のような流れになっているんです。

1. aに0を代入する。
2. ループ開始
3. 1回目でaに1を足す。これを新しくaとする。a＝1。
4. 2回目でaに1を足す。これを新しくaとする。a＝2。
5. 3回目でaに1を足す。これを新しくaとする。a＝3。
6. 4回目でaに1を足す。これを新しくaとする。a＝4。
7. 5回目でaに1を足す。これを新しくaとする。a＝5。
8. ループ終了。
9. 最終的なaの値をyに代入して値を返す。

さて、リストの1:5と言うのは実は何でも良くって、例えば3:7でも計算結果は同じになります。

function y = myloop2()
 a=0;
  for i = 3:7
   a = a + 1;
  endfor
 y = a; 

と言うのも、あくまで｢何回繰り返すのか?｣が大事であって、特に初期値がいくつで最終値がいくつであるのか、は関係無いからなんです。

ご覧の通り、プログラムmyloop2もmyloopと同じ計算結果を返してくれます。
当然、ループの最後の数を引数として与えて、任意に何回ループを繰り返すのか指定する事も出来ます。
例として、整数値nを入れると1からnまでの総和を返してくれるプログラムmysumを作ってみましょう。

function y = mysum(n)
 i = 0;
  for j = 1:n
   i = i + j;
  endfor
 y = i;

モンテカルロ･シミュレーション †

原則的にはif文とfor文の2つを使いこなせれば OCTAVE Workshopでのプログラムだけではなく、多言語でも殆ど全てのプログラムを書く事が出来ます。
他にも制御命令の構文は存在しますが*19、原理的にif文とfor文のちょっとしたヴァリエーションだったりします。
興味のある方はOCTAVE WorkshopのHelpのプルダウンメニューからHelp Contents⇒GNU Octave Manualへと進んで、12章辺りにOCTAVE Workshopで利用できる構文の説明がありますので、それを読んだら宜しいでしょう(もっともマニュアルは英語ですが)。

さて、ここからこのプログラム講座の最大の目玉、乱数を使ったモンテカルロ･シミュレーションのやり方を学んでいきましょう。

モンテカルロ･シミュレーションとは? †

突然ですが、現在僕らがインターネットに接続しているこれらパソコンは正確にはノイマン式計算機(コンピュータ)と呼びます。このノイマン式と言う冠ですが、これはこの形式のコンピュータの数学的ロジックを開発したユダヤ系の数学者、フォン･ノイマン博士から取られたものです。
この人は一種天才で、コンピュータの数学的基本理論の他、いわゆるゲーム理論でも名前が知られています。
ところで、元々コンピュータ自体も弾道計算をする為にアメリカの国家予算を使って作られたものであり、そう言った軍事機器を日常的に使用できる環境にいる、と言うのも何かヘンな気分がします。
さて、このコンピュータを利用した計算方法で非常に著名なのが、ここで紹介するモンテカルロ法と言うものです。この原理は意外と古くから理論面では知られていましたが、実際この方法をコンピュータで実現したのが、フォン･ノイマン博士を初めとする科学者のチームだったのです。そして、何の為の研究だったのか、と言うと、あの忌まわしい原爆製造の為のプロジェクト、マンハッタン計画での研究だったのです。 単純に言うと、原爆を製造する際、原子の振る舞いが予測できないといけません。ところが、原子は僕らの世界の視点とは違い、量子力学と言う力学の一分野に則って運動するとされています。そして量子力学とは平たく言うと、確率を組み入れた運動の記述なんです。
ところで、量子力学自体は一番簡単な水素原子を元にしてその理論基盤を発展させてきたのですが、原爆で使用する原子はウラニウムと言う水素原子より遥かに大きな原子です。そうすると、理論的な数式を用いた直接計算より遥かに計算過程が難しくなるのは何となく分かると思います。
そこで、フォン･ノイマンを初めとする科学者達が考えた方法論が、コンピュータ上で乱数を用いて、実際原子がどんな振る舞いをするのかシミュレーションしてみて、それを計算結果とする方法論です。このような乱数を用いたシミュレーション法をカジノで有名な都市にあやかってモンテカルロ･シミュレーションと呼ぶのです。
この手法はすさまじく強力で、現在では、軍事利用から転用されて様々な数学上/工学系の問題を解く場合にも用いられています。そして何が一番強力かと言うと、適切な状況設定さえしてやれば、数学や数式が苦手でもシミュレーションで暫定的な解が簡単に求まってしまう事なのです。
そして、第2次世界大戦中だったら一旦シミュレーションを始めたらそれこそ何日間も大型コンピュータを走らせ続けなければならなかったのですが、21世紀を迎えた現在、僕らは当時の大型コンピュータより遥かに高速のコンピュータを自宅で走らせているのです。つまり、僕らは自宅で簡単にモンテカルロ･シミュレーションが出来る大変ありがたい時代に住んでいるのですね。
軍事技術の競馬への転用をあまり面白く思わない向きもあるかもしれません。しかし、そもそもこのようなインターネット自体が軍事技術の転用ですし、ましてや競馬の肝、馬の品種改良でさえ元々軍事目的なんです。世界は軍事技術で溢れているのです。
ってなワケで、自宅でモンテカルロ･シミュレーションを平和利用してみませんか(笑)?

一様乱数 †

まずはモンテカルロ法のプログラムを組む前に、OCTAVE Workshopではどんな乱数を利用できるのか、そしてその乱数はどんな性質があるのか、ある程度知っていなければなりません。
原理的に、乱数とは、特定の確率分布で表される母集団から、無作為に抽出された数字の事を指します。つまり、確率分布により、得られる乱数の性質が違うんですね。
その中でもっとも利用頻度が高い乱数が一様乱数と言われるものです。通常単に乱数と言えば、この一様乱数を指すケースが多く、またもっとも基本的な乱数です。

この区間a〜bはどんな区間を選んでも構いません。通常は0〜1の範囲で設定されます。または0〜100も良く使われるでしょう。
なお、引数の行数、列数と言うのは、OCTAVE Workshopでは一様乱数を行列で返すように設計しているので、その行列の大きさを指定する為のものです。何も指定しなければ、1個の乱数を返してくれます。

上の例では0〜1の範囲での乱数を1個返しています。次の例では0〜1範囲での乱数を2×2行列で返しています。

コイン投げ †

さて、ここで、コイン投げのロジックを考えてみましょう。
通常、コインを投げて表が出る確率は50％、裏が出る確率も50％だと言われています。
そこで、これを一様乱数を使って表現してみましょう。考え方は次の通りです。

1. unifrnd(0,1)の結果が0〜1/2だったら表とする。
2. unifrnd(0,1)の結果が1/2〜1だったら裏とする。

そこでこれを利用して以下の様なコイン投げのプログラムcoinを書いてみましょう。

function y = coin()
 x = unifrnd(0, 1);
  if (x <= 1/2) 
   men = 1;
  else
   men = 0;
  endif
 y = men;

じゃんけんのプログラム †

同じ考え方を用いて、じゃんけんのプログラムを一様乱数を用いて書いてみましょう。
じゃんけんでグー/チョキ/パーを出す確率も等確率と考えられるので、

1. unifrnd(0,1)の結果が0〜1/3だったらグーとする。
2. unifrnd(0,1)の結果が1/3〜2/3だったらチョキとする。
3. unifrnd(0,1)の結果が2/3〜1だったらパーとする。

と取りあえず置いてみましょうか。
この考え方を用いたプログラムjankenは次のようになる筈です。

function y = janken()
 x = unifrnd(0, 1);
  if (x <= 1/3) 
   te = "gu";
  elseif (x <= 2/3)
   te = "choki";
  else
   te = "par";
  endif
 y = te;

簡単なシミュレーション †

上の二つのプログラムは一様乱数の使い方に慣れる為の練習でした。
一旦if構文と一様乱数の組み合わせ方に慣れれば、すぐに簡単なシミュレーションの為のプログラムは書けるようになります。
まずは1回のシミュレーション･プログラムを書く為の雛形を提示しましょう。

function y = mysimulation()
 count = 0; %  カウンタを0に戻す
            %%%
            %%%シミュレーション手順の記述
            %%%
 y = count; %  結果(カウンタの値)を出力

ではこの雛形を利用して、1枚のコインをn回投げた時に表が何回出るかをシミュレーションする為のプログラムmycointossを書いてみましょう。

function y = mycointoss(n)
 count = 0;
  for i = 1:n 
   x = unifrnd(0,1)
    if (x <= 1/2)
     coin = 1;
    else
     coin = 0;
    endif
    if (count == 1)
     count = count + 1;
    elseif
     count = count;
    endif
  endfor
 y = count;

上のプログラムmycointossを用いて10回コイン投げを試行した結果が以下のようになります。

プログラムの最適化 †

OCTAVE WorkshopやVBAやC言語等のプログラミング言語の優秀な所は、自分が定義したプログラムの内部で別に書いたプログラムを呼び出して実行できる辺りです。つまり、一つのプログラムで全てのルーティンを全部書く必要性がないのですね*21。
この発想はプログラムをなるべく短く効率良く書くには必要な発想です。これをプログラムの最適化と呼びます。
例えば前項で書いたプログラムmycointossの次(4行目〜9行目)の部分を見てみましょう。

   x = unifrnd(0,1)
    if (x <= 1/2)
     coin = 1;
    else
     coin = 0;
    endif

この部分はコイン投げで表が出るか裏が出るか判定する場所です。
ところが、この部分はそれ以前に定義したプログラムcoinで充分定義されている部分なんです。別な言い方をすれば、coinはcoinでプログラムをsaveしてたら、この部分ではそのプログラムcoinを呼び出してしまえばいいのです。
以上を鑑みて、プログラムmycointossを次のように書き換えます。

function y = mycointoss(n)
 count = 0;
  for i = 1:n
   x = coin();
    if (coin == 1)
     count = count + 1;
    else
     count = count;
    endif
  endfor
y = count;

これでもプログラムcoinが別に定義されていれば何の過不足もなくプログラムが実行できるのです。
また、上の修正されたプログラムはまだ改良が可能です。
と言うのも、｢コインを投げて表ならばカウントする｣作業はコインが表(1)だろうと裏(0)だろうと、そのまま結果を変数countに足せばいいだけなので、実はif文を利用する必要さえないのです。

function y = mycointoss(n)
 count = 0;
  for i = 1:n
   count = count + coin();
  endfor
y = count;

最初書いたプログラムは16行だったのに、今は既に6行までに圧縮されました。
このように、なるべく効率的に/少ない行数でプログラムを書くように考えればいいのです。

最適な乱数を設定する †

今までは一様乱数を利用してきました。
ところで、別の乱数を使えば実はプログラムmycointossは2行のプログラムに圧縮できるのです。

function y = mycointoss(n)
 y = binornd(n, 1/2);

ここで使用した乱数は二項乱数と呼ばれるものです。

binorndもunifrndと同じように行列で乱数を返す機能を持っているのですが、行数、列数の引数を省略すれば二項分布に従う正の整数を一個返してくれます。
このようにシミュレーションによって適切な乱数を選択すれば異様に短くプログラムを書く事も可能なのです*23。

本格的なモンテカルロ･シミュレーション †

今まで扱って来たような乱数の使い方とそれを用いた単純なシミュレーションの方法論さえ分かれば、いよいよ本格的なモンテカルロ法を行う事が可能になります。
ここで、OCTAVE Workshopに於けるモンテカルロ法のプログラムの雛形を提示したいと思います。

function y = mymontecarlo(n)
 count = 0;    %   カウンタを0に戻す
  for i = 1:n  
               %%% シミュレーション手順を記述し、
               %%% (シミュレーションをn回行う)
               %%% 結果をcountに足し算する
  endfor
  y = count/n; %   結果を出力

例えば、プログラムcoinを利用すれば

function y = mymontecarlo(n)
 count = 0;
  for i = 1:n
   count = count + coin();
  endfor
 y = count/n;

となって、コインをn回投げた場合のコインの表が出る確率を計算してくれます。
これは数学的に言うと、次の相対頻度の極限値としての母比率pの定義、

を地に行ってる計算なんです。
試しに今定義したmymontecarloで試行回数n＝10,000でコインの表が出る確率を計算してみましょう。

かなり0.5に近い値が出ます。
なお、このプログラムは二項乱数を利用して

 function y = mymontecarlo(n)
 y = binornd(n,1/2)/n;

と書く事も出来ます。

平均値に関するモンテカルロ･シミュレーション †

例えばここで、10回コインを投げたら平均で何回表がでるのか?考えてみましょう。
仮に10回コインを投げた時、表が出る回数は当然揺らぎます。ですから、単純に｢×回である｣という事は出来ません。
そこで発想を変えてみて、10回のコイン投げを1セットとして、それをn回繰り返してみるシミュレーションを考えてみればいいのです。
コイン投げのプログラムmycointossを利用すれば、平均値をモンテカルロ法で計算するプログラムcoinmontecarloは以下のように記述できます。

function y = coinmontecarlo(n)
 count = 0;
  for i = 1:n
   count = count + mycointoss(10);
  endfor
 y = count/n;

このプログラムcoinmontecarloを利用して、1,000セット、10,000セット、100,000セットでそれぞれモンテカルロ法を行って計算した結果が以下のようになりました。

理論上の平均値である、5に近くなって行ってるのが分かると思います。

二項乱数を利用した平均値を求めるモンテカルロ法の最適化 †

上で書いたプログラムcoinmontecarloは二項乱数を利用すればもっと短く、たったの二行で書く事が出来ます。

function y = coinmontecarlo(n)
 y = mean(binornd(10,1/2,1,n));

このプログラムの短さの秘訣は次のコマンドを併用した事によります。

そして、二項乱数binorndは二項分布に従う正の整数を乱数行列として返してくるコマンドでした。

つまり、1行n列なり、m行1列で出力すれば、それは乱数を要素としたベクトルvを出力するのと同義なんです。
この2つの関数の機能を組み合わせれば、前出のプログラムも最適化した状態で記述する事が出来るのです。

モンティ･ホールの問題 †

さて、ここでHRPTV5Cさんが示唆していたモンティ･ホールの問題に付いて考えてみたいと思います。
モンティ･ホールの問題とは確率論の中で有名な問題で、直感的な確率に対する思い込みと理論上の計算が喰い違う、と言ったちょっとパラドキシカルに見える問題です。
問題は次のようなものです。

モンティ・ホールの問題とは、テレビのバラエティ番組に起源する。

モンティ･ホールが司会者を務めるクイズ番組において、回答者である番組参加者は、
３つの扉の一つを選ぶ。その背後には、当たりである商品とはずれを意味する、
たとえばアヒルのおもちゃが隠されている。
 
回答者が扉Aを選ぶと、モンティはそれを伏せたまま、第二の扉、たとえば扉Cを
開く。
モンティは当たりの扉、ハズレの扉を知っており、必ず回答者が選ばなかった扉の内、
ハズレの方の一つを選ぶ。
そして、モンティは回答者に、もう一度、選択し直すことを提案する。回答者は、
扉Bに選択し直すことが有利であるか。

これはなかなか難問です。直感的には扉Aを選択しようと扉Bを選択しようと結果は変わらないように思えます。果たしてどうなんでしょうか?
この問題を解くにもモンテカルロ法が大変有用です。ちょっと考え方が難しく思えるんですが、次のように考えれば異様に簡単に整理できます。

1. まず、プログラムとしては扉を変更しないケースを中心として考える。
2. 0〜3までの値を取る乱数を利用して、回答者が選ぶ扉を設定する。具体的には、floor関数(切捨て)を使って、0、1、2のうちの一つの扉を選ぶものとする。
3. もう一つ0〜3までの値を取る乱数を利用して、当たりである商品が隠されている扉を設定する。これもfloor関数(切捨て)を使って、0、1、2のうちの一つの扉に当たりが隠されている。
4. どの道回答者は最初に選択した扉から動かないのが前提なので、回答者が選んだ扉と正解の扉が一致する確率が選択しなおさない回答者がゲームに勝つ確率になる。モンティが回答者が選ばなかった扉のどちらを開こうが、回答者は動かないので、モンテカルロ･シミュレーション上、モンティの行動を定式化する必要性はない。
5. 1−回答者が初めに選択した扉から動かない確率が、選択し直す事で勝つ確率になる。

以上の考え方を踏まえて、モンティ･ホールの問題を解くプログラムmontyhallを作ってみましょう。

function y = montyhall(n)
 staywin = 0;
  for i = 1:n
   choice = floor(unifrnd(0,3));
   answer = floor(unifrnd(0,3));
    if (choice == answer)
     staywin = staywin + 1;
    endif
  endfor
staying = staywin/n
y = (1-staywin)/n;

さて、モンテカルロ法(試行数3,000回)による解答はどうなるでしょうか?

実は計算によると、常に扉を選択しなおした方が有利なんですね。
これはなかなか面白いシミュレーションだと思います。

モンテカルロ法の様々な応用 †

モンテカルロ法は前述のような確率に絡んだ計算シミュレーションだけでなく、積分の近似計算に於いても用いられます。
数学的に言う積分と言うのは、図形の面積を求めたり体積を求めたりする計算法を指すんですが、実の事を言うと、これは確率計算と切っても切れないテクニックです。
と言うのも、専門的には測度論と言うのですが、確率とは広義ではある面積(体積)の全体の面積(体積)に対する割合の事を指します。よって、確率計算で使われるモンテカルロ法の積分に関する転用は数学的には実はかなり本質的な作業なんです。
ここではモンテカルロ法の積分計算への転用例を見てみたいと思います。
なお、雛形としては丸っきり同じで、

function y = mymontecarlo(n)
 count = 0;   %   カウンタを0へ戻す
  for i = 1:n
              %%% 砂粒をn個敷き詰め
              %%% ある領域に含まれる砂粒の
              %%% 数をcountに加算する
  endfor
 y = count/n  %   結果を出力する

となります。中間の捉え方だけが少々違うんですね。
これは図示すると以下の様な意味です。

ここで、面積S＝S1+S2だと言うのはお分かりだと思います。そして、Sの面積自体は分かっているものとします。
仮にこの領域Sに満遍なく乱数で発生させた砂を敷き詰めると、S1の面積に敷き詰められた砂粒の個数をn1、全体の砂の数をNとして

となり、確率の問題と同様に面積の問題も解く事が可能なのが分かるでしょう。
あとは、乱数で砂を実際に作り出してみればいいだけなんです。

円周率πを求めてみる †

一般的なモンテカルロ法の導入として良く紹介されるのが、円周率πの近似計算に対する適用です。これを実際、OCTAVE Workshopを用いてモンテカルロ法で求めてみましょう。考え方は以下の通りです。
まずは下の図のような、半径1の円の4分の1を考えます。この4分の1円は、1×1の正方形の中にスッポリと納まっているのは分かりますね?

さて、ここでまず適当に座標を設定します。つまり、正方形を(x,y)＝(0, 0)、(0, 1)、(1, 0)、(1, 1)、と4つの座標で囲まれている領域と考えて、この中に砂をまきます。この砂も(x, y)と言った座標上に落されるのです。
ここで三平方の定理、

で、砂の落ちた点の原点からの距離は厳密に計算する事ができます。
そして、円の半径は既に1と分かっているので、

1. 原点から砂の落ちた位置までの距離rが1以下であればカウントする。
2. 原点から砂の落ちた位置までの距離rが1より大きければカウントしない。

と言ったたった二つのルールで4分の1円内に落ちた砂の数をカウントできるのです。
そして、円の面積は半径×半径×πであるのは分かっているので、結果この4分の1円の四倍が円の面積になり、かつそれが求めたい円周率πの値でもあるのです。
以上を鑑みて、円周率πを求めるモンテカルロ･シミュレーションのプログラムpimontecarloを書いてみましょう。

function y = pimontecarlo(n)
 count = 0;
  for i = 1:n
   sand = unifrnd(0,1,1,2);
    if (sqrt(sand(1,1)^2 + sand(1,2)^2) <= 1)
     count = count + 1;
    endif
  endfor
 y = 4*count/n

このプログラムの要所は次の事です。

4行目:
sand = unifrnd(0,1,1,2);%ベクトルとして一様乱数を2個出力している。
                        %そのベクトルは列ベクトルで1行2列の大きさ。
                        %それがsandと言った変数の定義

5行目:
sand(1,1)^2+sand(1,2)^2 %ベクトルsandから要素を取り出している。
                        %1個目の要素をx座標、2個目の要素をy座標として
                        %見たたている。
                        %一般にベクトルvからi行目j列目の要素を取り出す場合、
                        %v(i,j)と記述する。

この2つがポイントです。
では、このモンテカルロ･シミュレーションの1万回試行の結果を見てみましょう。

3.14ピッタリとはいきませんが、かなり近い値が出ているのが分かると思います。
もっと試行回数を増やせば、どんどん3.14の値に近づいて行くのです。

なお、蛇足ですが、原理的には4分の1円を利用しないで円そのものを使ったプログラムも可能です。が、一般に砂を敷き詰める元の面積が大きければ大きいほど、敷き詰める砂の量が多くなる＝試行回数の増大を意味するので、ある限られた試行回数ではこの手のモンテカルロ法による近似積分計算の精度は悪くなります。
よって、なるべく必要最小限のコンパクトな面積を利用して近似計算のロジックを作った方がBetterだという事を申し添えておきます。

体積を求めてみよう�� †

例えば次のような図形の体積を考えてみます。

ここでxは0〜πの範囲、yは0〜2の範囲として、曲線はz＝sin xだとしましょう。
では、ここで曲線z＝sin x、0≦x≦π、0≦y≦2の範囲で囲まれた体積を求めよ、と言った問題を考えるとこれがなかなか難問なのに気付くと思います。
数学的に言うとこれは二重積分

の解であって、答えは一応計算上4となるんですが、こんな積分計算をせずとも、モンテカルロ法を利用すれば、近似解を求める事が可能なのです。

function V = taisekimontecarlo(n)
 count = 0;
  for i = 1:n
   sand = [unifrnd(0,pi) unifrnd(0,2) unifrnd(0,1)];
    if (sand(1,3) <= sin(sand(1,1)))
     count = count + 1;
    endif
  endfor
 V = count/n*pi*2*1;
  

この問題のポイントは、高さzは変数xだけに依存していて、yには全く影響を受けていない、と言ったところです。
そこで、プログラム上のテクニックとしては以下が上げられます。

4行目:
sand = [unifrnd(0,pi) unifrnd(0,2) unifrnd(0,1)];
%%一様乱数を3つ別々に生成して、それをベクトルsandとして纏めている。
%%この理由は、発生する座標の成分によって、定義域が違うから。
%%1つ目の乱数はx成分で、これは0〜πの範囲で乱数を生成。
%%2つ目の乱数はy成分で、これは0〜2の範囲で乱数を生成。
%%3つ目の乱数はz成分で、これは0〜1の範囲で乱数を生成。
%%なお、問題設定としてはz＝sin xで、かつ0≦x≦πなので、
%%zの最大の高さが1なのは数学的に自明。

5行目
sand(1,3) <= sin(sand(1,1))
%%前述の通り、zはxだけに依存してyには全く関係がない。
%%y方面からzを眺めると、まるで定数関数のように見える事からも
%%これは言える。
%%従って、乱数で生成されたzの値が、乱数で生成されたxを代入された
%%sinの値以下であれば、取りあえず題意は満たす事になる。

最終行:
V = count/n*pi*2*1
%%count/nの部分はいつもの通り。
%%これに全体の長方形の体積をかければ、求める答えとなる。
%%長方形の体積は縦×横×高さより、
%%π×2×1となる。 

上のプログラムtaisekimontecarloを利用して100,000回のモンテカルロ･シミュレーションを行った結果が以下の通りです。

理論値である4に極めて近い結果が出てるのがお分かりでしょう。

体積を求めてみよう�� †

モンテカルロ法を用いれば、次の図のようなわけの分かんない図形の体積さえ求める事が出来ます。

何か水分子の模型みたいですが(笑)、要するに中心が(0.3,0.3,0.7)と(0.7,0.7,0.7)の半径0.3の2つの球と中心が(0.5,0.5,0.5)で半径が0.4の球が合体しちゃってるんですよね。
球一つの体積なら問題なく数値計算できますが、このように互いに食い込んでいるような図形は、受験やテストの問題ならいざ知らず、あんまり解きたい問題ではありません(笑)。そこでこう言う場合も、｢概要だけ知りたい｣と言うのなら、コンピュータに丸投げしちゃって、モンテカルロ法で計算した方がマシなのです。

function y = myball(n)
 count = 0;
  for i = 1:n
   u = unifrnd(0,1,1,3);
   pt = (u(1,1)-0.3)^2 + (u(1,2)-0.3)^2 + (u(1,3)-0.7)^2;
    if (pt <= 0.3^2)
     count = count + 1;
    else
     pt = (u(1,1)-0.5)^2 + (u(1,2)-0.5)^2 + (u(1,3)-0.5)^2;
      if (pt <= 0.4^2)
       count = count + 1;
      else
       pt = (u(1,1)-0.7)^2 + (u(1,2)-0.7)^2 + (u(1,3)-0.7)^2;
        if (pt <= 0.3^2)
         count = count + 1;
        endif
      endif
    endif
  endfor
 y = count/n; 

もうだいぶモンテカルロ法の考え方に馴れたでしょうし、特に新しいテクニックも導入していないので、特に説明はしません。
が、要するに、3つの一様乱数の要素を持つベクトルuを定義し、それぞれの球の範囲内にあるかどうかif,else文を使って順次判定していってるだけの単純なプログラムです。
しかしながらこれでも、手計算で積分するより遥かに簡単に解の概要を知る事が出来るのです。
ちなみに1万回のシミュレーションによるモンテカルロ計算の結果は以下のようになりました。

その他の有用な乱数 †

既に二項乱数は紹介したんで、比較的使用頻度の高そうな二つの乱数プログラムを紹介して、OCTAVE Workshopによるモンテカルロ法入門は、取りあえず了としましょう。

OCTAVEとグラフィックス †

OCTAVE Workshopではplot関数を使ったグラフィック表示が可能です。実際はOCTAVEそのものにはグラフィックデヴァイスは無いのですが、その代わり、外部からGNUplotと言う描画系プログラムを呼び出してMATLAB互換のコマンドで使う形になっています。OCTAVEではGNUplotを別にコンピュータにインストールしなければ使えませんが、一方、OCTAVE Workshopではどちらもデフォルトでインストールする仕組みになっているのです。
ただし、本当の事を言うと、OCTAVE Workshopはまだベータヴァージョンで、かなりのバグを抱えており、特にこの描画系のコマンドplot関連で一番多くのバグを持っているようです。⇒参考:Octave Cygwinメモ
よって、ここでは不安定な描画/グラフィックス関係は扱わない事にします。いずれOCTAVE Workshopの安定版がリリースされた時、グラフィックス関連のトピックを扱いたいと思います。

今はここではGNUplotに関連するホームページの情報を纏めるだけに留めておきます。

OCTAVE参考リンク集 †

OCTAVE関連の情報は以下のサイトが参考になります。

コメント †

競馬@Wikiへ戻る