Find limit ...(極限値を見つける)

Maximaで極限値の計算をする場合は[Calculus]プルダウンメニューから[Find limit ...]メニューを選びます。



これによって現れるポップアップメニューは以下の通りです。

Limit of:極限値を求めたい数式を入力します。 when variable:上記の数式での変数を入力します。デフォルトはx。 goes to:変数をどこへ近づけたいのか入力します。デフォルトは0。なお、横の[Special]ボタンで他の場所も選択できます。 Pi:円周率へ近づけます。 E:自然対数の底(ネピア数)へ近づけます。 Infinity:正の無限大の方向へ近づけます。 - Infinity:負の無限大の方向へ近づけます。 from:関数の極限値は右から近づけた場合と左から近づけた場合と数値が食い違う場合があるんで、問題設定によって「どこから近づけるのか?」設定しなければならない場合があります。デフォルトはboth sides(両側から)です。 both sides:両側から極限値へ近づけます。 left:左側から極限値へ近づけます。 right:右側から極限値へ近づけます。 Use Taylor series:テイラー展開を用いて極限値を求めたい場合はここにチェックを入れます。

例題

極限値　を求めよ。

Maximaの[Limit]ポップアップのLimit of:に次のように入力します。



そして[OK]ボタンを押すと解が表示されます。



答えが−2だと言うことが分かります。


Maxima Manual: 極限

質問＜３６２８＞２００７／１０／１９from=BOXY「DO〜LOOPのBASICプログラム」

質問させてください。

100 FOR x=1 TO 100
110    FOR y=x TO 100
120       LET a=x
130       LET b=y
140       DO
150          LET r=MOD(a,b)
160 IF r=0 THEN EXIT DO
170 LET a=b
180 LET b=r
190 LOOP
200 IF b=1 THEN
210 LET z=SQR(x^2+y^2)
220 IF INT(z)=z THEN PRINT x,y,z
230 END IF
240 NEXT y
250 NEXT x
260 END

BASICプログラムですが、実行結果の表示とプログラムの意味が
理解できません。の表示になるのですが、
特に前半のDo〜LOOP構文の意味がいまいち理解できません。
アドバイスをお願いいたします。

★希望★完全解答★

---

いやあ、BASICですか。分かりませんね、こりゃ(苦笑)。
そもそもインデント無しのソースコードって読むの辛いです。どこがブロック構造なんだかちっとも分かりゃしない。現代的じゃない。
さすがFortran直系の言語です(苦笑)。
誰かが

今日のすべての言語はFortranかLispに基づいている。まずい方がFortranで、すばらしい方がLispだ。

とか言ってましたが、むべなるかな、と言うカンジです。
あ、別にBOXYさんのせいじゃありませんよ。こんなの教える、って決めた方が悪いと思います。時代遅れもはなはだしいです。

とか愚痴っててもしょーがないんで、取り敢えず走らせるだけ走らせてみましょうか。
今回、プログラムを走らせるBASICは十進BASICと呼ばれるモノです。Linux版があったので助かりました。
そしてその後、Schemeを使ってBASICのコードを翻訳してみます。その過程で「何が行われているか」分かるでしょう。僕の方も手探り状態なんで(笑)、まあ、一歩一歩進めてみましょうか。
なお、プログラミング言語Schemeに付いては大体の概略は質問＜3567＞、質問＜3590＞、質問＜3609＞辺りに書いていますんで、ザーっと目を通しておいてください。
ないしは初心者向けのサイトでScheme入門と言う易しい入門サイトもありますし、また、Windowsユーザーでしたら、タダで手に入り、遊びながらSchemeの基本が学べる魔法言語リリカル☆Lispと言うアドベンチャー・ゲーム風のチュートリアルソフトがあります(恐らくラストまで進むのに1日かからないでしょう)。

魔法言語リリカル☆Lisp画面写真

BASICを初めて習って大変なのに・・・・・とか思うかもしれませんが、絶対損はさせないのでやってみてください。

閑話休題:BASICは手続き型言語の入門用としては著名ですが、一方Schemeは関数型言語の入門用として著名な言語です。世の中にはたくさんの手続き型言語があって(例えばC言語、Pascal、Java等)、そっちの方が現時点で確かにポピュラーですが、現在は関数型言語の方がある意味、先端技術として注目を集めています(代表的なトコでHaskell、OCaml等)。ソフトウェアの開発効率としては関数型言語の方が優れてる、と言うのがその注目を集めてる理由なんだそうですが、一方、関数型言語は学ぶのが難しいと言われています。
しかしながら、ホントの事を言うと、関数型言語は学ぶのが難しい、と言うよりも、手続き型言語を徹底して学んだ人間にとっては、思考回路を関数型言語に合わせるのが難しい、と言うのが事実のようです。つまり、どうせやるならより初心者に近い時点で関数型言語の思考形式にある程度慣れた方が得策だ、と言う事です。こればっかりは知ってる方が知らないより損なのです。
と言うワケで、将来的な事を考えたら、この時点でSchemeを触った方がイイ、と言う事を明言しておきます。BASICに習熟すればするほど、関数型言語はどんどんどんどん分からなくなっていくでしょう。
一昔前だったらプログラミング言語を「買う」事で余計なコスト/リスクが生じて、こう言う「当たって砕けろ」式の学習法は勇気とお金のいる行動だったのですが、生憎、今は殆どのプログラミング言語は(Microsoft関係を除き)殆どタダで手に入るありがたい世の中になっています。分からなかったら分からなかったでうっちゃらかしておいて構わないので取りあえず触っておくだけ触っておけと言う事です。

では取りあえず、十進BASICを使って問題のソースコードを走らせてみましょうか。

	端末を使って十進BASICを起動する。
	十進BASICが立ち上がる。
	問題のソースコードを十進BASICにコピペする。
	[実行]ボタンを押してプログラムを走らせると、結果が表示される。

なるほどね。確かにの表示になっていますね。
さて、なかなか十進BASICは賢いようで、[実行]した後、ソースコードをインデント付きに書き換えてくれています。以下はその「インデント付き」のソースコードです。

100 FOR x=1 TO 100
110    FOR y=x TO 100
120       LET a=x
130       LET b=y
140       DO
150          LET r=MOD(a,b)
160          IF r=0 THEN EXIT DO
170          LET a=b
180          LET b=r
190       LOOP
200       IF b=1 THEN
210          LET z=SQR(x^2+y^2)
220          IF INT(z)=z THEN PRINT x,y,z
230       END IF
240    NEXT y
250 NEXT x
260 END

上のインデント付きのソースコードだとどこからどこまでがブロック（一区切りの処理）なのかとても分かり易くなっていますね。 行頭が揃っているトコロが1ブロックです。
何故ブロックが重要なのか?亀田は正直、全然BASICなんてやった事が無いですし、あまりハッキリとした事は言えないんですが、一般的に現代的な言語だとどれも、ある処理だけ取り出して単一のプログラムとして書くことが出来ます。これを専門的にはカプセル化と呼ぶらしいんですが、これがBASICには出来ないんです。非常に原始的と言わざるを得ない。
もちろん、BASICでもサブルーティンを設定して・・・・・と言うような事は出きるそうですが、もっとハッキリとしたGo To文なんかで明示しないと目標としたサブルーティンに飛べない。かつ、メインのプログラムに戻るときゃあどうすんだ?とか色々とメンド臭い事柄が生じます。従って、原則BASICの流儀としては一つのソースでは一つのプログラム、とならざるを得ない。ちょっと複雑な処理をしようとしたらたちまちソースが大きくなってしまって、たとえ自分がソース書いていても何やってんだか分からなくなりそうです。行数が32,768行なんかに肥大したらどうすんでしょうね(笑)?
一方、現代的なプログラミング言語の殆ど(例えばPerl、Python、Ruby等)はBASICの仲間であるインタプリタですが、こう言うメンド臭い作業は必要としません。同じ一枚のテキストファイルにプログラムを書き込む以上、どのプログラムからどのプログラムに飛ぼうが自動で見つけて処理してくれます。従ってGo To文なんて使う必要も無く、書く側から言うと「小さいプログラムを一つづつ書いて」全体として処理が一環していればいい、と言う利点があるんです。その方が書き間違いしても見つけやすいですしね。32,768行の大きな一つのソースを書くよりは8行程度の小さいプログラムが4,096個あった方がマシだ、と言う事です。そしてSchemeもこう言う特質を備えた言語です。
さて、では、まずは問題のDo〜LOOP構文(140〜190)を見ていきたいと思いますが、その前に一つプログラミング言語の約束事を押えておきましょう。
それはどんな言語(例えばC言語)でもそうですが、原則、

処理はソースの上から始まって下へと行われる

と言うことです。まあ、例外もひょっとしたらあるかもしれませんが、少なくとも亀田は下から上へ、ってのは知りませんね。あったとしてもそんなの使いづらくてしゃーないです(人間は横書きだったら上から下へと読んでいく、ってのが習慣だから、です)。
右から左へ・・・・・ってのも知りません(笑)。まあ、そんな言語作ったら是非ともムーディ勝山と名付けて欲しいトコなんですが(右から左へと受け流す言語・笑?)、いずれににせよ、今のトコ、そう言ったルールを採用している言語は無いと思います。
何が言いたいか、と言うと、原則Do〜LOOP構文(140〜190)の処理が終わらないとそれ以降の処理が為されない、と言う事です。つまり、200行目以降はDo〜LOOP構文(140〜190)の処理の結果を用いて「何か」を計算している。また、同様の理由で130行目までの処理が終わらないとDo〜LOOP構文(140〜190)の処理が始まらない、と言う事です。取りあえずBASICを知らない僕でもあくまでプログラミング言語の一般常識でそこまでは類推出きるのです。
さて、僕はBASICのDo〜LOOP構文は「繰り返し計算だろ」と言う予測は出来ましたが、ところが具体的な書式は知らない。そこで調べてみました。

◇ DO〜LOOP
　EXIT DO文が実行されるまで，DO行とLOOP行の間に書かれた文を繰り返す。

BASIC入門より

ああそうですか。なるほど。分かりました。
ではこの情報だけを元にDo〜LOOP構文(140〜190)だけをSchemeでカプセル化してみましょう。
ポイントとしてはa、bと言う数は何だか分からない。が、しかし、その数はDo〜LOOP構文(140〜190)以前のソースから渡される事だけは分かっています。従ってプログラム・・・・・名前はそのまんまdo-loopとしますが、do-loop内では定義する必要はない、と言う事です。
また繰り返し計算が何らかの条件を満たした場合、その計算結果を明示しないとなりません。その計算結果を利用して200行目以降のプログラムが処理してくれるわけですが、BASICの場合、EXIT DOだけ書けばいいようですが、Schemeの場合、明解な返り値を設定しなければならない、って事です。しかしこれは実はそんなに難しく無いです。
最初に言っておきますが、これはズバリ、再帰構造です。再帰構造に付いてはあとで説明するとして、では、Schemeで書いたプログラムdo-loopのソースをご紹介します。

(define do-loop ;プログラムdo-loopを定義

(lambda (a b) ;引数a,bを設定

(let ((r (modulo a b)) ;局所変数でrを計算

(if (zero? r) ;rが0だったら?

b ;bを返す
(do-loop b r))))) ;再帰呼び出し

恐らく、局所変数letに関しては、BASICもSchemeも同じようなものでしょう。ですからそれに付いては説明しません。また、if構文に関してもBASICでも似たような構造でしょう。その辺は上のソースを今まで培ってきたBASICの知識で読んでみてください。そしてそのあたりは特にポイントじゃないんです。
まず、一つ目のポイントとしては、BASICのソースにあるMOD、そしてSchemeで言うmoduloと言う命令です。
これは原則同じ命令を表していて、剰余を計算する命令です。例えば、BASICでは

MOD(3,2)→1

となるでしょうし、Schemeでも

(modulo 3 2)→1

となります。何故なら、3を2で割った余りは1だから、ですよね。単純に余りを求めているんです。
さて、BASICではrが0だったらEXIT DOしろ、と読めます。この部分ですよね。

160          IF r=0 THEN EXIT DO

この部分がSchemeで言うと以下の部分です。

(if (zero? r) ;rが0だったら?

b ;bを返す

SchemeがEXIT DOの代わりに何故bを返しているのか、と言うのは後で説明しますが、取りあえずここで分かるのは、aがbで割り切れる時、つまりrが0だった時にDo〜Loopを脱出しろ、って命令になっている。言い換えるとこれが脱出条件で、その時Do〜Loopが繰り返し計算を終えるのです。(つまり、これが無いと、プログラムはDo〜Loop内で延々と計算を繰り返す事となるんです)。
では問題はaがbで割り切れない場合はどうなるか、って事です。今、ここのブロックをカプセル化した時、a、bがどんな数かは問わない事としました。故にこの部分だけ見ると、a、bはどんな数を代入してもOKです。従って一般にaがbで割り切れるとは限らない。
さて、もう一度Schemeのソースに戻って良く見てほしいんですが、擬似コードとしては次のようになってるんです。

1. プログラムdo-loopを定義する


2. プログラムdo-loopの引数はaとbの二つの数とする


3. 局所変数としてrを定義する


4. rが0だったらbを返す


5. それ以外だったらプログラムdo-loopを引数b、rとして呼び出す

分かりますか?ポイントは5番です。今、プログラムdo-loopを作っていますが、プログラムdo-loop内でdo-loopを呼び出しています。つまり、自分の中で自分を呼び出していて、こう言うのを再帰的定義と呼ぶんです。
つまり、(do-loop a b)と言う形でプログラムを呼び出すと、rが0でない限り、中で(do-loop b r)が呼び出される。そして(do-loop b r)の中ではまた別のrが定義されて、それが0で無い限り・・・・とどっかで剰余計算の解が0にならない限り、延々とそう言う引数の書き換えだけが生じてdo-loopが呼び出され続けるんです。
これがBASICのソースで言うと、170行目と180行目のコードの意味です。

170          LET a=b
180          LET b=r

Schemeのように明示的な再帰呼び出しを示してはいませんが、この2行の命令の意味は、Schemeで言うトコの引数書き換えで、その状態のまま150行目に作業が戻ってるんです。従って、Do〜Loopの部分がまた計算されてrが0じゃなかったらまた引数が書き換えられて・・・・・・と、原理的にはSchemeの明示的な再帰呼び出しとまったく同じ作業を展開してるんです。
これが設問のDo〜LOOP構文(140〜190)の意味です。
なお、どうして返り値をbにしたのか、と言うと、元々のBASICのプログラムの200行目にその答えがあります。

200       IF b=1 THEN

先ほど、「プログラム言語は上から下に向かって実行される」と言いました。従って、Do〜LOOP構文の下に来て、なおかつDo〜LOOP構文に関わりのある引数が指定されている場所と言うのはここしかありません。aもrも他にはどこにも見当たらないのです。
と言うことは、一番欲しい数はbである、と言う事です。従って、計算結果としての返り値はbを指定しておけばいい、と言う事になります。

さて、ではついでに200行から240行までもSchemeでカプセル化してみましょう。ここではプログラム名をそのまま200-240とします。
プログラム200-240はちょっとしたBASICの欠点と、同じくSchemeの欠点を垣間見る事が出来るでしょう。お互いの長所、短所を比較するには面白いブロックだと思います。
ただし、プログラムとしてはそんなに難しくありません。基本の範疇だと思います。よって、ソースをチャチャッと書き下してみましょう。

(define 200-240 ;プログラム200-240を定義

(lambda (x y a b) ;4つの引数を設定

(let ((z (sqrt (+ (expt x 2) (expt y 2))))) ;zの計算

(if (and (= (do-loop a b) 1) ;条件節

(integer? z)) ;整数を確かめる
(begin (display x) ;beginで順次処理開始

(display " ") ;BASICのPRINT=display
(display y)
(display" ") ;文字列でスペース指定
(display z)
(newline) ;改行
(next-y x (+ y 1))) ;next-yを呼び出す
(next-y x (+ y 1))))))

基本的にはBASICでもお馴染みの単なるif-then構文ですね。
プログラム200-240は4つの引数x、y、a、bを持ちます。先ほどbだけが必要だ、と言ったのにaも引数として設定しているのはプログラム200-240内でdo-loopを呼び出す為です。do-loopはa、bの二つの引数が必要でした。
さて、Schemeの欠点その1。局所変数でzを書いていますが、まずその式が分かりづらいですね(苦笑)。これは演算子を最初に必ず置くと言う約束事、すなわち前置記法を採用しているSchemeの欠点です。もちろん前置記法には利点も多いんですが、その代わり、

なんて書かなきゃならない場合はちょっと大変です。少し考え込んじゃいますね(笑)。
まあ、丁寧に書いていけば決して中置記法(普通の数式の書き方)を前置記法に直すのは決して難しくないんですが、ちょっとした慣れが必要です。
なお、Schemeでは、sqrtと言う命令がBASICのSQRにあたり、また、(expt s t)と言う命令でを計算させます。
次はBASICの欠点です。BASICのソースコードでは

200       IF b=1 THEN
210          LET z=SQR(x^2+y^2)
220          IF INT(z)=z THEN PRINT x,y,z
230       END IF

となっていて、実は条件自体は

bが1でかつzが整数の時

と2つの条件を同時に考慮している筈なんですが、途中で局所変数によるzの定義が挿入されてたりしてあまりスマートではありません。最初にzを定義して、まとめて条件を記述してIF〜THEN構文に持ち込んだ方が感覚的にはキレイな筈なんですが、少々ゴチャゴチャしています。何故なんでしょう?
調べてみたら、BASICでは複数行にわたってしか書けないものを実行させたい場合、GOTO命令を使い行を飛ばす必要があるらしく、要するに一行で複数の条件を記述出来ないようです。従って、少々不格好なソースとならざるを得ません。一行には一つの条件、と言うキツい制約があるらしい。
その点、Schemeの方がソースとしては明解です。

(if (and (= (do-loop a b) 1) ;条件節

(integer? z)) ;整数を確かめる

とif節の中にandを使って上手く2つの条件節を纏めています。
なお、一つ目の条件節で先ほど作ったdo-loopを呼び出しています。BASICのオリジナルのソースでは、単にbを使って判定していましたが、Schemeヴァージョンでは関数丸ごと呼び出すようにしています。
と言うのも、bだけ記入しても何が何だか分かりませんし、明示的に返却されたbが必要だから、と言う理由に依ります。この条件内で先ほど書いた再帰定義プログラムdo-loopが実行されて、返り値が1かどうかここで判定されているわけです。プログラム同士が有機的に結合しているんですね。
もう一つがinteger?と言う命令で局所変数で定義されたzが整数かどうか判定しています。ここがBASICのオリジナルのソースでは220行目でのINT(z)=zと言う場所の部分なんですが、若干意味が違います。
元々BASICのINTと言う命令は、zを整数化せよ、と言う命令だそうです。すなわち、小数を四捨五入する命令なのでしょう。zを整数化したモノとzが同じだったら元々zは整数であり、数値が食い違えばzは整数ではない、と言う事です。この辺りはBASICでは結構まだるっこしい判定法を用いているようです。Schemeのinteger?に比べるとかなり回りくどい判定法ですね。
なお、Schemeではこのテの「判定用命令」は大体最後に「?」がくっついています。
Schemeの欠点その2。BASICのPRINTにあたるSchemeの命令がdisplayなんですが、原則displayは一つ、ないしは二つの引数しか表示出来ません。BASICの場合は一回PRINT命令を送れば複数の文字を出力できますが、displayの場合複数の文字を表示させるとなるとそれなりの行数の命令を書かないとなりません。従って、BASICでは

PRINT x, y, z

と一行で済むところをSchemeでは

(display x) ;xを出力
(display " ") ;スペースを出力
(display y) ;yを出力
(display" ") ;スペースを出力
(display z) ;zを出力
(newline) ;改行命令

と羅列しまくらないといけません。
一般にSchemeの入出力系列の命令は結構貧弱なんですね。と言うのも、教育用言語と言うことで必要最小限の入出力命令しか持ってないから、なんです。これがSchemeの欠点のうちの一つだと思います。
なお、beginと言う命令は、中に入っている命令を「順番に」実行せよ、と言う意味です。ここでx〜zまでの出力と改行までの命令を順次実行させています。そして大事なのは「次のyが何なのか」ここで決定しないといけないと言うことです。
もう一度プログラム200-240の擬似コードを見てみましょう。

1. プログラム200-240を定義する


2. 引数x、y、a、bを設定する


3. 局所変数を定義する


4. プログラムdo-loopの結果が1で、かつzが整数だった場合、x、スペース、y、スペース、z、改行の順で出力してプログラムnext-yを引数を一つ増やして呼び出す


5. そうじゃなくてもプログラムnext-yを引数を一つ増やして呼び出す

となっています。
ここでポイントとなってるのは、BASICのオリジナルのプログラムで見れば240行目のNEXT yにあたる部分で、ここでyの数を一個増やしてアタマの110行に戻っているんです。従って、ここでScheme版でも新しくnext-yと言う名前のプログラムを作成して呼び出す事、とします。また、条件節の結果が真であろうと偽であろうと、next-y自体は呼び出されなければならない、と言う事です(じゃないとプログラム全体が止まってしまいます)。
では次にプログラムnext-yを作ってみましょう。BASICオリジナルのソースの110行〜130行、そして240行目にあたります。
設計にあたって気をつけることは、

1. プログラムnext-yはx、yの二つの引数を持つ


2. yの値が100を越えたら引数xを一つ増やして、また引数yを一つ増やしたxの値をもって初期化する


3. 局所変数aをxとし、また局所変数bをyとする


4. プログラム200-240を呼び出す

と言う流れです。
ポイントとしては、プログラムnext-yからはプログラムdo-loopは呼び出さないで、敢えてプログラム200-240を直接呼び出している、と言う事です。と言うのも、プログラム200-240内でdo-loopを呼び出すようにしているので、それ以上呼び出す必要がないから、です。
また、問題設定から言うと、yが100に達したとき、xの数を1増やして、またその値をyの初期値としてセットしなおしています。この作業をプログラムnext-xを新しく作って呼び出す形にしています。
この2つのポイントを押えてプログラムnext-yを書くと以下のようになります。

(define next-y ;プログラムnext-yを定義する

(lambda (x y) ;引数をx、yとして設定

(if (> y 100) ;yが100を越えたら

(next-x (+ x 1) (+ x 1)) ;next-xを新しい初期値を引数として呼び出す
(let ((a x) ;局所変数aをxとする

(b y)) ;局所変数bをyとする
(200-240 x y a b))))) ;プログラム200-240を呼び出す

また新しいプログラムnext-xが出てきたので、これを作ります。これはオリジナルのBASICのソースの100行と250行にあたります。もうここまで来れば、xが一つづつ増えて行った末の終了条件さえハッキリ分かっていれば構いません。

(define next-x ;プログラムnext-xを定義

(lambda (x y) ;引数x、yを設定

(and (<= x 100) ;xが100以下で
(next-y x y)))) ;プログラムnext-yを呼び出す

andと言う命令は内側に含まれる命令が偽であった時、偽を返して終了します。従って、これによりxが100以下の時に限ってnext-yを呼び出し続けるのです。
最後にxとyに初期値を与えましょう。初期値は両者とも1でしたね。これをプログラムQ3628として定義します。

(define Q3628 ;プログラムQ3628を定義する

(lambda () ;引数は取らない

(next-x 1 1))) ;プログラムnext-xの引数を1、1として呼び出す

これでプログラムQ3628を呼び出した時点でプログラムnext-xが引数を1、1として呼び出され、それがまた別のプログラムを呼び出して・・・・・とまるでドミノ倒しのように全てが一つのプログラムとして動きます。
もう一回全部のソースコードを見ておきましょう。

(define Q3628

(lambda ()

(next-x 1 1)))

(define next-x

(lambda (x y)

(and (<= x 100)
(next-y x y))))

(define next-y

(lambda (x y) 　　

(if (> y 100)

(next-x (+ x 1) (+ x 1))
(let ((a x)

(b y))
(200-240 x y a b)))))

(define do-loop

(lambda (a b)

(let ((r (modulo a b)))

(if (zero? r)

b
(do-loop b r)))))

(define 200-240

(lambda (x y a b)

(let ((z (sqrt (+ (expt x 2) (expt y 2)))))

(if (and (= (do-loop a b) 1)

(integer? z))
(begin (display x)

(display " ")
(display y)
(display" ")
(display z)
(newline)
(next-y x (+ y 1)))
(next-y x (+ y 1))))))

そして全体的には次のような流れになっています。

Q3628→呼び出し→next-x→呼び出し→next-y→呼び出し→200-240→呼び出し→do-loop
　　　　　　　　　　　↑         ↓    ↑          ↓
　　　　　　　　　　　←←呼び出し←←     ←←呼び出し←←

このようにBASICのインデントに着目すれば、行頭が同じ位置にあるソースコード群をまとめてカプセル化する事が出来ます。また、プログラム実行の流れは「上から下へ」と言うのは事実なんですが、カプセル化するにあたっては特に時系列で書き連ねなければならない、って事はありません。気になった部分、好きな部分、ないしは簡単そうな部分からプログラムを書き始めて行って、最終的に整合性さえ取れていればいいのです。これが「現代的な」プログラムの書き方だそうです。「出来る部分から書く」のがワリに今では当たり前になってるんですね。
さて、ではSchemeでのソースの実行結果を見てみましょうか。




はい。全く同じ計算結果が表示されていますね。簡単でした。
取りあえず、今回のDo〜Loop構文のポイントは再帰構造です。そこだけ押えておけば取りあえずは難しくないでしょう。
また、このように「全く違う発想の」言語へ翻訳してみることによって、色々なプログラムの発想が得れるでしょう。勉強方法としてはメンド臭いように見えますが、意外とこの方が効率的に勉強できる、と言う事を付け加えておきます。

以上です。

---

さて、ここから余談です。
上の本編の方では再帰構造にこだわって解説しました。何故ならDo〜Loop構文が行ってる意味が分からない、と言う事だからです。
このテの問題の解説が難しいのは、数学的に意味が分からないのか、それともプログラムの部分で意味が分からないのか見えづらい、って事があるんです。どっちつかずなんですね。
それでDo〜Loop構文に絞って解説したんですが、一方、どうして剰余計算MOD(Schemeではmodulo)が出てくるのか意味が分からない、ってのも事実です。そこでまずはそれを解説したいと思います。
次に、BASICのソースコードをSchemeに翻訳してみせる事により、「プログラムとしては何をやってるのか?」明示する戦略を取ったんですが、しかし出来上がったSchemeのコードはあまりキレイではありません。
本編の方で

「開発効率の良さで関数型言語が注目を集めている」

と書きましたが、これじゃあんまり開発効率が良い、とは言えませんよね。
通常、開発効率の良さ、と言うのは

1. コードの読みやすさ


2. コードの量

で計られます。
読みやすさに関しては「慣れ」の問題が大きいのですが、「コードの量」ってのは別の指標を与えてくれます。同じ動作をするにあたって、書く量が増えるのは人間にとってたまったもんじゃありません。
どうしてC言語は動作は速いのに開発効率が悪いのか、と言う一つの理由がここにあります。何故なら「同じ事をやらせるのに」書く量が半端無く多いのです。ライブラリの呼び出しから始まって使われる変数の型の定義、その他諸々・・・・・・。コンピュータにとっては問題無い効率的な処理ですが（と言うよりC言語はコンピュータだけにとって効率的な命令になるように設計されている）「指令を出す」人間に取って大変なのです。人間にとって「効率的」じゃないんです （難しいのはコンピュータにとって効率的な事が人間にとっても効率的とは限らない辺りです！！！）。
何故今回Schemeで書いたコードが長くて開発効率が悪いように見えるのか?その理由は単純に一つです。それはBASICで書いたソースの構造をそのまま流用して書いたコードだからです。要するに、単純な手続き型言語の構造をそのまま関数型言語に当てはめてしまった。それがやたら冗長なコードになってしまった理由です。つまり、冗長なのはSchemeのせい、と言うより元々のBASICのソースが冗長なんです。
あくまで亀田個人の意見ですが、オリジナルのBASICのソースのように、内側にxのループがあり、その中にyのループがあり、そしてDo〜Loopの構文がある、なんてのは汚くて嫌です。これしかBASICでは書きようが無いのかもしれませんが、Schemeなら原則そんな書き方はしません。どうせやるんだったら、一回のループ構文で全てを表現した方がカッコいいのではないか、と思います。再帰構造を用いるならそっちの方が本来スマートな筈です。
この2点に絞って「余談」として、もっと効率的なSchemeのソースを提供してみたいと思います。

まずは最初のDo〜Loop構文の「数学的な意味」を考えてみましょう。と言うより、実は「数学的」って程のネタじゃないんですね。
多分、元々問題の発想としては、

、の範囲でが成り立つ整数x、y、zを求めよ

だと思っているでしょう?いわゆる3平方の定理ですよね。
それは基本的には間違っていないんです。しかしながら、上のような問題だと、わざわざDo〜Loop構文を使ってBASICのソースの内側に「何かを」忍ばせる必要が無いんです。
論より証拠。何故あそこでDo〜Loop構文が必要になるのか?次の表を見てみれば一発で分かると思います。これらが、の範囲でが成り立つ整数の組です。



はい。実は上の表のようになるんですね。単純に、でが成り立つ整数x、y、zの組み合わせを数えると‥・・・・実は59組も存在するんです。ちょっと見てみてください。
確かにBASICのソースを実行した時に示されている組も入っていますが、多すぎますね。何が違いかお分かりでしょうか?
そうですね、実は互いに素であるx、yの組み合わせ、って言う条件が付加されていないといけないんです。
例えば、ってのは成り立っています。同様にってのは成り立っています。しかしながら問題のBASICプログラムは3、4、5の組を計算ではじき出してくれますが、一方6、8、10の組は無視です。何故ならだからですよね。全部を4で割れば結果と変わりありません。つまり、6、8、10の組合わせは互いに素ではないのです。
互いが素である為には何が必要か?x、yの最大公約数が1である、ってのが条件となりますね。だから、返り値b＝1の時、と言う判定条件が付いたのです。
もうお分かりですか?実は問題のDo〜Loop構文は、平たく言うと、xとyの最大公約数を調べるプログラムなんです。 →参考：ユークリッドの互除法
一方、Schemeでは実は最大公約数を計算する命令、gcd(Greatest Common Diviser＝最大公約数の事)があります。これを上手く使ってプログラムQ3628を書き直してみたいと思います。
もう一つ、これは亀田の好みかもしれませんが、xを1から変化させていって100に到達させて、またyをxの値から変化させていって100に到達させる、と言う構造は好きくないです。僕だったら100から逆に辿っていきます。何故なら、一々yの値を設定しなおすのが面倒臭いから、ですよね。
条件としてはxが100から0に向かって変わっていき、0で全体の計算が終了する。また、yも100からxの値に向かって行ってx＝yになったときまた100に戻る、って考えた方が分かりやすいと思います。一々xの値が何なのか、チェックしてそれによってyの初期値が決まる、ってのはバカバカしいでしょう?後者だったらどの道100からスタートすれば良いのが自明です。
以上を考慮して、プログラムQ3628-newを作りますが、まずは大まかな再帰プログラムエンジン、Q3628-engineを記述します。ここでバラバラだったx、yの変化、そしてDo〜Loop構文を完全に一つの大きな再帰構造の中で記述する事になります。

(define Q3628-engine

(lambda (x y ls)

(if (zero? x)

(print-engine ls)
(let ((z (sqrt (+ (expt x 2) (expt y 2)))))

(Q3628-engine (if (= x y)

(- x 1)
x)
(if (= x y)

100
(- y 1))
(if (and (integer? z)

(= (gcd x y) 1))
(cons (list x y z) ls)
ls))))))

これが完全に一つの再帰構造で書かれたQ3628-engineです。
まずはポイントとしては今回は引数をx、y、lsの三つにしました。x、y、は問題のBASICソースにある通りですが、ここでは新しくリスト用の引数lsと言うのを設定しています。これは、条件にしたがったx、y、zの計算結果を次々と突っ込む為のスペースです。
全体の流れとしては次のようになっています。

1. プログラムQ3628-engineを定義する


2. x、y、lsの3つを引数として設定する


3. xが0の時、プログラムprint-engineを呼び出し、引数lsを渡す


4. それ以外の場合は局所変数zを計算し、Q3628-engineを呼び出す

たったこれだけ、です。非常にシンプルな構造になっています。
問題は再帰的に呼び出された4番のプロセスです。ここで呼び出されたQ3628-engineの引数を工夫してるんですね。
プログラムQ3628-engineは3つの引数を持っていました。この3つの引数が繰り返し計算中に書き換わっていくのですが、そこにif節を突っ込む事により繰り返し計算を効率的に表現しているんです(こう言う事ができるのが、手続き型言語と違った関数型言語の柔軟な所です)。
すなわち、再帰的に呼び出されたQ3628-engineのそれぞれの引数は、

• x:もしx＝yなら、xを1減らす。そうじゃなければxのまま。


• y:もしx＝yなら、yは100に戻る。そうじゃなければyを1減らす。


• ls:もし、zが整数で、かつ、xとyの最大公約数が1だったら、x、y、zのリストを作って、それを引数lsの先頭に突っ込む。そうじゃなければlsのまま。

註:consと言う命令はリストの最初に何かを突っ込む命令。また、listは任意の数のリストを生成する命令。両者とも質問＜3609＞参照。これらはSchemeの基本的な命令である。なお、リストとはBASICで言うと一種の配列にあたる。

と言うif構文によって制御されているのです。明解でしょ?かつオリジナルのBASICの制御構造を完璧に表している事が分かると思います。一々ループ構造に悩まないでスッキリ書けるのがSchemeの強みなのです。
なお、print-engineと言うプログラムはまだ書いていません。これは本編を読んでたら何となくどう言うプログラムになりそうかは想像が付くでしょう。その通りです(笑)。
取りあえずprint-engineも記述してみましょうか。

(define print-engine

(lambda (ls)

(if (null? ls)

(newline)
(begin (display (first (car ls)))

(display " ")
(display (second (car ls)))
(display " ")
(display (third (car ls)))
(newline)
(print-engine (cdr ls))))))

これも原則、質問＜3609＞に出てきたプログラムwrite-lsと同じものです。やっぱり再帰構造で書かれていますね。ある程度のリストをdisplayした後にリストの残り(cdr ls)にprint-engineを適用しています。
特に難しくない筈なんですが、新しく出てきたのは、まずfirst、second、thirdと言う命令です。これはつまり、リストの最初の部分、2番目の部分、3番目の部分、と言うのを表しています。意味的には簡単でしょう。
ところで、どうしてそうしてるのか、と言うと、プログラムQ3628-engineの返り値であるlsは、基本的には結局、次のようなリストを出力するから、です。

((3 4 5) (5 12 13) (7 24 25)・・・・・・)

つまり「リストのリスト」ですね。これをバラけさせてdisplayに渡さないといけない。
従って、例えば最初の出力3、4、5を得るためには、lsの最初のリストである(3 4 5)を取って、その中のfirst、second、thirdを出力した後改行、と言う命令にしないといけないのです。これがカラクリです。
ここで勘のイイ人は気づいたかもしれませんが、実はcarと言う命令と、firstと言う命令は全く同じモノです。また、リストの先頭を除去したリストを返す命令がcdrなので、実はsecondは(car (cdr リスト))と同じモノです。単に省略形ですね。
このように、どうして基本的なリスト操作命令であるcons、carやcdrがSchemeでは大事なのか、と言うと、この3つの命令を組み合わせるだけで色々なプログラム(命令)を簡単に書くことが可能だから、なんです。
まあいずれにせよ、プログラムprint-engineは質問＜3609＞のwrite-lsと原理的には全く同じコードなんで、そちらに詳しい説明を付けているのでそれを参照して下さい。
最後にプログラムQ3628-engineに引数の初期値を与えるプログラムQ3628-newを作って全て終了です。初期値はそれぞれ、xが100、yが100、そしてlsに空リストである()を与えて終わり、です。

(define Q3628-new

(lambda ()

(Q3628-engine 100 100 ())))

ご覧の通り、Q3628-newは特に引数は必要としません。Q3628-newを走らせれば、すぐに結果が返ってきます。
以下がQ3628-newを走らせた結果の画面写真です。



はい。全くオリジナルのBASICのソース実行と同じ結果が返ってきますね。そして、こちらの方がよりSchemeらしいコードの書き方です。
まあ、人によって感想は変わるかもしれませんが、BASICのDo〜Loop構文やFOR文の煩わしさよりよっぽどシンプル、かつ直感的にSchemeではコードを書くことが出来ます。これが「関数型言語」の強力なトコロですね。
また、BASICは初心者向けでとっつきが良い、と一般的には思われていますが、実は複雑な計算をさせようとすればするほど、扱いが極端に難しくなっていく言語だと思います。反面、Schemeは奥が深いのは確かですが、と同時にどんなに問題が複雑化してもあまり構文的には複雑にならない(と言うよりワンパターン)ので、むしろ初心者にとっては扱いが簡単な言語じゃないか、と思います。ハッキリ言ってBASICの方が難しいです。
まあ、こんな感じでSchemeのパワーを垣間見て下さい。また、BASICの問題をSchemeで書くことにより、かなり良質なプログラミングの練習（※）になると思います。

以上です。

※：これは必ずしもSchemeにとって良い練習、と言う意味ではなくってBASICにとって良い練習になるだろう、って意味です。もう一つの言語をBASICと同時に勉強して見比べた方が、逆にBASICへの理解が深まるでしょう。

---

追記:本編中でBasicの元ネタとなった世界で初めて作られた高級言語、Fortranに付いて書きましたが、そのFortranのオリジナルの作者、John Warner Backus 氏が2007年10月17日にお亡くなりになったようです。
ご冥福をお祈りいたします。


FORTRAN の父 John Warner Backus 氏 死去

Solve linear system ...(連立一次方程式を解く)

Maximaで連立一次方程式を解く場合は、[Equation]プルダウンメニューから[Solve linear system ...]を選びます。



これにより[Solve linear system]ポップアップが表れます。

Number of equations:方程式の数を入力します。例えば二元一次連立方程式の場合は2と入力し、三元一次連立方程式の場合は3と入力します。

方程式の数を入力し終わったら[OK]ボタンを押します(この例ではデフォルトで3を選択)。

Equation 1〜:解くべき方程式を入力していきます。 Variables:何の文字に付いて解くのか入力します。 なお、連立方程式の場合、当然解は複数存在するので、ここでそれらの文字を指定しますが、それらの文字は全てコンマ(,)で区切ります。

全て入力し終えたら[OK]ボタンをクリックします。

例題:

を解け。

Maximaに次のように入力します。



[OK]ボタンを押します。



これにより解は

だと言うことが分かります。


Maxima Manual: 方程式

質問＜３６０９＞AAZ「確率変数の分布関数と確率密度関数」

以下質問文です

X、Yの同時確率密度関数が

p(x,y)=k(x+y) (0≦x≦1,0≦y≦1)

=0 (その他)

のとき、次の問いに答えよ。

1. 定数kの値を求めよ。


2. X,Yの周辺確率密度関数、をそれぞれ求めよ。


3. X,Yの同時分布関数F(x,y)を求めよ。

★希望★完全解答★

---

ええと、ある意味凄く気持ち悪い問題だと思います(笑)。
まあ、実際はcqzypxさんが既に答えてくれているので、特に何ら書く事もないんですが、ある意味面白い問題だな、と思う部分もあるので、ちょっとしたアクロバットな方法論を紹介しておこう、と考えました。「数学じゃない」部分で解説するのもアリではないか、と(実際僕は数学らしい数学の話を書くことはめったにないのですが・笑)。

さて、ところで、実際問題、多重積分の問題としてはこの問題自体は大した問題ではありません。ただし、「2つの確率変数を扱った確率分布である」って前提が物凄く気持ち悪くなる原因なのです。
通常、多重定積分を行う場合、「積分する領域」を問題が示唆してくれたりするんですが、この問題の場合「それ」がありません。「確率分布である」と言う条件だけが提示されていて、それ以外は何もない。果して「確率分布である」と言う条件だけで積分領域がどう決まるのか、それがイメージしづらいのです。だから「気持ち悪い」。

確かに数学では「イメージする力」とか「想像力」が問われます。ところが、こういった問題の場合、その「イメージを掴む為の」とっかかりが見えない。そう言う理由で「多重積分の問題としては何て事が無い問題」が難問化するんですね。こりゃどーしようもねえ。

そこで。プログラムを書いて図示化する、って道を選びましょう。
普段の勉強に於いて、「どんな形であれ図示化しておく」練習をしておけば、いざと言うときこういう「掴み所が見えない」問題が出た時でも他の人よりは「想像力」を働かせる事が出来るでしょう。
また、実際、任意の二変数関数f(x,y)を用いた

の形の確率分布は「どーにでも捏造できる」モノです。ただし、このテのを専門的には正規化定数とか分配関数とか呼ぶんですが、一般にこのZは解析的に求まるかどうか分かりません。この問題の場合は当然「解析的に求まる」ワケですけど、あくまで「一般的に」その保証があるワケではないのです(ないしは凄く計算がシチメンド臭い、とか)。
こう言う場合、数理統計学者は一体どうしてるのか?と言う観点で、一種「最新技術」を紹介しておこうと思います。完全に「高校数学の」窓としては逸脱していますが(笑)、たまにはイイのではないか、と(笑)。また、理論面はどうあれ、高校生でも１回覚えておけば「バカの一つ覚え」で色々な確率分布を図示できるので、将来的な事を考えても有用な技術となるでしょう。今のうちに慣れておくに越した事ないです。

ところで、勘のイイ人だと、

「亀田と言えばコンピュータ・プログラム・・・・・・ああ、またモンテカルロ法でしょう?」

とか思うかもしれません。当たりです(笑)。ただし、今回は「ちょっと変わった」モンテカルロ法を行ってみます。
以前、モンテカルロ法を紹介した時に、

コンピュータを用いた乱数シミュレーションをモンテカルロ法と呼ぶ

と書きました。が、半分正解で半分違うんです。
実は「乱数を用いる」だけでなくって、「乱数をつくり出す」プログラムも広義では「モンテカルロ法」と呼びます。
つまり実際のトコ、「乱数絡み」のプログラムは全てモンテカルロ法、と呼んでも構わないのです。
この問題の場合、p(x,y)＝k(x+y)と言う「全くよくわからない確率分布」を相手にします。そして一般に、こういう「良く分からない」確率分布は無数に存在するんですね。
一方、そんな「良く分からない」確率分布でも、それが生成する「乱数」さえあれば、その確率分布がどういう形をしてるのか、平均は何なのか、分散はどんなモノなのか、それらの性質を知ることが出来ます。

「当たり前の事言ってるように聞こえるけど・・・・・・その確率分布そのものが"良く分からない"んじゃない?良く分からない確率分布から乱数を生成する?そんな事って出来るの?」

と思うでしょうが、実は出来るんです。
こういう「ちょっと特殊な、乱数を生成する為のモンテカルロ法」のうちの一つをメトロポリス法と呼びます。
今回はこのメトロポリス法を紹介したいと思います。専門書を買えば5,000円近くかかるネタなんですが(笑)、幸いこのサイトはタダなんで(笑)、ここでメトロポリス法を覚えたら絶対お得ですよ(笑)。さあ、お立ち会い(笑)!!!

方針

目的の1番目は、問題にある

p(x,y)=k(x+y)

と言う確率分布からそれに従う確率変数xとyの乱数列をメトロポリス法を用いて発生させる事です。
目的の2番目は、その乱数列xとyを用いて立体的な度数分布表(ヒストグラム)を作ります。これで「どういう形の確率分布なのか?」がハッキリと分かるでしょう。
なお、p(x,y)=k(x+y)に従う2つの乱数列を生成する為のプログラムはいつも通りDrSchemeを用いて記述しますが、それを図示化するには統計解析ソフトRを用います。と言うのも、DrSchemeはグラフィックス関係が弱い。その点、統計解析ソフトRは強力な作図機能を備えているので、データ(この場合はメトロポリス法で作った2つの乱数列)さえあれば、かなりラクにヒストグラムを作る事が出来ます。
一方、確かに全部の工程を統計解析ソフトRで行うことも可能ですが、正直、「プログラミング言語」としては統計解析ソフトRはあまりスマートではありません。少なくとも、個人的な意見では、DrSchemeの方が簡単で、かつ、キレイにプログラムを書く事が可能なので、頻雑さが無いのです。
よって、メンド臭く思うかもしれませんが、この2つのソフトの間で「ファイルをやり取り」した方が結果簡単だ、と言う事です。んで、思ったよりそれはメンド臭くないんで、心配しなくって大丈夫です。
なお、DrSchemeの基本的な使い方に付いては、ココやココに記述があるので、軽く目を通しておいてください。

メトロポリス法の基本

一応最初にメトロポリス法の概要を問題に従って説明したいと思います(こういう雛型提示を疑似コードと呼んだりします)。

1. 0以上1以下の一様乱数を使ってp(x,y)＝k(x+y)に従う確率変数の初期値、を決める。


2. 別の0以上1以下の一様乱数(random)を二つ使って、を計算式に従って作る。


3. 、とする。(この問題設定上、、でなければならず、そうでない場合は7番に移る。)


4. を計算する。


5. また新しくの一様乱数Rを生成する。


6. もし、なら生成した乱数を、として2番に戻って繰り返す。


7. それ以外は、生成した乱数を、として2番に戻って繰り返す。
   

これがメトロポリス法です。実際は、2番〜７番のステップを延々と繰り返せばxの乱数列とyの乱数列が生成されて、それらはp(x,y)＝k(x+y)に従う乱数列となるのです。
注意点をいくつか。
まずはメトロポリス法内で扱う一様乱数は結構数が多いので、どれがどれなのか混乱しないようにしましょう。
また、

「あれ?p(x,y)＝k(x+y)自体はどうなったの?これは生成される乱数列にはならないの?」

と言う疑問を持つ人もいるかもしれません。これは、最初は勘違いするんですが(僕も一時期してました・笑)、その通りです。乱数列となるのはあくまで、xとyの方であって、p(x,y)＝k(x+y)は「枠組を提供する(wを計算する)為だけに存在して」て、「結果自体」にはならないのです。
と言うのも、通常の数学的な枠組では、x、yは独立変数であって、f(x,y)がそれによって計算される「確固とした関数」(従属変数)になることは確かなんですが、一方、確率分布と言うのは、そう言う意味では「数学的な関数」ではないんです。むしろ、「確率変数」と言う言い方をしますが、こっちの方が、ある種「関数」そのもの(計算結果)なんですね。何の関数かって?平たく言えば「運」の関数なんですよ(笑)。
つまり、確率分布p(x,y)＝k(x+y)と言う数式自体は「数学的な関数」を表しているワケでも何でもなくって、「運の関数」である確率変数x、yの「出方」を制御する「枠組」なんですね。実は確率分布の数式自体は「数学的な関数として」の意味は持ってないんです。x、yが決まったからp(x,y)が一意に決まる、なんて構造はそもそも持っていません。
たぶんこの解説読んでも最初はピンと来ないでしょうが(笑)、まあいずれにせよ、「大事な結果」はxとyの方です。この事だけは覚えておいてください。
なおはメトロポリス法のプログラム上は、外部引数として与えるモノとし、固定した値は持たせない事とします。

まずは核(kernel)を作る

さて、取り敢えずは一番簡単な部分から作り始めましょう。
DrSchemeを起動して、上部のエディタに次のように入力します。

(define kernel

(lambda (x y)

(+ x y)))

これは数学的には次の関数を表しています。

これで問題の確率分布の表現は終わりです。
・・・・とか聞いて、

「え、でも、問題の確率分布はp(x,y)＝k(x+y)でしょ?kが消えちゃってるのはマズいんじゃないの?大体kを求めよ、ってのが題意でしょ?」

とか思う人もいるでしょうが、実はメトロポリス法に於いては、正規化定数の値は考慮しなくってイイのです。
「ええ?」って思うかもしれませんが、前出の疑似コードの4番を見てください。問題に従うと、

となって結局、正規化定数は約分されて消えてしまいます。
これはメトロポリス法の大きな特徴なんです。結果メトロポリス法に於いては全く正規化定数は関係無くなってしまうのです。
ちなみに、統計学的には用語が定まってないのですが、こう言った「確率分布の正規化定数以外の関数部分」をカーネル(核)と呼んだりする場合もあり、ここではそのカーネル、と言う単語を用います。
余談ですが、メトロポリス法を利用した実際の実装は殆んど無いのですが、例えば「解析的に積分が出来ない」確率分布にご存知正規分布があります。

この場合もカーネルは

と表現出来て、メトロポリス法を用いて正規分布に従う乱数(正規乱数)を作る事が出来ます。
メトロポリス法を覚えたら、是非とも実装してみてください。
なお、本題に戻ると、カーネルをDrSchemeに記述し終わったら、一旦ファイルを名付けて[保存]しましょう。



ここではファイル名をQ3609と名付けています。
そして、最終的に外部入力を促すQ3609と言う名前のプログラムを走らせて全ての計算をさせる、と言う意図があります。
次はそのQ3609と言うプログラムを書いてみましょう。

ダミーのプログラムQ3609を作る。

次にダミーのプログラムQ3609を作ります。
このプログラムは次の意図があります。

1. 主な計算自体を行うプログラムは別にmetropolisと言う名前を付けて作る。


2. 従ってQ3609は最小限の入力を人間(つまり僕等)から受け取ってそれをmetropolisに渡す役割、とする。


3. また、metropolisが必要とする引数の初期値を渡す。

この3つの役割をQ3609に持たせる事、とします。
なお、ここでQ3609を走らせる際に人間側が入力する引数は、何回繰り返し計算をさせるか指定するnと、との変位の幅を決めるの二つ、とします。
以上より、Q3609のコードは次のようになります。

(define Q3609

(lambda (delta n)

(metropolis (random) (random) delta n ()())))

なお、(random)と言うのは、Scheme上で0〜1の範囲の範囲で生成される一様乱数で、これでx、yの二つの変数の初期値を与えます。
また()ってのを空リストと呼びます。つまり、この二つに生成されたxの乱数とyの乱数を突っ込んで行くのです。これも乱数列の初期値、と言う事ですね(初期値なんで何も中に入っていない、って事です)。
つまり、ここで分かるのは、プログラムmetropolisを走らせる際、その書式は

(metropolis x y delta n xのリスト yのリスト)

となっている、って事です(別にこのプログラム自体を表向きには走らせませんが)。
そして、例えばn＝10,000を与えると、動作としては

0回目:(metropolis delta 10000 () ())
1回目:(metropolis delta 9999 () ())
2回目:(metropolis delta 9998 ( ) ( ))
3回目:(metropolis delta 9997 ( ) ( ))
・・・・・・

となる事を示唆しているんです。目論見通り、計算の回数nが減るに従って、リスト内の乱数の要素が増えていく事が分かりますし、それによって2つの乱数列を得る事が出来るのが分かると思います。
取り敢えず、ダミープログラムQ3609を先ほど作ったプログラムを先ほど作ったエディタ(ファイル)Q3609のkernelプログラムの下に改行して書き込み[保存]しましょう。



なお、ここでScheme上でのリストの基本操作命令であるconsに付いてちょっと解説します。
consと言うのは、リストの先頭に要素を追加する命令です。以下がその書式です。

(cons 要素 リスト)

これにより、リストの先頭に要素を追加する事が出来ます。
例えば、これはDrSchemeの下段のインタプリタで実際実験して欲しいんですが、次のような命令をプロンプトの後に書いてみて、どうなるか[実行]してみてください。

(cons 1 ())
(cons 2 '(1))
(cons 3 '(2 1))

はい、ご覧の通り、先頭に新しい要素を付け足したリストが作られますね。

註:なお、'はScheme用語でクオートと呼び、「これは単なる数字じゃないですよ」と言うサインを表しています。要するに、これで「リストです」と表現しているワケです。

このconsと言う命令が計算の本体であるプログラム、metropolisのある意味鍵を握ります。

メトロポリス法の実行エンジン、metropolisを作る

さて、いよいよ計算の本体、プログラムmetropolisを書いてみましょう。
プログラムmetoropolisも再帰と言う枠組で作ります。
ここでもう一度再帰の枠組を説明しておきましょう。

1. n＝0の時、metropolis法の実行結果を返す。


2. nが任意の自然数kの時、n＝k−1の時のmetropolis法の実行結果に対し、n＝kの時のmetropolis法の実行結果を適用する。

これだけ、です。
もう一度言いますが、これは数学的帰納法(ないしは数列/漸化式)の枠組とまるっきり同じだ、と言うことです。

1. P(0)は真である。


2. 任意の自然数kに対し、P(k)が真であれば、P(k+1)も真である。

全く同じ構造ですね。
つまり、こういう再帰、と言うプログラム法を学ぶ事により、数学的帰納法/数列/漸化式が苦手な人が結構多い(僕もそうですが・笑)んですが、「ああ、こういう考え方って使えるよな」と感じてほしい。そのテの数学的な考え方に慣れる為にはSchemeと言うプログラミング言語は絶好の武器なんですよ。高校生〜大学生の為の数学の自習向きのプログラミング言語、と言えます。
さて、再帰の枠組の1番に従って、「一体計算結果として何を返すのか?」考えます。取り敢えず暫定的にxの乱数列のリストとyの乱数列のリストの2つを要素としたリストを返すように設定してみましょうか。つまり「リストのリスト」を返すように設計してみましょう。
すなわち、

n＝0の時、(list xの乱数列のリスト yの乱数列のリスト)を返す

と取り敢えず置いてみます。なお、listと言う命令はリストを作る命令です。
あとは再帰の枠組の2番を考慮して、k番目のメトロポリス法をk−1番目のメトロポリス法の計算結果を用いて定義してやれば一丁上がりです。
以上を考慮して、プログラムmetropolisは次のようになります。

(define metropolis

(lambda (x y delta n lsx lsy)

(if (zero? n)

(list lsx lsy)
(let* ((xtrial (+ x (* delta (- (* 2 (random)) 1))))

(ytrial (+ y (* delta (- (* 2 (random)) 1))))
(w (/ (kernel xtrial ytrial) (kernel x y)))
(R (random)))
(metropolis (aor x xtrial w R)

(aor y ytrial w R)
delta
(- n 1)
(cons (aor x xtrial w R) lsx)
(cons (aor y ytrial w R) lsy))))))

意外に簡単なコードになったな、と思わないでしょうか?(まあ、そう感じなくても構いませんが・笑。しかし、他の言語だったらもっと複雑な形になったりします。)
取り敢えず、ここまでエディタ(Q3609ファイル)に書き込んで、そして[保存]しましょう。



なお、ある意味モンテカルロ法ってのはワンパターンなんで、一つ様式を覚えてしまえば他の問題にも簡単に転用出来ます。
が、ここで2つだけプログラムmetropolisに関してポイントを解説します。
一つ目は局所変数let*です。これは基本的には＜質問3590＞に出てきた局所変数letと同じモノです。これで定義された変数はプログラムmetropolisが1回実行されて終了するまでの範囲だけで有効なんです。つまり、1回metropolisが起動終了してしまえば、そこで定義された変数は廃棄されてしまいます。
「それが何か?」と思うかもしれませんが、こういう事です。つまり、k−1回目のメトロポリス法で使われた局所変数とk回目のメトロポリス法で使われる局所変数は同じ数ではない、と言う事です。これが意味するのは、「局所変数」と言うのは繰り返し計算の過程に於いて「どんどん書き換わっている」んですね。これが大事な性質なんです。
ここで定義されている局所変数は、数学的には次の4つです。

この四つはメトロポリス法の基本を鑑みれば分かると思いますが、「1回の計算が終わる度」書き換わらなければならない変数です。従ってlet*を使って局所変数として定義してるんですね。プログラム上の表記法と数学的な表記法は違いますが(例えばはmetropolis内ではxtrialと表記されています)、let*内の局所変数の定義と良く見比べて下さい。「全く同じ事を表現している」と言うのが分かるはずです。
なお、letとlet*の違い、と言うのは、たった一つで、letの場合は、局所変数内で定義された数を使って他の局所変数を定義出来ないのですが、let*の場合はそれが出来る、と言う部分です。
上の4つの式を見てもらいたいんですが、単純にwを定義するにあたって、別の局所変数とを参照しています。これはletに於いては許されていないんで、代わりにlet*を使っている、と。それだけなんです。
二つ目のポイントはプログラムaorです。k−1回目のメトロポリス法に「何かした結果」は、このプログラムaorを呼び出す事によって、

(metropolis (aor x xtrial w R) (aor y ytrial w R) delta (- n 1) (cons (aor x xtrial w R) lsx) (cons (aor y ytrial w R) lsy)

と記述することが出来るワケですが(consの使い方にも着目!)、実はこのプログラムaorはまだ定義していません(笑)。
よって、今からこのプログラムaorを記述します。

プログラムaorを作る

AORと聞くと、音楽が好きな古い世代の人たちは

「Adult Oriented Rock???」

とか思うでしょうが(古い話ですんません・笑)、全然関係ありません(笑)。いや、昔は(８０年代初頭くらい?)湘南辺りのサーファーとかAORってジャンルの音楽を良く聴いてたらしいです(笑)。
ええと、別にプログラムの名前をどう付けようが構わないんですが(笑)、ここでのaorの意味はaccept or rejectの省略形です。accept(受諾)かreject(棄却)か、って意味ですね。
メトロポリス法の基本の６番、７番にあたる仕事ですが、結局xであろうとyであろうと、acceptとrejectの仕事そのものはどっちの変数だろうと変わりません。よって、ここで汎用のプログラムaorを書いてしまおう、と言うのが狙いです。
ここではx、y、xtrial、ytrialに使える汎用の変数をz、ztrialと名付けて、次の事柄をプログラムaorに翻訳します。

1. かつ、の時はを使用。


2. それ以外はを使用。

これをif構文を利用して次のように書きます。

(define aor

(lambda (z ztrial w R)

(if (and (<= R w) (<= 0 ztrial 1))
ztrial
z)))

これも単なる基本的なif構文ですが、ポイントを二つほど。
まずは関数andですが、これは

(and 命題1 命題2 命題3 ....)

と複数の命題を並べ、全ての命題が真だった時、「真」を返します。そんなに難しくありませんね。
ここでは(<= R w)と言う命題と(<= 0 ztrial 1)と言う命題の2つが同時に真になるかどうか判定する役割を担っています。 二つ目のポイントは<=と言う命令です。これは「以下」かどうか調べる役割を担っています。 通常のプログラミング言語ですと、0以上１以下なんかを一行で記述するのは難しいのですが、生憎Schemeは前置記法を採用している為、(<= 0 ztrial 1)なんて書き方でを確かめる事も出来ますし、また、もっと複数の数の大小関係を調べるのもたった１行で書けるんです。
さて、エディタ(Q3609のファイル)にプログラムaorを書き込んで[保存]しましょう。



これで取り敢えず今回の問題に準拠したメトロポリス法の雛型は完成致しました。

メトロポリス法のソースコードを見てみる

今までの工程は初心者に対しては複雑に見えるでしょうから、一旦、プログラムQ3609の全体(ソースコード)をもう一回見てみましょう。

(define kernel

(lambda (x y)

(+ x y)))

(define Q3609

(lambda (delta n)

(metropolis (random) (random) delta n ()())))

(define metropolis

(lambda (x y delta n lsx lsy)

(if (zero? n)

(list lsx lsy)
(let* ((xtrial (+ x (* delta (- (* 2 (random)) 1))))

(ytrial (+ y (* delta (- (* 2 (random)) 1))))
(w (/ (kernel xtrial ytrial) (kernel x y)))
(R (random)))
(metropolis (aor x xtrial w R)

(aor y ytrial w R)
delta
(- n 1)
(cons (aor x xtrial w R) lsx)
(cons (aor y ytrial w R) lsy))))))

(define aor

(lambda (z ztrial w R)

(if (and (<= R w) (<= 0 ztrial 1))
ztrial
z)))

基本的にはプログラムQ3609は4つのパーツによって構成されています。
そしてこれらは、

Q3609→呼び出し→metropolis→呼び出し→kernel

↓
呼び出し
↓
aor

と言う依存関係になっています。
これら4つのパーツを同一ファイルに書き込む事によって、特にこちら側で何か設定しなくても、お互いが必要な時に必要なモノとして呼び出されて一つのプログラムのように動いてくれるのです。
では、実際にプログラムQ3609を走らせてみましょう。
[実行]ボタンを押してください。

Q3609を走らせてみる

Q3609の実行方法は、DrSchemeの下段のインタプリタのプロンプト以降に、Schemeの流儀に従って前置記法で次のように入力して[Enter]キーを叩きます。

(Q3609 delta 試行回数)

ここで、の値をどうするか、ですが、基本的には何でも構わないんですが、今回の問題のx,yの変位の幅「1」の大体三分の一から二分の一辺りが適当だと思います。すなわち0.3〜0.5くらいを与えれば良い。ここではくらいにしてみましょう。
また試行回数ですが、これも自由にして頂いて結構です。大体10,000回も走らせればイイのではないでしょうか?
取り敢えず亀田は100万回くらい走らせてみます。つまり、命令は

(Q3609 0.4 1000000)

としてみます。
以下の写真が実行結果ですね。



100万回もループさせるとさすがに実行時間がかかり、途中「止めときゃ良かった」とか思いましたが(笑)、一応計算は終了しました。
基本的に

((100万個のxの乱数列) (100万個のyの乱数列))

と言う形で出力されています。
が。色々乱数列を眺めて楽しむのもイイんですが、100万個×2列の乱数をただ眺めていたって何がどうなる、ってワケではありません。仮に1万個×2だとしても、人間が分析するには当たり前ですが、量が多すぎるのです。
そこで。こういう「メトロポリス法」で行った計算結果をファイルとして出力し、他の画像描画ソフトに渡してそこで作図させる作戦を取ります。つまり、今まで書いたQ3609のプログラムを若干手直しして改造します。
その前に「どう言ったファイルを作成させればいいのか?」ちょっと説明したいと思います。

CSVファイルと言う形式

ちなみに、基本的にパソコン上ではどんなデータでも原則文字だらけなんです。音声データだろうが、画像データだろうか、動画データだろうが、実は元々文字だらけのファイルなんですね。
「うそ!?」って思うかもしれませんが、原理的にはそうなんです。ただし、その「文字列」を記入したテキストファイルを「固有のフォーマット」に従って作成して、「適切な拡張子」を付加して定義する。例えば最近有名な音声データであるmp3ファイルなんかも原理的にはそうやって作られています。
さて、そんな中で純粋なテキストファイルとしてはこれまた色々なデータの記述方式があるんですが、その中でメジャーなモノにCSVファイルと呼ばれる形式があります。この形式のファイルは例えばMicrosoft Excelなんかでも表計算用データとして開く事も可能ですし、と言うことはMicrosoft Excelに読み込ませて作図させる事も可能です。
もっとも今回は統計解析ソフトRを用いる、と言う前提なんで、Microsoft Excelは用いません。大体、ちょっと多めのデータ(例えば1万件以上のデータ)を扱うとなるとMicrosoft Excelだと簡単にバグりますし、また描画も困難になるから、です。まあ、興味があったら、あとで自分で実験してみてください。
ところでCSVファイル形式、とは言っても実はそんなに作成が難しいファイル形式ではありません。以上の3点を満たせばCSVファイル、ってのは簡単に作成出来るんです。

1. ファイルの拡張子は.csvとする。


2. 列単位のデータの区切りはコンマ(,)を用いる


3. 行単位のデータの区切りは改行で表す

これだけ、です。意外に大した事ないでしょ(笑)?まあ、それだからこそ色々なソフトの間で「データの互換性を取るために」使われているワケですが。
具体例を見てみます。CSV形式では次のようなデータの並べ方をします。

,
,
,
,
......
,

このようにコンマで区切って対応するデータ同士を並べていけば良い。
これで乱数列をどうやってCSV形式に変換するか大体のトコ分かりました。
では実際に出力用のSchemeプログラムを書いてみましょう。

出力用プログラムwrite-lsを作る

取り敢えず、出力用プログラムのwrite-lsがどうなるか、書いてみましょう。

(define write-ls

(lambda (lsx lsy)

(if (or (null? lsx) (null? lsy));1、2参照

(newline);3参照
(and (display (car lsx));4、5参照

(display ",")
(display (car lsy))
(newline)
(write-ls (cdr lsx) (cdr lsy))))));6参照

実はこのプログラムwrite-ls自体も再帰構造で書かれています。
基本的なアイディアは以下の通りです。

1. xの乱数列、またはyの乱数列が無くなった時、改行して終了する。


2. それ以外では、xの乱数列の先頭を印字、コンマ(,)を印字、yの乱数列の先頭を印字、改行してwrite-lsをそれぞれの乱数列リストの残りに対して適用する

これだけ、です。
ではポイントをいくつか解説します。

1. OR関数:
   書式は
   

   
   (or 命題1 命題2 命題3 ...)
   

   
   と複数命題を並べる事が出来ます。命題のうちのどれかが真だった場合、真を返します。
   ここでは(null? lsx)と(null? lsy)と言う命題のうちのどちらかが真かどうか調べています。


2. null?関数:
   リストが空リストかどうか調べる関数です。空リストであった場合真を返します。
   書式は
   

   
   (null? 判定したいリスト)
   

   
   となります。
   

   
   例:(null? ())→真
   (null? '(1 2 3))→偽
   

   
   


3. newline関数:
   改行するための命令です。書式は
   

   
   (newline)
   

   
   です。


4. display関数:
   出力/表示の為の関数です。書式は
   

   
   (display 表示させたいモノ)
   

   
   です。「表示させたいモノ」は計算結果だろうと文字だろうと構いませんが、文字を表示する場合はダブルクオテーションマーク(")で挟みます。
   

   
   例:

   
   
   • (display (+ 2 3))→「5」を表示
   
   
   • (display "Hello!")→「Hello!」と表示
   
   



5. car関数:
   リスト基本操作関数のうちの一つ。car関数はリストの先頭要素を返します。
   書式は
   

   
   (car 先頭要素を返したいリスト)
   

   
   です。
   

   
   例:(car '(1 2 3))→「1」を返します。
   



6. cdr関数:
   リスト基本操作関数のうちの一つ。cdr関数はリストの先頭要素「以外の」リストを返します。
   書式は
   

   
   (cdr 先頭要素以外を返したいリスト)
   

   
   です。
   

   
   例:(cdr '(1 2 3))→(2 3)を返します。
   

これが出力用のプログラム、write-lsの全てです。
ではこれを用いてプログラムQ3609の改造を行います。
取り敢えずエディタ(Q3609のファイル)にプログラムwrite-lsを書き込んで[保存]しましょう。

プログラムQ3609をCSVテキストファイルを出力出来るように改造

まずは、もう一回全ソースコードを見て改造されたポイントを書いていきます。
なおセミコロン(;)以降はコメントと言ってプログラム本体の実行には何の影響もありません。

(define kernel

(lambda (x y)

(+ x y)))

(define Q3609

(lambda (fname delta n);引数fnameを追加

(with-output-to-file fname;新しい関数追加

(lambda ();ここの引数は無し

(metropolis (random) (random) delta n ()())))))

(define metropolis

(lambda (x y delta n lsx lsy)

(if (zero? n)

(write-ls lsx lsy);list関数の代わりに出力用関数write-lsを採用
(let* ((xtrial (+ x (* delta (- (* 2 (random)) 1))))

(ytrial (+ y (* delta (- (* 2 (random)) 1))))
(w (/ (kernel xtrial ytrial) (kernel x y)))
(R (random)))
(metropolis (aor x xtrial w R)

(aor y ytrial w R)
delta
(- n 1)
(cons (aor x xtrial w R) lsx)
(cons (aor y ytrial w R) lsy))))))

(define aor

(lambda (z ztrial w R)

(if (and (<= R w) (<= 0 ztrial 1))
ztrial
z)))

(define write-ls

(lambda (lsx lsy)

(if (or (null? lsx) (null? lsy))

(newline)
(and (display (car lsx))

(display ",")
(display (car lsy))
(newline)
(write-ls (cdr lsx) (cdr lsy))))))

基本的に改造するプログラムはQ3609とmetropolisの二箇所だけ、です。
まず、最初のポイントとしては、Q3690が外部から受け取る引数はfname、、試行回数nの3つとなり1つ増えています。なお、fnameってのはfilenameの省略形です。Q3609は出力するCSVファイルの名前を入力するように改造されたのです。
次のポイントは、最初書いたmetropolisでは最終的には(list lsx lsy)と言う命令によって「2つの乱数列リストのリスト」を返すように設計していました。これを先ほど作った新しいファイル出力用関数write-lsを用いて改造しています。今度はリストとして出力するのではなく、「CSVファイルとして」出力するのが目的だから、です。
最後のポイントがwith-output-to-file関数です。これがプログラムwrite-lsの出力をファイルに書き込む役割を担っています。書式は

(with-output-to-file ファイル名 (lambda () 作業内容))

となり、これをダミープログラムQ3609に組み込ませたわけです。

現時点では次のような作成したプログラム同士の呼び出しが行われています。

write-ls
↑
呼び出し
↑
Q3609→呼び出し→metropolis→呼び出し→kernel

↓
呼び出し
↓
aor

これでQ3609でCSVファイルを生成/出力が出来るようになりました。
では、一旦ファイルを[保存]して[実行]ボタンを押してみましょう。

Q3609を走らせてCSVファイルを生成する

では、Q3609を走らせて乱数データを記述したCSVファイルを実際作成してみましょう。
引数が一つ増えているので、Q3609の実行コマンドは以下のようになります。

(Q3609 出力ファイル名 delta n)

なお、注意事項としては、「出力ファイル名」は文字列データなので、ファイル名はダブルクオテーション(")で挟まなければならない、と言う約束事があります。
また、必ず拡張子(.csv)を名前の後に付けるのを忘れないようにしましょう。
ここでは出力ファイル名を"data.csv"と名付けて、、試行回数n＝100万回として次のようにDrScheme下段のインタプリタのプロンプトに入力して[Enter]キーを打ちます。

(Q3609 "data.csv" 0.4 1000000)

では結果を見てみましょう。



・・・・・・何も表示されませんね(笑)。
「失敗したか?」と思うかもしれませんが、いや、計算/ファイル出力は成功してるんで心配しないでください(笑)。キチンとプロンプト(>)が一つ増えています(これは計算が終わった事の証)。
また、「インタプリタ上に」計算結果を羅列するより早く動作を終えています。実はDrSchemeの場合、計算結果をインタプリタで表示する方が時間がかかるのです。そしてこれが意味するのは、実は先ほど実験した「100万回での試行実験」でも、計算そのものは比較的すぐ終わるんですが、「インタプリタ上への計算結果の羅列」の方にバカみたいに時間が取られている、って事です。意外なんですがテキストファイルへの出力の方が早く作業が終わるんですね。
では、生成したファイルをチェックしてみましょう。Windows版のDrSchemeの場合、デフォルトではファイルの生成先はユーザーフォルダ内の自分のフォルダ内のドキュメントフォルダ内に生成されていると思います。



パスで書くと

C:\Users\******\Documents\data.csv

となってる筈です(******はユーザー名)。
Linux版の場合は、homeディレクトリ内の自分のディレクトリにCSVファイルが生成されていると思います。



ではそのCSVファイル(data.csv)をマウスで右クリックして、テキストエディタ(Windows版ではメモ帳かWordpad)で開いてみましょう。

註:恐らくダブルクリックで開こうとすると、WindowsではMicrosoft Excelを使って開こうとする筈です。それでも構いませんが、今回亀田が試したような100万回での試行ですと、完全にMicrosoft Excelが扱えるデータ量の範疇を越えてしまうので、どのみち単なるテキストエディタで開いた方が無難です。

ご覧のように"data.csv"をテキストエディタで開くと、xの乱数列とyの乱数列がコンマ(,)で区切られた状態で2列に渡って記入されている事が分かります。成功です。
以上でDrSchemeでの作業は終わりです。
ここから統計解析ソフトRに作業が移ります。

統計解析ソフトRのダウンロードとインストール

Rもタダで手に入るオープンソースの統計解析用ソフトウェアです。強力な統計関数とグラフィック関係のデバイスが含まれているので、この機会にダウンロード/インストールしてみましょう。

ダウンロード&インストールの方法:

• Windows 版 R のインストール


• Mac 版 R のインストール


• Linux 版 R のインストール
  
  
  
  • Debian GNU/Linux (含む Knoppix)の場合 (インターネット経由でオリジナル版をインストールする場合)
    
    註:Ubuntu/Kubuntuもここに含まれます。ただし、KNOPPIX/MathはRをプリインストール済です。
  
  
  • Vine Linux 版インストール
  
  
  • openSUSE10.2 Linux 版インストール
  
  

※:Rの最新版は2.6.0です(2007年10月5日現在)。

それぞれのOS用の指示に従ってダウンロード/インストール/初期設定を行ってください。

Rの起動

Windows版の場合はダウンロード/インストールが終わるとデスクトップ上にRのアイコンが出来てる筈なので、それをダブルクリックすればRが立ち上がります。





Linux版の場合は、まずは端末(ターミナル)を起動します。





コマンドプロンプトにただRと大文字で打ち込むだけでRが端末上で起動します。

パッケージのインストール

Rの優秀なトコロは、数々の特殊な機能をまとめたプログラム群を追加インストール可能である、と言う部分です。これをパッケージと呼びます。
今回利用するパッケージはRcmdrと言うパッケージとgregmiscと呼ばれるパッケージの2種類です。
特にRcmdrと言うパッケージはコマンドラインバリバリのRにマウス操作の為のGUI環境を加える為のモノでこれは結構重宝します。詳しくはこのPDFに書いてあるので、参考にして下さい。マウスだけで基本的なRの操作が出来るようになります。

Windows版Rでのパッケージのインストール方法は、RGui上部にある[パッケージ]プルダウンメニューから[パッケージのインストール...]を選びます。



そうしたら、「どこのミラーサイトからダウンロードするのか?」尋ねてくる[CRAN mirror]ポップアップが出て来ます。



現時点では日本のRのダウンロード先は

• 福島県:福島県立会津大学 Japan(Aizu)


• 東京都:国立東京大学 Japan(Tokyo)


• 茨城県:国立筑波大学 Japan(Tsukuba)

の3箇所となっています。近場の大学のサイトをクリックして[OK]ボタンを押しましょう。



そうすると、インストール可能なパッケージのリスト、[Packages]ポップアップが出てきます。スクロールして、目的のパッケージを探し出してクリック、そして[OK]ボタンを押してください。



無事パッケージのダウンロードが終わります。

Linux版Rでのパッケージのインストール方法は、端末(ターミナル)で次のようなコマンドを入力します。


install.packages("ダウンロードしたいパッケージ")


例えば、今回Rcmdrとgregmiscが必要なワケですから、次のように端末のコマンドプロンプトに入力します。


install.packages("Rcmdr")
install.packages("gregmisc")


すると、R上で次のようにパッケージのダウンロードが始まります。

註:Ubuntu/Kubuntuの場合は、R経由でインストールするより、Synapticパッケージ・マネージャないしはAdeptパッケージ・マネージャ経由でインストールした方が確実かもしれません。

Rcmdrの起動方法

Rcmdr(Rコマンダー)はWindows版では、まずはRGui上部にある[パッケージ]プルダウンメニューから[パッケージの読み込み...]を選びます。



そうすると[1つを選択してください]ポップアップが出てきますので、パッケージの一覧から[Rcmdr]を探し出してクリックして[OK]ボタンを押します。



それにより、Rコマンダーが以下のように起動します。



以降は、基本的に亀田はLinuxユーザーなので、Linux版の画面写真を元に解説していきますが、Rコマンダー及びRはLinux版でもWindows版でも全く同じように操作出来るので、ご心配は無用です。

Rcmdr(Rコマンダー)はLinux版では、コマンドラインで起動します。
端末(ターミナル)で次のように打って、Rcmdrを起動します。


library(Rcmdr)


そうするとRコマンダーが起動します。



これを使えば基本的なRの操作を全てマウスで行う事が出来ます。
ではいよいよDrSchemeのメトロポリス法で作成したCSVファイルをRコマンダーを使ってRに読み込んでみましょう。

Rコマンダーでデータのインポート

RコマンダーでデータをインポートするにはRコマンダー上部の[データ]プルダウンメニューから[データのインポート]→[テキストファイルまたはクリップボードから]を選びます。



そうすると[テキストファイルまたはクリップボードからデータを読み込む]ポップアップが表れます。



ここでは次のように設定をします。

	データセット名を入力:好きな名前を入力して構いません。ここではデフォルトのまま"Dataset"とします。 ファイル内に変数名あり:ファイル内に変数名は付けていないので、このチェックボックスは外します。 フィールドの区切り記号:CSVファイルはフィールドをコンマ(,)で区切っているのでカンマを選びます。 以上の設定が終わったら[OK]ボタンを押します。
	[Open]ポップアップが開くので、先ほどDrSchemeで作ったCSVファイル(data.csv)を捜し出して指定し、[Open]ボタンを押します。
	そうすると、data.csvがRに読み込まれます。

では、ホントにdata.csvがRに無事読み込まれたかどうか確認してみましょう。
Rコマンダーの上部に[データセットを表示]と言うボタンがあります。それをクリックしてみます。



そうすると上の図のようなポップアップが表れます。キチンと読み込まれていますね。こうなったらもう鬼の首を取ったも同然です(笑)。
ちなみにxの乱数列はV1と言う列に格納されていて、yの乱数列はV2と言う列に格納されているのを確認しておいて下さい。
ではポップアップを閉じてください。次の作業へと移ります。

R上でのデータの扱い方

Rでのデータは基本的にはベクトルとして扱います。
今ここに、DatasetのV1列にxの乱数列が存在し、そしてV2列にyの乱数列が存在することが分かっています。取り敢えずこれらをxを言うデータベクトルとyと言うデータベクトルにそれぞれ格納しなければなりません。
ここで初めてRコマンダー上でコマンドラインを用います。Rコマンダーの上段のスクリプトウィンドウに次のように記述して下さい。


x <- c(Dataset\$V1)
y <- c(Dataset\$V2)


これはxとyにDatasetのV1とV2をベクトルとして代入せよ、と言った様な意味です。
今回は細かく説明しませんが、基本的にRでは、代入をするときは＝記号の代わりに<-を用います。また、ベクトルをc()と言う書き方をします。そして何々と言うデータセットの何々と言うデータ、と表現する場合、上のようにドル記号(\$)で「〜の」を表現します。

そして、上図のようにコマンドをスクリプトウィンドウに書いた後、それらのコマンドをマウスで選択し、右クリックボタンを押した後、[実行]を選びます。



そうすると、下段の[出力ウィンドウ]に「命令を実行したよ」と言う意味で入力→実行したコマンドが表記されます。
気になるんでしたら、xないしはyの文字だけ選択して同じように右クリックで[実行]してみてください。xデータベクトルの中身やyデータベクトルの中身が下の[出力ウィンドウ]にズラズラ出力される筈です。

Rコマンダーでgregmiscパッケージを読み込ませる

まずは3D描画をする前に下準備です。
先ほどダウンロードしておいたgregmiscパッケージをRコマンダー上で有効にします。
Rコマンダー上部に[ツール]プルダウンメニューがあるので、それから[パッケージのロード]を選んで下さい。



そうすると、下の図のように[パッケージのロード]ポップアップが表れるので、そのパッケージリストの中からgregmiscを捜し出して、選択、[OK]ボタンを押してください。



これでヒストグラムの3D描画の為の準備は完了です。

3Dのヒストグラムを描く

ここからは舟尾暢男さんのR-Tips辺りからコードをパクって来ます(笑)。
次のコードをコピーして、Rコマンダーの上部の[スクリプトウィンドウ]にペーストし、そしてそのコードをマウスで選択→右クリックで実行して下さい。

  h2d <- hist2d(x, y, show=FALSE, same.scale=TRUE, nbins=c(20,30))
  persp( h2d\$x, h2d\$y, h2d\$counts,
         ticktype="detailed", theta=60, phi=30,
         expand=0.5, shade=0.5, col="cyan", ltheta=-30)

そうしたら3Dのヒストグラムが生成される筈です。



はい、こういう感じの度数分布になってるんですね。これがメトロポリス法で作成した乱数列2つを用いたヒストグラムです。(まだデコボコしてますが、極限では平面になりそうだ、って事まで分かります!)
最初に「多重積分の問題としては非常に簡単だ」と書きました。ただし、積分領域がイメージしづらいだろう、と。
んで、多分AAZさんも実際に色々な積分領域を試したんでしょうが、「どれもピンと来なかった」ってのがホントのトコだと思います。
実はメトロポリス法に依る実験ですと、要するに、と言う「凄くオーソドックスな」積分領域で構わないのです。そしてその多重積分が確率分布の性質、を満たすとすると、kは1にならざるを得ないのです。お分かりでしょうか?
多分直感的にはx+yの分布は「一様分布みたいになるんじゃないか?」と思ってたんでしょう。自分の直観と相反する計算結果だったのが戸惑いの原因、「気持ち悪さの原因」だったのではないか?
しかしながら、メトロポリス法でk(x+y)の分布を見る限り、「一様分布」と言う直感は否定されて、実はx=1、y=1に近づけば近づく程、「その目の出現の確率」が高くなっていくんです。逆にx=0、y=0に近づく程、「その目の出現の確率」は低くなっていくのです。面白いですね(笑)。
これがcqzypxさんが書かれた回答のコンピュータ実験的な図示化なのです。

周辺分布を見てみる

では次に周辺分布を見てみましょう。とは言っても考え方は3Dのヒストグラムより簡単です。
単純に言うと、の周辺分布とは、データベクトルxだけのヒストグラムの事を指します。
また、も同じく、データベクトルyだけのヒストグラムですね。
Rコマンダーで単純なヒストグラムを作るには、上部の[グラフ]プルダウンメニューから[ヒストグラム]を選べばイイだけです。



そうすると[ヒストグラム]ポップアップが表れます。



[変数]に表示されているV1(xの乱数列)、V2(yの乱数列)から任意の変数を選んで下さい。
また、[軸の尺度]もどれを選んでも構いません。ここではデフォルトの頻度の計算を選んでみましょう。

	xの周辺分布(ヒストグラム)
	yの周辺分布(ヒストグラム)

まあ、大して変わり映えのしないグラフが二つ出来ただけ、ですが(笑)。
ところで、周辺分布の「頻度の高さ」と言うのは、一般に3Dのヒストグラムより高くなります。
と言うのも当然で、「周辺分布を作る為の積分操作」と言うのは、例えばxの周辺分布を作るには、あるxの値に従ってy軸方向に散らばった「全てのデータ」を足しあげる、って事です。
また、yの周辺分布に付いてもそうですよね。x軸方向に散らばったあるyに従ったデータを「全て足しあわせて」いるのです。
こういう「積分操作の意味」と言うのも図解の方が分かり易いと思います。
念の為に密度表示でもグラフを作ってみましょうか。

	xの周辺分布(密度表示)
	yの周辺分布(密度表示)

これらも極限値の状態で滑らかになった場合、cqzypxさんが書かれた回答に近くなりそうだな、と言うことが分かります。特に切片辺りの値、そしてその最大値が大体一致しそうな事が分かりますね。

こんなカンジで現代の統計学者も「良く分からない確率分布」「計算が大変そうな確率分布」を見かけた場合、ここで行ったようなメトロポリス法でそれらの確率分布に従う乱数を生成して、それらの性質を研究しているそうです。
色々な確率分布に従う乱数を是非ともメトロポリス法で作成して、実験してみてください。


R-Tips