RAND関数による擬似乱数の生成
○魚住 龍史 浜田 知久馬東京理科大学大学院 工学研究科 経営工学専攻 Generating pseudo-random numbers using RAND function
Ryuji Uozumi Chikuma Hamada
Department of Management Science, Graduate School of Engineering, Tokyo University of Science
要旨
ある臨床試験デザインに対して，予め見積もった症例数の妥当性を評価する 1 つの方法として，乱数 を用いたモンテカルロシミュレーションによる評価が挙げられる．SAS では，擬似乱数を生成させるた めの関数として，V9 から RAND 関数が正規版として導入された．本発表では，RAND 関数で利用でき る確率分布を整理した上で，それぞれの確率分布に従う擬似乱数を生成する方法について取りあげる．
キーワード：RAND，乱数，Mersenne Twister，CALL STREAMINIT，確率分布，正規分布，テーブル分
布，モンテカルロ法，逆関数法
1 はじめに
乱数とは，ある特定の確率分布を持つ数列である．モンテカルロ法とは，乱数を利用する技法の総称であ り，シミュレーションや確率計算においてランダムな要素を含むものはモンテカルロ法の対象である [14, 16]． 例えば，医薬品開発において，ある臨床試験デザインに対する症例数を統計学的に見積もることが求めら れる．このとき，見積もった症例数は検出力が名義水準を満たすことが求められ，検出力を評価するための 1 つの方法として，擬似乱数を生成させて実施する，モンテカルロシミュレーションによる評価が挙げられ る．本ユーザー総会においても，症例数設計法に対して，モンテカルロシミュレーションによる性能評価が 多く行われている [9]．また，より複雑な臨床試験デザイン[8] に対する症例数設計を行う場合，モンテカルロ シミュレーションは非常に有用である．SAS では，擬似乱数を生成させるための関数として，RAND 関数が提供されている．本稿では，RAND
関数で採用されている擬似乱数生成アルゴリズムについて述べる．そして，RAND 関数で利用できる確 率分布を整理した上で，それぞれの確率分布に従う擬似乱数を生成する方法について取りあげる．なお，
本稿で取りあげる確率分布は単変量の場合 [2, 3, 4]
とし，多変量の確率分布については Gentle (2003) ある いはSAS Institute Inc. (2008) を参照されたい [1, 7]．
2 RAND 関数による擬似乱数の生成
RAND 関数は V8.2 で評価版の機能として提供され，V9 から正規版として提供されている．RAND 関数が 提供される前は，RANUNI 関数や RANNOR 関数で擬似乱数生成の関数が提供されていた．RAND 関数では， Mersenne Twister と呼ばれるアルゴリズムに基づいて擬似乱数が生成される [6]．このアルゴリズムは，旧来の RANUNI 関数などで採用されている乗算合同法に基づく方法と比較して，以下の点で優れている． 周期性 旧来の RANUNI 関数などが生成する乱数の周期は2311である一方，RAND 関数による乱数は
1
2
19937
と極めて長い周期をとる． 統計的性質 生成された乱数は，623 次元超立方体の中に均等に分布する．これは，乱数を組み合わせて ベクトル形式で使用したとしても，その乱数に偏りがほとんど見られないことを意味する． 生成速度 旧来の RANUNI 関数などの関数と比べて，生成速度は同程度以上の速さである． RAND 関数による擬似乱数の生成は，データステップ内で実施する．RAND 関数の構文を図 1 に示す． data random; call streaminit(seed); /* シードの指定 */ do i=1to10000; /* 生成させる擬似乱数のオブザベーション数 */x=rand('distribution', param-1, ..., param-k); output; end; run; 図 1 : RAND 関数の構文 RAND 関数は，第 1 引数で確率分布を指定し，第 2 引数以降で確率分布に対応する数のパラメータを指定 する．図 1 のプログラムでは，データステップ内で RAND 関数を用いることによって，10000 個の擬似乱数 を生成することができる．また，擬似乱数の満たすべき条件として，再現性が挙げられる．RAND 関数に加 え，CALL ルーチン CALL STREAMINIT を記述することで，シードを指定した再現性のある擬似乱数を生成 することが可能である．なお，指定できるシードの範囲は
2
31
1
未満の整数である．2.1 RAND 関数で指定できる確率分布
RAND 関数において，第 1 引数で指定できる確率分布のうち，離散型の確率分布について確率関数及び RAND 関数の記述方法を表 1 に示す．さらに，RAND 関数で指定できる確率分布のうち，連続型の確率分布 について確率密度関数及び RAND 関数における記述方法を表 2 に示す．
表 1 : RAND 関数で利用できる離散型の確率分布 確率分布 RAND 関数における確率関数 RAND 関数における記述 二項分布 x n x
p
p
x
n
x
f
(
1
)
)
(
,x
0
,
1
,
2
,...,
n
x = RAND('BINOMIAL', p , n ) ベルヌーイ分布f
(
x
)
p
x(
1
p
)
1x,x
0
,
1
x = RAND('BERNOULLI', p ) ポアソン分布!
)
exp(
)
(
x
x
f
x
,x
0
,
1
,
2
,...
x = RAND('POISSON', ) 超幾何分布
n
N
x
n
R
N
x
R
x
f
(
)
,)
,
min(
)
,
0
max(
n
R
N
x
n
R
x = RAND('HYPER', N , R , n ) 負の二項分布 x kp
p
k
x
k
x
f
(
1
)
1
1
)
(
,x
0
,
1
,
2
,...
x = RAND('NEGBINOMIAL', p , k ) 幾何分布f
(
x
)
p
(
1
p
)
x,x
0
,
1
,
2
,...
x = RAND('GEOMETRIC', p ) テーブル分布† ip
i
f
(
)
, x1,2,...,n
n i i p n f 1 1 ) 1 ( x = RAND('TABLE', p , 1 p , ...) 2 †順序カテゴリカルデータを生成 表 2 : RAND 関数で利用できる連続型の確率分布 確率分布 RAND 関数における確率密度関数 RAND 関数における記述 正規分布 2 2 2 ) ( exp 2 1 ) ( x x f ,
x
x = RAND('NORMAL', , ) 対数正規分布 2 ) (log exp 2 1 ) ( 2 x x x f , x0 x = RAND('LOGNORMAL') t 分布 2 1 2 1 2 2 1 ) ( x x f ,
x
x = RAND('T', ) コーシー分布 2 1 1 1 ) ( x x f ,
x
x = RAND('CAUCHY') F 分布 2 2 1 1 2 2 2 2 1 2 1 2 1 2 1 1 2 1 ) ( 2 2 2 ) ( x x f , x x = RAND('F', 1, 2) ガンマ分布 ) exp( ) ( 1 ) ( x 1 x a x f a , x0 x = RAND('GAMMA', a )アーラン分布 ) exp( ) ( 1 ) ( x 1 x a x f a , x0 x = RAND('ERLANG', a‡) ‡ a : 整数 カイ二乗分布 2 exp 2 2 ) ( 2 1 2 x x x f , x0 x = RAND('CHISQUARE', ) 指数分布 f(x)exp(x), x0 x = RAND('EXPONENTIAL') ワイブル分布 x x x f( ) 1exp , x0 x = RAND('WEIBULL', , ) ベータ分布 1 1 ) 1 ( ) ( ) ( ) ( ) ( xa xb b a b a x f , 0x1 x = RAND('BETA', a , b ) 一様分布 f(x)1, 0x1 x = RAND('UNIFORM') 三角分布 ) 1 ( 1 ) 1 ( 2 ) 0 ( 2 ) ( x h h x h x h x x f x = RAND('TRIANGLE', h ) SAS プログラム 理論式 data random; call streaminit(20130718); do i=1to10000;
do case=1; x=rand('tabled',0.4,0.3,0.2);output;end; do case=2; x=rand('tabled',0.1,0.1,0.2,0.2,0.3);output;end; do case=3; x=rand('tabled',0.7,0.4,0.3);output; end; end; run; 1 1
n i i p ip
i
f
(
)
, i1,2,...,n,
n i i p n f 1 1 ) 1 ( 1 1
n i i p ip
i
f
(
)
, i1,2,...,j1,
1 1 1 ) ( j i i p j f 乱数のヒストグラム 図 2 : テーブル分布に従う擬似乱数の生成SAS プログラム 理論分布 data random;
call streaminit(20130718); do i=1to10000;
do case=1; x=rand('normal',0,1);output;end; do case=2; x=rand('normal',2,1);output;end; do case=3; x=rand('normal',2,0.5);output;end;
end; run; 乱数のヒストグラム 図 3 : 正規分布に従う擬似乱数の生成 ここで，離散型の確率分布に従う擬似乱数の例として，テーブル分布に従う擬似乱数 ~ ( , ,..., ) 2 1 p pn p Table X ， 連続型の確率分布に従う擬似乱数の例として，正規分布に従う擬似乱数X ~N(
,
2)を考える．それぞれの 擬似乱数を生成させる SAS プログラム，理論分布・理論式，及び生成させた擬似乱数の分布を図 2，図 3 に 示す．なお，ヒストグラムの作成には SGPLOT プロシジャを用いた [10, 11, 12, 13]．理論分布・理論式を対照に すると，求める擬似乱数を生成できたことがわかる．また，正規分布のパラメータを指定せずに実行すると， 標準正規分布X
~ N
(
0
,
1
)
に従う擬似乱数が生成される．さらに，テーブル分布のパラメータの指定を1
1
n i ip
の状態の下，第 (n + 1) パラメータを指定せずに実行すると，
n i ip
n
f
11
)
1
(
を求めた上で，
n i i n p p p p Table X 1 2 1, ,..., ,1 ~ に従う擬似乱数が生成される． なお，RAND 関数で指定できる確率分布のうち，確率 (密度) 関数を求めるための PDF 関数を簡略化した パラメータで定義されているものがある．例えば，各関数で定義されている指数分布及びワイブル分布の確 率密度関数を表 3 に示す．表 3 の場合，RAND 関数でパラメータ の指数分布に従う擬似乱数を生成させる 1 つの方法として，それぞれの確率分布間の関係を活用することが挙げられる．一部の確率分布の間には， パラメータに特定の値を代入した特別な場合，変数変換，漸近近似などによって，関連付けられるものがあ る [5]．表 3 : 指数分布及びワイブル分布の確率密度関数 確率分布 確率密度関数 RAND 関数 指数分布 PDF('EXPONENTIAL', ) x x f( ) 1exp , x0 x = RAND('EXPONENTIAL') ) exp( ) (x x f , x0 ワイブル分布 PDF('WEIBULL', , ) x x x f( ) 1exp , x0 x = RAND('WEIBULL', , ) x x x f( ) 1exp , x0 1 f(x) 1exp x 図 4 : RAND 関数で指定できる確率分布の関係 図 4 に，RAND 関数で指定できる確率分布の関係を示す．表 3 の場合で考えると，指数分布はワイブル分 布の特別な場合であるといえる．よって，RAND 関数の確率分布として，表 3 に示したワイブル分布のパラ メータを指定することにより，パラメータ の指数分布に従う擬似乱数を生成することが可能である．その 他の確率分布についても，確率分布間の関係を活用して，関連する異なった分布から擬似乱数を生成するこ
とができる．なお，図 4 における一部の確率分布は，表 1 及び表 2 で示した確率 (密度) 関数とは異なった定 義の確率分布で考えている [5]．
2.2 逆関数法による擬似乱数の生成
RAND 関数では指定できない確率分布に従う擬似乱数の生成が必要な場合がある．このとき，解決策の 1 つとして，逆関数法を用いることが挙げられる．いくつかある手法のうち，考え方がシンプルな方法である といえる．例として，対数ロジスティック分布を取りあげると，RAND 関数において，対数ロジスティック 分布は第 1 引数に指定できる確率分布として含まれていない．そこで，逆関数法を用いて，ワイブル分布及 び対数ロジスティック分布に従う擬似乱数を生成する方法を示す． 逆関数法 [15] は，指定する確率分布の累積分布関数F
(x
)
の逆関数}
)
(
:
inf{
)
(
1y
x
F
x
y
F
,0
y
1
が存在する場合，図 5 より，標準一様分布U
~ U
(
0
,
1
)
からの変数変換)
(
1U
F
X
として擬似乱数を生成することができる．なお，XがF(x)に従うことは，)
(
)
(
1x
F
y
x
y
F
より，)
(
)}
(
{
}
)
(
{
)
(
X
x
P
F
1U
x
P
U
F
x
F
x
P
でわかる． 図 5 : 逆関数法の原理 ここで，RAND 関数によって生成させた一様分布に従う擬似乱数を用いて，逆関数法によりワイブル分布 に従う擬似乱数を生成させる．さらに，RAND 関数における確率分布としてワイブル分布を指定させて生成 させた擬似乱数と比較する．それぞれの方法で生成させた擬似乱数のヒストグラムを図 6 に示す．これを見 ると，逆関数法によって生成させた擬似乱数は，RAND 関数で生成させた擬似乱数と同様の分布であること が確認できる． また，RAND 関数においてサポートされていない確率分布に従う擬似乱数の生成が必要な場合がある．例 として，RAND 関数の確率分布において，対数ロジスティック分布 [5, 17] は第 1 引数に指定できる確率分布と して含まれていない．RAND 関数で指定できる確率分布を用いて，対数ロジスティック分布のような確率分布に従う擬似乱数を生成させる方法の 1 つとして，逆関数法を用いることが挙げられる．そこで，逆関数法 を用いて，対数ロジスティック分布に従う擬似乱数を生成する方法を，ワイブル分布に従う擬似乱数を生成 する方法とともに表 4 に示す． 表 4 : 各関数で定義されている指数分布及びワイブル分布の確率密度関数 ワイブル分布 対数ロジスティック分布 確率密度関数 x x x f( ) 1exp , x0
x x x f 1 ) ( 1 , x0 累積分布関数 x x F( ) 1 exp , x0
x x F 1 1 1 ) ( , x0 累積分布関数の 逆関数
1/ 1 ) 1 log( ) (x U F
/ 1 1 1 1 ) ( U U x F 確率変数
1/ logU X
/ 1 1 1 U U X 図 6 : それぞれの方法で生成させたワイブル分布に従う擬似乱数の分布3 まとめ
本稿では，SAS で提供されている，RAND 関数で採用されている擬似乱数生成アルゴリズムについて， 旧来の RANUNI 関数などで採用されている方法と比べ，周期の長さ，生成速度，そして統計的性質の点から 優れていることを述べた．次に，RAND 関数による擬似乱数生成方法について，RAND 関数で指定できる 確率分布及びそのパラメータを詳述した．また，RAND 関数において，PDF 関数とは異なったパラメー タで定義されている確率分布があることを取りあげた．さらに，対数ロジスティック分布のように， RAND 関数においてサポートされていない確率分布に従う擬似乱数を生成させる方法として，逆関数法 について取りあげた．以上より，医薬品開発で見積もった症例数に対する検出力の評価のように，モンテカルロ法による擬似乱数を用いたシミュレーションが求められる場合，RAND 関数によって生成させ た擬似乱数を用いることは有用であるといえる．
参考文献
[1] Gentle JE. Random Number Generation and Monte Carlo Methods, Second Edition. Springer-Verlag; 2003. [2] Johnson NL, Kotz S, Balakrishnan N. Continuous Univariate Distributions, Vol.1, Second Edition. New York: John
Wiley and Sons; 1994.
[3] Johnson NL, Kotz S, Balakrishnan N. Continuous Univariate Distributions, Vol.2, Second Edition. New York: John Wiley and Sons; 1995.
[4] Johnson NL, Kemp AW, Kotz S. Univariate Discrete Distributions, Third Edition. New York: John Wiley and Sons; 2005.
[5] Leemis LM, McQueston JT. Univariate Distribution Relationships. The American Statistician 2008; 62: 45–53. [6] Matsumoto M, Nishimura T. Mersenne Twister: A 623–Dimensionally Equidistributed Uniform Pseudo-Random
Number Generator. ACM Transactions on Modeling and Computer Simulation 1998; 8: 3–30. [7] SAS Institute Inc. SAS/IML(R) 9.2 User’s Guide, Cary, NC, USA: SAS Institute Inc; 2008.
[8] Uozumi R, Hamada C. An Adaptive Subpopulation Selection Design Using Time-to-Event Endpoints. Abstracts for the XXVIth International Biometric Conference, 26-31 August, 2012, Kobe, Japan. International Biometric Society (http://secretariat.ne.jp/ibc2012/programme.html). ISBN 978-0-9821919-2-7.
[9] 魚住龍史, 水澤純基, 浜田知久馬. 生存時間解析における Lakatos の症例数設計法の有用性の評価. SAS ユ ーザー総会 論文集 2009; 19–28.
[10] 魚住龍史, 浜田知久馬. SG (Statistical Graphics) Procedures による Kaplan-Meier プロットの作成. SAS ユ ーザー総会 論文集 2011, 185–199.
[11] 魚 住 龍 史 , 浜 田 知 久 馬 . が ん 臨 床 試 験 に お け る 腫 瘍 縮 小 効 果 の 検 討 に 有 用 な グ ラ フ の 作 成 –SGPLOT プロシジャの最新機能を活用–. SAS ユーザー総会 論文集 2012, 151–165.
[12] 高浪洋平. SG プロシジャと GTL によるグラフの作成と ODS PDF による統合解析帳票の作成 ～TQT 試 験における活用事例～. SAS ユーザー総会 論文集 2011, 201–219.
[13] 高浪洋平. SAS と HTML アプリケーションによる CDISC ADaM 形式の解析用データセットを用いた有害 事象の解析帳票・グラフ簡易作成ツールの開発事例. SAS ユーザー総会 論文集 2012, 185–205. [14] 津田孝夫. モンテカルロ法とシミュレーション. 培風館, 1995. [15] 東京大学教養学部統計学教室. 自然科学の統計学. 東京大学出版会, 1992. [16] 伏見正則. 乱数. 東京大学出版会, 1989. [17] 蓑谷千凰彦. 確率統計ハンドブック. 朝倉書店, 2003.
