2014-10-31

二変量正規分布発生のプログラム

統計

昨日の記事のアルゴリズムを書いてみました。

<a href="http://nishiru3.hatenablog.com/entry/2014/10/30/123522">二変量正規乱数発生のアルゴリズム - nishiru3の日記</a>
二変量正規乱数発生のアルゴリズム - nishiru3の日記

きれいではありません。
マジックナンバーも入っております。
が、ご容赦ください。

このままだと、二変数の標準正規分布になります。
相関係数rhoを1.0とすれば、直線に並びます。

#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;
const int pi = 3.141592;
void makeGauss (double* mu,double* sigma, double rho,int imax
		,double *x, double *y) {
    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の計算
	x1 = sqrt(-2.0*log(num1))*cos(2*pi*num2);
	// 線形変換によるz1の計算
	z1  = mu[0] + sigma[0]*x1;
	// ボックスミュラー法によるx1の計算
	x2  = sqrt(-2.0*log(num1))*sin(2*pi*num2);
	// 線形変換によるz2の計算
	z2  = mu[1]+rho*sigma[1]/sigma[0]*(z1-mu[0])
	    +sqrt((1.0-rho*rho)*sigma[1]*sigma[1])*x2;
	cout << z1 << "," << z2 << endl;
	x[i] = z1; y[i] = z2;
    }
}

int main(void) {
    // ボックスミューラー法による正規分布の作成
    int imax = 1000;
    double mu[2],sigma[2],rho;
    double x[1000],y[1000];
    mu[0] = 0.; mu[1] = 0.;
    sigma[0] = 1.; sigma[1]=1.;
    rho = 0.0;
    // 平均ベクトル、分散、相関係数、発生個数
    makeGauss(mu,sigma,rho,imax,x,y);
}

2014-10-30

二変量正規乱数発生のアルゴリズム

統計

前回の記事から、二変量の正規乱数の発生方法が欲しいなと思いましたので、エッセンスだけ取りまとめました。

<a href="http://nishiru3.hatenablog.com/entry/2014/10/24/090340">ボックスミュラーのアルゴリズム - nishiru3の日記</a>
ボックスミュラーのアルゴリズム - nishiru3の日記

検索をかけたら、非常に良い資料が出てきました。
下記URLのうち、アルゴリズムに関係ある部分だけをピックアップしております。
http://www012.upp.so-net.ne.jp/doi/sas/simulation/multi_co/bivariate_normal.pdf

二変量の場合、平均値と分散だけではなく、共分散も必要になります。xとyの関係も必要ということです。

$\begin{eqnarray} \left( \begin{array}{cc} x \\ y \end{array} \right)\sim N \left( \begin{array}{cc} \left( \begin{array}{cc} \mu_x \\ \mu_y \end{array} \right) , \left( \begin{array}{cc} {\sigma_x}^2 & \sigma_{yx} \\ \sigma_{xy} & {\sigma_y}^2 \end{array} \right) \end{array} \right) \end{eqnarray}$

まず下記の条件で $x$ の正規乱数を発生させます。
$\begin{eqnarray} x \sim N\left(\mu_x,{\sigma_x}^2\right) \end{eqnarray}$
ここで、 $x$ の平均を $\mu_x$ 、分散を ${\sigma_x}^2$ とします。
これは、 $x$ に関する正規乱数になります。発生方法はボックスミューラー法を用いることにします。
さて、 $y$ についても同様に正規乱数を発生させる必要がありますが、そのままだと、yも独立に正規乱数を発生させることになります。

そこで、yの発生方法は以下のように考えます。
$\begin{eqnarray} y \sim N\left(\mu_y+\rho \frac{\sigma_y}{\sigma_x}(x-\mu_x),(1-\rho^2){\sigma_y}\right) \end{eqnarray}$
ここで、当初 ${\sigma_{xy}=\rho \sigma_x \sigma_y}$ と考えます。 $\rho$ はx、yの相関係数です。

ちなみに、x,yの平均値を0、相関係数を $\rho$ =0とすると、標準正規分布となるようです。

多変量であっても、各変量の相関係数を設定してあげれば、同様に発生させることができると思います。ただ、ちょっと煩雑になるかもしれません。

2014-10-29

二変数正規分布

統計

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

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

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

2014-10-28

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

統計

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

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

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

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

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);	
}

2014-10-27

Perlのサブルーチン

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";

2014-10-26

Perlのサブルーチン

Perl

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

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

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

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

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

#!/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);