RAND 関数による擬似乱数の生成 魚住龍史浜田知久馬 東京理科大学大学院工学研究科経営工学専攻 Geeratig pseudo-radom umbers usig RAND fuctio Ryuji Uozumi Chikuma Hamada Departmet of Maagemet Sciece Graduate School of Egieerig Tokyo Uiversity of Sciece 要旨 ある臨床試験デザインに対して 予め見積もった症例数の妥当性を評価する つの方法として 乱数を用いたモンテカルロシミュレーションによる評価が挙げられる.SAS では 擬似乱数を生成させるための関数として V9 から RAND 関数が正規版として導入された. 本発表では RAND 関数で利用できる確率分布を整理した上で それぞれの確率分布に従う擬似乱数を生成する方法について取りあげる. キーワード :RAND 乱数 Mersee TwisterCALL STREAMINIT 確率分布 正規分布 テーブル分 布 モンテカルロ法 逆関数法 はじめに 乱数とは ある特定の確率分布を持つ数列である. モンテカルロ法とは 乱数を利用する技法の総称であり シミュレーションや確率計算においてランダムな要素を含むものはモンテカルロ法の対象である [4 6]. 例えば 医薬品開発において ある臨床試験デザインに対する症例数を統計学的に見積もることが求められる. このとき 見積もった症例数は検出力が名義水準を満たすことが求められ 検出力を評価するための つの方法として 擬似乱数を生成させて実施する モンテカルロシミュレーションによる評価が挙げられる. 本ユーザー総会においても 症例数設計法に対して モンテカルロシミュレーションによる性能評価が多く行われている [9]. また より複雑な臨床試験デザイン [8] に対する症例数設計を行う場合 モンテカルロシミュレーションは非常に有用である. SAS では 擬似乱数を生成させるための関数として RAND 関数が提供されている. 本稿では RAND 関数で採用されている擬似乱数生成アルゴリズムについて述べる. そして RAND 関数で利用できる確率分布を整理した上で それぞれの確率分布に従う擬似乱数を生成する方法について取りあげる. なお 本稿で取りあげる確率分布は単変量の場合 [ 3 4] とし いは SAS Istitute Ic. (008) を参照されたい [ 7]. 多変量の確率分布については Getle (003) ある
RAND 関数による擬似乱数の生成 RAND 関数は V8. で評価版の機能として提供され V9 から正規版として提供されている.RAND 関数が提供される前は RANUNI 関数や RANNOR 関数で擬似乱数生成の関数が提供されていた.RAND 関数では Mersee Twister と呼ばれるアルゴリズムに基づいて擬似乱数が生成される [6]. このアルゴリズムは 旧来の RANUNI 関数などで採用されている乗算合同法に基づく方法と比較して 以下の点で優れている. 周期性旧来の RANUNI 関数などが生成する乱数の周期は 3 である一方 RAND 関数による乱数は 9937 と極めて長い周期をとる. 統計的性質 生成された乱数は 63 次元超立方体の中に均等に分布する. これは 乱数を組み合わせて ベクトル形式で使用したとしても その乱数に偏りがほとんど見られないことを意味する. 生成速度旧来の RANUNI 関数などの関数と比べて 生成速度は同程度以上の速さである. RAND 関数による擬似乱数の生成は データステップ内で実施する.RAND 関数の構文を図 に示す. data radom; call streamiit(seed); /* シードの指定 */ do i= to 0000; /* 生成させる擬似乱数のオブザベーション数 */ =rad('distributio' param-... param-k); output; ed; ru; 図 : RAND 関数の構文 RAND 関数は 第 引数で確率分布を指定し 第 引数以降で確率分布に対応する数のパラメータを指定 する. 図 のプログラムでは データステップ内で RAND 関数を用いることによって 0000 個の擬似乱数 を生成することができる. また 擬似乱数の満たすべき条件として 再現性が挙げられる.RAND 関数に加 え CALL ルーチン CALL STREAMINIT を記述することで シードを指定した再現性のある擬似乱数を生成 することが可能である. なお 指定できるシードの範囲は 3 未満の整数である.. RAND 関数で指定できる確率分布 RAND 関数において 第 引数で指定できる確率分布のうち 離散型の確率分布について確率関数及び RAND 関数の記述方法を表 に示す. さらに RAND 関数で指定できる確率分布のうち 連続型の確率分布について確率密度関数及び RAND 関数における記述方法を表 に示す.
表 : RAND 関数で利用できる離散型の確率分布確率分布 RAND 関数における確率関数 RAND 関数における記述二項分布 = RAND('BINOMIAL' p ) p ( p) 0... ベルヌーイ分布 f ( ) p ( p) 0 = RAND('BERNOULLI' p ) ポアソン分布超幾何分布負の二項分布幾何分布 ep( ) 0...! R N R N ma( 0 R N) mi( R) k ) k k ( p p 0... p( p 0... ) = RAND('POISSON' ) = RAND('HYPER' N R ) = RAND('NEGBINOMIAL' p k ) = RAND('GEOMETRIC' p ) テーブル分布 f ( i) p... f ( ) i p i i = RAND('TABLE' p p...) 順序カテゴリカルデータを生成 表 : RAND 関数で利用できる連続型の確率分布 確率分布 RAND 関数における確率密度関数 RAND 関数における記述 正規分布 対数正規分布 t 分布 ( ) ep (log ) ep 0 = RAND('NORMAL' ) = RAND('LOGNORMAL') = RAND('T' ) コーシー分布 分布 ガンマ分布 f ( ) ( ) a ep( ) 0 ( a) = RAND('CAUCHY') = RAND('' ) = RAND('GAMMA' a )
アーラン分布 a ep( ) 0 ( a) = RAND('ERLANG' a ) a : 整数 カイ二乗分布 ep 0 = RAND('CHISQUARE' ) 指数分布 ep( ) 0 = RAND('EXPONENTIAL') ワイブル分布 ep 0 = RAND('WEIBULL' ) ベータ分布 ( a b) a b ( ) 0 ( a) ( b) = RAND('BETA' a b ) 一様分布 0 = RAND('UNIORM') 三角分布 h ( ) h (0 h) ( h ) = RAND('TRIANGLE' h ) SAS プログラム data radom; call streamiit(03078); do i= to 0000; do case=; =rad('tabled'0.40.30.);output;ed; do case=; =rad('tabled'0.0.0.0.0.3);output;ed; do case=3; =rad('tabled'0.70.40.3);output; ed; ed; ru; 乱数のヒストグラム 理論式 i i p i p i i f ( i) p i i... f ( ) p i j f ( i) p i i... j f ( j) p i i 図 : テーブル分布に従う擬似乱数の生成
SAS プログラム data radom; call streamiit(03078); do i= to 0000; do case=; =rad('ormal'0);output;ed; do case=; =rad('ormal');output;ed; do case=3; =rad('ormal'0.5);output;ed; ed; ru; 乱数のヒストグラム 理論分布 図 3 : 正規分布に従う擬似乱数の生成 ここで 離散型の確率分布に従う擬似乱数の例として テーブル分布に従う擬似乱数 X Table( p p... p ) ~ 連続型の確率分布に従う擬似乱数の例として 正規分布に従う擬似乱数 X ~ N( ) を考える. それぞれの 擬似乱数を生成させる SAS プログラム 理論分布 理論式 及び生成させた擬似乱数の分布を図 図 3 に 示す. なお ヒストグラムの作成には SGPLOT プロシジャを用いた [0 3]. 理論分布 理論式を対照に すると 求める擬似乱数を生成できたことがわかる. また 正規分布のパラメータを指定せずに実行すると 標準正規分布 X ~ N(0 ) に従う擬似乱数が生成される. さらに テーブル分布のパラメータの指定を i pi の状態の下 第 ( + ) パラメータを指定せずに実行すると f ( ) p i を求めた上で X ~ Table p p... p p i に従う擬似乱数が生成される. i なお RAND 関数で指定できる確率分布のうち 確率 ( 密度 ) 関数を求めるための PD 関数を簡略化した パラメータで定義されているものがある. 例えば 各関数で定義されている指数分布及びワイブル分布の確 率密度関数を表 3 に示す. 表 3 の場合 RAND 関数でパラメータ の指数分布に従う擬似乱数を生成させる つの方法として それぞれの確率分布間の関係を活用することが挙げられる. 一部の確率分布の間には パラメータに特定の値を代入した特別な場合 変数変換 漸近近似などによって 関連付けられるものがあ る [5]. i
表 3 : 指数分布及びワイブル分布の確率密度関数 確率分布 確率密度関数 RAND 関数 指数分布 PD('EXPONENTIAL' ) = RAND('EXPONENTIAL') ep( ) 0 ep 0 ワイブル分布 PD('WEIBULL' ) = RAND('WEIBULL' ) ep 0 ep 0 ep 図 4 : RAND 関数で指定できる確率分布の関係 図 4 に RAND 関数で指定できる確率分布の関係を示す. 表 3 の場合で考えると 指数分布はワイブル分布の特別な場合であるといえる. よって RAND 関数の確率分布として 表 3 に示したワイブル分布のパラメータを指定することにより パラメータ の指数分布に従う擬似乱数を生成することが可能である. その他の確率分布についても 確率分布間の関係を活用して 関連する異なった分布から擬似乱数を生成するこ
とができる. なお 図 4 における一部の確率分布は 表 及び表 で示した確率 ( 密度 ) 関数とは異なった定 義の確率分布で考えている [5].. 逆関数法による擬似乱数の生成 RAND 関数では指定できない確率分布に従う擬似乱数の生成が必要な場合がある. このとき 解決策の つとして 逆関数法を用いることが挙げられる. いくつかある手法のうち 考え方がシンプルな方法である といえる. 例として 対数ロジスティック分布を取りあげると RAND 関数において 対数ロジスティック 分布は第 引数に指定できる確率分布として含まれていない. そこで 逆関数法を用いて ワイブル分布及 び対数ロジスティック分布に従う擬似乱数を生成する方法を示す. 逆関数法 [5] は 指定する確率分布の累積分布関数 () の逆関数 ( y) if{ : ( ) y} 0y が存在する場合 図 5 より 標準一様分布 U ~U(0 ) からの変数変換 X ( U) として擬似乱数を生成することができる. なお X が () に従うことは より でわかる. ( y) P( X ) P{ y ( ) ( U) } P{ U ( )} ( ) 図 5 : 逆関数法の原理 ここで RAND 関数によって生成させた一様分布に従う擬似乱数を用いて 逆関数法によりワイブル分布に従う擬似乱数を生成させる. さらに RAND 関数における確率分布としてワイブル分布を指定させて生成させた擬似乱数と比較する. それぞれの方法で生成させた擬似乱数のヒストグラムを図 6 に示す. これを見ると 逆関数法によって生成させた擬似乱数は RAND 関数で生成させた擬似乱数と同様の分布であることが確認できる. また RAND 関数においてサポートされていない確率分布に従う擬似乱数の生成が必要な場合がある. 例 として RAND 関数の確率分布において 対数ロジスティック分布 [5 7] は第 引数に指定できる確率分布と して含まれていない.RAND 関数で指定できる確率分布を用いて 対数ロジスティック分布のような確率分
布に従う擬似乱数を生成させる方法の つとして 逆関数法を用いることが挙げられる. そこで 逆関数法 を用いて 対数ロジスティック分布に従う擬似乱数を生成する方法を ワイブル分布に従う擬似乱数を生成 する方法とともに表 4 に示す. 表 4 : 各関数で定義されている指数分布及びワイブル分布の確率密度関数 確率密度関数累積分布関数累積分布関数の逆関数確率変数 ワイブル分布 ep 0 ( ) ep 0 / log( ) ( ) U / X logu 対数ロジスティック分布 0 ( ) 0 U ( ) U U X U / / 図 6 : それぞれの方法で生成させたワイブル分布に従う擬似乱数の分布 3 まとめ 本稿では SAS で提供されている RAND 関数で採用されている擬似乱数生成アルゴリズムについて 旧来の RANUNI 関数などで採用されている方法と比べ 周期の長さ 生成速度 そして統計的性質の点から優れていることを述べた. 次に RAND 関数による擬似乱数生成方法について RAND 関数で指定できる確率分布及びそのパラメータを詳述した. また RAND 関数において PD 関数とは異なったパラメータで定義されている確率分布があることを取りあげた. さらに 対数ロジスティック分布のように RAND 関数においてサポートされていない確率分布に従う擬似乱数を生成させる方法として 逆関数法について取りあげた. 以上より 医薬品開発で見積もった症例数に対する検出力の評価のように モン
テカルロ法による擬似乱数を用いたシミュレーションが求められる場合 RAND 関数によって生成させ た擬似乱数を用いることは有用であるといえる. 参考文献 [] Getle JE. Radom Number Geeratio ad Mote Carlo Methods Secod Editio. Spriger-Verlag; 003. [] Johso NL Kotz S Balakrisha N. Cotiuous Uivariate Distributios Vol. Secod Editio. New York: Joh Wiley ad Sos; 994. [3] Johso NL Kotz S Balakrisha N. Cotiuous Uivariate Distributios Vol. Secod Editio. New York: Joh Wiley ad Sos; 995. [4] Johso NL Kemp AW Kotz S. Uivariate Discrete Distributios Third Editio. New York: Joh Wiley ad Sos; 005. [5] Leemis LM McQuesto JT. Uivariate Distributio Relatioships. The America Statisticia 008; 6: 45 53. [6] Matsumoto M Nishimura T. Mersee Twister: A 63 Dimesioally Equidistributed Uiform Pseudo-Radom Number Geerator. ACM Trasactios o Modelig ad Computer Simulatio 998; 8: 3 30. [7] SAS Istitute Ic. SAS/IML(R) 9. User s Guide Cary NC USA: SAS Istitute Ic; 008. [8] Uozumi R Hamada C. A Adaptive Subpopulatio Selectio Desig Usig Time-to-Evet Edpoits. Abstracts for the XXVIth Iteratioal Biometric Coferece 6-3 August 0 Kobe Japa. Iteratioal Biometric Society (http://secretariat.e.jp/ibc0/programme.html). ISBN 978-0-9899--7. [9] 魚住龍史 水澤純基 浜田知久馬. 生存時間解析における Lakatos の症例数設計法の有用性の評価. SAS ユーザー総会論文集 009; 9 8. [0] 魚住龍史 浜田知久馬. SG (Statistical Graphics) Procedures による Kapla-Meier プロットの作成. SAS ユーザー総会論文集 0 85 99. [] 魚住龍史 浜田知久馬. がん臨床試験における腫瘍縮小効果の検討に有用なグラフの作成 SGPLOT プロシジャの最新機能を活用. SAS ユーザー総会論文集 0 5 65. [] 高浪洋平. SG プロシジャと GTL によるグラフの作成と ODS PD による統合解析帳票の作成 ~TQT 試験における活用事例 ~. SAS ユーザー総会論文集 0 0 9. [3] 高浪洋平. SAS と HTML アプリケーションによる CDISC ADaM 形式の解析用データセットを用いた有害事象の解析帳票 グラフ簡易作成ツールの開発事例. SAS ユーザー総会論文集 0 85 05. [4] 津田孝夫. モンテカルロ法とシミュレーション. 培風館 995. [5] 東京大学教養学部統計学教室. 自然科学の統計学. 東京大学出版会 99. [6] 伏見正則. 乱数. 東京大学出版会 989. [7] 蓑谷千凰彦. 確率統計ハンドブック. 朝倉書店 003. 連絡先 E-mail : uozumi@ms.kagu.tus.ac.jp