2014/08/01
データいろいろ
 >  Octava で中心極限定理
事情で少し MATLAB を使うことになったが自分の PC にはない。文法をなるべく忘れないよう、互換性のある Octava で適当にテストしてみた。当初は GUI がある FreeMat を試したが、自分の環境(2014/07/28 を参照)では基本的な関数の hist を実行すると落ちてしまった。

↓ から Windows 用の Octava をダウンロードする。最新の 3.6.4 をクリック。

http://sourceforge.net/projects/octave/files/Octave%20Windows%20binaries/


↓ インストーラが 71.5 MB ある。時間がかかったが無事ダウンロードできた。


↓ インストーラは昔 foobar2000 で見かけたのと同じ形式。大半の画面で設定項目があるが、全てデフォルトで進み、問題なくインストールできた。








↓ 普通のスタートメニューができるのでクリックすると、いきなりコマンドプロンプトだけが出る。Octava に GUI がないと知らない人は、エラーか何かと勘違いするかも。



以下、中心極限定理の簡単なシミュレーションのコードと動画。まず一様乱数で1〜10個の確率変数を取って平均を出し、そのヒストグラムを描くと次第に正規分布に近づく。確率変数が2個だと三角形型の分布。環境によって動画を再生できないと思うので、リンクも置く。動画埋め込みの詳細は 2014/07/16 を参照。最後まで再生したら初期画面に戻ってしまう。そのまま止めたいが、今のところ方法が見つからない。

あと、正確にシミュレーションというなら平均を出す所で再び rand を使わず、直前に作った乱数から選ぶべき。R の sample 関数に当たるものが分かったら書き直す。

■ Ex.1
https://kenpg.up.seesaa.net/image/20140801_oct_1.flv
% 一様乱数
n = 100000 ;
u = linspace(0, 1, 100) ;
for Nsum = 1 : 10
h1 = hist(rand(1, n), u) ;
h2 = hist(sum(rand(Nsum, n)) / Nsum, u) ;
plot(u, h1, u, h2)
pause(2)
end



↓ 一様乱数どうしを掛け合わせると -log(x) の分布になる。詳細は下記を参照。この分布から、先と同様に1〜10個の確率変数を取るとやっぱり正規分布に近づく。

■ wrong, rogue and log : 一様乱数の積は一様乱数にはならない
http://kashino.exblog.jp/12885005/

■ Ex.2
https://kenpg.up.seesaa.net/image/20140801_oct_2.flv
% 一様乱数×一様乱数
n = 100000 ;
u = linspace(0, 1, 100) ;
for Nsum = 1 : 10
h1 = hist(rand(1, n) .* rand(1, n), u) ;
h2 = hist(sum(rand(Nsum, n) .* rand(Nsum, n)) / Nsum, u) ;
plot(u, h1, u, h2)
pause(2)
end



↓ 一様乱数に一様乱数の累乗をすると指数関数に似た右上がりの分布になるが、詳細は不明。ゼロ付近は何か違う。ここから1〜10個の確率変数を取ると、正規分布に近づくようで、でも微妙に右に寄っている気もする。これは今後の宿題。

■ Ex.3
https://kenpg.up.seesaa.net/image/20140801_oct_3.flv
% 一様乱数 ^ 一様乱数
n = 100000 ;
u = linspace(0, 1, 100) ;
for Nsum = 1 : 10
h1 = hist(rand(1, n) .^ rand(1, n), u) ;
h2 = hist(sum(rand(Nsum, n) .^ rand(Nsum, n)) / Nsum, u) ;
plot(u, h1, u, h2)
pause(2)
end



↓ Ex.2 の平方根を取ると、左に歪んだ釣鐘型になる。ここから1〜10個の確率変数を取ると、すごく尖った正規分布になる。

■ Ex.4
https://kenpg.up.seesaa.net/image/20140801_oct_4.flv
% 一様乱数×一様乱数の平方根
n = 100000 ;
u = linspace(0, 1, 100) ;
for Nsum = 1 : 10
h1 = hist(sqrt(rand(1, n) .* rand(1, n)), u) ;
h2 = hist(sqrt(sum(rand(Nsum, n) .* rand(Nsum, n))) / Nsum, u) ;
plot(u, h1, u, h2)
pause(2)
end

<< Word 2007 でラベル付き矢印(1)
ベクトルを円周上に並べる >>