3-1

3-1解答

3-1

変数変換

$\begin{displaymath}\begin{cases}y_1 &= \sqrt{-2\ln x_1} \, \cos{2\pi x_2} \\ y_2 &= \sqrt{-2\ln x_1} \, \sin{2\pi x_2} \\ \end{cases}\end{displaymath}$

(11)

を用いて

$\displaystyle \int_0^{1}dx_1 \int_0^1 dx_2$

を計算して

$\displaystyle \Int dy \, \frac{1}{\sqrt{2\pi}} e^{-y^2/2} =1$

(12)

を証明せよ。この変数変換をも用いて、

の一様乱数からガウス分布に確率分布が従う確率変数を発生させる方法をBox-Muller法という。

3-1解答

積分の変数変換に於いて、変換後と変換前の微小積分要素の間には以下の関係が成り立つ。

$\displaystyle dx_1 dx_2 \dots dx_n= \vert J\vert dy_1 dy_2\dots dy_n$

(13)

ここで $\vert J\vert$ はJacobianで

$\displaystyle \vert J\vert = \left\vert \frac{\partial (x_1,x_2,\dots, x_n)}{\p... ...ots \\ \del{x_1}{y_n} & \del{x_2}{y_n} & \ldots & \del{x_n}{y_n} \end{vmatrix}$

(14)

と書ける。変数変換Eq.(11)より

$\begin{displaymath} y_1^2 + y_2^2 = -2 \ln x_1, \quad \begin{cases} 2 y_1 = -2 ... ...x_1}{y_1} = -x_1 y_1 \\ \del{x_1}{y_2} = -x_1 y_2 \end{cases}\end{displaymath}$

$\begin{displaymath} \dfrac{y_2}{y_1} = \tan{2\pi x_2}, \quad \begin{cases} -\df... ...2} &= \dfrac{\cos^2{2\pi x_2}}{2\pi} \dfrac{1}{y_1} \end{cases}\end{displaymath}$

であるから、Jacobianは

$\displaystyle \vert J\vert$	$\displaystyle = (-x_1 y_1)\left( \dfrac{\cos^2{2\pi x_2}}{2\pi} \dfrac{1}{y_1} ... ...pi x_2}}{2\pi} \left[1+\left(\frac{y_2}{y_1}\right)^2\right] =-\frac{x_1}{2\pi}$
	$\displaystyle =-\frac{e^{-(y_1^2+y_2^2)/2}}{2\pi} =-\left( \frac{1}{\sqrt{2\pi} }e^{-y_1^2/2}\right) \left( \frac{1}{\sqrt{2\pi} }e^{-y_2^2/2}\right)$

となる。よって、

$% latex2html id marker 2290 $\displaystyle 1= \int_0^1 dx_1\int_0^1 dx_2 = \left... ...2}\right)\, ; \quad \therefore\, \Int dy \, \frac{1}{\sqrt{2\pi}} e^{-y^2/2} =1$$

となり、証明された。

実際にBox-Muller法を用いて確率変数を生成し、分布のヒストグラムを書いてみると以下のようになる。

**図 4:** Box-Muller法による正規分布の生成。
$\includegraphics[width=14truecm,scale=1.1]{Box-Muller_Gauss.eps}$

例えばC言語でプログラムすると、以下のように書ける（試行数100000回）。

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>

#define PI           (M_PI)
#define MAX          (5.0)
#define MIN          (-MAX)
#define DELTA        (0.05)
#define HEIKIN(i)    ( MIN+DELTA*(i)+DELTA/2.0 )
#define NUMBER       ((int)( (MAX-MIN)/DELTA ))
#define SHIKOU       (100000)

double get_random(void) ;
void calc_y(double*,double*) ;

int main(void)
{
  FILE *fptr1, *fptr2 ;

  double y1, y2 ;
  int i, shikou ;
  int Y1[NUMBER] = {0}, Y2[NUMBER] = {0} ;

  for( shikou=0; shikou<SHIKOU; shikou++){
    calc_y( &y1, &y2) ;

    i = 0 ;
    do{
      if( y1 >= MIN+DELTA*i  && y1 < MIN+DELTA*(i+1) ) Y1[i]++ ;
      else if( y2 >= MIN+DELTA*i  && y2 < MIN+DELTA*(i+1) ) Y2[i]++ ;
      i++ ;
    }while( i<NUMBER ) ;
  }

  fptr1 = fopen("Y1_n.dat", "w") ;
  fptr2 = fopen("Y2_n.dat", "w") ;

  for( i=0; i<NUMBER; i++){
    fprintf( fptr1, "%lf   %lf\n", HEIKIN(i), ((double)Y1[i])/(SHIKOU*DELTA)) ;
    fprintf( fptr2, "%lf   %lf\n", HEIKIN(i), ((double)Y2[i])/(SHIKOU*DELTA)) ;
  }
  fclose( fptr1) ;
  fclose( fptr2) ;

  return 0 ;
}

double get_random()
{
  static int flag=0;

  if(flag==0){
    srand(time(NULL)) ;
    flag=1 ;
  }
  return 1.0*rand()/(RAND_MAX+1.) -0.0 ;
}

void calc_y(double *y1, double *y2)
{
  double x1, x2 ;

  x1 = get_random() ;
  x2 = get_random() ;

  *y1 = sqrt( -2*log(x1) ) * cos(2*PI*x2) ;
  *y2 = sqrt( -2*log(x1) ) * sin(2*PI*x2) ;
}

著者: 茅根裕司 chinone_at_astr.tohoku.ac.jp