「ILE RPG で「乱数によるシミュレーション」」の続きです。
前回は乱数の発生させ方のいろいろを紹介しましたが、今回は、前回収まりきらなかったポアソン乱数と、そうした乱数を利用した「大数の法則」と「中心極限定理」のシミュレーション例です。
Program 9-4 (ポアソン乱数)
* Program 9-4 (ポアソン乱数) H DFTACTGRP(*NO) BNDDIR('QC2LE') * D pornd PR 10i 0 D a 8f * D exprnd PR 8f D a 8f * D srand PR extproc('srand') D var 10u 0 value * D rand PR 10i 0 extproc('rand') * D log PR 8f extproc('log') D var 8f value * D time PR 20i 0 extproc('time') D var 20I 0 * D a S 10i 0 D i S 10i 0 D numrnd S 10i 0 D seed1 S 8f inz(1) D seed3 S 8f inz(3) D NULL S 20i 0 inz(x'00') * /free srand(%uns(time(NULL))) ; numrnd = 200 ; for i = 1 to numrnd ; a = pornd(seed3) ; dsply %char(%dech(a:15:5)) ; endfor ; return ; /end-free * P pornd B D pornd PI 10i 0 D t 8f * D t1 S 8f D k S 10i 0 inz(0) * /free t1 = exprnd(seed1) ; dow (t1<t) ; k = k+1 ; t1 = t1 + exprnd(seed1) ; enddo ; return k ; /end-free P pornd E * P exprnd B D exprnd PI 8f D a 8f * D b S 8f D c S 8f * /free b = rand()/32768.0 ; if b = 0 ; b = 1/32768.0 ; endif ; c = -(1/a)*log(b) ; return c ; /end-free P exprnd E
Program 9-5 (大数の法則)
* Program 9-5 (大数の法則) H DFTACTGRP(*NO) BNDDIR('QC2LE') * D mean1 PR 8f D n 10i 0 * D rand PR 10i 0 extproc('rand') * D i S 10i 0 D n S 10i 0 D numrnd S 10i 0 D a S 8f * /free numrnd = 200 ; n = 5 ; for i = 1 to numrnd ; a = mean1(n) ; dsply %char(%dech(a:15:5)) ; endfor ; return ; /end-free * P mean1 B D mean1 PI 8f D n 10i 0 * D a S 8f D b S 8f D c S 8f D i S 10i 0 * /free a = 0 ; for i = 1 to n ; b = rand()/32768.0 ; a = a+b ; endfor ; c = a/n ; return c ; /end-free P mean1 E
Program 9-6 (中心極限定理)
* Program 9-6 (中心極限定理) H DFTACTGRP(*NO) BNDDIR('QC2LE') * D mean1 PR 8f D n 10i 0 * D mean2 PR 8f D n 10i 0 * D rand PR 10i 0 extproc('rand') * D sqrt PR 8f extproc('sqrt') D var 8f value * D i S 10i 0 D n S 10i 0 D numrnd S 10i 0 D a S 8f * /free numrnd = 200 ; n = 5 ; for i = 1 to numrnd ; a = mean2(n) ; dsply %char(%dech(a:15:5)) ; endfor ; return ; /end-free * P mean1 B D mean1 PI 8f D n 10i 0 * D a S 8f D b S 8f D c S 8f D i S 10i 0 * /free a = 0 ; for i = 1 to n ; b = rand()/32768.0 ; a = a+b ; endfor ; c = a/n ; return c ; /end-free P mean1 E * P mean2 B D mean2 PI 8f D n 10i 0 * D a S 8f D b S 8f * /free a = mean1(n) ; b = sqrt(12.0*n)*(a-0.5) ; return b ; /end-free P mean2 E