nishiru3の日記

備忘録です。ネットのゴミ。

二変数正規分布

二変数以上のものも必要になったので備忘録です。

一変数の場合は、平均と分散があれば大丈夫ですが、
二変数の場合は、各変数の平均ベクトルと分散共分散行列が必要になります。

yが変数ベクトル、\muが平均ベクトル、\Sigmaが分散共分散ベクトルです。
グラフは、平均(0,0)分散共分散行列が対角行列で成分が各1とした場合です。

f(y;\mu,\Sigma)=\frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}}\exp \left\{-\frac{1}{2}(y-\mu)^T \Sigma^{-1} ({y}-{\mu})\right\}

グラフはGraphRを使わせていただきました。

f:id:nishiru3:20141029123101p:plain

ボックスミューラーのプログラム

ボックスミューラーのプログラムです。
アルゴリズムは下記の記事です。


ボックスミュラーのアルゴリズム - nishiru3の日記

平均、分散、個数を入力し、
発生させ、標準出力するようにしています。

ボックスミューラーで標準正規分布を作成し、線形変換より、
指定した平均と分散になるようにしています。

出力値のヒストグラムを作れば、妥当性をチェックできると思います。

C/C++で書いてます。
本当であれば、ポインタを返すようにすれば良いですが・・・。

#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;
const int pi = 3.141592;
void boxmuller (double mu,double sigma,int imax) {
    double num1,num2,x1,x2,z1,z2;
    for (int i = 0; i < imax; i++){
	// 0~1の一様乱数生成
	num1 = (double)rand()/(double)RAND_MAX;
	num2 = (double)rand()/(double)RAND_MAX;
	// ボックスミュラー法
	x1 = sqrt(-2.0*log(num1))*cos(2*pi*num2);
	x2 = sqrt(-2.0*log(num1))*sin(2*pi*num2);
	// 線形変換
	z1  = mu + sigma*x1;
	z2  = mu + sigma*x2;
	cout << z1 << "," << z2 << endl; 
    }
}
int main(void) {
    // ボックスミューラー法による正規分布の作成
    // 平均値、分散、個数
    boxmuller(0,1,1000);	
}

Perlのサブルーチン

Perlのリファレンス渡しについてです。

Perl入学式2013の第4回資料の一部を元に考えました。

https://github.com/perl-entrance-org/workshop-2013-04/blob/master/slide.md

例題として、ベクトルの内積を求めるサブルーチンを考えます。

ポイントは、サブルーチンに配列を渡すとき、

  • 配列をそのまま引数で渡すと、中身が展開されてしまう。
  • 展開されないようにするために、リファレンスで渡す。

のように理解しました。
まあ、合っているかどうかわかりませんが・・・。

Perl入学式の第4回では、ここら辺を確認したいと思います。

use strict;
use warnings;

my @a = qw/ 1 2 /;
my @b  = qw/ 2 3 /;
# 内積
sub dot {
    my ($v1, $v2) = @_;
    my $result = 0 ;
    # ベクトルの次元は2
    for my $i (0..1) {
	$result += @$v1[$i]*@$v2[$i];
    }
    return $result;
}
# 各配列を区別したい場合はリファレンス渡しにする
my $out = dot(\@a, \@b);
print "$out"."\n";

Perlのサブルーチン

Perlのサブルーチンのメモです。

引数は「@_」に入るので、サブルーチンの中で、引数の受け渡しが必要です。
もちろん省略は可能のようです。

サブルーチンは以下の形がひな形のようです。

sub ブルーチン名 {
    #処理の記述
}

ちょっとまだよくわかってませんが、そのうちPerl入学式第4回を受講すれば、
ある程度整理されるかと思います。(ちょっとした予習です)

#!/usr/bin/env perl
use strict;
use warnings;
#足し算のサブルーチン
sub add {
    # 引数は「@_」
    my @nums = @_;
    print "@nums\n";
    my $sum = 0.0;
    # 引数の個数分足し合わせる
    for my $num (@nums){
	$sum += $num;
    }
    print "合計=$sum\n";
}
# サブルーチン呼び出し
add(1,2,3,4,5,6,7,8,9,10);

Perlの読みたい本は、以下です。
読んで理解したら、ようやくヒヨッコかなと思います。

初めてのPerl 第6版

初めてのPerl 第6版

プログラミングPerl〈VOLUME1〉

プログラミングPerl〈VOLUME1〉

モダンPerl入門 (CodeZine BOOKS)

モダンPerl入門 (CodeZine BOOKS)

Fortran

Fortranの本です。
出てすぐ買いました。

当時、Windowsでのf90のインストールおよびその使い方がわからなくて、
むさぼるように読みました。

本書籍では、G95が使われていますが、今はバギーですね。
今は、商用をのぞけばgfortran一択です。

この本で良いところは、数値計算でのプログラムの書き方について、言及しているところです。
それまでの私が読んできた本は、「文法はこれこれです」「サンプルプログラムはこれです」が多かったので、

が書いてあって、独学できる本です。

一方、構造体については記述されていないので、別途、リファレンスを読む必要があります。

紹介でした。

数値計算のためのFortran90/95プログラミング入門

数値計算のためのFortran90/95プログラミング入門

ボックスミュラーのアルゴリズム

仕事で、乱数生成が必要になったので、アルゴリズムだけ備忘録を残します。
pythonであれば、numpyを使うのですが、業務によってはライブラリを使わない場合もあります。

参考にしたのは、下記の文献です。

自然科学の統計学 (基礎統計学)

自然科学の統計学 (基礎統計学)

このシリーズは3冊あるようで、3冊とも購入しました。
最近の書籍には珍しく?数式だったり証明が記述されており、
中身は充実してると思います。初版から20年以上たっていますが、
いまだに売れているということからも、この本の評価がわかると思います。
(もしかしたら、大学の教科書として使われていることが多いのかもしれませんが)

ちなみに、「人文・社会科学の統計学」については、経済学における統計学の使われ方も多少記述されており、
「ほうほう」と思わせてくれるものもありました。

人文・社会科学の統計学 (基礎統計学)

人文・社会科学の統計学 (基礎統計学)

以下、アルゴリズムです。
プログラムを書く場合は、一様乱数を呼び出せば、どの言語でもいけると思います。

一様乱数U1、U2から
X_1=\sqrt{-2\log{U_1}}\cos\left(2\pi U_2\right)
X_1=\sqrt{-2\log{U_1}}\sin\left(2\pi U_2\right)
により得られるX_1とX_2は互いに独立に標準正規分布N(0,1)に従う。
ちなみにlogは自然対数log_eである。

一般の正規分布N(\mu,\sigma^2)に従うZの場合は、
上記のXに対して、下記の線形変換で求める。

Z=\mu+\sigma X

Perlメモ

スクリプトを置いた箇所のファイル一覧を出すプログラムです。
下記の書籍を大いに使っております。

もっと自在にサーバを使い倒す 業務に役立つPerl (Software Design plus)

もっと自在にサーバを使い倒す 業務に役立つPerl (Software Design plus)

書籍の良い点を上げておきます。

  • 薄いこと(良くまとまっており、持ち運びに便利)
  • プログラムの一行一行に解説があること。
  • お作法についてちゃんと記述されていること。
  • リファレンスについて記述されていること。
  • オブジェクト指向について多少記述されていること。

加えて下記の本は、やさしさが溢れております。
リファレンスについては触れられておりませんが、
ちょっとしたテキスト処理なら、できるようになると思います。

新版Perl言語プログラミングレッスン入門編

新版Perl言語プログラミングレッスン入門編

use strict;
use warnings;
use Data::Dumper;
# ディレクトリの取得
my $dir = "./";
# ディレクトリハンドルをオープン
opendir my $dh, $dir
    or die qq/Can't/;
#「.」、「..」以外のファイルを取得
my $files = [];
while ( my $file = readdir $dh) {
    next if $file eq '.' || $file eq '..';
    push $files, $file;
}
print Dumper $files;
print "@$files"."\n";
closedir $dh;