シミュレーションデータ作成プログラム

by

IMG_1486.JPG
※去年の春にいった鳥取砂丘、気持ちいい

ARモデルを用いて、数値シミュレーション用のデータを自動生成するプログラムを作りました。

一応簡単なアルゴリズムとしては、ARモデルでノイズを期待値0、分散1のガウス分布に従って発生させ、1000時点毎に変化点があるようなデータとするために、ARモデルで発生させた数値に、変化量の初期値を5として、変化点毎に(変化量)-0.5していくような値を足していくようにしてます。

実行方法は、

gcc -lm ファイル名.c

でコンパイルした後、

./実行ファイル名 発生データ数 ARパラメータ1 ARパラメータ2

とすると、標準出力されます。

ファイルで用いる場合は、リダイレクトなどでテキストに落としてください。


ソースはこのようになってますので、変化点のアルゴリズムを弄ることで、発生データを自由にかえることができます。

ノイズ発生の関数の、期待値と分散の値を変えることでも、ノイズを変えることができます。

たとえばこのプログラムを、

./a.out 10000 0.6 0.5

で発生させると、
data_create_ar.gif
のようになります。

数値シミュレーションで変化点らしさを調べたいときとかに使います。
ちゃんと変化点らしさをデータの変化量を反映してるか、とかですね。

では、ソースを載せときます。

// The data of the numerical simulation is created by AR Model
// Written by Ryosuke Matsumoto
// in 2007/10/13

#include
#include
#include

double gauss_rand(void)
{
int i;
double r, ex, sigma;

ex = 0.0;
sigma = 1.0;

while(1)
{
r = (double)rand()/RAND_MAX;
if(1.0/sqrt(2*M_PI*sigma)*exp(-(r-ex)*(r-ex)/2.0/sigma/sigma) >= (double)rand()/RAND_MAX)
break;
}

return r;
}

int main(int argc,char **argv)
{

int k,t,s,max_points;
double *x,*z,*a,*e,*u;

// Check the arguments
if (argc != 4)
{
fprintf(stderr,”Usage: %s POINTS A1_PARAM A2_PARAM \n”,argv[0]);
exit(1);
}

k = 2;
max_points = atoi(argv[1]);

if ((x = (double *)malloc((max_points+1)*sizeof(double))) == NULL) {
fprintf(stderr,”Unable to malloc memory – fatal!\n”);
exit(-1);
}
if ((z = (double *)malloc((max_points+1)*sizeof(double))) == NULL) {
fprintf(stderr,”Unable to malloc memory – fatal!\n”);
exit(-1);
}
if ((a = (double *)malloc((k+1)*sizeof(double))) == NULL) {
fprintf(stderr,”Unable to malloc memory – fatal!\n”);
exit(-1);
}
if ((e = (double *)malloc((max_points+1)*sizeof(double))) == NULL) {
fprintf(stderr,”Unable to malloc memory – fatal!\n”);
exit(-1);
}
if ((u = (double *)malloc((max_points+1)*sizeof(double))) == NULL) {
fprintf(stderr,”Unable to malloc memory – fatal!\n”);
exit(-1);
}

t = 1;
s = 1;
a[0] = 0;
a[1] = atof(argv[2]);
a[2] = atof(argv[3]);
// Change Points Initial Value
u[0] = 5.0;

while (t <= max_points) { if (t > 2)
{
u[t] = u[t-1];
e[t] = gauss_rand();
z[t] = a[1] * z[t-1] – a[2] * z[t-2] + e[t];
x[t] = z[t] + u[t];

// Change Points Create Rull
if (t == 1000 * s)
{
u[t] += u[0] – 0.5;
u[0] -= 0.5;
s++;
}
}
else
{
u[t] = u[t-1];
e[t] = gauss_rand();
x[t] = e[t];
}

printf(“%lf\n”,x[t]);
t++;
}

if (x != NULL)
free(x);
if (z != NULL)
free(z);
if (a != NULL)
free(a);
if (e != NULL)
free(e);
if (u != NULL)
free(u);

return 0;
}