フィルタ

カルマンフィルタの基礎

カルマンフィルタの基礎

特に仕事で必要だった訳ではなく、気になっていたので、買ってみました。
ソースも書いてありますが、Matlabです。
「Octave」「Scilab」であればあまり変更せずに使えると思います。

時系列解析入門

時系列解析入門

カルマンフィルタもそうですが、パーティクルフィルタ（or モンテカルロフィルタ）について記述されています。
ソースコードはWebに公開されています。

Perlの勉強

Perl入学式の資料を引用しています。

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

文字列の検索です。
マッチするかどうかは「=~//」でいけるようです。

my $str = 'python,perl,ruby,php,awk,sed';
if ($str =~ /perl/) {
    print "'$str'は'perl'を含みます.";
}

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

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

二変量正規乱数発生のアルゴリズム - 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);
}

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

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

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

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

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

まず下記の条件での正規乱数を発生させます。

ここで、の平均を、分散をとします。
これは、に関する正規乱数になります。発生方法はボックスミューラー法を用いることにします。
さて、についても同様に正規乱数を発生させる必要がありますが、そのままだと、yも独立に正規乱数を発生させることになります。

そこで、yの発生方法は以下のように考えます。

ここで、当初と考えます。はx、yの相関係数です。

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

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

二変数正規分布

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

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

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

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

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

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

ボックスミュラーのアルゴリズム - 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";