(57)【特許請求の範囲】
【請求項１】
コンピュータ上でユアン及びベントラーの方法Ｉに基づいて多変量非正規分布に従う乱数を発生する乱数発生方法において、
コンピュータを用いてｎ次元経験分布にｎ次元多変量非正規分布をあてはめるステップと、
コンピュータを用いて乱数を発生するステップとを有し、
前記あてはめるステップは、前記経験分布の３次及び４次モーメントに関するあてはめについて、次の関係式（１）及び（２）を用いる
多変量非正規分布に従う乱数発生方法。
【数１】
0003708928.tif 000195
【数２】
0003708928.tif 000196
ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄｎは、ｎ次のデュプリケーション行列であり、Ｄｎ+は、Ｄｎのムーア・ペンローズの一般化逆行列であり、
【数３】
0003708928.tif 000197
は、クロネッカー積であり、Ｅiiは、ｅiを第ｉ列単位としてｅiｅi'である。
ここで、前記ユアン及びベントラーの方法Ｉは、次の通りである。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。
ξjに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。
階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔij）は、行列Σ＝（σij）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。
このとき、次の式（３）で与えられる確率ベクトルＸ＝（ｘ1，・・・，ｘn）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ1，・・・，ξm）´である。
【数４】
0003708928.tif 000198
【請求項２】
確率変数ξ1，・・・，ξm及びνには、ピアソン分布系に従う乱数を使う請求項１に記載の多変量非正規分布に従う乱数発生方法。
【請求項３】
確率変数ξ1，・・・，ξm及びνには、少なくとも２以上の型のピアソン分布を用いる請求項２に記載の多変量非正規分布に従う乱数発生方法。
【請求項４】
前記あてはめるステップは、式（４）及び（５）でそれぞれ与えられるｎ次元の経験分布（確率ベクトル（Ｘ１，・・・，Ｘｎ）´とする。）の３次及び４次モーメントに関して、式（６）の値を最小にするパラメータＴ，ζ＝（ζ1，・・・，ζｍ），γ，κ＝（κ1，・・・，κｍ），βの少なくとも１つを求める請求項１乃至３のいずれか１項に記載の多変量非正規分布に従う乱数発生方法。
【数５】
0003708928.tif 000199
【数６】
0003708928.tif 000200
【数７】
0003708928.tif 000201
ただし、ｆijk（Ｔ，ζ，γ，κ，β）及びｆijkl（Ｔ，ζ，γ，κ，β）は、３次モーメントＥ（ｘiｘjｘk）及び４次モーメントＥ（ｘiｘjｘkｘｌ）にそれぞれ対応する関係式（１）及び（２）による表現とする。また、ｗijk及びｗijklは、所定の重みである。
【請求項５】
ｗijk＝１，ｗijkl＝１である請求項４に記載の多変量非正規分布に従う乱数発生方法。
【請求項６】
ｎ＝２であって、前記あてはめるステップは、式（７）及び（８）でそれぞれ与えられる２次元経験分布（確率ベクトル（Ｘ１，Ｘ２）´とする。）の３次及び４次モーメントと、式（９）及び（１０）で与えられる３次モーメントＥ（ｘiｘjｘk）及び４次モーメントＥ（ｘiｘjｘkｘｌ）にそれぞれ対応する表現とについて、式（１１）で与えられる値を最小にするようにパラメータＴ，ζ＝（ζ1，ζ２），γ，κ＝（κ1，κ２），βの少なくとも１つを求める請求項４又は５に記載の多変量非正規分布に従う乱数発生方法。
【数８】
0003708928.tif 000202
【数９】
0003708928.tif 000203
【数１０】
0003708928.tif 000204
【数１１】
0003708928.tif 000205
【数１２】
0003708928.tif 000206
【請求項７】
ｗi＝１である請求項６に記載の多変量非正規分布に従う乱数発生方法。
【請求項８】
パラメータＴ，ζ＝（ζ１，・・・，ζｍ），γ，κ＝（κ１，・・・，κｍ），βの少なくとも１つを最尤法により推定する請求項１乃至７のいずれか１項に記載の多変量非正規分布に従う乱数発生方法。
【請求項９】
式（１２）で与えられるκの値を用いて、表１に示すκの値と分布の型との対応関係の少なくとも１つを用いて前記確率変数の属する型を決定する請求項１乃至８のいずれか１項に記載の多変量非正規分布に従う乱数発生方法。
【数１３】
0003708928.tif 000207
【表１】
0003708928.tif 000208
ただし、β１とβ２は、それぞれ歪度の２乗と尖度であり、式（１２）及び表１におけるκは、前記パラメータκ＝（κ1，・・・，κｍ）とは異なる。
【請求項１０】
コンピュータを用いて、式（１３）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、
コンピュータを用いて、棄却法により乱数を発生するステップと、
によってピアソンＩＶ型分布に従う乱数を発生させる請求項１乃至９のいずれか１項に記載の多変量非正規分布に従う乱数発生方法。
【数１４】
0003708928.tif 000209
【請求項１１】
前記乱数を発生するステップは、コンピュータを用いて発生した、加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数に基づく請求項１乃至１０のいずれか１項に記載の乱数発生方法。
【請求項１２】
コンピュータを用いて、ユアン及びベントラーの方法Ｉに基づいて、次の関係式（３１）及び（３２）を用いてｎ次元経験分布にｎ次元多変量非正規分布をあてはめるステップと、
コンピュータを用いて、ｎ次元経験分布がＸ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ≧ｎである。）個の確率変数ξ１〜ξｍからなる確率ベクトル（ξ１，・・・，ξｍ）´である。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。また、ξjに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。）を満たす非正規分布を有すると仮定して、前記Ｘについて乱数を計算するステップと、
コンピュータを用いて、各確率変数ｘ１〜ｘｎについての空間を、所定間隔Δｈで分割し、このｎ次元空間を（Δｈ）ｎの超立方体へ分割するステップと、
コンピュータを用いて、各データ・ベクトルの属する区間（Δｈ）ｎ内に存在する乱数の個数を乱数の総数で除して、当該区間（Δｈ）ｎの確率Ｐｉ（θ）を計算するステップと、
コンピュータを用いて、前記確率Ｐｉ（θ）をｎ次元体積（Δｈ）ｎで除して、前記データ・ベクトルの属する区間（Δｈ）ｎの尤度ｆｉ（θ）を計算するステップと、
コンピュータを用いて、前記尤度ｆｉ（θ）の積Πｉ＝１Ｎｆｉ（θ）又は対数尤度の和Σi=1Nｌｏｇ  ｆi（θ）が最大になるようなパラメータθを計算するステップと、
を有する多変量非正規分布のパラメータ推定方法。
【数３５】
0003708928.tif 000210
【数３６】
0003708928.tif 000211
ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄｎは、ｎ次のデュプリケーション行列であり、Ｄｎ+は、Ｄｎのムーア・ペンローズの一般化逆行列であり、
【数３７】
0003708928.tif 000212
は、クロネッカー積であり、Ｅiiは、ｅiを第ｉ列単位としてｅiｅi'である。
ここで、前記ユアン及びベントラーの方法Ｉは、次の通りである。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。
ξjに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。
階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔij）は、行列Σ＝（σij）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。
このとき、次の式（３３）で与えられる確率ベクトルＸ＝（ｘ1，・・・，ｘn）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ1，・・・，ξm）´である。
【数３８】
0003708928.tif 000213
【請求項１３】
確率変数ξ1，・・・，ξm及びνには、ピアソン分布系に従う乱数を使う請求項１２に記載の多変量非正規分布のパラメータ推定方法。
【請求項１４】
確率変数ξ1，・・・，ξm及びνには、少なくとも２以上の型のピアソン分布を用いる請求項１３に記載の多変量非正規分布のパラメータ推定方法。
【請求項１５】
コンピュータを用いて、式（３４）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、
コンピュータを用いて、棄却法により乱数を発生するステップと、
によってピアソンＩＶ型分布に従う乱数を発生させる請求項１２乃至１４のいずれか１項に記載の多変量非正規分布のパラメータ推定方法。
【数３９】
0003708928.tif 000214
【請求項１６】
コンピュータを用いて加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する請求項１２乃至１５のいずれか１項に記載の多変量非正規分布のパラメータ推定方法。
【請求項１７】
請求項１乃至１６のいずれか１項に記載された手順を実現する、コンピュータにより実行可能なコンピュータ・プログラム。
【請求項１８】
請求項１７記載のコンピュータ・プログラムを記録した記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、多変量非正規分布に従う乱数発生方法及び多変量非正規分布のパラメータ推定方法、ピアソンＩＶ型分布に従う乱数発生方法及びコンピュータ・プログラムに関する。
【０００２】
【従来の技術】
従来、多変量非正規分布を生成する方法は、数えるしか存在していない。このような分布は、多変量解析と統計モデリングの分野において、切望されている。
【０００３】
フライシュマン（Fleishmann）（１９７８）、タジカマラ（Tadikamalla）（１９８０）、ベイルとマブレリ（Vale and Mavrelli）（１９８３）は、多変量非正規分布を生成する方法を提案した。これらの方法は、シミュレーションした母集団の確率ベクトルの４次のモーメント行列を計算することが難しいという欠点がある。ユアン及びベントラー（１９９７）は、４次のモーメント行列を考慮し、特定の周辺の歪度と尖度を有する多変量非正規分布を生成する手法を提案した。
【０００４】
本明細書では、ユアン等が方法Ｉと称したものは、次の通りである。ξ1，
・・・，ξmは、独立な確率変数であり、Ｅ（・）を期待値として、
Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj及びξ＝（ξ1，・・・，ξ1）´である。νは、ξに独立な確率変数であり、
Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。また、Ｔ＝（ｔij）は、階数ｐのｐ行ｍ列の非確率行列であり、ＴＴ´＝Σとｍ≧ｐを満たす。すると、次の式（３５）で与えられる確率ベクトルは、一般に、非楕円分布に従い、Ｃｏｖ（Ｘ）＝Σを満たす。ここで、Ｃｏｖ（・）は、確率ベクトルＸ＝（ｘ1，・・・，ｘn）´の分散共分散行列である。
【０００５】
【数４０】
0003708928.tif 000002
ｘiの周辺の歪度（skew）と尖度（kurt）は、それぞれ次の式（３６）及び
（３７）によって与えられる。
【０００６】
【数４１】
0003708928.tif 000003
【数４２】
0003708928.tif 000004
Ｘは、ξ〜Ｎm（０，Ｉ）及びＰｒ（ν≧０）＝１（ただし、Ｎmはｍ次元正規分布、Ｐｒは確率を表す。）が満たされるとき、楕円分布に従う（ファン、コッツ及びング（Fang, Kotz and Ng）（１９９０））。ν＝１とし、
ｚ1，・・・，ｚpを独立な標準正規変数とし、
ξ＝｜ｚ1｜−Ｅ（｜ｚ1｜），ξj＝ｚj，ｊ＝２，・・・，ｐとすると、式（３５）におけるＸは、アザリニ及びバレ（Azzalini and Valle）（１９９６）によって定義された歪んだ正規分布に従う。
【０００７】
この方法には、母集団の４次モーメント行列を容易に計算できるという利点がある。Ｓを式（３５）のＸからの大きさｎの標本に基づいた標本共分散とする。
ｖｅｃｈ（・）を対称行列から重複しない要素を取り出してベクトルに変換する演算子とし、ｓ＝ｖｅｃｈ（Ｓ）とする。すると、ｓの漸近的分散共分散行列はΓ／ｎによって与えられる。ここで、
Γ＝ｖａｒ（ｖｅｃｈ（ＸＸ´））であり、次の式（３８）によって与えられる。ｖａｒは、分散である。
【０００８】
【数４３】
0003708928.tif 000005
ただし、ｔjは、Ｔ＝（ｔij）の第ｊ列ベクトルであり、
【数４４】
0003708928.tif 000006
は、クロネッカー積である。Ｄpは、ｐ次のデュプリケーション行列であり、Ｄp+は、Ｄpのムーア・ペンローズの一般化逆行列である。
【０００９】
この後、ユアン等は、検定統計の理論を研究した（ユアン及びベントラー（１９９９ａ，１９９９ｂ，２０００））。ユアン等は、確率ベクトルを生成するにはフライシュマン（１９７８）の方法を推奨している。しかし、この分布は、あまり広い歪度及び尖度を有していない。これに対して、ピアソン分布系は、次のような特徴を有する。
【００１０】
（１）ピアソン分布系は、様々な歪度と尖度を有する広いクラスの分布を表現することができる。ピアソン分布系は、ガンマ分布、ベータ分布、ｔ分布等のいくつかの良く知られた分布を含む。
【００１１】
（２）ピアソン分布系は、４つのモーメントによって特徴付けられる。この特徴は、ユアン等の方法を適用するためには十分である。
【００１２】
（３）ＩＶ型を除いたピアソン分布系に従う乱数は、正規分布とガンマ分布に従う乱数を用いて生成される。
【００１３】
（４）最近、実現することが難しかったピアソンＩＶ型分布を含む殆ど全ての種類を用いる実用的な手法が開発された（永原（２００２））。
【００１４】
ここで、ピアソン分布系を導入する。ケイ・ピアソン（K. Pearson）（１８９５）は、確率密度関数ｐに関する次の式（３９）の微分方程式に従い、ピアソン分布系を定義した。
【００１５】
【数４５】
0003708928.tif 000007
ピアソン第ＩＶ型分布は、正規化定数を計算することが難しいので、実現することが難しかったが、論文（永原（１９９９））において、計算に用いる再帰的公式を導出し、様々な問題を解決した。この後、非ガウスフィルタに用いるピアソンＩＶ型分布を含む分布系の殆ど全ての種類を用いる実用的な手法を開発した（永原（２００２））。これらの研究により、ピアソン分布系が実用的に利用できるようになった。
【００１６】
ここで、ピアソン分布系に関して参照される研究者と文献を掲げる（ピアソン（１９１４）、ピアソン及びハートリー（１９５４）、ピアソン（１９６３）、ホードリ−（１９６８）、エルダートン及びジョンソン（１９６９）、オード（１９７２）、パリシュ（１９８３）、ジョンソン、コッツ及びバラカリシュナン（１９９４）、スチュアート及びオード（１９９４））。
【００１７】
参考文献としては、次のものが挙げられる。
【００１８】
［１］アザリニ，エイ．及びバレ，エイ．ディー．（Azzalini, A. and Valle, A. D.）（１９９６）「多変量の歪んだ正規分布（The multivariate skew-normal distribution）」、バイオメトリカ（Biometrika）、８３巻、７１５〜７２６頁
［２］ドブロイ，エル．（Devroye. L.）（１９８４）「対数凹密度を有する確率変数を生成する簡単なアルゴリズム（A simple algorithm for generating random variates with a log-concave density）」、コンピューティング（Computing）、３３巻、２４７〜２５７頁
［３］ドブロイ，エル．（Devroye. L.）（１９８６）「非一様な確率変数の生成（Non-Uniform Random Variate Generation）」、スプリンガー（Springer）、ニューヨーク（New York）
［４］エルダートン，ダブリュ．ピー．及びジョンソン，エヌ．エル．（Elderton, W. P. and Johnson, N. P.）（１９６９）「頻度曲線の系（Systems of Frequency Curves）」、ケンブリッジ大学出版（Cambridge University Press）
［５］ファング  ケイ．ティー，コッツ，エス及びング，ケイ．ダブリュ．（Fang K. −T, Kotz, S. and Ng, K. W.）（１９９０）「対称多変量及び関連する分布（Symmetric Multivariate and Related Distributions）」、チャップマン及びホール（Chapman & Hall）、ロンドン（London）
［６］フライシュマン，エイ．（Fleishmann, A.）（１９７８）「非正規分布のシミュレーション方法（A method for simulating non-normal distributions）」、サイコメトリカ（Psychometrika）、４３巻、５２１〜５３２頁
［７］ホードリー，エイ．ビー．（Hoadley, A. B.）（１９６８）「左端と最初の３つのモーメントが知られた歪んだ密度を近似するためのピアソン分布の使用（Use of the Pearson densities for approximating a skew density whose left terminal and first three moments are known）」、バイオメトリカ（Biometrika）、５５巻、３号、５５９〜５６３頁
［８］ジョンソン，エヌ．エル，コッツ，エス．及びバラクリシュナン，エヌ.（Johnson, N. L, Kotz, S. and Balakrishnan, N.）（１９９４）「連続単変量分布１（Continuous univariate distribution）」、第２版、ジョン・ワイリー（John Wiley）、
［９］マグナス，ジェイ．アール．及びニューデッカー，エイチ．（Magnus, J. R. and Neudecker, H.）（１９８８）「統計及び計量経済学における応用を有する行列微分法（Matrix Differential Calculus with Application in Statistics and Econometrics）」、ワイリー（Wiley）、ニューヨーク（New York）
［１０］マクグラス，イー．ジェイ．及びアービング，ディー．シー（McGrath, E. J. and Irving, D. C.）（１９７３）「効果的なモンテカルロ・シミュレーションの技術 第２巻：選択した確率分布についての乱数の生成（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）
［１１］ムーアヘッド，アール．ジェー．（Muirhead, R. J.）（１９８２）「多変量統計理論の諸問題（Aspscts of Multivariate Statistical Theory）」、ワイリー（Wiley）、ニューヨーク（New York）
［１２］永原，ワイ．（Nagahara, Y.）（１９９９）「ピアソンＩＶ型分布の確率密度関数及び分布関数とパラメータの最尤推定（The PDF and CF of Pearson type IV distributions and the ML estimation of the parameters）」、スタティスティクス・アンド・プロバビリティー・レターズ（Statistics & Probability Letters）、４３巻、２５１〜２６４頁
［１３］永原，ワイ（Nagahara, Y）（２００２）「ピアソン分布系に基づく非ガウスフィルタ及びスムーザー（Non-Gaussian filter and smoother based on the Pearson distribution system）」、ジャーナル・オブ・タイム・シリーズ・アナリシス（Journal of Time Series Analysis）、近刊
［１４］オード，ジェイ．ケイ．（Ord, J. K.）（１９７２）「頻度分布の族（Family of Frequency Distributions）」、グリフィン（Griffin）出版
［１５］パリシュ，アール．（Parrish, R.）（１９８３）「ピアソン分布についてのメンバー選択とパラメータ推定への積分アプローチ（On an integrated approach to member selection and parameter estimation for Pearson distributions）」、コンピュテイショナル・スタティスティクス及びデータ解析（Computational Statistics & Data Analysis）、１巻、２３９〜２５５頁
［１６］ピアソン，イー．エス．（Pearson, E. S.）（１９６３）「モーメントを用いた確率分布近似で生じるいくつかの問題（Some problems arising in approximating to probability distributions, using moments）」、バイオメトリカ（Biometrika）、５０巻、９５〜１１１頁
［１７］ピアソン，ケイ（Pearson, K.）（１８９５）「均質な物質における歪んだ分散についての研究（Memoir on skew variation in homogeneous material）」、フィロソフィカル・トランザクションズ・オブ・ザ・ローヤル・ササイアティー（Phil. Trans. Roy. Soc.）、エイ（Ａ）、１８６巻、３４３〜４１４頁
［１８］ピアソン，ケイ（Pearson, K.）（１９１４）「統計及び計量生物学のための数表（Tables for Statistics and Biometrics）」、１巻、ケンブリッジ大学出版（Cambridge University Press）
［１９］ピアソン，イー．エス．及びハートリー，エイチ．オー．（Pearson, E. S. and Hartrey, H. O. ）（１９５４）「統計のための計量生物学数表（Biometrika Tables for Statistics）」、１巻、ケンブリッジ大学出版（Cambridge University Press）
［２０］ショラック，ジー．アール．（Shorack, G. R.）（２０００）「統計のための確率（Probability for Statistics）」、スプリンガー（Springer）、ニューヨーク（New York）
［２１］スチュアート，エイ及びオード，ジェイ．ケイ．（Stuart, A. and Ord, J. K.）（１９９４）「ケンダルによる統計の進んだ理論（Kendall's Acvanced Theory of Statistics）」、１巻、分布理論（Distribution Theory）、６版、エドワード・アーノルド（Edward Arnold）
［２２］タジカマラ，ピー．アール．（Tadikamalla, P. R.）（１９８０）「非正規分布のシミュレーションについて（On simulating non-normal distributions）」、サイコメトリカ（Psychometrika）、４５巻、２７３〜２７９頁
［２３］ユアン，ケイ．エイチ．及びベントラー，ピー．エム．（Yuan K, -H. and Bentler P. M.）（１９９７）「特定の周辺の歪度及び尖度を有する多変量分布の生成（Generating multivariate distributions with specified marginal skewness and kurtosis）」、ソフトスタット‘９７−統計ソフトウェアにおける発展（ダブリュ．バンディラ及びエフ．ファウルバウム編集）内（in SoftStat'97-Advances in Statistical Software 6, （W. Bandilla and F. Faulbaum, Eds.））、３８５〜３９１頁、ルシアス・アンド・ルシアス（Lucius & Lucius）、シュトゥットガルト（Stuttgart）、ドイツ
［２４］ユアン，ケイ．エイチ．及びベントラー，ピー．エム．（Yuan K, -H. and Bentler P. M.）（１９９９ａ）「非正規分布における２つのクラスの下での共分散構造解析における正規理論及び付随する検定統計（On normal theory and associated test statistics in covariance structure analysis under two classes of nonnormal distributions）」、スタティスティカ・シニカ（Statist. Sinica）、９巻、８３１〜８５３頁
［２５］ユアン，ケイ．エイチ．及びベントラー，ピー．エム．（Yuan K, -H. and Bentler P. M.）（１９９９ｂ）「いくつかの非正規分布の下での共分散構造解析における正規理論ＭＬＥの漸近的分布について（On asymptotic distributions of normal theory MLE in covariance structure analysis under some nonnormal distributions）」、スタティスティクス・アンド・プロバビリティー・レターズ（Statistics & Probability Letters）、４２巻、１０７〜１１３頁
［２６］ユアン，ケイ．エイチ．及びベントラー，ピー．エム．（Yuan K, -H. and Bentler P. M.）（２０００）「非正規分布のいくつかのクラスにおける相関係数の推定（Inferences on correlation coefficients in some classes of nonnormal distributions）」、ジャーナル・オブ・マルチヴァリエイト・アナリシス（Journal of Multivariate Analysis）、７２巻、２３０〜２４８頁
［２７］ベイル，ディー及びモーレイ，ブイ．エイ．（Vale, D. and Maurelli, V. A.）（１９８３）「多変量非正規分布シミュレーション（Simulating Multivariate Nonormal Distributions）」、サイコメトリカ（Psycometrika）、４８巻、４６５〜４７１頁
【発明が解決しようとする課題】
ところで、検定統計のみならず統計モデリングの見地からも、多変量非正規分
布を導出し、経験分布を近似する分布を構築することがさらに求められている。
また、少数の標本からも多変量非正規分布を安定して推定することが求められて
いる。
【００１９】
本発明は、前述の実情に鑑みて提案されるものであって、ピアソン分布系及び一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンＩＶ型分布に従う乱数発生方法及びコンピュータ・プログラムを提供することを目的とする。
【００２０】
【課題を解決するための手段】
前述の課題を解決するために、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータ上でユアン及びベントラーの方法Ｉに基づいて多変量非正規分布に従う乱数を発生する乱数発生方法において、コンピュータを用いてｎ次元経験分布にｎ次元多変量非正規分布をあてはめるステップと、コンピュータを用いて乱数を発生するステップとを有し、前記あてはめるステップは、前記経験分布の３次及び４次モーメントに関するあてはめについて、次の関係式（４０）及び（４１）を用いる。
【００２１】
【数４６】
0003708928.tif 000008
【数４７】
0003708928.tif 000009
ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄｎは、ｎ次のデュプリケーション行列であり、Ｄｎ+は、Ｄｎのムーア・ペンローズの一般化逆行列であり、
【数４８】
0003708928.tif 000010
は、クロネッカー積であり、Ｅiiは、ｅiを第ｉ列単位としてｅiｅi'である。
【００２２】
ここで、前記ユアン及びベントラーの方法Ｉは、次の通りである。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。
【００２３】
ξjに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。
【００２４】
階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔij）は、行列Σ＝（σij）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。
【００２５】
このとき、次の式（４２）で与えられる確率ベクトルＸ＝（ｘ1，・・・，ｘn）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ1，・・・，ξm）´である。
【００２６】
【数４９】
0003708928.tif 000011
好ましくは、確率変数ξ1，・・・，ξm及びνには、ピアソン分布系に従う乱数を使う。
【００２７】
好ましくは、確率変数ξ1，・・・，ξm及びνには、少なくとも２以上の型のピアソン分布を用いる。
【００２８】
好ましくは、前記あてはめるステップは、式（４３）及び（４４）でそれぞれ与えられるｎ次元の経験分布（確率ベクトル（Ｘ１，・・・，Ｘｎ）´とする。
）の３次及び４次モーメントに関して、式（４５）の値を最小にするパラメータＴ，ζ＝（ζ1，・・・，ζｍ），γ，κ＝（κ1，・・・，κｍ），βの少なくとも１つを求める。
【００２９】
【数５０】
0003708928.tif 000012
【数５１】
0003708928.tif 000013
【数５２】
0003708928.tif 000014
ただし、ｆijk（Ｔ，ζ，γ，κ，β）及びｆijkl（Ｔ，ζ，γ，κ，β）は、３次モーメントＥ（ｘiｘjｘk）及び４次モーメントＥ（ｘiｘjｘkｘｌ）にそれぞれ対応する関係式（４０）及び（４１）による表現とする。また、ｗijk及びｗijklは、所定の重みである。
ここで、記号
【数５３】
0003708928.tif 000015
は、
【数５４】
0003708928.tif 000016
（ただし、ｉ≦ｊ≦ｋ）の意味で、１≦ｉ≦ｎ，１≦ｊ≦ｎ，１≦ｋ≦ｎ，ｉ≦ｊ≦ｋの各条件を満たす全てのｉ，ｊ，ｋの組み合わせについて式に代入し、それらの総和を取ったものである。同様に、記号
【数５５】
0003708928.tif 000017
は、
【数５６】
0003708928.tif 000018
（ただし、ｉ≦ｊ≦ｋ≦ｌ）の意味で、１≦ｉ≦ｎ，１≦ｊ≦ｎ，１≦ｋ≦ｎ，１≦ｌ≦ｎ，ｉ≦ｊ≦ｋ≦ｌの各条件を満たす全てのｉ，ｊ，ｋ，ｌの組み合わせについて式に代入し、それらの総和を取ったものである。
また、前記ｍｉｊｋ＝Ｅ（ＸiＸjＸk）及びｍｉｊｋｌ＝
Ｅ（ＸiＸjＸkＸl）の意味について説明する。データの分布
（Ｘ１，・・・，Ｘｎ）´は、次の式（Ａ１）のようにＮ個のデータ・ベクトル
【数５７】
0003708928.tif 000019
（ただし、
【数５８】
0003708928.tif 000020
は各々のデータ（実現値）を示す。１≦ｑ≦N）から構成されている。
【数５９】
0003708928.tif 000021
前記ｍｉｊｋ及びｍｉｊｋｌは、このようなデータ・ベクトルの要素となるデータを用いて次の式（Ａ２）及び（Ａ３）によって計算される。
【数６０】
0003708928.tif 000022
【数６１】
0003708928.tif 000023
【００３０】
好ましくは、ｗijk＝１，ｗijkl＝１である。
【００３１】
好ましくは、ｎ＝２であって、前記あてはめるステップは、式（４６）及び
（４７）でそれぞれ与えられる２次元経験分布（確率ベクトル（Ｘ１，Ｘ２）´とする。）の３次及び４次モーメントと、式（４８）及び（４９）で与えられる３次モーメントＥ（ｘiｘjｘk）及び４次モーメントＥ（ｘiｘjｘkｘｌ）にそれぞれ対応する表現とについて、式（５０）で与えられる値を最小にするようにパラメータＴ，ζ＝（ζ1，ζ２），γ，κ＝（κ1，κ２），βの少なくとも１つを求める。
【００３２】
【数６２】
0003708928.tif 000024
【数６３】
0003708928.tif 000025
【数６４】
0003708928.tif 000026
【数６５】
0003708928.tif 000027
【数６６】
0003708928.tif 000028
好ましくは、ｗi＝１である。
【００３３】
好ましくは、パラメータＴ，ζ＝（ζ１，・・・，ζｍ），γ，κ＝
（κ１，・・・，κｍ），βの少なくとも１つを最尤法により推定する。
【００３４】
好ましくは、式（５１）で与えられるκの値を用いて、表８に示すκの値と分布の型との対応関係の少なくとも１つを用いて前記確率変数の属する型を決定する。
【００３５】
【数６７】
0003708928.tif 000029
【表８】
0003708928.tif 000030
ただし、β１とβ２は、それぞれ歪度の２乗と尖度であり、式（５１）及び表８におけるκは、前記パラメータκ＝（κ1，・・・，κｍ）とは異なる。
【００３６】
好ましくは、前記決定した型に応じて、次の表９及び表１０に示す生成法の少なくとも１つを用いて乱数Ｚを発生する。
【００３７】
【表９】
0003708928.tif 000031
【表１０】
0003708928.tif 000032
好ましくは、コンピュータを用いて、式（５２）に示すピアソンＩＶ型分布の
正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、
コンピュータを用いて、棄却法により乱数を発生するステップと、によってピア
ソンＩＶ型分布に従う乱数を発生させる。
【００３８】
【数６８】
0003708928.tif 000033
好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、式（５３）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンＩＶ型分布に従う乱数を発生させる。
【００３９】
【数６９】
0003708928.tif 000034
好ましくは、前記乱数を発生するステップは、コンピュータを用いて発生した、加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数に基づく。
【００４０】
好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、ｎ個のデータ
【数７０】
0003708928.tif 000035
からなるデータ・ベクトルＸｉ´について、経験分布のデータ｛Ｘ１´，・・・，ＸＮ´｝を取得するステップと、コンピュータを用いて、前記データ｛Ｘ１´，・・・，ＸＮ´｝を標準化して｛Ｘ１，・・・，ＸＮ｝とするステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ１，・・・，ＸＮ｝に基づいて、分散共分散行列Σを計算するステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ１，・・・，ＸＮ｝に基づいて、３次モーメントｍｉｊｋ（１≦ｉ≦ｊ≦ｋ≦ｎ）を計算するステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ１，・・・，ＸＮ｝に基づいて、４次モーメントｍｉｊｋｌ（１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）を計算するステップと、コンピュータを用いて、前記分散共分散行列Σから行列Ｔを計算するステップと、前記データ・ベクトルが従う分布が、Ｘ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ＞ｎである。）個の確率変数ξ１〜ξｍからなる確率ベクトル（ξ１，・・・，ξｍ）´である。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。また、ξjに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。）を満たす非正規分布を有すると仮定して、コンピュータを用いて、３次モーメントｍｉｊｋ（１≦ｉ≦ｊ≦ｋ≦ｎ）及び４次モーメントｍｉｊｋｌ（１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）とｆｉｊｋ（Ｔ，ζ，γ，κ，β），ｆｉｊｋｌ（Ｔ，ζ，γ，κ，β）（ここでζ＝（ζ１，・・・，ζｍ），κ＝（κ１，・・・，κｍ）である。）との差を損失とする損失関数を導入し、ＴＴ´＝Σ（Ｔ´はＴの転置行列である。）という条件の下で、全体として評価するリスク関数を最小化するようにパラメータＴ，（ζ１，・・・，ζｍ），γ，（κ１，・・・，κｍ），βを決定するステップと、コンピュータを用いて、前記決定されたパラメータ（ζ１，・・・，ζｍ），γ，（κ１，・・・，κｍ），βに基づいて、式（５４）及び表１１により（β１とβ２は、それぞれ歪度の２乗と尖度である。）、前記確率ベクトル（ξ１，・・・，ξｍ）´及び確率変数νがピアソン分布のどの型に属するかを決定するステップと、
【数７１】
0003708928.tif 000036
【表１１】
0003708928.tif 000037
コンピュータを用いて乱数を発生し、この乱数に基づいて前記確率ベクトル（ξ１，・・・，ξｍ）´及び確率変数νについて乱数を計算するステップと、コンピュータを用いて、Ｘ＝νＴξに基づいて、標準化されたＸの乱数を計算するステップと、コンピュータを用いて、標準化されたＸの乱数を標準化前の乱数に変換するステップと、を有する。
【００４１】
ただし、式（５４）及び表１１において、β１とβ２は、それぞれ歪度の２乗と尖度であり、κは、前記パラメータκ＝（κ1，・・・，κｍ）とは異なる。
【００４２】
好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、ｎ個のデータ
【数７２】
0003708928.tif 000038
からなるデータ・ベクトルＸｉ´について、経験分布のデータ
｛Ｘ１´，・・・，ＸＮ´｝を取得するステップと、コンピュータを用いて、前記データ｛Ｘ１´，・・・，ＸＮ´｝を標準化して｛Ｘ１，・・・，ＸＮ｝とするステップと、コンピュータを用いて、前記標準化されたデータ
｛Ｘ１，・・・，ＸＮ｝に基づいて、分散共分散行列Σを計算するステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ１，・・・，ＸＮ｝に基づいて、３次モーメントｍｉｊｋ（１≦ｉ≦ｊ≦ｋ≦ｎ）を計算するステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ１，・・・，ＸＮ｝に基づいて、４次モーメントｍｉｊｋｌ（１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）を計算するステップと、コンピュータを用いて、前記分散共分散行列ΣからＴ＝Σ１／２（行列の平方根）を満たす行列Ｔを計算するステップと、前記データ・ベクトルが従う分布が、Ｘ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ＝ｎである。）個の確率変数ξ１〜ξｍからなる確率ベクトル（ξ１，・・・，ξｍ）´である。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，
Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。また、ξjに独立な確率変数νは、パラメータγ及びβについて、
Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。）を満たす非正規分布を有すると仮定して、コンピュータを用いて、３次モーメントｍｉｊｋ（１≦ｉ≦ｊ≦ｋ≦ｎ）及び４次モーメントｍｉｊｋｌ（１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）とｆｉｊｋ（ζ，γ，κ，β），ｆｉｊｋｌ（ζ，γ，κ，β）（ここでζ＝（ζ１，・・・，ζｍ），κ＝（κ１，・・・，κｍ）である。）との差を損失とする損失関数を導入し、全体として評価するリスク関数を最小化するようにパラメータ（ζ１，・・・，ζｍ），γ，（κ１，・・・，κｍ），βを決定するステップと、コンピュータを用いて、前記決定されたパラメータ
（ζ１，・・・，ζｍ），γ，（κ１，・・・，κｍ），βに基づいて、式（５５）及び表１２により（β１とβ２は、それぞれ歪度の２乗と尖度である。）、前記確率ベクトル（ξ１，・・・，ξｍ）´及び確率変数νがピアソン分布のどの型に属するかを決定するステップと、
【数７３】
0003708928.tif 000039
【表１２】
0003708928.tif 000040
コンピュータを用いて乱数を発生し、この乱数に基づいて前記確率ベクトル
（ξ１，・・・，ξｍ）´及び確率変数νについて乱数を計算するステップと、コンピュータを用いて、Ｘ＝νＴξに基づいて、標準化されたＸの乱数を計算するステップと、コンピュータを用いて、標準化されたＸを標準化前の乱数に変換するステップと、を有する。
ここで、経験分布のデータ｛Ｘ１´，・・・，ＸＮ´｝とは、次の意味である。
【数７４】
0003708928.tif 000041
ここでは、Ｘｒ´（１≦ｒ≦Ｎ）はデータ・ベクトルを表し、先のデータの分布（Ｘ１，・・・，ＸＮ）´での表記Ｘとは異なる。標準化されたデータ｛Ｘ１，・・・，ＸＮ｝も同様の意味である。
【００４３】
ただし、式（５５）及び表１２において、β１とβ２は、それぞれ歪度の２乗と尖度であり、κは、前記パラメータκ＝（κ1，・・・κｍ）とは異なる。
【００４４】
好ましくは、前記損失関数（Ｌ）及びリスク関数（Ｒ）は、式（５６）及び（５７）、式（５８）及び（５９）、又は式（６０）及び（６１）のいずれかの組で与えられる請求項１３又は１４に記載の多変量非正規分布のパラメータ推定方法。
【００４５】
【数７５】
0003708928.tif 000042
【数７６】
0003708928.tif 000043
【数７７】
0003708928.tif 000044
【数７８】
0003708928.tif 000045
【数７９】
0003708928.tif 000046
【数８０】
0003708928.tif 000047
好ましくは、前記決定した型に応じて、次の表１３及び表１４に示す生成法の少なくとも１つを用いて乱数Ｚを発生する。
【００４６】
【表１３】
0003708928.tif 000048
【表１４】
0003708928.tif 000049
好ましくは、コンピュータを用いて、式（６２）に示すピアソンＩＶ型分布の
正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、
コンピュータを用いて、棄却法により乱数を発生するステップと、によってピア
ソンＩＶ型分布に従う乱数を発生させる。
【００４７】
【数８１】
0003708928.tif 000050
好ましくは、前記乱数を発生するステップは、加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。
【００４８】
本発明に係る多変量非正規分布のパラメータ推定方法は、コンピュータを用いて、式（６３）に示すように、所定の確率分布に従う確率変数ν、階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ、所定の確率分布に従う確率ベクトルξ＝（ξ1，・・・，ξm）´の積として与えられる確率ベクトルＸ＝（ｘ1，・・・，ｘn）´で与えられる多次元非正規分布を計算するステップと、
【数８２】
0003708928.tif 000051
最尤法によりパラメータを推定するステップと、を有する。
【００４９】
好ましくは、前記最尤法によりパラメータを推定するステップは、コンピュータを用いて、ｎ次元空間を分割するステップと、コンピュータを用いて、経験分布のデータ｛Ｘ1，・・・，ＸN｝について、Ｘiが属する各分割に属する乱数の数を乱数の総数で除して、次の式（６４）で与えられる当該分割の確率を求めるステップと、
【数８３】
0003708928.tif 000052
コンピュータを用いて、前記分割の確率を当該分割のＮ次元体積で除して前記Ｘiが属する、次の式（６５）で与えられる分割の尤度ｆi（θ）を求めるステップと、
【数８４】
0003708928.tif 000053
コンピュータを用いて、尤度の積Πi=1Nｆi（θ）又は対数尤度の和Σi=1Nｌｏｇ  ｆi（θ）が最大になるようなパラメータθを推定するステップと、を有する。
【００５０】
好ましくは、確率変数ξ1，・・・，ξm及びνには、ピアソン分布系に従う乱数を使う。
【００５１】
好ましくは、確率変数ξ1，・・・，ξm及びνには、少なくとも２以上の型のピアソン分布を用いる。
【００５２】
好ましくは、コンピュータを用いて、式（６６）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンＩＶ型分布に従う乱数を発生させる。
【００５３】
【数８５】
0003708928.tif 000054
好ましくは、コンピュータを用いて、ユアン及びベントラーの方法Ｉに基づいて、次の関係式（６７）及び（６８）を用いてｎ次元経験分布にｎ次元多変量非正規分布をあてはめるステップを有し、このあてはめるステップで求めたパラメータの近傍について前記パラメータθを推定する。
【００５４】
【数８６】
0003708928.tif 000055
【数８７】
0003708928.tif 000056
ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄｎは、ｎ次のデュプリケーション行列であり、Ｄｎ+は、Ｄｎのムーア・ペンローズの一般化逆行列であり、
【数８８】
0003708928.tif 000057
は、クロネッカー積であり、Ｅiiは、ｅiを第ｉ列単位としてｅiｅi'である。
【００５５】
ここで、前記ユアン及びベントラーの方法Ｉは、次の通りである。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。
【００５６】
ξjに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。
【００５７】
階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔij）は、行列Σ＝（σij）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。
【００５８】
このとき、次の式（６９）で与えられる確率ベクトルＸ＝（ｘ1，・・・，ｘn）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ1，・・・，ξm）´である。
【００５９】
【数８９】
0003708928.tif 000058
好ましくは、前記パラメータθは、パラメータＴ，ζ＝（ζ1，・・・，ζｍ），γ，κ＝（κ1，・・・，κｍ），βの少なくとも１つを表す。
【００６０】
好ましくは、コンピュータを用いて加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。
【００６１】
好ましくは、本発明に係る多変量非正規分布のパラメータ推定方法は、コンピュータを用いて、ユアン及びベントラーの方法Ｉに基づいて、次の関係式（７０）及び（７１）を用いてｎ次元経験分布にｎ次元多変量非正規分布をあてはめるステップと、コンピュータを用いて、ｎ次元経験分布がＸ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ≧ｎである。）個の確率変数ξ１〜ξｍからなる確率ベクトル（ξ１，・・・，ξｍ）´である。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，
Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。また、ξjに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。）を満たす非正規分布を有すると仮定して、前記Ｘについて乱数を計算するステップと、コンピュータを用いて、各確率変数ｘ１〜ｘｎについての空間を、所定間隔Δｈで分割し、このｎ次元空間を（Δｈ）ｎの超立方体へ分割するステップと、コンピュータを用いて、各データ・ベクトルの属する区間（Δｈ）ｎ内に存在する乱数の個数を乱数の総数で除して、当該区間（Δｈ）ｎの確率Ｐｉ（θ）を計算するステップと、コンピュータを用いて、前記確率Ｐｉ（θ）をｎ次元体積（Δｈ）ｎで除して、前記データ・ベクトルの属する区間（Δｈ）ｎの尤度ｆｉ（θ）を計算するステップと、コンピュータを用いて、前記尤度ｆｉ（θ）の積Πｉ＝１Ｎｆｉ（θ）又は対数尤度の和Σi=1Nｌｏｇ  ｆi（θ）が最大になるようなパラメータθを計算するステップと、を有する。
【００６２】
【数９０】
0003708928.tif 000059
【数９１】
0003708928.tif 000060
ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄｎは、ｎ次のデュプリケーション行列であり、Ｄｎ+は、Ｄｎのムーア・ペンローズの一般化逆行列であり、
【数９２】
0003708928.tif 000061
は、クロネッカー積であり、Ｅiiは、ｅiを第ｉ列単位としてｅiｅi'である。
【００６３】
ここで、前記ユアン及びベントラーの方法Ｉは、次の通りである。独立な確率変数ξ1，・・・，ξmは、パラメータζj及びκjについて、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj（１≦ｊ≦ｍ）を満たす。
【００６４】
ξjに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。
【００６５】
階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔij）は、行列Σ＝（σij）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。
【００６６】
このとき、次の式（７２）で与えられる確率ベクトルＸ＝（ｘ1，・・・，ｘn）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ1，・・・，ξm）´である。
【００６７】
【数９３】
0003708928.tif 000062
好ましくは、確率変数ξ1，・・・，ξm及びνには、ピアソン分布系に従う乱数を使う。
【００６８】
好ましくは、確率変数ξ1，・・・，ξm及びνには、少なくとも２以上の型のピアソン分布を用いる。
【００６９】
好ましくは、コンピュータを用いて、式（７３）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンＩＶ型分布に従う乱数を発生させる。
【００７０】
【数９４】
0003708928.tif 000063
好ましくは、コンピュータを用いて加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。
【００７１】
本発明に係る予想損失額の計算方法は、前記手順に従い、ポートフォリオを構成するｎ種類の資産のｎ次元損益分布にｎ次元多変量非正規分布をあてはめ、シミュレーションによって予想損失額を計算する。
【００７２】
本発明に係る金融商品のシミュレーション方法は前記手順を用い、前記多変量非正規分布に従う乱数に基づいて金融商品の特性をコンピュータを用いてシミュレーションする。
【００７３】
本発明に係る天候・保険デリバティブの設計又は価格付け方法は、前記手順を用いて天候・保険デリバティブを設計又は価格付けする。
【００７４】
本発明に係るアセット・アロケーションの方法は、前記手順を用いてアセット・アロケーションを行う。
【００７５】
本発明に係るイオン注入のシミュレーション方法は、前記手順を用い、前記多変量非正規分布に従う乱数に基づいてイオン注入によるイオンの分布をコンピュータを用いてシミュレーションする。
【００７６】
本発明に係る能力評価の方法は、前記手順を用いて能力評価を行う。
【００７７】
本発明に係るコンピュータにより実行可能なコンピュータ・プログラムは、前記手順を実現するものである。好ましくは、前記コンピュータ・プログラムは、光ディスク、磁気ディスク、光磁気ディスク又は磁気テープの少なくとも１つによって提供される。
【００７８】
【発明の実施の形態】
以下、本発明に係る多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンＩＶ型分布に従う乱数発生方法及びコンピュータ・プログラムの実施の形態並びにリスク・マネジメントにおける予想損失額の計算方法について、図面を参照して詳細に説明する。
【００７９】
以下の本実施の形態は、大きくは次の構成からなっている。
【００８０】
１．ピアソン分布系及び乱数生成法
１．１  ピアソン分布系
１．２  ピアソン分布ＩＶ型に従う乱数の生成
１．３  ピアソン分布系に従う乱数の生成
２．  ピアソン分布系を用いた多変量非正規分布の生成
２．１  多変量正規分布の確率ベクトル
２．２  シミュレーション
２．３  応用のためのあてはめプロシージャ
３．  最尤法
４．  実施例
５．  変形例
５．１  金融商品のシミュレーション方法
５．２  半導体へのイオン注入
５．３  アセット・アロケーション
５．４  能力評価
本実施の形態は、様々な歪度と尖度を有する広いクラスの分布を表すことができるピアソン分布系を用い、多変量非正規分布をシミュレーションする方法を提案する。第１に、実現が困難であったピアソンＩＶ型分布を含むピアソン分布系に従う乱数を生成するプロシージャを導く。第２に、ピアソン分布系を適用し、ユアンとベントラー（１９９７）によって開発された方法を用い、多変量非正規分布に従う乱数を生成する。第３に、本願発明者によって求められた３次のモーメント行列を用い、種々の標本分布について３次と４次のモーメント行列をあてはめ、近似的な多変量非正規分布を構築することができる。第４に、最尤法によるパラメータ推定方法を提案する。このパラメータ推定方法によって、少数のデータからモーメント法によっては実現できない正確さで、多変量非正規分布を推定することができることを示す。そして、これらの方法が、様々なビジネスに応用できることを、リスク・マネジメントにおける予想損失額の計算について具体的に説明する。最後に、変形例として、金融商品のシミュレーション方法、半導体へのイオン注入、アセット・アロケーション及び能力評価について簡単に触れる。
【００８１】
１  ピアソン分布系及び乱数生成法
１．１  ピアソン分布系
ピアソン分布系の型は、次のパラメータκによって分類される。パラメータκは、次の式（７４）によって定義される。
【００８２】
【数９５】
0003708928.tif 000064
ただし、β1とβ2は、それぞれ歪度の２乗と尖度である。
【００８３】
ピアソン分布は、κ＜０，０＜κ＜１及び１＜κについて、それぞれＩ型、ＩＶ型及びＶＩ型と呼ばれる。これらの３つの場合は、「主要型」と呼ばれる（エルダートン及びジョンソン（Elderton and Johnson）（１９６９））。ＩＩＩ型（κ＝±∞）は、Ｉ型とＶＩ型の境界にある。Ｖ型（κ＝１）は、第ＩＶ型と第ＶＩ型の境界にある。κ＝０の場合については、ＶＩＩ型、ＩＩ型及び正規分布が含まれる。
【００８４】
ピアソン分布系においては、型が指定されると、この分布の３つ又は４つのパラメータは、平均、分散、歪度及び尖度から決定される。歪度の２乗と尖度の間の典型的な関係は、多数の教科書に示されている（ジョンソン、コッツ及びバラカリシュナン（Johnson, Kotz and Balakrishnan）（１９９４）、スチュアート及びオード（Stuart and Ord）（１９９４）等）。
【００８５】
まず、表１５〜１８に示すように、ピアソン分布系は、最も適切な表現によると、８つの型に分類され、密度関数、及び平均、分散、歪度及び尖度による陽な形が与えられる。
【００８６】
次に、モーメントから各型のパラメータを得るための変換式とプロシージャを示す。各型の変換式を表１５〜１８に示す。変換式によって、全てのパラメータを順に容易に計算することができる。
【００８７】
【表１５】
0003708928.tif 000065
【表１６】
0003708928.tif 000066
表１６において、ピアソンＩ型分布について、ｐ，ｑ＞１のとき単峰型、０＜
ｑ≦１＜ｐのときＪ型、０＜ｐ≦１＜ｑのときＪ型、０＜ｐ，ｑ≦１のときＵ型
になる。なお、ｐ＝ｑ＝１のときには一様分布になる。また、ピアソンＩＩ型分
布について、ｐ＞１のとき単峰型、０＜ｐ＜１のときＵ型になる。なお、ｐ＝１
のときには一様分布になる。
【００８８】
【表１７】
0003708928.tif 000067
【表１８】
0003708928.tif 000068
１．２  ピアソン分布ＩＶ型に従う乱数の生成
本実施の形態においては、ピアソン分布ＩＶ型の確率密度関数は、次の式（７５）によって定義される。ただし、Γ（・）はガンマ関数、Ｂ（・）はベータ関数である。
【００８９】
【数９６】
0003708928.tif 000069
ドブロイ（Devrye）（１９８６）によると、確率変数XがピアソンＩＶ型分布であってｂ＞１のとき、次の式（７６）により分布するとすると、
ａｒｃｔａｎ（ｘ）は、式（７７）によって定義される対数凹密度を有する。
【００９０】
【数９７】
0003708928.tif 000070
【数９８】
0003708928.tif 000071
ドブロイ（１９８４）は、対数凹密度についての棄却アルゴリズム（指数版）が、この分布に従う乱数を生成するのに適用できることを示唆した。しかし、その正規化定数は、著書において問題として残されていた（ドブロイ（１９８６）の練習２．７）。論文（永原（１９９９））では、この問題を解析的に解き、その正規化定数が容易に計算できることを示した（式（７８））。ピアソンＩＶ型分布に従う乱数の生成のアルゴリズムを次に示す。Ｃは、ピアソンＩＶ型分布の正規化定数を示す。
【００９１】
【数９９】
0003708928.tif 000072
0003708928.tif 000073
このアルゴリズムは、正規化定数の解析解又はこの展開を計算するステップと、棄却法により乱数を発生するステップを含んでいる。
【００９２】
ここで、乱数のシーズについて説明する。ある確率分布に従う乱数を発生させるにあたっては、一様乱数がその基本になっている。
【００９３】
原理的には、一様乱数が与えられれば、逆関数法によって、様々な確率分布に従う乱数が発生できる。本発明においても、この一様乱数をもとにして、ピアソン分布系の確率分布に従う乱数の発生のために使う指数分布、ガンマ分布、正規分布に従う乱数を発生している。この一様乱数を発生させる方法には、大別すると物理乱数、擬似乱数及び準乱数及び低ディスクレパンシー列が含まれる。準乱数及び低ディスクレパンシー列とは、擬似乱数とは違って、ランダムな性質はもっていないが、一様に分布する点列である。
【００９４】
擬似乱数には、合同法、加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法（Generalized Feedback Shift-register; ＧＦＳＲ）法及びメルセンヌツイスターが含まれる。準乱数には、リチトメヤー（Ritchtmyer）列、ヘイセルグローブ（Haselgrove）列、ハマーズレイ（Hammersley）列、低ディスクレパンシー（Low Discrepancy）列である、ホルトン（Halton）列、一般化ホルトン（Halton）列、ソボル（Sobol）列、改良ソボル（Sobol）列、フォーレ（Faure）列、一般化フォーレ（Faure）列、一般化ニーダーライター（Niederreiter）列が含まれる。近年、低ディスクレパンシー列は、多く研究されてきて、一つの定義として定着してきている。
【００９５】
なお、メルセンヌツイスターの参考文献には、マツモト，エム．及びニシムラ，ティー．（Matsumoto, M. and Nishimura, T）「メルセンヌツイスター：  ６２３次元等分布一様擬似乱数発生（Mersenne twister: A 623-dimensionally equidistributed uniform pseudorandom number generator）」、１９９８年、エイシーエム・トランザクションズ・オン・モデリング・アンド・コンピューター・シミュレーション（ACM Transaction on Modeling and Computer Simulation）、８巻、３〜３０頁がある。
【００９６】
また、ソボル列の参考文献には、プレス，ダブリュ．，ツコルスキー，エス．，ヴェタリング，ダブリュ．，フラネリー，ビー．（Press, W., Teukolsky, S., Vetterling, W., Flannery, B）「フォートラン７７の数値処理法（Numerical Recipes in Fortran 77）」、１９９２年、第２版、ケンブリッジ大学出版（Cambridge Univ. Press.）がある。
【００９７】
この中で、一般的な方法は、擬似乱数の合同法である。ところが、最近では、合同法よりはるかに長い周期を持つ、ＧＦＳＲ法やメルセンヌツイスターが開発されてきている。さらに、準乱数や低ディスクレパンシー列が研究されていて、実用にも供されてきている。さて、これらの周期のより長い擬似乱数（例えばメルセンヌツイスター）や準乱数（低ディスクレパンシー列）を本実施の形態に適用したところ、本実施の形態の実施に重要な乱数の性質を大幅に改善することができた。従来の合同法では実現できなかった正確性で、目標とする確率分布の平均、分散、歪度、尖度を再現した乱数を得ることができた。
【００９８】
本実施の形態の一番の重要な点は、非常に多数の乱数で、もはや解析的には表現できないような複雑な性質を持った確率分布を表現することである。この目的のためには、従来の合同法よりさらに正確な確率分布を再現できる乱数が望まれ、特に本実施の形態の特徴である、歪度、尖度が正規分布と違う分布（ピアソンＩＶ型等）については、擬似乱数では合同法より長い周期のもの（例えばメルセンヌツイスター）、又は準乱数を使うことによって大幅に改善された、より実用的な正確性が達成された。
【００９９】
以下に１次元の例を述べる。まず、ピアソンＩＶ型に関する例では、乱数の個数が多くなると、合同法に比べ、メルセンヌツイスター、準乱数（ここでは、ソボル列を使用する。ソボル列は、低ディスクレパンシー列でもある。）を使った方が、平均、分散、歪度、尖度すべてにわたって目標とする値（表の一段目の数）に大幅に近くなっている。
【０１００】
例えば、表２１（合同法）、表２５（メルセンヌツイスター）、表２９（ソボル列）において、１０００００個の乱数で１０００回シミュレーションした結果、平均について、合同法では０．０１０２３、メルセンヌツイスターでは０．０００１８、準乱数では、−０．００００１となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは４桁、準乱数では、５桁までと大幅に改善している。同様に、分散について、合同法では１．０３４３３、メルセンヌツイスターでは１．００００９、準乱数では、０．９９９９９となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは５桁、準乱数では、６桁までと大幅に改善している。
【０１０１】
また、本実施の形態で重要な歪度と尖度についてはさらに改善している。同じ表で、歪度について目標値は０．７５で、合同法では０．７９６５１、メルセンヌツイスターでは０．７５０４１、準乱数では、０．７５００７となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは４桁、準乱数では、５桁までと大幅に改善している。尖度について目標値は５．０で、合同法では５．３６３３４、メルセンヌツイスターでは５．０００５７、準乱数では、５．００２５７となり、合同法での小数点以下１桁目の正確さが、メルセンヌツイスターでは４桁、準乱数でも、３桁までと大幅に改善している。また、メルセンヌツイスター及び準乱数では、乱数の数が増えれば、それに従って目標の値にそれぞれ収束している点が合同法より優れており、この性質は、本実施の形態において乱数が持つべき最も重要な性質である。
【０１０２】
ここで、モンテカルロ・シミュレーション技法として対照変量（負の相関）法、制御変量法、モーメント適合法、回帰分析法、非線形回帰分析法、加重抽出法、層化抽出法、ラテン・ハイパーキューブ法、マルチンゲール分散減少法、条件付モンテカルロ法のような分散減少法が容易に適用できる。これら技法は、乱数の発生法に依らず、計算を改良する方法である。
【０１０３】
表１９〜３０において、数は、乱数の個数を表し、値は、１０００回の平均である。各括弧内の値は各値の標準偏差である。この結果、アルゴリズムを用いて妥当な乱数を得た。
【０１０４】
【表１９】
0003708928.tif 000074
【表２０】
0003708928.tif 000075
【表２１】
0003708928.tif 000076
【表２２】
0003708928.tif 000077
【表２３】
0003708928.tif 000078
【表２４】
0003708928.tif 000079
【表２５】
0003708928.tif 000080
【表２６】
0003708928.tif 000081
【表２７】
0003708928.tif 000082
【表２８】
0003708928.tif 000083
【表２９】
0003708928.tif 000084
【表３０】
0003708928.tif 000085
１．３  ピアソン分布系に従う乱数の生成
ＩＶ型を除くピアソン分布系については、正規分布とガンマ分布に従う乱数を用いることで、容易に乱数を生成できることが知られている（マクグラス及びアービング（McGrath and Irving）（１９７３））。本実施の形態では、表１５〜１８のような定義を実現する実用的なプロシージャを開発する。これらの結果は、表３１と表３２に示されている。典型的な例は、表３３〜５３に示されている。数は、乱数の個数を示し、値は１０００回の平均である。各括弧内の値は、各値の標準偏差である。この結果、アルゴリズムを用いて妥当な乱数を得る。
【０１０５】
ピアソンＩＶ型以外の分布系については、合同法に比べ、メルセンヌツイスター、準乱数（ここでは、ソボル列を使用。ソボル列は低ディスクレパンシー列でもある。）を使った方が、平均、分散では、ほとんどすべてで大幅に改善した。重要な歪度、尖度については、幅の広い領域を受け持つピアソンＶＩ型及びピアソンＩ型では、１０００００個の乱数で１０００回シミュレーションの結果で、まず、ピアソンＶＩ型で、表３３（合同法）、表４０（メルセンヌツイスター）及び表４４（ソボル列）において、歪度について目標値は１．０で、合同法では０．９９３２３、メルセンヌツイスターでは０．９９９７５、準乱数では０．９９９６８となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数でも、同様に４桁までと大幅に改善している。同じ表で、尖度について目標値は４．８で、合同法では４．７３３８６、メルセンヌツイスターでは４．７９８６３、準乱数では４．７９８３６となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは３桁、準乱数でも、同様に３桁までと改善している。
【０１０６】
次に、ピアソンＩ型について、表３４（合同法）、表４１（メルセンヌツイスター）及び表４８（ソボル列）において、歪度について目標値は１．０で、合同法では１．００６１０、メルセンヌツイスターでは０．９９９４５、準乱数では１．００００４となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数は５桁までと大幅に改善している。同じ表で、尖度について目標値は４．０で、合同法では４．００５８０、メルセンヌツイスターでは３．９９７３７、準乱数では４．０００３７となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは同様の３桁、準乱数では、４桁までと改善している。
【０１０７】
ピアソンＩＩ型では、表３５（合同法）、表４２（メルセンヌツイスター）及び表４９（ソボル列）において、歪度について目標値は０．０で、合同法では０．０００１２、メルセンヌツイスターでは−０．００００５、準乱数では０．００００５となり、合同法での小数点以下４桁目の正確さが、メルセンヌツイスターでは５桁、準乱数でも、同様に５桁までと大幅に改善している。尖度について目標値は２．５で、同じ表で、合同法では２．４９４９０、メルセンヌツイスターでは２．４９９９２、準乱数では２．５００１４となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは５桁、準乱数でも、４桁までと改善している。
【０１０８】
正規分布では、表３６（合同法）、表４３（メルセンヌツイスター）及び表５０（ソボル列）において、歪度について目標値は０．０で、合同法では０．００４３９、メルセンヌツイスターでは０．０００２３、準乱数では０．０００００となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数では６桁までと大幅に改善している。尖度について目標値は３．０で、合同法では２．９８１３０、メルセンヌツイスターでは３．０００１８、準乱数では３．００００１となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは４桁、準乱数でも、５桁までと改善している。
【０１０９】
ピアソンＶＩＩ型では、表３７（合同法）、表４４（メルセンヌツイスター）、表５１（ソボル列）において、歪度について目標値は０．０で、合同法では０．００２８２、メルセンヌツイスターでは０．０００６２、準乱数では−０．０００４７となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数でも４桁までと大幅に改善している。尖度について目標値は４．０で、合同法では３．９６６３９、メルセンヌツイスターでは３．９９６９０、準乱数では３．９９８９４となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは３桁、準乱数でも、３桁までと改善している。
【０１１０】
ピアソンＩＩＩ型では、表３８（合同法）、表４５（メルセンヌツイスター）、表５２（ソボル列）において、歪度について目標値は１．０で、合同法では１．００１２８、メルセンヌツイスターでは０．９９９５１、準乱数では０．９９９９６となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数では５桁までと大幅に改善している。尖度について目標値は４．５で、合同法では４．４５７１１、メルセンヌツイスターでは４．４９８６４、準乱数では４．４９９７８となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは３桁、準乱数では、４桁までと改善している。
【０１１１】
ピアソンＶ型では、表３９（合同法）、表４６（メルセンヌツイスター）、表５３（ソボル列）において、歪度について目標値は１．０で、合同法では１．００４５６、メルセンヌツイスターでは０．９９８４２、準乱数では０．９９９９２となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは同様の３桁、準乱数では５桁までと大幅に改善している。尖度について目標値は４．９７０で、合同法では５．１５５０５、メルセンヌツイスターでは４．９５８２９、準乱数では４．９６９２５となり、合同法での小数点以下１桁目の正確さが、メルセンヌツイスターでは２桁、準乱数では、４桁までと大幅に改善している。また、以上分布についても、ピアソンＩＶ型のときと同様に、メルセンヌツイスター及び準乱数では、乱数の数が増えれば、それに従って目標の値にそれぞれ収束している。この点が合同法より優れており、この性質は本発明において乱数が持つべき最も重要な性質である。
【０１１２】
【表３１】
0003708928.tif 000086
【表３２】
0003708928.tif 000087
【表３３】
0003708928.tif 000088
【表３４】
0003708928.tif 000089
【表３５】
0003708928.tif 000090
【表３６】
0003708928.tif 000091
【表３７】
0003708928.tif 000092
【表３８】
0003708928.tif 000093
【表３９】
0003708928.tif 000094
【表４０】
0003708928.tif 000095
【表４１】
0003708928.tif 000096
【表４２】
0003708928.tif 000097
【表４３】
0003708928.tif 000098
【表４４】
0003708928.tif 000099
【表４５】
0003708928.tif 000100
【表４６】
0003708928.tif 000101
【表４７】
0003708928.tif 000102
【表４８】
0003708928.tif 000103
【表４９】
0003708928.tif 000104
【表５０】
0003708928.tif 000105
【表５１】
0003708928.tif 000106
【表５２】
0003708928.tif 000107
【表５３】
0003708928.tif 000108
２  ピアソン分布系を用いた多変量非正規分布の生成
２．１  多変量非正規分布の確率ベクトル
ユアン等（１９９７）によると、ユアン等が方法Iと称したものは、次の通りである。ξ1，・・・，ξmは、独立な確率変数であり、Ｅ（ξj）＝０，Ｅ（ξj2）＝１，Ｅ（ξj3）＝ζj，Ｅ（ξj4）＝κj及びξ＝
（ξ1，・・・，ξm）´である。νは、ξに独立な確率変数であり、Ｅ（ν）＝０，Ｅ（ν2）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。ここで、この分布は、νが１のとき、ピアソン分布系の確率密度関数によって表されることが特記される。また、Ｔ＝（ｔij）は、階数ｐの非確率ｐ行ｍ列の行列であり、ＴＴ´＝Σ及びｍ≧ｐを満たす。すると、次の式（７９）で与えられる確率ベクトルは、一般に、Ｃｏｖ（Ｘ）＝Σを満たす非楕円分布に従うようになる。ここで、Ｃｏｖ（・）は、ベクトルＸ＝（ｘ1，・・・，ｘn）´の分散共分散行列である。
【０１１３】
【数１００】
0003708928.tif 000109
ｘiの周辺の歪度及び尖度は、それぞれ次の式（８０）及び（８１）によって与えられる。
【０１１４】
【数１０１】
0003708928.tif 000110
【数１０２】
0003708928.tif 000111
そして、４次のモーメント行列Γ＝ｖａｒ（ｖｅｃｈ（ＸＸ´））は、次の式（８２）によって与えられる。
【０１１５】
【数１０３】
0003708928.tif 000112
ただし、ｔjは、Ｔ＝（ｔij）の第ｊ列ベクトルであり、
【数１０４】
0003708928.tif 000113
は、クロネッカー積である。Ｄpは、ｐ次のデュプリケーション行列であり、
Ｄp+は、Ｄpのムーア・ペンローズの一般化逆行列である。
【０１１６】
統計モデリングの見地から、３次のモーメント行列が必要である。本願発明者は、３次のモーメント行列について、次の定理を得ることができた。
【０１１７】
定理３．１（３次のモーメント行列）３次のモーメント行列Λは、次の式（８３）によって与えられる。
【０１１８】
【数１０５】
0003708928.tif 000114
ここで、ｅiは、第ｉ列単位であり、Ｅii＝ｅiｅi'である。
【０１１９】
証明  ３次のモーメント行列Λは、Ｙ＝Ｔξとして、次の式（８４）のようになる。
【０１２０】
【数１０６】
0003708928.tif 000115
ＹＹ´＝Ｔξξ´Ｔ´及びｖｅｃｈ（ＹＹ´）＝Ｄp+ｖｅｃ（ＹＹ´）であるので、
【数１０７】
0003708928.tif 000116
となる。ここで、記号
【数１０８】
0003708928.tif 000117
は、行列のクロネッカー積を表す。以下でも同様である。
【０１２１】
したがって、次の式（８６）を得る。
【０１２２】
【数１０９】
0003708928.tif 000118
次の式（８７）は、明らかである。
【０１２３】
【数１１０】
0003708928.tif 000119
この定理は、式（８４），（８６）及び（８７）から得られる。
【０１２４】
証明終わり
ユアン及びベントラー（１９９９ａ，１９９９ｂ及び２０００）は、検定統計理論を研究した。ユアン等は、次の式（８８）を示唆した。
【０１２５】
【数１１１】
0003708928.tif 000120
ただし、Π＝∂ｆ／∂ｖｅｃｈ´（Σ）の階数ｓの行列として、Ω＝ΠΓΠ´であり、Ｓは、標本分散共分散行列である。また、ユアン等は、以前、Πｖｅｃｈ（ｔjｔj´）＝０，ｊ＝１，・・・，ｍを仮定すると、次の式
（８９）が成立することを示唆している。
【０１２６】
【数１１２】
0003708928.tif 000121
ただし、
【数１１３】
0003708928.tif 000122
である。
【０１２７】
本実施の形態は、様々な３次と４次のモーメント行列を用い、多変量非正規分布に従う乱数を得ることを目的としている。Ｔ，ξ，νを変化させることによって、異なる３次と４次のモーメント行列を有する様々な多変量非正規分布を得ることができる。まず、Ｘ1，Ｘ2，ξ1及びξ2を有するｐ＝ｍ＝２の典型的な場合についてシミュレーションを行う。次に、３次と４次のモーメントを経験データの１つにあてはめるプロシージャを提案する。
【０１２８】
なお、ピアソン分布系のほかには、4次までのモーメントで表せるジョンソン分布系やラムダ分布などがνやξの分布として採用できる。
【０１２９】
２．２  シミュレーション
本実施の形態では、ｘ1，ｘ2，ξ1及びξ2を有するｐ＝ｍ＝２の場合についてシミュレーションする。このシミュレーションにおいて、３次と４次のモーメント行列の要素は、次の通りである。
【０１３０】
［３次モーメント］
【数１１４】
0003708928.tif 000123
ただし、ｆ1，ｆ2，ｆ3，ｆ4は、Ｘ＝νＴξで発生される確率ベクトルＸ＝（ｘ１，ｘ２）´のそれぞれ３次モーメントＥ（ｘ13），Ｅ（ｘ23），Ｅ（ｘ12ｘ2），Ｅ（ｘ1ｘ22）の表現である。
【０１３１】
［４次モーメント］
【数１１５】
0003708928.tif 000124
ただし、ｆ5，ｆ6，ｆ7，ｆ8，ｆ9は、Ｘ＝νＴξで発生される確率ベクトルＸ＝（ｘ１，ｘ２）´のそれぞれ４次モーメントＥ（ｘ14），Ｅ（ｘ24），Ｅ（ｘ13ｘ2），Ｅ（ｘ12ｘ22），Ｅ（ｘ1ｘ23）の表現である。ここで、σ１１＝σ２２＝１であり、また、ｔ１１２＝ｔ２２２＝１−ｔ１２２＝１−ｔ２１２である。
【０１３２】
ユアンとベントラー（２０００）の例を採用し、ρ＝０．５及びｔ11＝ｔ22＝０．２５８８，ｔ12＝ｔ21＝０．９６５９とする。結果は、表５４〜８０に示されている。特定の歪度及び尖度を有する、各分布のξ1，ξ2，νから乱数を得た後、この乱数を平均と分散がそれぞれ０と１になるように標準化する。続いて、式Ｘ＝νＴξによって、ｘ1とｘ2のための乱数Ｘ１，Ｘ２を得る。各表において、Ｅ（・）は、これらの乱数を用いて計算される。括弧内は、３次と４次のモーメント行列の式を用いて計算された理論値である。これらの結果によると、１０，０００，０００点と１，０００，０００点の場合、この方法による乱数は、理論的なモーメント、すなわち分散共分散、３次と４次のモーメント行列を有する目的の分布を正確に表すことができる。
【０１３３】
ここで、異なるピアソン分布の型の組み合わせも可能であることは特記されるべきであることである。例えば、金融リスク、オペレーショナル・リスク、保険リスク、災害リスク、事業リスクは、一般に異なったピアソン分布の型で記述される。本実施の形態によると、異なったピアソン分布の組み合わせも可能なので、このような異なった種類のリスクの組み合わせも統合して管理することができる。例えば、各種の資産を含んだポートフォリオの予想損失額を統合して計算することができる。
【０１３４】
なお、金融での最大損失額の計算では、極値分布（ガンベル分布、フレッシュ分布、ワイブル分布）や一般化パレート分布・対数正規分布が使われるが、これらの分布もこの発明の方法で、νやξの分布として採用できる。このことによって、多変量化や、よりデータを正確に表す分布を推定できる。
【０１３５】
ピアソンＩＶ型の場合、２次元においてξ１、ξ２、νすべてに、ピアソンＩＶ型を用いた例を表５４（合同法）、表６３（メルセンヌツイスター）、表７２（準乱数  ソボル列）に示す。以下、ｘ１、ｘ２は確率変数、Ｘ１，Ｘ２は乱数を表す。
【０１３６】
まず、目標とする分布の特性はξ１（歪度０．７５、尖度５．０）、ξ２（歪度１．０、尖度５．５）、ν（歪度０．５０、尖度４．０）で、それから生成される多変量非正規分布のモーメント特性の理論値は、Ｅ（ｘ１３）＝０．４５７０８、Ｅ（ｘ２３）＝０．３４６６０、Ｅ（ｘ１２ｘ２）＝０．１４４９９、
Ｅ（ｘ１ｘ２２）＝０．１２２８９、Ｅ（ｘ１４）＝２０．７３８６５、
Ｅ（ｘ２４）＝１９．００６７８、Ｅ（ｘ１３ｘ２）＝８．４６５１５、
Ｅ（ｘ１２ｘ２２）＝７．１２３９０、Ｅ（ｘ１ｘ２３）＝８．０３２２１である。
【０１３７】
例えば、１０，０００，０００個のシミュレーションの場合、表５４（合同法）では、それぞれ、ξ１（歪度０．８３０８２、尖度５．２８９２８）、
ξ２（歪度１．０５９０６、尖度６．１９４９５）、ν（歪度０．５８８０７、尖度４．１７７７４）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ１３）＝０．５７５０７（０．４５７０８）、
Ｅ（Ｘ２３）＝０．４５４９９（０．３４６６０）、Ｅ（Ｘ１２Ｘ２）＝０．１８７１５（０．１４４９９）、Ｅ（Ｘ１Ｘ２２）＝０．１６３８３（０．１２２８９）、Ｅ（Ｘ１４）＝２４．１４４９３（２０．７３８６５）、
Ｅ（Ｘ２４）＝２０．９４８９７（１９．００６７８）、Ｅ（Ｘ１３Ｘ２）＝９．５３０５４（８．４６５１５）、Ｅ（Ｘ１２Ｘ２２）＝７．７１１６８（７．１２３９０）、Ｅ（Ｘ１Ｘ２３）＝８．７５５００（８．０３２２１）である。
【０１３８】
括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ１、ξ２、νの実際のシミュレーション時の歪度、尖度（例えば、上の例ではξ１（歪度０．８３０８２、尖度５．２８９２８）、ξ２（歪度１．０５９０６、尖度６．１９４９５）、ν（歪度０．５８８０７、尖度４．１７７７４））から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との開きは大きい。
【０１３９】
ところが、これをメルセンヌツイスターでシミュレーションすると、表６３（メルセンヌツイスター）にあるように、それぞれ、ξ１（歪度０．７５２７９、尖度５．０２１２３）、ξ２（歪度１．００００４、尖度５．５０８３１）、ν（歪度０．５０４１４、尖度４．００７８０）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ１３）＝０．４５１８５（０．４５７０８）、Ｅ（Ｘ２３）＝０．３５０１９（０．３４６６０）、
Ｅ（Ｘ１２Ｘ２）＝０．１４２０２（０．１４４９９）、
Ｅ（Ｘ１Ｘ２２）＝０．１２４０４（０．１２２８９）、Ｅ（Ｘ１４）＝２０．０３２１０（２０．７３８６５）、Ｅ（Ｘ２４）＝１９．１６３４１（１９．００６７８）、Ｅ（Ｘ１３Ｘ２）＝８．１６２４９（８．４６５１５）、
Ｅ（Ｘ１２Ｘ２２）＝７．０２１２７（７．１２３９０）、Ｅ（Ｘ１Ｘ２３）＝８．０５１１７（８．０３２２１）である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。
【０１４０】
また、準乱数（ソボル列）では表７２（ソボル列）にあるように、それぞれ、ξ１（歪度０．７４９６１、尖度４．９９６９７）、ξ２（歪度０．９９９８９、尖度５．４９７１６）、ν（歪度０．４９９８４、尖度３．９９６６１）で、それから生成される多変量非正規分布のモーメント特性の計測値は、
Ｅ（Ｘ１３）＝０．４６５１３（０．４５７０８）、Ｅ（Ｘ２３）＝０．３４７２４（０．３４６６０）、Ｅ（Ｘ１２Ｘ２）＝０．１５１９５（０．１４４９９）、Ｅ（Ｘ１Ｘ２２）＝０．１３１４７（０．１２２８９）、
Ｅ（Ｘ１４）＝２０．９６１４０（２０．７３８６５）、Ｅ（Ｘ２４）＝１８．９６４９３（１９．００６７８）、Ｅ（Ｘ１３Ｘ２）＝８．６２４４４（８．４６５１５）、Ｅ（Ｘ１２Ｘ２２）＝７．２３７７９（７．１２３９０）、
Ｅ（Ｘ１Ｘ２３）＝８．１５８６６（８．０３２２１）である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。
【０１４１】
ピアソンＶＩ型の場合、２次元においてξ１、ξ２、νすべてに、ピアソンＶＩ型を用いた例を表５７（合同法）、表６６（メルセンヌツイスター）、表７５（準乱数  ソボル列）に示す。まず、目標とする分布の特性はξ１（歪度１．０、尖度４．８）、ξ２（歪度１．２５、尖度５．５）、ν（歪度０．７５、尖度４．０）で、それから生成される多変量非正規分布のモーメント特性の理論値は、Ｅ（ｘ１３）＝０．８５７８３、Ｅ（ｘ２３）＝０．６９２１１、
Ｅ（ｘ１２ｘ２）＝０．２７４８８、Ｅ（ｘ１ｘ２２）＝０．２４１７４、
Ｅ（ｘ１４）＝２０．７３５０６、Ｅ（ｘ２４）＝１８．３１０４５、
Ｅ（ｘ１３ｘ２）＝８．４５１７６、Ｅ（ｘ１２ｘ２２）＝７．０７３９１、
Ｅ（ｘ１ｘ２３）＝７．８４５６３である。
【０１４２】
例えば、１０，０００，０００個のシミュレーションの場合、表５７（合同法）では、それぞれ、ξ１（歪度０．９９１２６、尖度４．６４４９５）、ξ２（歪度１．２２６１５、尖度５．２２３８６）、ν（歪度０．７５３０４、尖度３．９２１２３）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ１３）＝０．８４３０２（０．８５７８３）、Ｅ（Ｘ２３）＝０．６７３４２（０．６９２１１）、Ｅ（Ｘ１２Ｘ２）＝０．２６８６５（０．２７４８８）、Ｅ（Ｘ１Ｘ２２）＝０．２４０２５（０．２４１７４）、Ｅ（Ｘ１４）＝１９．２００２１（２０．７３５０６）、Ｅ（Ｘ２４）＝１７．１４７０５（１８．３１０４５）、Ｅ（Ｘ１３Ｘ２）＝７．９６７９４（８．４５１７６）、Ｅ（Ｘ１２Ｘ２２）＝６．７６２１５（７．０７３９１）、Ｅ（Ｘ１Ｘ２３）＝７．４１６７９（７．８４５６３）である。
【０１４３】
括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ１、ξ２、νの実際のシミュレーション時の歪度、尖度（例えば、上の例ではξ１（歪度０．９９１２６、尖度４．６４４９５）、ξ２（歪度１．２２６１５、尖度５．２２３８６）、ν（歪度０．７５３０４、尖度３．９２１２３））から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との間に開きがある。
【０１４４】
ところが、これをメルセンヌツイスターでシミュレーションすると、表６６（メルセンヌツイスター）にあるように、それぞれ、ξ１（歪度０．９９９６８、尖度４．８００６２）、ξ２（歪度１．２５０４２、尖度５．４９７０３）、ν（歪度０．７５０２２、尖度３．９９８７７）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ１３）＝０．８６１８５（０．８５７８３）、Ｅ（Ｘ２３）＝０．６７７４６（０．６９２１１）、
Ｅ（Ｘ１２Ｘ２）＝０．２７７１３（０．２７４８８）、
Ｅ（Ｘ１Ｘ２２）＝０．２３８６８（０．２４１７４）、Ｅ（Ｘ１４）＝２０．８７６６８（２０．７３５０６）、Ｅ（Ｘ２４）＝１８．０１４６４（１８．３１０４５）、Ｅ（Ｘ１３Ｘ２）＝８．５９７６０（８．４５１７６）、
Ｅ（Ｘ１２Ｘ２２）＝７．１０８３７（７．０７３９１）、Ｅ（Ｘ１Ｘ２３）＝７．７５２９８（７．８４５６３）である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。特に４次のモーメントがより近づいている。
【０１４５】
また、準乱数（ソボル列）では表７５（ソボル列）にあるように、それぞれ、それぞれ、ξ１（歪度１．０００２５、尖度４．８０８３２）、ξ２（歪度１．２４９９０、尖度５．４９８６２）、ν（歪度０．７５０５２、尖度４．００５３３）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ１３）＝０．８５８３４（０．８５７８３）、
Ｅ（Ｘ２３）＝０．６７５８６（０．６９２１１）、Ｅ（Ｘ１２Ｘ２）＝０．２６４８６（０．２７４８８）、Ｅ（Ｘ１Ｘ２２）＝０．２２９５７（０．２４１７４）、Ｅ（Ｘ１４）＝２０．７５１１１（２０．７３５０６）、
Ｅ（Ｘ２４）＝１８．３５８９８（１８．３１０４５）、Ｅ（Ｘ１３Ｘ２）＝８．３４４００（８．４５１７６）、Ｅ（Ｘ１２Ｘ２２）＝６．９４１４３（７．０７３９１）、Ｅ（Ｘ１Ｘ２３）＝７．７０２９１（７．８４５６３）である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。特に４次のモーメントがより近づいている。
【０１４６】
ピアソンＩ型の場合、２次元においてξ１、ξ２、νすべてに、ピアソンＩ型を用いた例を表６０（合同法）、表６９（メルセンヌツイスター）、表７８（準乱数  ソボル列）に示す。まず、目標とする分布の特性はξ１（歪度１．０、尖度４．０）、ξ２（歪度１．２５、尖度５．０）、ν（歪度０．７５、尖度３．０）で、それから生成される多変量非正規分布のモーメント特性の理論値は、
Ｅ（ｘ１３）＝０．８５７８３、Ｅ（ｘ２３）＝０．６９２１１、
Ｅ（ｘ１２ｘ２）＝０．２７４８８、Ｅ（ｘ１ｘ２２）＝０．２４１７４、
Ｅ（ｘ１４）＝１４．２３４９０、Ｅ（ｘ２４）＝１１．６３７１０、
Ｅ（ｘ１３ｘ２）＝５．９４８８１、Ｅ（ｘ１２ｘ２２）＝５．０６１７３、
Ｅ（ｘ１ｘ２３）＝５．２９９３９である。
【０１４７】
例えば、１０，０００，０００個のシミュレーションの場合、表６０（合同法）では、それぞれ、ξ１（歪度１．０２８０２、尖度４．０４７６０）、
ξ２（歪度１．２９３２３、尖度５．１５２６４）、ν（歪度０．７６２０８、尖度２．９９８７１）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ１３）＝０．９０４９１（０．８５７８３）、
Ｅ（Ｘ２３）＝０．７２６３４（０．６９２１１）、Ｅ（Ｘ１２Ｘ２）＝０．２９２４４（０．２７４８８）、Ｅ（Ｘ１Ｘ２２）＝０．２５６４２（０．２４１７４）、Ｅ（Ｘ１４）＝１４．５７５８７（１４．２３４９０）、
Ｅ（Ｘ２４）＝１１．８３１０５（１１．６３７１０）、Ｅ（Ｘ１３Ｘ２）＝６．０７５６４（５．９４８８１）、Ｅ（Ｘ１２Ｘ２２）＝５．１２７２８（５．０６１７３）、Ｅ（Ｘ１Ｘ２３）＝５．３７７８４（５．２９９３９）である。
【０１４８】
括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ１、ξ２、νの実際のシミュレーション時の歪度、尖度（例えば、上の例ではξ１（歪度１．０２８０２、尖度４．０４７６０）、ξ２（歪度１．２９３２３、尖度５．１５２６４）、ν（歪度０．７６２０８、尖度２．９９８７１））から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との間に開きがある。
【０１４９】
ところが、これをメルセンヌツイスターでシミュレーションすると、表６９（メルセンヌツイスター）にあるように、ξ１（歪度１．０００３９、尖度４．００４０２）、ξ２（歪度１．２４９１４、尖度４．９９５６７）、ν（歪度０．７５００９、尖度３．０００６１）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ１３）＝０．８６２１５（０．８５７８３）、Ｅ（Ｘ２３）＝０．６９７９１（０．６９２１１）、
Ｅ（Ｘ１２Ｘ２）＝０．２７６９６（０．２７４８８）、Ｅ（Ｘ１Ｘ２２）＝０．２４３８１（０．２４１７４）、Ｅ（Ｘ１４）＝１４．２８１２８（１４．２３４９０）、Ｅ（Ｘ２４）＝１１．６６９１０（１１．６３７１０）、
Ｅ（Ｘ１３Ｘ２）＝５．９３１９９（５．９４８８１）、Ｅ（Ｘ１２Ｘ２２）＝５．０５６０２（５．０６１７３）、Ｅ（Ｘ１Ｘ２３）＝５．２９９１６（５．２９９３９）である。この結果からすべての目標のモーメント特性の理論値と非常に近くなっていることがわかる。
【０１５０】
また、準乱数（ソボル列）では表７８（ソボル列）にあるように、ξ１（歪度０．９９９７３、尖度３．９９８７８）、ξ２（歪度１．２５０８２、尖度５．００６４０）、ν（歪度０．７４９５４、尖度２．９９７６６）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ１３）＝０．８５１５１  （０．８５７８３）、Ｅ（Ｘ２３）＝０．６９４２５（０．６９２１１）、Ｅ（Ｘ１２Ｘ２）＝０．２７２７１（０．２７４８８）、
Ｅ（Ｘ１Ｘ２２）＝０．２４３５８（０．２４１７４）、Ｅ（Ｘ１４）＝１４．１２４５４（１４．２３４９０）、Ｅ（Ｘ２４）＝１１．６３９４０（１１．６３７１０）、Ｅ（Ｘ１３Ｘ２）＝５．８７２８１（５．９４８８１）、
Ｅ（Ｘ１２Ｘ２２）＝５．０２３８１（５．０６１７３）、Ｅ（Ｘ１Ｘ２３）＝５．２８３３５（５．２９９３９）である。この結果からすべての目標のモーメント特性の理論値と非常に近くなっていることがわかる。
【０１５１】
【表５４】
0003708928.tif 000125
【表５５】
0003708928.tif 000126
【表５６】
0003708928.tif 000127
【表５７】
0003708928.tif 000128
【表５８】
0003708928.tif 000129
【表５９】
0003708928.tif 000130
【表６０】
0003708928.tif 000131
【表６１】
0003708928.tif 000132
【表６２】
0003708928.tif 000133
【表６３】
0003708928.tif 000134
【表６４】
0003708928.tif 000135
【表６５】
0003708928.tif 000136
【表６６】
0003708928.tif 000137
【表６７】
0003708928.tif 000138
【表６８】
0003708928.tif 000139
【表６９】
0003708928.tif 000140
【表７０】
0003708928.tif 000141
【表７１】
0003708928.tif 000142
【表７２】
0003708928.tif 000143
【表７３】
0003708928.tif 000144
【表７４】
0003708928.tif 000145
【表７５】
0003708928.tif 000146
【表７６】
0003708928.tif 000147
【表７７】
0003708928.tif 000148
【表７８】
0003708928.tif 000149
【表７９】
0003708928.tif 000150
【表８０】
0003708928.tif 000151
２．３  応用のためのあてはめプロシージャ
統計モデリングの見地から、経験分布の１つに３次と４次のモーメント行列をあてはめるプロシージャが必要である。各要素と標本のモーメントの式の間の偏りの２乗の和を最小にすることによって、Ｔ，ζ，γ，κ及びβを計算するプロシージャを提案する。すなわち、ＴＴ´＝Σという条件の下、次の式（９３）によって定義される関数の値を最小化することによって、適切なＴ，ζ，γ，κ及びβを得ることができる。
【０１５２】
ここで、ｍ＝ｐのときＴ＝Σ１／２（行列の平方根）であり、固定されている。行列の平方根の計算は、実対称行列の固有値問題に帰着され、おもな方法として、ヤコビ法、ハウスホルダー法、二分法、逆反復法、二分法・同時逆反復法、キューエル（ＱＬ）法、ディバイド・アンド・コンクウェア（Divide and
Conquer）法がある。
【０１５３】
【数１１６】
0003708928.tif 000152
ただし、ｆiは、３次と４次のモーメントの各要素の式であり、例えば、ｐ＝２について、ｆ1＝Ｅ（ｘ13）＝γｔ113ζ1＋γｔ123ζ2，ｋ＝９であり、ｍiは、標本分布のモーメントである。ｐ＝ｍ＝２の場合、Ｔは固定されている（ユアン及びベントラー（２０００））。例えば、標本モーメントｍiは、表８１に従う。式（９３）の値を最小にすることによって、ζ1＝０．６５９６７，ζ2＝１．０３６２４，γ＝０．６３５５３，κ1＝５．２７４７５，κ2＝６．６９１７６，β＝４．０２６８０というパラメータを得る。これらのパラメータによって得られた行列要素は、表８１に示されている。ここで、ｗiは、重みである。
【０１５４】
【表８１】
0003708928.tif 000153
３  最尤法
前述のモーメント法にかわる方法として最尤法がある。モーメント法では、３次と４次モーメントのような高次モーメントについては不安定になることがあるが、最尤法ではこの点が改善される可能性がある。
【０１５５】
最尤法では、生成した乱数によって密度関数を数値的に近似する。例として、特定の多変量非正規分布を１０，０００，０００点の乱数により非常に正確にシミュレーションする。最尤法のプロシージャは、以下の通りである。
【０１５６】
（１）最尤法に用いるパラメータの初期値をモーメント法から採用する。
【０１５７】
（２）多変量非正規分布から例えば１０，０００，０００点の乱数を得る。
【０１５８】
（３）各次元をステップΔｈによって分割に分ける。超立方体となる各分割内にある点の個数を数え、全ての点の個数で各分割内にある点の個数を除する。
【０１５９】
（４）次のように対数尤度を求める。
【０１６０】
｛Ｘ1，・・・，ＸN｝は、ｎ次元データの集合（Ｘｉは、ｎ次元データ・ベクトルである。）を示す。特定の多変量非正規分布のパラメータθを有するＸiの分割の確率Ｐi（θ）は、次の式（９４）によって定義される。ここで、パラメータθは、ユアンとベントラーの方法Ｉで用いられるパラメータＴ，ζ＝（ζ1，・・・，ζｍ），γ，κ＝（κ1，・・・，κｍ）、βの少なくとも１つを表すものとする。
【０１６１】
【数１１７】
0003708928.tif 000154
したがって、パラメータθを有する近似した尤度ｆi（θ）は、次の式
（９５）によって定義される。
【０１６２】
【数１１８】
0003708928.tif 000155
対数尤度の和は、次の式（９６）によって計算される。
【０１６３】
【数１１９】
0003708928.tif 000156
そして、この対数尤度の和を最大にすることにより、最尤推定量
【数１４６】
0003708928.tif 000157
を得る。このようにして、最尤推定パラメータが得られた。なお、尤度の積を最大にしても同様である。
【０１６４】
多変量非正規分布推定量
【数１４７】
0003708928.tif 000158
を、対数尤度の和を最大化することにより推定する。この方法を予想損失額に適用するため、１年の取引日に対応する２５０点を採用する。２５０点のモーメントを表８２に示し、最尤法による推定分布のモーメントを表８３〜８５に示す。
【０１６５】
２５０点は、合同法を用いて得られた、目的とする分布の特性が、ξ１（歪度０．７５、尖度５．０）、ξ２（歪度１．００、尖度５．５）、ν（歪度０．５、尖度４．０）を持つ（先のピアソンＩＶ型の多変量シミュレーションで使われた目標値）乱数である。２５０点からなる分布の特性の計測値は表８２に示されており、目的とする分布の特性とはかなり違ったものになっている。
【０１６６】
次に、多数の乱数からなる（ここでは、１０，０００，０００点）分布による最尤法の推定結果を示す。なお、分散共分散のみを経験分散共分散０．４１８１９と同じに設定するため、ｔ11＝０．２１４８，ｔ12＝０．９７６７に設定する。まず、表８３は合同法による推定結果であるが、２５０点の計測値よりも、目標とされた分布の特性値（表５４）に近づいていることがわかる。このことから、少数のデータの分布の推定には、モーメント法より、最尤法がすぐれていることがわかる。なお、最大対数尤度は、−５６２．７３になった。特に、金融のリスク管理などでは、一年のデータは２５０点程度であり，これから、多変量非正規分布を用い、良好な推定をするためには、最尤法を用いらなければならないといえよう。
【０１６７】
さらに、表８４（メルセンヌツイスター）、表８５（準乱数  ソボル列）を示す。前の多変量非正規分布に従う乱数のシミュレーション結果から、これらの方法は合同法に比べ、非常に目標とする分布に近い乱数を発生することがわかっており、これは最尤法を実行するための重要な性質である。結果として、両者の方法も、より目標値に近い値を推定している。なお、メルセンヌツイスターについての最大対数尤度は−５６１．３６、ソボル列についての最大対数尤度は−５６０．３６となった。
【０１６８】
本発明では、最尤法について主に連続型確率変数の場合の密度関数を乱数から直接計算している方法を中心にしているが、離散型確率変数の場合は密度関数の代わりに確率関数となる。また、密度関数や確率関数を積分した分布関数として表現する方法も可能である。
【０１６９】
本実施の形態は、一般性を失わないように標準化された平均０、分散１の設定（Ｅ［ｘｉ］＝０，Ｖａｒ［ｘｉ］＝σｉｉ＝１（Ｅ［ｘｉ２］＝σｉｉ＝１），１≦ｉ≦ｎ）で行われているが、一般の場合でも標準化の逆変換を行えば容易に拡張でき、また、最尤法を使えば、さらに、その標準化のための平均と標準偏差もパラメータとして、直接、生のデータ（標準化などの前処理をしていないという意味）から分布を推定することが容易に可能である。
【０１７０】
【表８２】
0003708928.tif 000159
【表８３】
0003708928.tif 000160
【表８４】
0003708928.tif 000161
【表８５】
0003708928.tif 000162
なお、さらに、分布の１次から４次までの一部または全部のモーメントが解析的に容易にもとまらない場合や、発散などで存在しない場合、あるいは密度関数でなく特性関数で定義される分布などの場合でも、その分布の乱数を発生させてそれからモーメントを計算する方法（必ず計算可能）や、最尤法を行うにあたって分布を特徴付けるパラメータやその他のパラメータ（両パラメータを含めてθと表す。）を適当な間隔で適当に与え、その分布の乱数を用いて幅広く尤度を計算する方法を使えば、そのような分布もこの発明のνやξの分布として採用できる。そのような分布として例えば、一般の場合、２次のモーメントが発散し、密度関数が陽表的に書けないが特性関数によって表現できる安定分布などがある。ゆえに、本発明の最尤法は、νやξの分布としては、乱数が発生できれば１次から４次までの一部又は全部のモーメントが求まらない分布でも採用可能となり、ユアン・ベントラーの方法からさらに拡張された方法となっている。
【０１７１】
４  実施例
図１は、２次元を例とした予想損失額計算を示すフローチャートである。ここでは、日本株式リターンと外国株式リターンを２次元の変数に取る。この予想損失額計算の一連の手順は、コンピュータにおいて所定のコンピュータ・プログラムを実行することによって実現される。
【０１７２】
ステップＳ１においては、日本株式リターンＸ1´と外国株式リターンＸ2´を標準化する。
【０１７３】
標準化は、次のように経験分布の平均を０、分散を１とすることにより行う。日本株式リターンと外国株式リターンのある実現値を持つ確率変数Ｘ1´及びＸ2´について、それぞれの平均をＥ（Ｘ1´）及びＥ（Ｘ2´）、分散をＶ（Ｘ1´）及びＶ（Ｘ2´）とする。このとき、標準化日本株式リターンＸ1と標準化外国株式リターンＸ2を、次の式（９７）と（９８）により与える。
【０１７４】
【数１２０】
0003708928.tif 000163
【数１２１】
0003708928.tif 000164
ステップＳ２においては、標準化日本株式リターンＸ1と標準化外国株式リターンＸ2から分散共分散行列Σを計算する。
【０１７５】
ステップＳ３においては、標準化日本株式リターンＸ1と標準化外国株式リターンＸ2について３次モーメントを計算する。３次モーメントは、次の式（９９）に示すように、ｍ1〜ｍ4までの４種類がある。具体的には、Ｅ（・）は、Ｘ１，Ｘ２のデータ（実現値）から計算される。
【０１７６】
【数１２２】
0003708928.tif 000165
ステップＳ４においては、標準化日本株式リターンＸ1と標準化外国株式リターンＸ2について４次モーメントを計算する。４次モーメントは、次の式（１００）に示すように、ｍ5〜ｍ9までの５種類がある。
【０１７７】
【数１２３】
0003708928.tif 000166
ステップＳ５においては、ステップＳ２において計算した分散共分散行列Σについて、Σ＝ＴＴ´を満たす行列Ｔ＝（ｔij）を計算する。いまξの次元とＸｉの次元がともに２次元なので、行列Ｔは、Ｔ＝Σ１／２（行列の平方根）と固定される。
【０１７８】
具体的には、次の式（１０１）の値を最小にするパラメータζ＝（ζ1，ζ2），γ，κ＝（κ1，κ2），βを決定する。
【０１７９】
【数１２４】
0003708928.tif 000167
ここで、重みｗiによって、モーメントに軽重をつけることができる。本実施例では、簡単のためにｗi＝１とする。
【０１８０】
式（１０１）においては、３次及び４次モーメントの表現として、次の式（１０２）と（１０３）の表現を用いる。行列Ｔは、Ｔ＝Σ１／２（行列の平方根）と計算されているので、式（１０１）に最小値を与えるパラメータζ＝（ζ1，ζ2），γ，κ＝（κ1，κ2），βを求める。
【０１８１】
［３次モーメント］
【数１２５】
0003708928.tif 000168
ただし、ｆ1，ｆ2，ｆ3，ｆ4は、Ｘ＝νＴξで発生される確率ベクトルＸ＝（ｘ１，ｘ２）´のそれぞれ３次モーメントＥ（ｘ13），Ｅ（ｘ23），Ｅ（ｘ12ｘ2），Ｅ（ｘ1ｘ22）の表現である。
【０１８２】
［４次モーメント］
【数１２６】
0003708928.tif 000169
ただし、ｆ5，ｆ6，ｆ7，ｆ8，ｆ9は、Ｘ＝νＴξで発生される確率ベクトルＸ＝（ｘ１，ｘ２）´のそれぞれ４次モーメントＥ（ｘ14），Ｅ（ｘ24），Ｅ（ｘ13ｘ2），Ｅ（ｘ12ｘ22），Ｅ（ｘ1ｘ23）の表現である。ここで、σ１１＝σ２２＝１であり、また、ｔ１１２＝ｔ２２２＝１−ｔ１２２＝１−ｔ２１２である。
【０１８３】
ステップＳ６においては、確率変数ξ1, ξ2，νがピアソン分布のどの型に属するか決定する。
【０１８４】
確率変数ξ1は、Ｅ（ξ1）＝０，Ｅ（ξ12）＝１，Ｅ（ξ13）＝ζ1，
Ｅ（ξ14）＝κ1を満たす。同様に、確率変数ξ2は、Ｅ（ξ2）＝０，
Ｅ（ξ22）＝１，Ｅ（ξ23）＝ζ2，Ｅ（ξ24）＝κ2を満たす。確率変数νは、Ｅ（ν）＝０，Ｅ（ν２）＝１，Ｅ（ν3）＝γ，Ｅ（ν4）＝βを満たす。規格化された確率変数ξ1, ξ2，νの歪度はそれぞれζ1，ζ2，γであり、尖度はそれぞれκ1，κ2，βである。
【０１８５】
ここで、次の式（１０４）と表８６を用い、確率変数ξ1，ξ2，νがピアソン分布のどの型になるか決定する。式（１０４）におけるβ1とβ2は、それぞれ歪度の２乗と尖度である。
【０１８６】
【数１２７】
0003708928.tif 000170
【表８６】
0003708928.tif 000171
ステップＳ７においては、確率変数ξ1，ξ2，νのそれぞれについて、ステップＳ６で求めた型のピアソン分布に従う乱数を発生する。
【０１８７】
具体的に、ピアソンＩＶ型以外の場合には、次の表８７と表８８に示すプロシージャによって乱数を発生する。
【０１８８】
【表８７】
0003708928.tif 000172
【表８８】
0003708928.tif 000173
ピアソンＩＶ型の場合は、次のアルゴリズムによって乱数を発生する。
【０１８９】
【数１２８】
0003708928.tif 000174
0003708928.tif 000175
このアルゴリズムは、ピアソンＩＶ型分布の正規化定数を解析解又はこの展開（式（１０５））に基づいて計算するステップと、棄却法により乱数を発生するステップを含む。
【０１９０】
ステップＳ８においては、ステップＳ７において発生した乱数ξ1，ξ2，νを次の式（１０６）によって、標準化日本株式リターン従う乱数ｘ1と標準化外国株式リターンに従う乱数ｘ2に変換する。
【０１９１】
【数１２９】
0003708928.tif 000176
ステップＳ９においては、ステップＳ８で得られた乱数ｘ1，ｘ2を次の式（１０７）と（１０８）によって、標準化前の日本株式リターンに従う乱数ｘ1´と標準化前の外国株式リターンに従う乱数ｘ2´に変換する。ここで、分散Ｖ（・）と平均Ｅ（・）には、ステップＳ１で標準化した際の値を用いる。
【０１９２】
【数１３０】
0003708928.tif 000177
【数１３１】
0003708928.tif 000178
ステップＳ１０においては、ステップＳ９で得られた乱数ｘ1´と乱数ｘ2´を用いて、予想損失額（ＶＡＲ）を計算する。
【０１９３】
日本株式と外国株式の資産の残高をＡ1，Ａ2とすると、これらの重みを用いて日本株式と外国株式を併せたリターンに従う乱数を得る。次の式（１０９）と（１１０）にそれぞれ日本株式と外国株式の資産の重みを示し、日本株式と外国株式を併せたリターンに従う乱数Ｗを式（１１１）に示す。
【０１９４】
【数１３２】
0003708928.tif 000179
【数１３３】
0003708928.tif 000180
【数１３４】
0003708928.tif 000181
乱数Ｗを用いたモンテカルロ法によるシミュレーションによって、日本株式と外国株式を併せたリターン分布をマイナス無限大から積分した領域が、全体の領域の１％になったところのリターンを用い、このリターン（通常、マイナスである。）の絶対値に日本株式と外国株式を併せた資産の残高を乗じたものが標準的な予想損失額である。
【０１９５】
具体的には、日本株式と外国株式を併せたリターン分布に従う乱数Ｗを複数発生させ、低いものから並べて下から１％目（通常、銀行規制などで使われている基準）のところのリターンを用いる。なお、予想損失額に用いるリターンは１％に限らず、所望の値を用いることができる。
【０１９６】
本実施例においては、日本株式と外国株式の２種類の資産を有する２次元の場合について説明したが、ｎ種類の資産を有する一般のｎ次元の場合には、容易に拡張することができる。以上の方法は、損失額自体の分布によるリスク管理にも適用できる。
【０１９７】
この場合、ステップＳ５においては、式（１１２）において、
ｆijk（Ｔ，ζ，γ，κ，β）とｆijkl（Ｔ，ζ，γ，κ，β）には、３次と４次モーメントの表現として式（１１３）と（１１４）の関係を用いる。また、
ｍijkとｍijklは、式（１１５）と（１１６）によって定義される。具体的には、Ｅ（・）は、Ｘｉのデータ（実現値）から計算される。
【０１９８】
【数１３５】
0003708928.tif 000182
【数１３６】
0003708928.tif 000183
【数１３７】
0003708928.tif 000184
【数１３８】
0003708928.tif 000185
【数１３９】
0003708928.tif 000186
そして、式（１１２）を最小にするようなパラメータζ＝（ζ1，・・・，ζｍ），γ，κ＝（κ1，・・・，κｍ），βを求める。
【０１９９】
ステップＳ６においては、これらのパラメータから、確率変数ν，ξ1〜ξｍがどの型のピアソン分布か決定し、ステップＳ７においてその型の乱数を発生する。ステップＳ８においては、ステップＳ７において発生した乱数を乱数ｘ1，ｘ2・・・，ｘnに変換し、さらにステップＳ９においては、乱数ｘ1´，ｘ2´・・・，ｘn´に変換する。
【０２００】
ステップＳ１０においては、ｎ種類の資産の残高をＡ1，Ａ2，・・・，Ａnとすると、これらの重みを乗じたｎ種類の資産を併せたリターンに従う乱数を得る。次の式（１１７）にｎ種類の資産の重みを示し、式（１１８）にｎ種類の資産を併せたリターン分布に従う乱数Ｗを示す。
【０２０１】
【数１４０】
0003708928.tif 000187
【数１４１】
0003708928.tif 000188
このような乱数Ｗを用いて、ｎ種類の資産を有する場合にも２次元の場合と同様にモンテカルロ法によるシミュレーションによって予想損失額を計算することができる。すなわち、複数の乱数Ｗを発生させ、低いものから並べて下から１％目のところのリターンを取り出し、このリターンの絶対値にｎ種類の資産の残高を乗じたものを予想損失額とする。
【０２０２】
さらに、本実施例においては、モーメント法に代わって最尤法によって推定したパラメータを用いて予想損失額を計算することができる。前述したように、株式リターンを用いる場合、年間取引日は２５０日と少数である。最尤法によると、例えば２５０点のような少数の標本点であっても、パラメータを正確に算出することができる。
【０２０３】
図２は、最尤法を用いて予想損失額を計算する一連の手順を示すフローチャートである。この一連の手順においては、ステップＳ５とＳ６の間にステップＳ１０１が加わった他は、図１に示した一連の手順と同様であるので、ステップＳ１０１を除いて説明を省略するものとする。
【０２０４】
ステップＳ１０１において、ステップＳ５においてモーメント法により計算したパラメータを初期値として、最尤法によりパラメータζ1，ζ2，γ，κ1，κ2，βの推定値を計算する。
【０２０５】
具体的に、ステップＳ１０１においては、次の手順によりパラメータの推定値を計算する。なお、以下このステップＳ１０１においては、パラメータＴ，ζ1，ζ2，γ，κ1，κ2，βの少なくとも１つをθで表記することにする。
【０２０６】
（１）ステップＳ５において、モーメント法から計算したパラメータを初期値として採用する。
【０２０７】
（２）多変量非正規分布から例えば１０，０００，０００点の乱数を得る。乱数の発生は、ステップＳ６とＳ７と同様にして行う。
【０２０８】
（３）各次元をステップΔｈによって分割に分ける。超立方体となる各分割内にある点の個数を数え、全ての点の個数で各分割内にある点の個数を除する。
【０２０９】
（４）次のように対数尤度を求める。
【０２１０】
｛Ｘ1，・・・，ＸN｝は、ｎ次元データの集合（Ｘｉは、ｎ次元データ・ベクトルである。）を示す。特定の多変量非正規分布のパラメータθを有するＸiの分割の確率Ｐi（θ）は、次の式（１１９）によって定義される。ここで、パラメータθは、ユアンとベントラーの方法Ｉで用いられるパラメータＴ，ζ＝（ζ1，・・・，ζｍ），γ，κ＝（κ1，・・・，κｍ）、βの少なくとも１つを表すものとする。
【０２１１】
【数１４２】
0003708928.tif 000189
したがって、パラメータθを有する近似した尤度ｆi（θ）は、次の式（１２０）によって定義される。
【０２１２】
【数１４３】
0003708928.tif 000190
対数尤度の和は、次の式（１２１）によって計算される。
【０２１３】
【数１４４】
0003708928.tif 000191
そして、この対数尤度の和を最大にすることにより、最尤推定量
【数１４８】
0003708928.tif 000192
を得る。この最尤推定量をパラメータとして用いる。なお、尤度の積を最大にしても同様である。
【０２１４】
本実施の形態によると、従来のリスクをより正確に測定出来るだけでなく、不動産やオルタネティブ投資などの歪んだリターン分布を含んだポートフォリオなどのリスク管理も可能になる。また、従来、リスク管理で使われていた、リターン（収益率）ベースの方法だけでなく、損失額自体の分布を直接評価できるようになり、現物だけでなく、オプションなど様々なデリバティブを含むポートフォリオのリスク管理や、保険・天候・地震などのその他の様々なリスク（金融リスク、オペレーショナル・リスク、保険リスク、災害リスク、事業リスクなど）も統合したリスク管理が可能となる。ここでは、リスク計測対象の損失額自体の分布や収益率の分布など、リスクを計測する際の評価のための分布を総称して損益分布と呼ぶ。
【０２１５】
本実施の形態においては、ステップＳ２、Ｓ３及びＳ４で求めた３次及び４次モーメントｍ１〜ｍ９を用いてステップＳ５でパラメータζ，γ，κ，βを求め、そのパラメータを初期値として、ステップＳ１０１において前記初期値の近傍に属する最大対数尤度を与えるパラメータを求めた。しかしながら、パラメータの初期値を定めるステップＳ２〜Ｓ５の一部または全部を省略することもできる。この場合、パラメータの初期値は未定として、ステップＳ１０１において、パラメータ空間における離散的な格子点について幅広い領域で対数尤度を計算し、最大対数尤度（又は最大尤度）を与えるパラメータを求める。このような方法は、格子探索法（グリッドサーチ）と称される。
【０２１６】
この方法によれば、１次から４次までのモーメントの一部又は全部が存在しない確率分布であっても、乱数さえ発生できればξやνの分布として採用可能となる。そのような場合、最尤法のための、分布を特徴付けるようなパラメータやその他あらゆる必要なパラメータを広くθと記す。
【０２１７】
５  変形例
５．１  金融商品のシミュレーション方法
次に、本実施の形態の第１の変形例として、本実施の形態の乱数を使った多変量非正規分布の特性を用いた金融商品のシミュレーション方法を説明する。
【０２１８】
この方法では、次のアルゴリズムによって発生されるピアソンＩＶ型分布に従う乱数を使ったモンテカルロ法によるシミュレーションを用いることができる。このアルゴリズムは、ピアソンＩＶ型分布の正規化定数の解析解又はこの展開（式（１２２））を計算するステップと、棄却法により乱数を発生するステップとを含む。
【０２１９】
【数１４５】
0003708928.tif 000193
0003708928.tif 000194
この方法は、多資産の高次の相関を考慮した乱数を発生させ、実際のマーケットの特性を正確に表すリターンを発生することができる。例えば、金融商品のパフォーマンス測定などに用いられる。
【０２２０】
この方法は、前記多変量非正規分布（ξやνに、ピアソン分布系を使うことに限らず一般的な確率分布を使うことを含む。）を用いることにより、金融の分野における分布に特徴的ないわゆるファット・テールを正確に再現することができる。また、様々なデリバティブの対象データを本実施の形態の方法により、正確な確率分布（ξやνに、ピアソン分布系を使うことに限らず一般的な確率分布を使うことを含む。）をあてはめ（単変量分布の場合も、多変量の特殊な場合として同様に扱える。）、その分布に従うモンテカルロ・シミュレーションをすることで価格が計算できる。デリバティブの設計には、多地点での気象データや様々な指数の相関が必要で、その分析にも本実施の形態が応用される。
【０２２１】
例えば、天候・保険デリバティブやリアル・オプションのプライシングや設計などに応用できる。特に、プライシングに多変量非正規分布を使う例としては、コリレーティング・トリガー・オプション（複数の気象要素を選んで、それぞれに補償金の支払基準を設定するもので、オプション購入者は、契約終了時にその中で一番有利なものを選択・行使できるオプション。）のような複合オプションや複数地点オプションなどがある。
【０２２２】
５．２  半導体へのイオン注入
次に、本実施の形態の第２の変形例として、半導体へのイオン注入への応用を説明する。イオン注入への応用では、半導体に注入されるイオンの分布をモンテカルロ法によってシミュレーションする。
【０２２３】
イオンの空間分布や電荷特性は、一般に多変量非正規分布となるが、ピアソン分布を含む前記多変量非正規分布に従う乱数を用いた方法により正確に再現することができる。
【０２２４】
このように本実施の形態の多変量非正規分布を適用することにより、従来の解析的なモデルを、本発明の多変量非正規分布の乱数表現によって実際の分布により近い正確なモデルに改良できる。
【０２２５】
また、現在、イオン注入後のイオン分布を知るための物理モデルのモンテカルロ・シミュレーションは、非常に長時間かかっているが、本実施の形態の多変量非正規分布の乱数とパラメータ推定法によって、より少数の物理シミュレーションの結果から短時間にその分布の全体を正確に推定可能になり、大幅に現在のモンテカルロ・シミュレーションの短時間化を行える。
【０２２６】
５．３  アセット・アロケーション
本実施の形態の方法によって、多変量非正規分布に従う資産間のアセット・アロケーションを、従来の平均・分散による最適化だけでなく歪度も考慮できる最適化手法と組み合わせて実行する。この歪度まで考慮できる最適化手法については、参考文献に今野浩「理財工学ＩＩ」、日科技連、１９９８年、コンノ，エイチ．，スズキ，ティー.及びコバヤシ，ディー．（Konno,H and Suzuki, T. and Kobayashi, D.）（１９９８）「平均リスク歪度ポートフォリオ問題を解くための分岐と境界アルゴリズム（A Branch and Bound Algorithm for Solving Mean-Risk-Skewness Portfolio Models）」、最適化方法及びソフトウェア（Optimization Methods and Software）、１０巻、２９７〜３１７頁。コンノ，エイチ．，シラカワ，エイチ．，及びヤマザキ，エイチ．，（Konno, H., Shirakawa, H., and Yamazaki, H.）（１９９３）「平均絶対偏差ポートフォリオ最適化モデル（A Mean-Absolute Deviation Skewness Portfolio Optimization Model）」、オペレーションズリサーチ年報（Annals of Operations Research）、４５巻、２０５〜２２０頁、コンノ，エイチ．及びスズキ，ケイ.（Konno,H and Suzuki, K.）（１９９５）「平均分散歪度ポートフォリオ最適化モデル（A Mean-Variance-Skewness Portfolio Optimization Model）」、日本オペレーションズリサーチ学会誌（Journal of the Operations Research Society of Japan）、３８巻、１７３〜１８７頁がある。さらに、尖度まで考慮した最適化を使えば、拡張も容易である。
【０２２７】
具体的には、本実施の形態の方法により、実際の資産のリターンや価格の分布を推定し、ポートフォリオとしての平均、分散、歪度、尖度を導出する。そして、平均、分散だけでなく、歪度まで考慮できる最適化手法によって、最適な資産配分を計算する。なお、歪度だけでなく尖度まで考慮できる最適化手法によって、最適な資産配分を計算することができる。
【０２２８】
この方法により、不動産やオルタネティブ投資などの従来の正規分布の近似では無理がある資産も含めて、正確で有効なアセット・アロケーションが可能となる。また、年金コンサルテーションや資産運用への応用としてパフォーマンスが改善するアセット・アロケーションや、多変量非正規分布の資産モデル構築も可能である。
【０２２９】
５．４  能力評価
例えば知能検査など能力の測定結果は、正規分布を形成しないことが知られている。例えば知能検査の場合には、１次元の場合にはピアソン分布が使われている。能力評価の統計的分析については、次の参考文献がある。バート，シー．（１９６３）「知能の分布は正規分布に従うか？（Is Intelligence Distributed Normally ?）」、英国統計心理学雑誌（The British Journal of Statistical Psychology）、１６（XVI）巻、第２部。
【０２３０】
異なった能力検査や試験の結果に、本発明の実施によって、正確な多変量非正規分布が推定され、能力評価を正確に実施することができる。このことから、入試判定や人事アセスメント（選抜基準の設定や選抜有効性の検証）などをより正確化することができる。また、パーセンタイルなど様々な基準により、総合的な能力値を定義することができる。
【０２３１】
さらに、大学入学試験の結果と入学者の卒業成績、就職試験の結果と入社後の評価など、一部の周辺分布が大量のデータであるとき、本実施の形態の方法により、残りの少数のデータの分布を最尤法により推定することもできる。
【０２３２】
産業上の利用可能性
本発明によると、ピアソン分布系や一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布の推定方法、及びこれらの手順を格納した媒体を提供することができる。
【０２３３】
具体的には、多変量非正規分布に従う乱数発生方法及び多変量非正規分布のパラメータ推定方法は、様々な分野における統計モデリングに有効である。具体的には、統計ソフトウェアのパッケージ、リスク・マネジメントにおける予想損失額の計算への応用、金融商品のシミュレーション、アセット・アロケーションへの応用、天候・保険デリバティブやリアル・オプションのプライシングと設計、半導体へのイオン注入のシミュレーション、教育や人事アセスメントにおける能力評価に有効である。
【０２３４】
【発明の効果】
本発明によると、ピアソン分布系及び一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンＩＶ型分布に従う乱数発生方法及びコンピュータ・プログラムを提供することができる。
【図面の簡単な説明】
【図１】図１は、予想損失額を計算する一連の手順を示すフローチャートである。
【図２】図２は、最尤法を用いて予想損失額を計算する一連の手順を示すフローチャートである。
【図１】
0003708928.tif 000215
【図２】
0003708928.tif 000216