TOP > 国内特許検索 > 多変量非正規分布に従う乱数発生方法及びそのパラメータの推定方法

多変量非正規分布に従う乱数発生方法及びそのパラメータの推定方法

国内特許コード P140011047
整理番号 2002-P29
掲載日 2014年10月30日
出願番号 特願2003-082816
公開番号 特開2004-005510
登録番号 特許第3708928号
出願日 平成15年3月25日(2003.3.25)
公開日 平成16年1月8日(2004.1.8)
登録日 平成17年8月12日(2005.8.12)
優先権データ
  • 特願2002-093355 (2002.3.28) JP
  • 特願2002-093356 (2002.3.28) JP
発明者
  • 永原 裕一
出願人
  • 学校法人明治大学
発明の名称 多変量非正規分布に従う乱数発生方法及びそのパラメータの推定方法
発明の概要 【課題】ピアソン分布系及び一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築する。
【解決手段】コンピュータ上でユアン及びベントラーの方法Iに基づいて多変量非正規分布に従う乱数を発生する乱数発生方法であって、コンピュータを用いてn次元経験分布にn次元多変量非正規分布をあてはめるステップと、コンピュータを用いて乱数を発生するステップとを有し、前記あてはめるステップは、3次及び4次モーメントについての所定の関係式を用いて、前記経験分布の3次及び4次モーメントに関するあてはめを行う。また、この乱数発生法を用いて最尤法によってパラメータを推定する。さらに、前記乱数発生方法及びパラメータ推定方法を金融分野、半導体へのイオン注入等へのシミュレーションに応用する。
【選択図】図2
従来技術、競合技術の概要


従来、多変量非正規分布を生成する方法は、数えるしか存在していない。このような分布は、多変量解析と統計モデリングの分野において、切望されている。



フライシュマン(Fleishmann)(1978)、タジカマラ(Tadikamalla)(1980)、ベイルとマブレリ(Vale and Mavrelli)(1983)は、多変量非正規分布を生成する方法を提案した。これらの方法は、シミュレーションした母集団の確率ベクトルの4次のモーメント行列を計算することが難しいという欠点がある。ユアン及びベントラー(1997)は、4次のモーメント行列を考慮し、特定の周辺の歪度と尖度を有する多変量非正規分布を生成する手法を提案した。



本明細書では、ユアン等が方法Iと称したものは、次の通りである。ξ1
・・・,ξmは、独立な確率変数であり、E(・)を期待値として、
E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj及びξ=(ξ1,・・・,ξ1)´である。νは、ξに独立な確率変数であり、
E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。また、T=(tij)は、階数pのp行m列の非確率行列であり、TT´=Σとm≧pを満たす。すると、次の式(35)で与えられる確率ベクトルは、一般に、非楕円分布に従い、Cov(X)=Σを満たす。ここで、Cov(・)は、確率ベクトルX=(x1,・・・,xn)´の分散共分散行列である。



【数40】


iの周辺の歪度(skew)と尖度(kurt)は、それぞれ次の式(36)及び
(37)によって与えられる。



【数41】


【数42】


Xは、ξ~Nm(0,I)及びPr(ν≧0)=1(ただし、Nmはm次元正規分布、Prは確率を表す。)が満たされるとき、楕円分布に従う(ファン、コッツ及びング(Fang, Kotz and Ng)(1990))。ν=1とし、
1,・・・,zpを独立な標準正規変数とし、
ξ=|z1|-E(|z1|),ξj=zj,j=2,・・・,pとすると、式(35)におけるXは、アザリニ及びバレ(Azzalini and Valle)(1996)によって定義された歪んだ正規分布に従う。



この方法には、母集団の4次モーメント行列を容易に計算できるという利点がある。Sを式(35)のXからの大きさnの標本に基づいた標本共分散とする。
vech(・)を対称行列から重複しない要素を取り出してベクトルに変換する演算子とし、s=vech(S)とする。すると、sの漸近的分散共分散行列はΓ/nによって与えられる。ここで、
Γ=var(vech(XX´))であり、次の式(38)によって与えられる。varは、分散である。



【数43】


ただし、tjは、T=(tij)の第j列ベクトルであり、
【数44】


は、クロネッカー積である。Dpは、p次のデュプリケーション行列であり、Dp+は、Dpのムーア・ペンローズの一般化逆行列である。



この後、ユアン等は、検定統計の理論を研究した(ユアン及びベントラー(1999a,1999b,2000))。ユアン等は、確率ベクトルを生成するにはフライシュマン(1978)の方法を推奨している。しかし、この分布は、あまり広い歪度及び尖度を有していない。これに対して、ピアソン分布系は、次のような特徴を有する。



(1)ピアソン分布系は、様々な歪度と尖度を有する広いクラスの分布を表現することができる。ピアソン分布系は、ガンマ分布、ベータ分布、t分布等のいくつかの良く知られた分布を含む。



(2)ピアソン分布系は、4つのモーメントによって特徴付けられる。この特徴は、ユアン等の方法を適用するためには十分である。



(3)IV型を除いたピアソン分布系に従う乱数は、正規分布とガンマ分布に従う乱数を用いて生成される。



(4)最近、実現することが難しかったピアソンIV型分布を含む殆ど全ての種類を用いる実用的な手法が開発された(永原(2002))。



ここで、ピアソン分布系を導入する。ケイ・ピアソン(K. Pearson)(1895)は、確率密度関数pに関する次の式(39)の微分方程式に従い、ピアソン分布系を定義した。



【数45】


ピアソン第IV型分布は、正規化定数を計算することが難しいので、実現することが難しかったが、論文(永原(1999))において、計算に用いる再帰的公式を導出し、様々な問題を解決した。この後、非ガウスフィルタに用いるピアソンIV型分布を含む分布系の殆ど全ての種類を用いる実用的な手法を開発した(永原(2002))。これらの研究により、ピアソン分布系が実用的に利用できるようになった。



ここで、ピアソン分布系に関して参照される研究者と文献を掲げる(ピアソン(1914)、ピアソン及びハートリー(1954)、ピアソン(1963)、ホードリ-(1968)、エルダートン及びジョンソン(1969)、オード(1972)、パリシュ(1983)、ジョンソン、コッツ及びバラカリシュナン(1994)、スチュアート及びオード(1994))。



参考文献としては、次のものが挙げられる。



[1]アザリニ,エイ.及びバレ,エイ.ディー.(Azzalini, A. and Valle, A. D.)(1996)「多変量の歪んだ正規分布(The multivariate skew-normal distribution)」、バイオメトリカ(Biometrika)、83巻、715~726頁
[2]ドブロイ,エル.(Devroye. L.)(1984)「対数凹密度を有する確率変数を生成する簡単なアルゴリズム(A simple algorithm for generating random variates with a log-concave density)」、コンピューティング(Computing)、33巻、247~257頁
[3]ドブロイ,エル.(Devroye. L.)(1986)「非一様な確率変数の生成(Non-Uniform Random Variate Generation)」、スプリンガー(Springer)、ニューヨーク(New York)
[4]エルダートン,ダブリュ.ピー.及びジョンソン,エヌ.エル.(Elderton, W. P. and Johnson, N. P.)(1969)「頻度曲線の系(Systems of Frequency Curves)」、ケンブリッジ大学出版(Cambridge University Press)
[5]ファング ケイ.ティー,コッツ,エス及びング,ケイ.ダブリュ.(Fang K. -T, Kotz, S. and Ng, K. W.)(1990)「対称多変量及び関連する分布(Symmetric Multivariate and Related Distributions)」、チャップマン及びホール(Chapman & Hall)、ロンドン(London)
[6]フライシュマン,エイ.(Fleishmann, A.)(1978)「非正規分布のシミュレーション方法(A method for simulating non-normal distributions)」、サイコメトリカ(Psychometrika)、43巻、521~532頁
[7]ホードリー,エイ.ビー.(Hoadley, A. B.)(1968)「左端と最初の3つのモーメントが知られた歪んだ密度を近似するためのピアソン分布の使用(Use of the Pearson densities for approximating a skew density whose left terminal and first three moments are known)」、バイオメトリカ(Biometrika)、55巻、3号、559~563頁
[8]ジョンソン,エヌ.エル,コッツ,エス.及びバラクリシュナン,エヌ.(Johnson, N. L, Kotz, S. and Balakrishnan, N.)(1994)「連続単変量分布1(Continuous univariate distribution)」、第2版、ジョン・ワイリー(John Wiley)、
[9]マグナス,ジェイ.アール.及びニューデッカー,エイチ.(Magnus, J. R. and Neudecker, H.)(1988)「統計及び計量経済学における応用を有する行列微分法(Matrix Differential Calculus with Application in Statistics and Econometrics)」、ワイリー(Wiley)、ニューヨーク(New York)
[10]マクグラス,イー.ジェイ.及びアービング,ディー.シー(McGrath, E. J. and Irving, D. C.)(1973)「効果的なモンテカルロ・シミュレーションの技術 第2巻:選択した確率分布についての乱数の生成(Techniques for efficient Monte Carlo simulation: volume 2: random number generation for selected probability distributions)」、テクニカル・リポート(Technical Report SAI-72-590-LJ)、サイエンス・アプリケーションズ社(Science Applications, Inc.),ラヨラ(La Jolla)、カリフォルニア(CA)
[11]ムーアヘッド,アール.ジェー.(Muirhead, R. J.)(1982)「多変量統計理論の諸問題(Aspscts of Multivariate Statistical Theory)」、ワイリー(Wiley)、ニューヨーク(New York)
[12]永原,ワイ.(Nagahara, Y.)(1999)「ピアソンIV型分布の確率密度関数及び分布関数とパラメータの最尤推定(The PDF and CF of Pearson type IV distributions and the ML estimation of the parameters)」、スタティスティクス・アンド・プロバビリティー・レターズ(Statistics & Probability Letters)、43巻、251~264頁
[13]永原,ワイ(Nagahara, Y)(2002)「ピアソン分布系に基づく非ガウスフィルタ及びスムーザー(Non-Gaussian filter and smoother based on the Pearson distribution system)」、ジャーナル・オブ・タイム・シリーズ・アナリシス(Journal of Time Series Analysis)、近刊
[14]オード,ジェイ.ケイ.(Ord, J. K.)(1972)「頻度分布の族(Family of Frequency Distributions)」、グリフィン(Griffin)出版
[15]パリシュ,アール.(Parrish, R.)(1983)「ピアソン分布についてのメンバー選択とパラメータ推定への積分アプローチ(On an integrated approach to member selection and parameter estimation for Pearson distributions)」、コンピュテイショナル・スタティスティクス及びデータ解析(Computational Statistics & Data Analysis)、1巻、239~255頁
[16]ピアソン,イー.エス.(Pearson, E. S.)(1963)「モーメントを用いた確率分布近似で生じるいくつかの問題(Some problems arising in approximating to probability distributions, using moments)」、バイオメトリカ(Biometrika)、50巻、95~111頁
[17]ピアソン,ケイ(Pearson, K.)(1895)「均質な物質における歪んだ分散についての研究(Memoir on skew variation in homogeneous material)」、フィロソフィカル・トランザクションズ・オブ・ザ・ローヤル・ササイアティー(Phil. Trans. Roy. Soc.)、エイ(A)、186巻、343~414頁
[18]ピアソン,ケイ(Pearson, K.)(1914)「統計及び計量生物学のための数表(Tables for Statistics and Biometrics)」、1巻、ケンブリッジ大学出版(Cambridge University Press)
[19]ピアソン,イー.エス.及びハートリー,エイチ.オー.(Pearson, E. S. and Hartrey, H. O. )(1954)「統計のための計量生物学数表(Biometrika Tables for Statistics)」、1巻、ケンブリッジ大学出版(Cambridge University Press)
[20]ショラック,ジー.アール.(Shorack, G. R.)(2000)「統計のための確率(Probability for Statistics)」、スプリンガー(Springer)、ニューヨーク(New York)
[21]スチュアート,エイ及びオード,ジェイ.ケイ.(Stuart, A. and Ord, J. K.)(1994)「ケンダルによる統計の進んだ理論(Kendall's Acvanced Theory of Statistics)」、1巻、分布理論(Distribution Theory)、6版、エドワード・アーノルド(Edward Arnold)
[22]タジカマラ,ピー.アール.(Tadikamalla, P. R.)(1980)「非正規分布のシミュレーションについて(On simulating non-normal distributions)」、サイコメトリカ(Psychometrika)、45巻、273~279頁
[23]ユアン,ケイ.エイチ.及びベントラー,ピー.エム.(Yuan K, -H. and Bentler P. M.)(1997)「特定の周辺の歪度及び尖度を有する多変量分布の生成(Generating multivariate distributions with specified marginal skewness and kurtosis)」、ソフトスタット‘97-統計ソフトウェアにおける発展(ダブリュ.バンディラ及びエフ.ファウルバウム編集)内(in SoftStat'97-Advances in Statistical Software 6, (W. Bandilla and F. Faulbaum, Eds.))、385~391頁、ルシアス・アンド・ルシアス(Lucius & Lucius)、シュトゥットガルト(Stuttgart)、ドイツ
[24]ユアン,ケイ.エイチ.及びベントラー,ピー.エム.(Yuan K, -H. and Bentler P. M.)(1999a)「非正規分布における2つのクラスの下での共分散構造解析における正規理論及び付随する検定統計(On normal theory and associated test statistics in covariance structure analysis under two classes of nonnormal distributions)」、スタティスティカ・シニカ(Statist. Sinica)、9巻、831~853頁
[25]ユアン,ケイ.エイチ.及びベントラー,ピー.エム.(Yuan K, -H. and Bentler P. M.)(1999b)「いくつかの非正規分布の下での共分散構造解析における正規理論MLEの漸近的分布について(On asymptotic distributions of normal theory MLE in covariance structure analysis under some nonnormal distributions)」、スタティスティクス・アンド・プロバビリティー・レターズ(Statistics & Probability Letters)、42巻、107~113頁
[26]ユアン,ケイ.エイチ.及びベントラー,ピー.エム.(Yuan K, -H. and Bentler P. M.)(2000)「非正規分布のいくつかのクラスにおける相関係数の推定(Inferences on correlation coefficients in some classes of nonnormal distributions)」、ジャーナル・オブ・マルチヴァリエイト・アナリシス(Journal of Multivariate Analysis)、72巻、230~248頁
[27]ベイル,ディー及びモーレイ,ブイ.エイ.(Vale, D. and Maurelli, V. A.)(1983)「多変量非正規分布シミュレーション(Simulating Multivariate Nonormal Distributions)」、サイコメトリカ(Psycometrika)、48巻、465~471頁

産業上の利用分野


本発明は、多変量非正規分布に従う乱数発生方法及び多変量非正規分布のパラメータ推定方法、ピアソンIV型分布に従う乱数発生方法及びコンピュータ・プログラムに関する。

特許請求の範囲 【請求項1】
コンピュータ上でユアン及びベントラーの方法Iに基づいて多変量非正規分布に従う乱数を発生する乱数発生方法において、
コンピュータを用いてn次元経験分布にn次元多変量非正規分布をあてはめるステップと、
コンピュータを用いて乱数を発生するステップとを有し、
前記あてはめるステップは、前記経験分布の3次及び4次モーメントに関するあてはめについて、次の関係式(1)及び(2)を用いる
多変量非正規分布に従う乱数発生方法。
【数1】


【数2】


ただし、E(・)は、期待値であり(以下同じ)、vech(・)は、対称行列から重複しない行列要素を取り出したベクトルであり、Dは、n次のデュプリケーション行列であり、D
+は、Dのムーア・ペンローズの一般化逆行列であり、
【数3】


は、クロネッカー積であり、Eiiは、eiを第i列単位としてeii
'である。
ここで、前記ユアン及びベントラーの方法Iは、次の通りである。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。
ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。
階数nのn行m列(m≧n)の非確率行列T=(tij)は、行列Σ=(σij)について、TT´=Σを満たす。ただし、行列T´はTの転置行列である。
このとき、次の式(3)で与えられる確率ベクトルX=(x1,・・・,xn)´は、Cov(X)=Σを満たす。ただし、Cov(・)は、ベクトルの分散共分散行列であり、ξ=(ξ1,・・・,ξm)´である。
【数4】



【請求項2】
確率変数ξ1,・・・,ξm及びνには、ピアソン分布系に従う乱数を使う請求項1に記載の多変量非正規分布に従う乱数発生方法。

【請求項3】
確率変数ξ1,・・・,ξm及びνには、少なくとも2以上の型のピアソン分布を用いる請求項2に記載の多変量非正規分布に従う乱数発生方法。

【請求項4】
前記あてはめるステップは、式(4)及び(5)でそれぞれ与えられるn次元の経験分布(確率ベクトル(X,・・・,X)´とする。)の3次及び4次モーメントに関して、式(6)の値を最小にするパラメータT,ζ=(ζ1,・・・,ζ),γ,κ=(κ1,・・・,κ),βの少なくとも1つを求める請求項1乃至3のいずれか1項に記載の多変量非正規分布に従う乱数発生方法。
【数5】


【数6】


【数7】


ただし、fijk(T,ζ,γ,κ,β)及びfijkl(T,ζ,γ,κ,β)は、3次モーメントE(xijk)及び4次モーメントE(xijk)にそれぞれ対応する関係式(1)及び(2)による表現とする。また、wijk及びwijklは、所定の重みである。

【請求項5】
ijk=1,wijkl=1である請求項4に記載の多変量非正規分布に従う乱数発生方法。

【請求項6】
n=2であって、前記あてはめるステップは、式(7)及び(8)でそれぞれ与えられる2次元経験分布(確率ベクトル(X,X)´とする。)の3次及び4次モーメントと、式(9)及び(10)で与えられる3次モーメントE(xijk)及び4次モーメントE(xijk)にそれぞれ対応する表現とについて、式(11)で与えられる値を最小にするようにパラメータT,ζ=(ζ1,ζ),γ,κ=(κ1,κ),βの少なくとも1つを求める請求項4又は5に記載の多変量非正規分布に従う乱数発生方法
【数8】


【数9】


【数10】


【数11】


【数12】



【請求項7】
i=1である請求項6に記載の多変量非正規分布に従う乱数発生方法。

【請求項8】
パラメータT,ζ=(ζ,・・・,ζ),γ,κ=(κ,・・・,κ),βの少なくとも1つを最尤法により推定する請求項1乃至7のいずれか1項に記載の多変量非正規分布に従う乱数発生方法。

【請求項9】
式(12)で与えられるκの値を用いて、表1に示すκの値と分布の型との対応関係の少なくとも1つを用いて前記確率変数の属する型を決定する請求項1乃至8のいずれか1項に記載の多変量非正規分布に従う乱数発生方法。
【数13】


【表1】


ただし、βとβは、それぞれ歪度の2乗と尖度であり、式(12)及び表1におけるκは、前記パラメータκ=(κ1,・・・,κ)とは異なる。

【請求項10】
コンピュータを用いて、式(13)に示すピアソンIV型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、
コンピュータを用いて、棄却法により乱数を発生するステップと、
によってピアソンIV型分布に従う乱数を発生させる請求項1乃至のいずれか1項に記載の多変量非正規分布に従う乱数発生方法。
【数14】



【請求項11】
前記乱数を発生するステップは、コンピュータを用いて発生した、加算生成法、M系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数に基づく請求項1乃至1のいずれか1項に記載の乱数発生方法。

【請求項12】
コンピュータを用いて、ユアン及びベントラーの方法Iに基づいて、次の関係式(31)及び(32)を用いてn次元経験分布にn次元多変量非正規分布をあてはめるステップと、
コンピュータを用いて、n次元経験分布がX=νTξ(ここで、ξは、m(ここで、m≧nである。)個の確率変数ξ~ξからなる確率ベクトル(ξ,・・・,ξ)´である。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。また、ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。)を満たす非正規分布を有すると仮定して、前記Xについて乱数を計算するステップと、
コンピュータを用いて、各確率変数x~xについての空間を、所定間隔Δhで分割し、このn次元空間を(Δh)の超立方体へ分割するステップと、
コンピュータを用いて、各データ・ベクトルの属する区間(Δh)内に存在する乱数の個数を乱数の総数で除して、当該区間(Δh)の確率P(θ)を計算するステップと、
コンピュータを用いて、前記確率P(θ)をn次元体積(Δh)で除して、前記データ・ベクトルの属する区間(Δh)の尤度f(θ)を計算するステップと、
コンピュータを用いて、前記尤度f(θ)の積Πi=1(θ)又は対数尤度の和Σi=1Nlog fi(θ)が最大になるようなパラメータθを計算するステップと、
を有する多変量非正規分布のパラメータ推定方法。
【数35】


【数36】


ただし、E(・)は、期待値であり(以下同じ)、vech(・)は、対称行列から重複しない行列要素を取り出したベクトルであり、Dは、n次のデュプリケーション行列であり、D
+は、Dのムーア・ペンローズの一般化逆行列であり、
【数37】


は、クロネッカー積であり、Eiiは、eiを第i列単位としてeii
'である。
ここで、前記ユアン及びベントラーの方法Iは、次の通りである。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。
ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。
階数nのn行m列(m≧n)の非確率行列T=(tij)は、行列Σ=(σij)について、TT´=Σを満たす。ただし、行列T´はTの転置行列である。
このとき、次の式(33)で与えられる確率ベクトルX=(x1,・・・,xn)´は、Cov(X)=Σを満たす。ただし、Cov(・)は、ベクトルの分散共分散行列であり、ξ=(ξ1,・・・,ξm)´である。
【数38】



【請求項13】
確率変数ξ1,・・・,ξm及びνには、ピアソン分布系に従う乱数を使う請求項12に記載の多変量非正規分布のパラメータ推定方法。

【請求項14】
確率変数ξ1,・・・,ξm及びνには、少なくとも2以上の型のピアソン分布を用いる請求項13に記載の多変量非正規分布のパラメータ推定方法。

【請求項15】
コンピュータを用いて、式(34)に示すピアソンIV型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、
コンピュータを用いて、棄却法により乱数を発生するステップと、
によってピアソンIV型分布に従う乱数を発生させる請求項12乃至14のいずれか1項に記載の多変量非正規分布のパラメータ推定方法。
【数39】



【請求項16】
コンピュータを用いて加算生成法、M系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する請求項12乃至15のいずれか1項に記載の多変量非正規分布のパラメータ推定方法。

【請求項17】
請求項1乃至16のいずれか1項に記載された手順を実現する、コンピュータにより実行可能なコンピュータ・プログラム。

【請求項18】
請求項17記載のコンピュータ・プログラムを記録した記録媒体
国際特許分類(IPC)
Fターム
画像

※ 画像をクリックすると拡大します。

JP2003082816thum.jpg
出願権利状態 登録
掲載中の発明について更に詳しい内容の説明を御希望の際は、お気軽にお問い合せください。


PAGE TOP

close
close
close
close
close
close
close