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

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

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