相互相関係数

二変数の相関係数を相互相関係数と言います。
下記の式で表現されるようです。

データはxyの二つの時系列データがある場合を考えています。

\begin{eqnarray}
{r}_{xy}(h)= 
\frac{
\Sigma_{t=1}^{T-h}(x_t-\overline{x})(y_{t+h}-\overline{y})}
{\displaystyle\sqrt{\Sigma_{t=1}^{T-h}(x_t-\overline{x})^2}
\displaystyle\sqrt{\Sigma_{t=1}^{T}(y_{t+h}-\overline{y})^2}
}
\end{eqnarray}

その他の変数の意味は下記をご参照ください。


自己相関係数 - nishiru3の日記

自己相関や相互相関は、自己回帰モデルで使われるので、
重要な概念です。また、何度も言いますがスペクトル解析でも出てきます。

確率統計については、理系、文系問わずに必要な学問だなぁと改めて感じていますし、多くの分野で使われるということは、それだけ応用が利くということです。

自己相関係数のプログラム

自己相関係数を求めるためのプログラムです。

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
using namespace std;
void self_corr (double* x, int imax,double* r){
    // 平均値の算定
    double x_ave = 0.0;
    for(int i = 0; i < imax; i++){
	x_ave+=x[i];
    }
    x_ave/=(double)imax;
    // 平均値を差し引く
    double *y = new double[imax];
    for(int i = 0; i < imax; i++){
	y[i] = x[i]-x_ave;
    }
    // 自己相関係数のうち、分母denomの算定
    double denom=0.0;		// 分母denom
    for(int i = 0; i< imax; i++){
	denom=denom+y[i]*y[i];
    }
    // 自己相関係数のうち、分子numerおよび自己相関係数rの算定
    for(int h = 0; h <imax; h++){
	double numer = 0.0;	// 分子
	for(int i = 0; i < imax-h; i++){
	    numer=numer+y[i]*y[i+h];
	}
	r[h] = numer/denom;	// 自己相関係数
    }
}
//
int main(void) {
    // 入力ファイル
    fstream ifs("in.txt");
    // 時系列の個数
    int imax = 7306;
    // 時系列データxと自己相関係数rの配列確保
    double* x = new double[imax];
    double* r = new double[imax];
    // 入力ファイルからのデータ入力
    for(int i = 0; i < imax; i++) {
	ifs >> x[i];
    }
    // 自己相関係数の算定
    self_corr ( x, imax, r);
    // 自己相関係数の出力
    for(int i = 0; i < imax; i++) {
	cout << r[i] << endl;
    }
}

気象庁のホームページから、日平均気温を引っ張ってきて、
自己相関係数を求めてみました。


気象庁|過去の気象データ検索

日平均気温であれば、365日周期で、相関係数が1に近づきます。
(当たり前ですね。一年前との気温は概ね同じはず。)
ただし、ラグが大きくなればなるほど、周期はそのままに、相関係数は小さくなります。
このように当たり前のことであれば、自己相関を取るまでもなく、その特徴がわかるかと思います。

さらに踏み込むと、二変数に対して相関係数をとる方法に相互相関係数があります。
たとえば、雨と河川流量とかです。なんとなく、雨が降ると流量は増えそうだなと、直感でわかりますが、その時間差は?と聞かれると、不明です。

次回はぜひそれについて。

f:id:nishiru3:20141105184403j:plain

自己相関係数

乱流の話で出てきたのが、はじめてだったと思います。
あと、気温とか、雨量とかの分析でも使われます。

気温なんかは日変動量や月変動量等の理屈はわかっているので、あまり意味がないですが、
練習問題としてはわかりやすいです。

スペクトル解析でも当然のように出てきます。
非常に重要です。

以下、自己相関係数の定義です。

r(h) = \frac{\Sigma_{t=1}^{T-h}(x_t-\overline{x})(x_{t+h}-\overline{x})}{\Sigma_{t=1}^{T}(x_t-\overline{x})^2}

Tは時系列の個数、hはラグ、xはデータ、\overline{x}はデータの平均値です。h=0の場合、当然ですがr(0)=1です。ラグが0の場合は、その値と一致するからです。

フィルタ

カルマンフィルタの基礎

カルマンフィルタの基礎

カルマンフィルタの基礎

特に仕事で必要だった訳ではなく、気になっていたので、買ってみました。
ソースも書いてありますが、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の関係も必要ということです。

\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とすると、標準正規分布となるようです。

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