JP4350062B2 - 非正規分布に従う乱数発生方法及びそのパラメータの推定方法 - Google Patents

Description

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

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

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

本明細書では、ユアン等が方法Ｉと称したものは、次の通りである。ξ₁，・・・，ξ_mは、独立な確率変数であり、Ｅ（・）を期待値として、Ｅ（ξ_j）＝０，Ｅ（ξ_j ²）＝１，Ｅ（ξ_j ³）＝ζ_j，Ｅ（ξ_j ⁴）＝κ_j及びξ＝（ξ₁，・・・，ξ₁）´である。νは、ξに独立な確率変数であり、Ｅ（ν）＝０，Ｅ（ν²）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。また、Ｔ＝（ｔ_ij）は、階数ｐのｐ行ｍ列の非確率行列であり、ＴＴ´＝Σとｍ≧ｐを満たす。すると、次の式（３５）で与えられる確率ベクトルは、一般に、非楕円分布に従い、Ｃｏｖ（Ｘ）＝Σを満たす。ここで、Ｃｏｖ（・）は、確率ベクトルＸ＝（ｘ₁，・・・，ｘ_n）´の分散共分散行列である。

ｘ_iの周辺の歪度（skew）と尖度（kurt）は、それぞれ次の式（３６）及び（３７）によって与えられる。

Ｘは、ξ〜Ｎ_m（０，Ｉ）及びＰｒ（ν≧０）＝１（ただし、Ｎ_mはｍ次元正規分布、Ｐｒは確率を表す。）が満たされるとき、楕円分布に従う（ファン、コッツ及びング（Fang, Kotz and Ng）（１９９０））。ν＝１とし、ｚ₁，・・・，ｚ_pを独立な標準正規変数とし、ξ＝｜ｚ₁｜−Ｅ（｜ｚ₁｜），ξ_j＝ｚ_j，ｊ＝２，・・・，ｐとすると、式（３５）におけるＸは、アザリニ及びバレ（Azzalini and Valle）（１９９６）によって定義された歪んだ正規分布に従う。

この方法には、母集団の４次モーメント行列を容易に計算できるという利点がある。Ｓを式（３５）のＸからの大きさｎの標本に基づいた標本共分散とする。ｖｅｃｈ（・）を対称行列から重複しない要素を取り出してベクトルに変換する演算子とし、ｓ＝ｖｅｃｈ（Ｓ）とする。すると、ｓの漸近的分散共分散行列はΓ／ｎによって与えられる。ここで、Γ＝ｖａｒ（ｖｅｃｈ（ＸＸ´））であり、次の式（３８）によって与えられる。ｖａｒは、分散である。

ただし、ｔ_jは、Ｔ＝（ｔ_ij）の第ｊ列ベクトルであり、

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

この後、ユアン等は、検定統計の理論を研究した（ユアン及びベントラー（１９９９ａ，１９９９ｂ，２０００））。ユアン等は、確率ベクトルを生成するにはフライシュマン（１９７８）の方法を推奨している。しかし、この分布は、あまり広い歪度及び尖度を有していない。これに対して、ピアソン分布系は、次のような特徴を有する。

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

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

（３）ＩＶ型を除いたピアソン分布系に従う乱数は、正規分布とガンマ分布に従う乱数を用いて生成される。

（４）最近、実現することが難しかったピアソンＩＶ型分布を含む殆ど全ての種類を用いる実用的な手法が開発された（永原（２００２））。

ここで、ピアソン分布系を導入する。ケイ・ピアソン（K. Pearson）（１８９５）は、確率密度関数ｐに関する次の式（３９）の微分方程式に従い、ピアソン分布系を定義した。

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

ここで、ピアソン分布系に関して参照される研究者と文献を掲げる（ピアソン（１９１４）、ピアソン及びハートリー（１９５４）、ピアソン（１９６３）、ホードリ−（１９６８）、エルダートン及びジョンソン（１９６９）、オード（１９７２）、パリシュ（１９８３）、ジョンソン、コッツ及びバラカリシュナン（１９９４）、スチュアート及びオード（１９９４））。

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

［１］アザリニ，エイ．及びバレ，エイ．ディー．（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）、４８巻、４６５〜４７１頁

ところで、検定統計のみならず統計モデリングの見地からも、多変量非正規分布を導出し、経験分布を近似する分布を構築することがさらに求められている。また、少数の標本からも多変量非正規分布を安定して推定することが求められている。

本発明は、前述の実情に鑑みて提案されるものであって、ピアソン分布系及び一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンＩＶ型分布に従う乱数発生方法及びコンピュータ・プログラムを提供することを目的とする。

前述の課題を解決するために、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータ上でユアン及びベントラーの方法Ｉに基づいて多変量非正規分布に従う乱数を発生する乱数発生方法において、コンピュータを用いてｎ次元経験分布にｎ次元多変量非正規分布をあてはめるステップと、コンピュータを用いて乱数を発生するステップとを有し、前記あてはめるステップは、前記経験分布の３次及び４次モーメントに関するあてはめについて、次の関係式（４０）及び（４１）を用いる。

ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄ_ｎは、ｎ次のデュプリケーション行列であり、Ｄ_ｎ ⁺は、Ｄ_ｎのムーア・ペンローズの一般化逆行列であり、

は、クロネッカー積であり、Ｅ_iiは、ｅ_iを第ｉ列単位としてｅ_iｅ_i ^'である。

ここで、前記ユアン及びベントラーの方法Ｉは、次の通りである。独立な確率変数ξ₁，・・・，ξ_mは、パラメータζ_j及びκ_jについて、Ｅ（ξ_j）＝０，Ｅ（ξ_j ²）＝１，Ｅ（ξ_j ³）＝ζ_j，Ｅ（ξ_j ⁴）＝κ_j（１≦ｊ≦ｍ）を満たす。

ξ_jに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν²）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。

階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔ_ij）は、行列Σ＝（σ_ij）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。

このとき、次の式（４２）で与えられる確率ベクトルＸ＝（ｘ₁，・・・，ｘ_n）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ₁，・・・，ξ_m）´である。

好ましくは、確率変数ξ₁，・・・，ξ_m及びνには、ピアソン分布系に従う乱数を使う。

好ましくは、確率変数ξ₁，・・・，ξ_m及びνには、少なくとも２以上の型のピアソン分布を用いる。

好ましくは、前記あてはめるステップは、式（４３）及び（４４）でそれぞれ与えられるｎ次元の経験分布（確率ベクトル（Ｘ_１，・・・，Ｘ_ｎ）´とする。）の３次及び４次モーメントに関して、式（４５）の値を最小にするパラメータＴ，ζ＝（ζ₁，・・・，ζ_ｍ），γ，κ＝（κ₁，・・・，κ_ｍ），βの少なくとも１つを求める。

ただし、ｆ_ijk（Ｔ，ζ，γ，κ，β）及びｆ_ijkl（Ｔ，ζ，γ，κ，β）は、３次モーメントＥ（ｘ_iｘ_jｘ_k）及び４次モーメントＥ（ｘ_iｘ_jｘ_kｘ_ｌ）にそれぞれ対応する関係式（４０）及び（４１）による表現とする。また、ｗ_ijk及びｗ_ijklは、所定の重みである。
ここで、記号

は、

（ただし、ｉ≦ｊ≦ｋ）の意味で、１≦ｉ≦ｎ，１≦ｊ≦ｎ，１≦ｋ≦ｎ，ｉ≦ｊ≦ｋの各条件を満たす全てのｉ，ｊ，ｋの組み合わせについて式に代入し、それらの総和を取ったものである。同様に、記号

は、

（ただし、ｉ≦ｊ≦ｋ≦ｌ）の意味で、１≦ｉ≦ｎ，１≦ｊ≦ｎ，１≦ｋ≦ｎ，１≦ｌ≦ｎ，ｉ≦ｊ≦ｋ≦ｌの各条件を満たす全てのｉ，ｊ，ｋ，ｌの組み合わせについて式に代入し、それらの総和を取ったものである。
また、前記ｍ_ｉｊｋ＝Ｅ（Ｘ_iＸ_jＸ_k）及びｍ_ｉｊｋｌ＝Ｅ（Ｘ_iＸ_jＸ_kＸ_l）の意味について説明する。データの分布（Ｘ_１，・・・，Ｘ_ｎ）´は、次の式（Ａ１）のようにＮ個のデータ・ベクトル

（ただし、

は各々のデータ（実現値）を示す。１≦ｑ≦N）から構成されている。

前記ｍ_ｉｊｋ及びｍ_ｉｊｋｌは、このようなデータ・ベクトルの要素となるデータを用いて次の式（Ａ２）及び（Ａ３）によって計算される。

好ましくは、ｗ_ijk＝１，ｗ_ijkl＝１である。

好ましくは、ｎ＝２であって、前記あてはめるステップは、式（４６）及び（４７）でそれぞれ与えられる２次元経験分布（確率ベクトル（Ｘ_１，Ｘ_２）´とする。）の３次及び４次モーメントと、式（４８）及び（４９）で与えられる３次モーメントＥ（ｘ_iｘ_jｘ_k）及び４次モーメントＥ（ｘ_iｘ_jｘ_kｘ_ｌ）にそれぞれ対応する表現とについて、式（５０）で与えられる値を最小にするようにパラメータＴ，ζ＝（ζ₁，ζ_２），γ，κ＝（κ₁，κ_２），βの少なくとも１つを求める。

好ましくは、ｗ_i＝１である。

好ましくは、パラメータＴ，ζ＝（ζ_１，・・・，ζ_ｍ），γ，κ＝（κ_１，・・・，κ_ｍ），βの少なくとも１つを最尤法により推定する。

好ましくは、式（５１）で与えられるκの値を用いて、表８に示すκの値と分布の型との対応関係の少なくとも１つを用いて前記確率変数の属する型を決定する。

ただし、β_１とβ_２は、それぞれ歪度の２乗と尖度であり、式（５１）及び表８におけるκは、前記パラメータκ＝（κ₁，・・・，κ_ｍ）とは異なる。

好ましくは、前記決定した型に応じて、次の表９及び表１０に示す生成法の少なくとも１つを用いて乱数Ｚを発生する。

好ましくは、コンピュータを用いて、式（５２）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンＩＶ型分布に従う乱数を発生させる。

好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、式（５３）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンＩＶ型分布に従う乱数を発生させる。

好ましくは、前記乱数を発生するステップは、コンピュータを用いて発生した、加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数に基づく。

好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、ｎ個のデータ

からなるデータ・ベクトルＸ_ｉ´について、経験分布のデータ｛Ｘ_１´，・・・，Ｘ_Ｎ´｝を取得するステップと、コンピュータを用いて、前記データ｛Ｘ_１´，・・・，Ｘ_Ｎ´｝を標準化して｛Ｘ_１，・・・，Ｘ_Ｎ｝とするステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ_１，・・・，Ｘ_Ｎ｝に基づいて、分散共分散行列Σを計算するステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ_１，・・・，Ｘ_Ｎ｝に基づいて、３次モーメントｍ_ｉｊｋ（１≦ｉ≦ｊ≦ｋ≦ｎ）を計算するステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ_１，・・・，Ｘ_Ｎ｝に基づいて、４次モーメントｍ_ｉｊｋｌ（１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）を計算するステップと、コンピュータを用いて、前記分散共分散行列Σから行列Ｔを計算するステップと、前記データ・ベクトルが従う分布が、Ｘ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ＞ｎである。）個の確率変数ξ_１〜ξ_ｍからなる確率ベクトル（ξ_１，・・・，ξ_ｍ）´である。独立な確率変数ξ₁，・・・，ξ_mは、パラメータζ_j及びκ_jについて、Ｅ（ξ_j）＝０，Ｅ（ξ_j ²）＝１，Ｅ（ξ_j ³）＝ζ_j，Ｅ（ξ_j ⁴）＝κ_j（１≦ｊ≦ｍ）を満たす。また、ξ_jに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν²）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。）を満たす非正規分布を有すると仮定して、コンピュータを用いて、３次モーメントｍ_ｉｊｋ（１≦ｉ≦ｊ≦ｋ≦ｎ）及び４次モーメントｍ_ｉｊｋｌ（１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）とｆ_ｉｊｋ（Ｔ，ζ，γ，κ，β），ｆ_ｉｊｋｌ（Ｔ，ζ，γ，κ，β）（ここでζ＝（ζ_１，・・・，ζ_ｍ），κ＝（κ_１，・・・，κ_ｍ）である。）との差を損失とする損失関数を導入し、ＴＴ´＝Σ（Ｔ´はＴの転置行列である。）という条件の下で、全体として評価するリスク関数を最小化するようにパラメータＴ，（ζ_１，・・・，ζ_ｍ），γ，（κ_１，・・・，κ_ｍ），βを決定するステップと、コンピュータを用いて、前記決定されたパラメータ（ζ_１，・・・，ζ_ｍ），γ，（κ_１，・・・，κ_ｍ），βに基づいて、式（５４）及び表１１により（β_１とβ_２は、それぞれ歪度の２乗と尖度である。）、前記確率ベクトル（ξ_１，・・・，ξ_ｍ）´及び確率変数νがピアソン分布のどの型に属するかを決定するステップと、

コンピュータを用いて乱数を発生し、この乱数に基づいて前記確率ベクトル（ξ_１，・・・，ξ_ｍ）´及び確率変数νについて乱数を計算するステップと、
コンピュータを用いて、Ｘ＝νＴξに基づいて、標準化されたＸの乱数を計算するステップと、コンピュータを用いて、標準化されたＸの乱数を標準化前の乱数に変換するステップと、を有する。

ただし、式（５４）及び表１１において、β_１とβ_２は、それぞれ歪度の２乗と尖度であり、κは、前記パラメータκ＝（κ₁，・・・，κ_ｍ）とは異なる。

好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、ｎ個のデータ

からなるデータ・ベクトルＸ_ｉ´について、経験分布のデータ｛Ｘ_１´，・・・，Ｘ_Ｎ´｝を取得するステップと、コンピュータを用いて、前記データ｛Ｘ_１´，・・・，Ｘ_Ｎ´｝を標準化して｛Ｘ_１，・・・，Ｘ_Ｎ｝とするステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ_１，・・・，Ｘ_Ｎ｝に基づいて、分散共分散行列Σを計算するステップと、
コンピュータを用いて、前記標準化されたデータ｛Ｘ_１，・・・，Ｘ_Ｎ｝に基づいて、３次モーメントｍ_ｉｊｋ（１≦ｉ≦ｊ≦ｋ≦ｎ）を計算するステップと、コンピュータを用いて、前記標準化されたデータ｛Ｘ_１，・・・，Ｘ_Ｎ｝に基づいて、４次モーメントｍ_ｉｊｋｌ（１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）を計算するステップと、コンピュータを用いて、前記分散共分散行列ΣからＴ＝Σ^１／２（行列の平方根）を満たす行列Ｔを計算するステップと、前記データ・ベクトルが従う分布が、Ｘ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ＝ｎである。）個の確率変数ξ_１〜ξ_ｍからなる確率ベクトル（ξ_１，・・・，ξ_ｍ）´である。独立な確率変数ξ₁，・・・，ξ_mは、パラメータζ_j及びκ_jについて、Ｅ（ξ_j）＝０，Ｅ（ξ_j ²）＝１，Ｅ（ξ_j ³）＝ζ_j，Ｅ（ξ_j ⁴）＝κ_j（１≦ｊ≦ｍ）を満たす。また、ξ_jに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν²）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。）を満たす非正規分布を有すると仮定して、コンピュータを用いて、３次モーメントｍ_ｉｊｋ（１≦ｉ≦ｊ≦ｋ≦ｎ）及び４次モーメントｍ_ｉｊｋｌ（１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）とｆ_ｉｊｋ（ζ，γ，κ，β），ｆ_ｉｊｋｌ（ζ，γ，κ，β）（ここでζ＝（ζ_１，・・・，ζ_ｍ），κ＝（κ_１，・・・，κ_ｍ）である。）との差を損失とする損失関数を導入し、全体として評価するリスク関数を最小化するようにパラメータ（ζ_１，・・・，ζ_ｍ），γ，（κ_１，・・・，κ_ｍ），βを決定するステップと、コンピュータを用いて、前記決定されたパラメータ（ζ_１，・・・，ζ_ｍ），γ，（κ_１，・・・，κ_ｍ），βに基づいて、式（５５）及び表１２により（β_１とβ_２は、それぞれ歪度の２乗と尖度である。）、前記確率ベクトル（ξ_１，・・・，ξ_ｍ）´及び確率変数νがピアソン分布のどの型に属するかを決定するステップと、

コンピュータを用いて乱数を発生し、この乱数に基づいて前記確率ベクトル（ξ_１，・・・，ξ_ｍ）´及び確率変数νについて乱数を計算するステップと、コンピュータを用いて、Ｘ＝νＴξに基づいて、標準化されたＸの乱数を計算するステップと、コンピュータを用いて、標準化されたＸを標準化前の乱数に変換するステップと、を有する。
ここで、経験分布のデータ｛Ｘ_１´，・・・，Ｘ_Ｎ´｝とは、次の意味である
。

ここでは、Ｘ_ｒ´（１≦ｒ≦Ｎ）はデータ・ベクトルを表し、先のデータの分布（Ｘ_１，・・・，Ｘ_Ｎ）´での表記Ｘとは異なる。標準化されたデータ｛Ｘ_１，・・・，Ｘ_Ｎ｝も同様の意味である。

ただし、式（５５）及び表１２において、β_１とβ_２は、それぞれ歪度の２乗と尖度であり、κは、前記パラメータκ＝（κ₁，・・・κ_ｍ）とは異なる。

好ましくは、前記損失関数（Ｌ）及びリスク関数（Ｒ）は、式（５６）及び（５７）、式（５８）及び（５９）、又は式（６０）及び（６１）のいずれかの組で与えられる請求項１３又は１４に記載の多変量非正規分布のパラメータ推定方法。

好ましくは、前記決定した型に応じて、次の表１３及び表１４に示す生成法の少なくとも１つを用いて乱数Ｚを発生する。

好ましくは、コンピュータを用いて、式（６２）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンＩＶ型分布に従う乱数を発生させる。

好ましくは、前記乱数を発生するステップは、加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。

本発明に係る多変量非正規分布のパラメータ推定方法は、コンピュータを用いて、式（６３）に示すように、所定の確率分布に従う確率変数ν、階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ、所定の確率分布に従う確率ベクトルξ＝（ξ₁，・・・，ξ_m）´の積として与えられる確率ベクトルＸ＝（ｘ₁，・・・，ｘ_n）´で与えられる多次元非正規分布を計算するステップと、

最尤法によりパラメータを推定するステップと、を有する。

好ましくは、前記最尤法によりパラメータを推定するステップは、コンピュータを用いて、ｎ次元空間を分割するステップと、コンピュータを用いて、経験分布のデータ｛Ｘ₁，・・・，Ｘ_N｝について、Ｘ_iが属する各分割に属する乱数の数を乱数の総数で除して、次の式（６４）で与えられる当該分割の確率を求めるステップと、

コンピュータを用いて、前記分割の確率を当該分割のＮ次元体積で除して前記Ｘ_iが属する、次の式（６５）で与えられる分割の尤度ｆ_i（θ）を求めるステップと、

コンピュータを用いて、尤度の積Π_i=1 ^Nｆ_i（θ）又は対数尤度の和Σ_i=1 ^Nｌｏｇ ｆ_i（θ）が最大になるようなパラメータθを推定するステップと、を有する。

好ましくは、確率変数ξ₁，・・・，ξ_m及びνには、ピアソン分布系に従う乱数を使う。

好ましくは、確率変数ξ₁，・・・，ξ_m及びνには、少なくとも２以上の型のピアソン分布を用いる。

好ましくは、コンピュータを用いて、式（６６）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンＩＶ型分布に従う乱数を発生させる。

好ましくは、コンピュータを用いて、ユアン及びベントラーの方法Ｉに基づいて、次の関係式（６７）及び（６８）を用いてｎ次元経験分布にｎ次元多変量非正規分布をあてはめるステップを有し、このあてはめるステップで求めたパラメータの近傍について前記パラメータθを推定する。

ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄ_ｎは、ｎ次のデュプリケーション行列であり、Ｄ_ｎ ⁺は、Ｄ_ｎのムーア・ペンローズの一般化逆行列であり、

は、クロネッカー積であり、Ｅ_iiは、ｅ_iを第ｉ列単位としてｅ_iｅ_i ^'である。

ここで、前記ユアン及びベントラーの方法Ｉは、次の通りである。独立な確率変数ξ₁，・・・，ξ_mは、パラメータζ_j及びκ_jについて、Ｅ（ξ_j）＝０，Ｅ（ξ_j ²）＝１，Ｅ（ξ_j ³）＝ζ_j，Ｅ（ξ_j ⁴）＝κ_j（１≦ｊ≦ｍ）を満たす。

ξ_jに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν²）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。

階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔ_ij）は、行列Σ＝（σ_ij）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。

このとき、次の式（６９）で与えられる確率ベクトルＸ＝（ｘ₁，・・・，ｘ_n）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ₁，・・・，ξ_m）´である。

好ましくは、前記パラメータθは、パラメータＴ，ζ＝（ζ₁，・・・，ζ_ｍ），γ，κ＝（κ₁，・・・，κ_ｍ），βの少なくとも１つを表す。

好ましくは、コンピュータを用いて加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。

好ましくは、本発明に係る多変量非正規分布のパラメータ推定方法は、コンピュータを用いて、ユアン及びベントラーの方法Ｉに基づいて、次の関係式（７０）及び（７１）を用いてｎ次元経験分布にｎ次元多変量非正規分布をあてはめるステップと、コンピュータを用いて、ｎ次元経験分布がＸ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ≧ｎである。）個の確率変数ξ_１〜ξ_ｍからなる確率ベクトル（ξ_１，・・・，ξ_ｍ）´である。独立な確率変数ξ₁，・・・，ξ_mは、パラメータζ_j及びκ_jについて、Ｅ（ξ_j）＝０，Ｅ（ξ_j ²）＝１，Ｅ（ξ_j ³）＝ζ_j，Ｅ（ξ_j ⁴）＝κ_j（１≦ｊ≦ｍ）を満たす。また、ξ_jに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν²）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。）を満たす非正規分布を有すると仮定して、前記Ｘについて乱数を計算するステップと、コンピュータを用いて、各確率変数ｘ_１〜ｘ_ｎについての空間を、所定間隔Δｈで分割し、このｎ次元空間を（Δｈ）^ｎの超立方体へ分割するステップと、コンピュータを用いて、各データ・ベクトルの属する区間（Δｈ）^ｎ内に存在する乱数の個数を乱数の総数で除して、当該区間（Δｈ）^ｎの確率Ｐ_ｉ（θ）を計算するステップと、コンピュータを用いて、前記確率Ｐ_ｉ（θ）をｎ次元体積（Δｈ）^ｎで除して、前記データ・ベクトルの属する区間（Δｈ）^ｎの尤度ｆ_ｉ（θ）を計算するステップと、コンピュータを用いて、前記尤度ｆ_ｉ（θ）の積Π_ｉ＝１ ^Ｎｆ_ｉ（θ）又は対数尤度の和Σ_i=1 ^Nｌｏｇ ｆ_i（θ）が最大になるようなパラメータθを計算するステップと、を有する。

ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄ_ｎは、ｎ次のデュプリケーション行列であり、Ｄ_ｎ ⁺は、Ｄ_ｎのムーア・ペンローズの一般化逆行列であり、

は、クロネッカー積であり、Ｅ_iiは、ｅ_iを第ｉ列単位としてｅ_iｅ_i ^'である。

ここで、前記ユアン及びベントラーの方法Ｉは、次の通りである。独立な確率変数ξ₁，・・・，ξ_mは、パラメータζ_j及びκ_jについて、Ｅ（ξ_j）＝０，Ｅ（ξ_j ²）＝１，Ｅ（ξ_j ³）＝ζ_j，Ｅ（ξ_j ⁴）＝κ_j（１≦ｊ≦ｍ）を満たす。

ξ_jに独立な確率変数νは、パラメータγ及びβについて、Ｅ（ν）＝０，Ｅ（ν²）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。

階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔ_ij）は、行列Σ＝（σ_ij）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。

このとき、次の式（７２）で与えられる確率ベクトルＸ＝（ｘ₁，・・・，ｘ_n）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ₁，・・・，ξ_m）´である。

好ましくは、確率変数ξ₁，・・・，ξ_m及びνには、ピアソン分布系に従う乱数を使う。

好ましくは、確率変数ξ₁，・・・，ξ_m及びνには、少なくとも２以上の型のピアソン分布を用いる。

好ましくは、コンピュータを用いて、式（７３）に示すピアソンＩＶ型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンＩＶ型分布に従う乱数を発生させる。

好ましくは、コンピュータを用いて加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。

本発明に係る予想損失額の計算方法は、前記手順に従い、ポートフォリオを構成するｎ種類の資産のｎ次元損益分布にｎ次元多変量非正規分布をあてはめ、シミュレーションによって予想損失額を計算する。

本発明に係る金融商品のシミュレーション方法は前記手順を用い、前記多変量非正規分布に従う乱数に基づいて金融商品の特性をコンピュータを用いてシミュレーションする。

本発明に係る天候・保険デリバティブの設計又は価格付け方法は、前記手順を用いて天候・保険デリバティブを設計又は価格付けする。

本発明に係るアセット・アロケーションの方法は、前記手順を用いてアセット・アロケーションを行う。

本発明に係るイオン注入のシミュレーション方法は、前記手順を用い、前記多変量非正規分布に従う乱数に基づいてイオン注入によるイオンの分布をコンピュータを用いてシミュレーションする。

本発明に係る能力評価の方法は、前記手順を用いて能力評価を行う。

本発明に係るコンピュータにより実行可能なコンピュータ・プログラムは、前記手順を実現するものである。好ましくは、前記コンピュータ・プログラムは、光ディスク、磁気ディスク、光磁気ディスク又は磁気テープの少なくとも１つによって提供される。

本発明によると、ピアソン分布系及び一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンＩＶ型分布に従う乱数発生方法及びコンピュータ・プログラムを提供することができる。

以下、本発明に係る多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンＩＶ型分布に従う乱数発生方法及びコンピュータ・プログラムの実施の形態並びにリスク・マネジメントにおける予想損失額の計算方法について、図面を参照して詳細に説明する。

以下の本実施の形態は、大きくは次の構成からなっている。

１．ピアソン分布系及び乱数生成法
１．１ ピアソン分布系
１．２ ピアソン分布ＩＶ型に従う乱数の生成
１．３ ピアソン分布系に従う乱数の生成
２． ピアソン分布系を用いた多変量非正規分布の生成
２．１ 多変量正規分布の確率ベクトル
２．２ シミュレーション
２．３ 応用のためのあてはめプロシージャ
３． 最尤法
４． 実施例
５． 変形例
５．１ 金融商品のシミュレーション方法
５．２ 半導体へのイオン注入
５．３ アセット・アロケーション
５．４ 能力評価
本実施の形態は、様々な歪度と尖度を有する広いクラスの分布を表すことができるピアソン分布系を用い、多変量非正規分布をシミュレーションする方法を提案する。第１に、実現が困難であったピアソンＩＶ型分布を含むピアソン分布系に従う乱数を生成するプロシージャを導く。第２に、ピアソン分布系を適用し、ユアンとベントラー（１９９７）によって開発された方法を用い、多変量非正規分布に従う乱数を生成する。第３に、本願発明者によって求められた３次のモーメント行列を用い、種々の標本分布について３次と４次のモーメント行列をあてはめ、近似的な多変量非正規分布を構築することができる。第４に、最尤法によるパラメータ推定方法を提案する。このパラメータ推定方法によって、少数のデータからモーメント法によっては実現できない正確さで、多変量非正規分布を推定することができることを示す。そして、これらの方法が、様々なビジネスに応用できることを、リスク・マネジメントにおける予想損失額の計算について具体的に説明する。最後に、変形例として、金融商品のシミュレーション方法、半導体へのイオン注入、アセット・アロケーション及び能力評価について簡単に触れる。

１ ピアソン分布系及び乱数生成法
１．１ ピアソン分布系
ピアソン分布系の型は、次のパラメータκによって分類される。パラメータκは、次の式（７４）によって定義される。

ただし、β₁とβ₂は、それぞれ歪度の２乗と尖度である。

ピアソン分布は、κ＜０，０＜κ＜１及び１＜κについて、それぞれＩ型、ＩＶ型及びＶＩ型と呼ばれる。これらの３つの場合は、「主要型」と呼ばれる（エルダートン及びジョンソン（Elderton and Johnson）（１９６９））。ＩＩＩ型（κ＝±∞）は、Ｉ型とＶＩ型の境界にある。Ｖ型（κ＝１）は、第ＩＶ型と第ＶＩ型の境界にある。κ＝０の場合については、ＶＩＩ型、ＩＩ型及び正規分布が含まれる。

ピアソン分布系においては、型が指定されると、この分布の３つ又は４つのパラメータは、平均、分散、歪度及び尖度から決定される。歪度の２乗と尖度の間の典型的な関係は、多数の教科書に示されている（ジョンソン、コッツ及びバラカリシュナン（Johnson, Kotz and Balakrishnan）（１９９４）、スチュアート及びオード（Stuart and Ord）（１９９４）等）。

まず、表１５〜１８に示すように、ピアソン分布系は、最も適切な表現によると、８つの型に分類され、密度関数、及び平均、分散、歪度及び尖度による陽な形が与えられる。

次に、モーメントから各型のパラメータを得るための変換式とプロシージャを示す。各型の変換式を表１５〜１８に示す。変換式によって、全てのパラメータを順に容易に計算することができる。

表１６において、ピアソンＩ型分布について、ｐ，ｑ＞１のとき単峰型、０＜ｑ≦１＜ｐのときＪ型、０＜ｐ≦１＜ｑのときＪ型、０＜ｐ，ｑ≦１のときＵ型になる。なお、ｐ＝ｑ＝１のときには一様分布になる。また、ピアソンＩＩ型分布について、ｐ＞１のとき単峰型、０＜ｐ＜１のときＵ型になる。なお、ｐ＝１のときには一様分布になる。

１．２ ピアソン分布ＩＶ型に従う乱数の生成
本実施の形態においては、ピアソン分布ＩＶ型の確率密度関数は、次の式（７５）によって定義される。ただし、Γ（・）はガンマ関数、Ｂ（・）はベータ関数である。

ドブロイ（Devrye）（１９８６）によると、確率変数XがピアソンＩＶ型分布であってｂ＞１のとき、次の式（７６）により分布するとすると、ａｒｃｔａｎ（ｘ）は、式（７７）によって定義される対数凹密度を有する。

ドブロイ（１９８４）は、対数凹密度についての棄却アルゴリズム（指数版）が、この分布に従う乱数を生成するのに適用できることを示唆した。しかし、その正規化定数は、著書において問題として残されていた（ドブロイ（１９８６）の練習２．７）。論文（永原（１９９９））では、この問題を解析的に解き、その正規化定数が容易に計算できることを示した（式（７８））。ピアソンＩＶ型分布に従う乱数の生成のアルゴリズムを次に示す。Ｃは、ピアソンＩＶ型分布の正規化定数を示す。

このアルゴリズムは、正規化定数の解析解又はこの展開を計算するステップと、棄却法により乱数を発生するステップを含んでいる。

ここで、乱数のシーズについて説明する。ある確率分布に従う乱数を発生させるにあたっては、一様乱数がその基本になっている。

原理的には、一様乱数が与えられれば、逆関数法によって、様々な確率分布に従う乱数が発生できる。本発明においても、この一様乱数をもとにして、ピアソン分布系の確率分布に従う乱数の発生のために使う指数分布、ガンマ分布、正規分布に従う乱数を発生している。この一様乱数を発生させる方法には、大別すると物理乱数、擬似乱数及び準乱数及び低ディスクレパンシー列が含まれる。準乱数及び低ディスクレパンシー列とは、擬似乱数とは違って、ランダムな性質はもっていないが、一様に分布する点列である。

擬似乱数には、合同法、加算生成法、Ｍ系列、一般化フィードバック・シフトレジスタ法（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.）がある。

この中で、一般的な方法は、擬似乱数の合同法である。ところが、最近では、合同法よりはるかに長い周期を持つ、ＧＦＳＲ法やメルセンヌツイスターが開発されてきている。さらに、準乱数や低ディスクレパンシー列が研究されていて、実用にも供されてきている。さて、これらの周期のより長い擬似乱数（例えばメルセンヌツイスター）や準乱数（低ディスクレパンシー列）を本実施の形態に適用したところ、本実施の形態の実施に重要な乱数の性質を大幅に改善することができた。従来の合同法では実現できなかった正確性で、目標とする確率分布の平均、分散、歪度、尖度を再現した乱数を得ることができた。

本実施の形態の一番の重要な点は、非常に多数の乱数で、もはや解析的には表現できないような複雑な性質を持った確率分布を表現することである。この目的のためには、従来の合同法よりさらに正確な確率分布を再現できる乱数が望まれ、特に本実施の形態の特徴である、歪度、尖度が正規分布と違う分布（ピアソンＩＶ型等）については、擬似乱数では合同法より長い周期のもの（例えばメルセンヌツイスター）、又は準乱数を使うことによって大幅に改善された、より実用的な正確性が達成された。

以下に１次元の例を述べる。まず、ピアソンＩＶ型に関する例では、乱数の個数が多くなると、合同法に比べ、メルセンヌツイスター、準乱数（ここでは、ソボル列を使用する。ソボル列は、低ディスクレパンシー列でもある。）を使った方が、平均、分散、歪度、尖度すべてにわたって目標とする値（表の一段目の数）に大幅に近くなっている。

例えば、表２１（合同法）、表２５（メルセンヌツイスター）、表２９（ソボル列）において、１０００００個の乱数で１０００回シミュレーションした結果、平均について、合同法では０．０１０２３、メルセンヌツイスターでは０．０００１８、準乱数では、−０．００００１となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは４桁、準乱数では、５桁までと大幅に改善している。同様に、分散について、合同法では１．０３４３３、メルセンヌツイスターでは１．００００９、準乱数では、０．９９９９９となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは５桁、準乱数では、６桁までと大幅に改善している。

また、本実施の形態で重要な歪度と尖度についてはさらに改善している。同じ表で、歪度について目標値は０．７５で、合同法では０．７９６５１、メルセンヌツイスターでは０．７５０４１、準乱数では、０．７５００７となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは４桁、準乱数では、５桁までと大幅に改善している。尖度について目標値は５．０で、合同法では５．３６３３４、メルセンヌツイスターでは５．０００５７、準乱数では、５．００２５７となり、合同法での小数点以下１桁目の正確さが、メルセンヌツイスターでは４桁、準乱数でも、３桁までと大幅に改善している。また、メルセンヌツイスター及び準乱数では、乱数の数が増えれば、それに従って目標の値にそれぞれ収束している点が合同法より優れており、この性質は、本実施の形態において乱数が持つべき最も重要な性質である。

ここで、モンテカルロ・シミュレーション技法として対照変量（負の相関）法、制御変量法、モーメント適合法、回帰分析法、非線形回帰分析法、加重抽出法、層化抽出法、ラテン・ハイパーキューブ法、マルチンゲール分散減少法、条件付モンテカルロ法のような分散減少法が容易に適用できる。これら技法は、乱数の発生法に依らず、計算を改良する方法である。

表１９〜３０において、数は、乱数の個数を表し、値は、１０００回の平均である。各括弧内の値は各値の標準偏差である。この結果、アルゴリズムを用いて妥当な乱数を得た。

１．３ ピアソン分布系に従う乱数の生成
ＩＶ型を除くピアソン分布系については、正規分布とガンマ分布に従う乱数を用いることで、容易に乱数を生成できることが知られている（マクグラス及びアービング（McGrath and Irving）（１９７３））。本実施の形態では、表１５〜１８のような定義を実現する実用的なプロシージャを開発する。これらの結果は、表３１と表３２に示されている。典型的な例は、表３３〜５３に示されている。数は、乱数の個数を示し、値は１０００回の平均である。各括弧内の値は、各値の標準偏差である。この結果、アルゴリズムを用いて妥当な乱数を得る。

ピアソンＩＶ型以外の分布系については、合同法に比べ、メルセンヌツイスター、準乱数（ここでは、ソボル列を使用。ソボル列は低ディスクレパンシー列でもある。）を使った方が、平均、分散では、ほとんどすべてで大幅に改善した。重要な歪度、尖度については、幅の広い領域を受け持つピアソンＶＩ型及びピアソンＩ型では、１０００００個の乱数で１０００回シミュレーションの結果で、まず、ピアソンＶＩ型で、表３３（合同法）、表４０（メルセンヌツイスター）及び表４４（ソボル列）において、歪度について目標値は１．０で、合同法では０．９９３２３、メルセンヌツイスターでは０．９９９７５、準乱数では０．９９９６８となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数でも、同様に４桁までと大幅に改善している。同じ表で、尖度について目標値は４．８で、合同法では４．７３３８６、メルセンヌツイスターでは４．７９８６３、準乱数では４．７９８３６となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは３桁、準乱数でも、同様に３桁までと改善している。

次に、ピアソンＩ型について、表３４（合同法）、表４１（メルセンヌツイスター）及び表４８（ソボル列）において、歪度について目標値は１．０で、合同法では１．００６１０、メルセンヌツイスターでは０．９９９４５、準乱数では１．００００４となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数は５桁までと大幅に改善している。同じ表で、尖度について目標値は４．０で、合同法では４．００５８０、メルセンヌツイスターでは３．９９７３７、準乱数では４．０００３７となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは同様の３桁、準乱数では、４桁までと改善している。

ピアソンＩＩ型では、表３５（合同法）、表４２（メルセンヌツイスター）及び表４９（ソボル列）において、歪度について目標値は０．０で、合同法では０．０００１２、メルセンヌツイスターでは−０．００００５、準乱数では０．００００５となり、合同法での小数点以下４桁目の正確さが、メルセンヌツイスターでは５桁、準乱数でも、同様に５桁までと大幅に改善している。尖度について目標値は２．５で、同じ表で、合同法では２．４９４９０、メルセンヌツイスターでは２．４９９９２、準乱数では２．５００１４となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは５桁、準乱数でも、４桁までと改善している。

正規分布では、表３６（合同法）、表４３（メルセンヌツイスター）及び表５０（ソボル列）において、歪度について目標値は０．０で、合同法では０．００４３９、メルセンヌツイスターでは０．０００２３、準乱数では０．０００００となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数では６桁までと大幅に改善している。尖度について目標値は３．０で、合同法では２．９８１３０、メルセンヌツイスターでは３．０００１８、準乱数では３．００００１となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは４桁、準乱数でも、５桁までと改善している。

ピアソンＶＩＩ型では、表３７（合同法）、表４４（メルセンヌツイスター）、表５１（ソボル列）において、歪度について目標値は０．０で、合同法では０．００２８２、メルセンヌツイスターでは０．０００６２、準乱数では−０．０００４７となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数でも４桁までと大幅に改善している。尖度について目標値は４．０で、合同法では３．９６６３９、メルセンヌツイスターでは３．９９６９０、準乱数では３．９９８９４となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは３桁、準乱数でも、３桁までと改善している。

ピアソンＩＩＩ型では、表３８（合同法）、表４５（メルセンヌツイスター）、表５２（ソボル列）において、歪度について目標値は１．０で、合同法では１．００１２８、メルセンヌツイスターでは０．９９９５１、準乱数では０．９９９９６となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは４桁、準乱数では５桁までと大幅に改善している。尖度について目標値は４．５で、合同法では４．４５７１１、メルセンヌツイスターでは４．４９８６４、準乱数では４．４９９７８となり、合同法での小数点以下２桁目の正確さが、メルセンヌツイスターでは３桁、準乱数では、４桁までと改善している。

ピアソンＶ型では、表３９（合同法）、表４６（メルセンヌツイスター）、表５３（ソボル列）において、歪度について目標値は１．０で、合同法では１．００４５６、メルセンヌツイスターでは０．９９８４２、準乱数では０．９９９９２となり、合同法での小数点以下３桁目の正確さが、メルセンヌツイスターでは同様の３桁、準乱数では５桁までと大幅に改善している。尖度について目標値は４．９７０で、合同法では５．１５５０５、メルセンヌツイスターでは４．９５８２９、準乱数では４．９６９２５となり、合同法での小数点以下１桁目の正確さが、メルセンヌツイスターでは２桁、準乱数では、４桁までと大幅に改善している。また、以上分布についても、ピアソンＩＶ型のときと同様に、メルセンヌツイスター及び準乱数では、乱数の数が増えれば、それに従って目標の値にそれぞれ収束している。この点が合同法より優れており、この性質は本発明において乱数が持つべき最も重要な性質である。

２ ピアソン分布系を用いた多変量非正規分布の生成
２．１ 多変量非正規分布の確率ベクトル
ユアン等（１９９７）によると、ユアン等が方法Iと称したものは、次の通りである。ξ₁，・・・，ξ_mは、独立な確率変数であり、Ｅ（ξ_j）＝０，Ｅ（ξ_j ²）＝１，Ｅ（ξ_j ³）＝ζ_j，Ｅ（ξ_j ⁴）＝κ_j及びξ＝（ξ₁，・・・，ξ_m）´である。νは、ξに独立な確率変数であり、Ｅ（ν）＝０，Ｅ（ν²）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。ここで、この分布は、νが１のとき、ピアソン分布系の確率密度関数によって表されることが特記される。また、Ｔ＝（ｔ_ij）は、階数ｐの非確率ｐ行ｍ列の行列であり、ＴＴ´＝Σ及びｍ≧ｐを満たす。すると、次の式（７９）で与えられる確率ベクトルは、一般に、Ｃｏｖ（Ｘ）＝Σを満たす非楕円分布に従うようになる。ここで、Ｃｏｖ（・）は、ベクトルＸ＝（ｘ₁，・・・，ｘ_n）´の分散共分散行列である。

ｘ_iの周辺の歪度及び尖度は、それぞれ次の式（８０）及び（８１）によって与えられる。

そして、４次のモーメント行列Γ＝ｖａｒ（ｖｅｃｈ（ＸＸ´））は、次の式（８２）によって与えられる。

ただし、ｔ_jは、Ｔ＝（ｔ_ij）の第ｊ列ベクトルであり、

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

統計モデリングの見地から、３次のモーメント行列が必要である。本願発明者は、３次のモーメント行列について、次の定理を得ることができた。

定理３．１（３次のモーメント行列）３次のモーメント行列Λは、次の式（８３）によって与えられる。

ここで、ｅ_iは、第ｉ列単位であり、Ｅ_ii＝ｅ_iｅ_i ^'である。

証明 ３次のモーメント行列Λは、Ｙ＝Ｔξとして、次の式（８４）のようになる。

ＹＹ´＝Ｔξξ´Ｔ´及びｖｅｃｈ（ＹＹ´）＝Ｄ_p ⁺ｖｅｃ（ＹＹ´）である
ので、

となる。ここで、記号

は、行列のクロネッカー積を表す。以下でも同様である。

したがって、次の式（８６）を得る。

次の式（８７）は、明らかである。

この定理は、式（８４），（８６）及び（８７）から得られる。

証明終わり
ユアン及びベントラー（１９９９ａ，１９９９ｂ及び２０００）は、検定統計理論を研究した。ユアン等は、次の式（８８）を示唆した。

ただし、Π＝∂ｆ／∂ｖｅｃｈ´（Σ）の階数ｓの行列として、Ω＝ΠΓΠ´であり、Ｓは、標本分散共分散行列である。また、ユアン等は、以前、Πｖｅｃｈ（ｔ_jｔ_j´）＝０，ｊ＝１，・・・，ｍを仮定すると、次の式（８９）が成立することを示唆している。

ただし、

である。

本実施の形態は、様々な３次と４次のモーメント行列を用い、多変量非正規分布に従う乱数を得ることを目的としている。Ｔ，ξ，νを変化させることによって、異なる３次と４次のモーメント行列を有する様々な多変量非正規分布を得ることができる。まず、Ｘ₁，Ｘ₂，ξ₁及びξ₂を有するｐ＝ｍ＝２の典型的な場合についてシミュレーションを行う。次に、３次と４次のモーメントを経験データの１つにあてはめるプロシージャを提案する。

なお、ピアソン分布系のほかには、4次までのモーメントで表せるジョンソン分布系やラムダ分布などがνやξの分布として採用できる。

２．２ シミュレーション
本実施の形態では、ｘ₁，ｘ₂，ξ₁及びξ₂を有するｐ＝ｍ＝２の場合についてシミュレーションする。このシミュレーションにおいて、３次と４次のモーメント行列の要素は、次の通りである。

［３次モーメント］

ただし、ｆ₁，ｆ₂，ｆ₃，ｆ₄は、Ｘ＝νＴξで発生される確率ベクトルＸ＝（ｘ_１，ｘ_２）´のそれぞれ３次モーメントＥ（ｘ₁ ³），Ｅ（ｘ₂ ³），Ｅ（ｘ₁ ²ｘ₂），Ｅ（ｘ₁ｘ₂ ²）の表現である。

［４次モーメント］

ただし、ｆ₅，ｆ₆，ｆ₇，ｆ₈，ｆ₉は、Ｘ＝νＴξで発生される確率ベクトルＸ＝（ｘ_１，ｘ_２）´のそれぞれ４次モーメントＥ（ｘ₁ ⁴），Ｅ（ｘ₂ ⁴），Ｅ（ｘ₁ ³ｘ₂），Ｅ（ｘ₁ ²ｘ₂ ²），Ｅ（ｘ₁ｘ₂ ³）の表現である。ここで、σ_１１＝σ_２２＝１であり、また、ｔ_１１ ^２＝ｔ_２２ ^２＝１−ｔ_１２ ^２＝１−ｔ_２１ ^２である。

ユアンとベントラー（２０００）の例を採用し、ρ＝０．５及びｔ₁₁＝ｔ₂₂＝０．２５８８，ｔ₁₂＝ｔ₂₁＝０．９６５９とする。結果は、表５４〜８０に示されている。特定の歪度及び尖度を有する、各分布のξ₁，ξ₂，νから乱数を得た後、この乱数を平均と分散がそれぞれ０と１になるように標準化する。続いて、式Ｘ＝νＴξによって、ｘ₁とｘ₂のための乱数Ｘ_１，Ｘ_２を得る。各表において、Ｅ（・）は、これらの乱数を用いて計算される。括弧内は、３次と４次のモーメント行列の式を用いて計算された理論値である。これらの結果によると、１０，０００，０００点と１，０００，０００点の場合、この方法による乱数は、理論的なモーメント、すなわち分散共分散、３次と４次のモーメント行列を有する目的の分布を正確に表すことができる。

ここで、異なるピアソン分布の型の組み合わせも可能であることは特記されるべきであることである。例えば、金融リスク、オペレーショナル・リスク、保険リスク、災害リスク、事業リスクは、一般に異なったピアソン分布の型で記述される。本実施の形態によると、異なったピアソン分布の組み合わせも可能なので、このような異なった種類のリスクの組み合わせも統合して管理することができる。例えば、各種の資産を含んだポートフォリオの予想損失額を統合して計算することができる。

なお、金融での最大損失額の計算では、極値分布（ガンベル分布、フレッシュ分布、ワイブル分布）や一般化パレート分布・対数正規分布が使われるが、これらの分布もこの発明の方法で、νやξの分布として採用できる。このことによって、多変量化や、よりデータを正確に表す分布を推定できる。

ピアソンＩＶ型の場合、２次元においてξ_１、ξ_２、νすべてに、ピアソンＩＶ型を用いた例を表５４（合同法）、表６３（メルセンヌツイスター）、表７２（準乱数 ソボル列）に示す。以下、ｘ_１、ｘ_２は確率変数、Ｘ_１，Ｘ_２は乱数を表す。

まず、目標とする分布の特性はξ_１（歪度０．７５、尖度５．０）、ξ_２（歪度１．０、尖度５．５）、ν（歪度０．５０、尖度４．０）で、それから生成される多変量非正規分布のモーメント特性の理論値は、Ｅ（ｘ_１ ^３）＝０．４５７０８、Ｅ（ｘ_２ ^３）＝０．３４６６０、Ｅ（ｘ_１ ^２ｘ_２）＝０．１４４９９、Ｅ（ｘ_１ｘ_２ ^２）＝０．１２２８９、Ｅ（ｘ_１ ^４）＝２０．７３８６５、Ｅ（ｘ_２ ^４）＝１９．００６７８、Ｅ（ｘ_１ ^３ｘ_２）＝８．４６５１５、Ｅ（ｘ_１ ^２ｘ_２ ^２）＝７．１２３９０、Ｅ（ｘ_１ｘ_２ ^３）＝８．０３２２１である。

例えば、１０，０００，０００個のシミュレーションの場合、表５４（合同法）では、それぞれ、ξ_１（歪度０．８３０８２、尖度５．２８９２８）、ξ_２（歪度１．０５９０６、尖度６．１９４９５）、ν（歪度０．５８８０７、尖度４．１７７７４）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．５７５０７（０．４５７０８）、Ｅ（Ｘ_２ ^３）＝０．４５４９９（０．３４６６０）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．１８７１５（０．１４４９９）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．１６３８３（０．１２２８９）、Ｅ（Ｘ_１ ^４）＝２４．１４４９３（２０．７３８６５）、Ｅ（Ｘ_２ ^４）＝２０．９４８９７（１９．００６７８）、Ｅ（Ｘ_１ ^３Ｘ_２）＝９．５３０５４（８．４６５１５）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝７．７１１６８（７．１２３９０）、Ｅ（Ｘ_１Ｘ_２ ^３）＝８．７５５００（８．０３２２１）である。

括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ_１、ξ_２、νの実際のシミュレーション時の歪度、尖度（例えば、上の例ではξ_１（歪度０．８３０８２、尖度５．２８９２８）、ξ_２（歪度１．０５９０６、尖度６．１９４９５）、ν（歪度０．５８８０７、尖度４．１７７７４））から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との開きは大きい。

ところが、これをメルセンヌツイスターでシミュレーションすると、表６３（メルセンヌツイスター）にあるように、それぞれ、ξ_１（歪度０．７５２７９、尖度５．０２１２３）、ξ_２（歪度１．００００４、尖度５．５０８３１）、ν（歪度０．５０４１４、尖度４．００７８０）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．４５１８５（０．４５７０８）、Ｅ（Ｘ_２ ^３）＝０．３５０１９（０．３４６６０）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．１４２０２（０．１４４９９）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．１２４０４（０．１２２８９）、Ｅ（Ｘ_１ ^４）＝２０．０３２１０（２０．７３８６５）、Ｅ（Ｘ_２ ^４）＝１９．１６３４１（１９．００６７８）、Ｅ（Ｘ_１ ^３Ｘ_２）＝８．１６２４９（８．４６５１５）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝７．０２１２７（７．１２３９０）、Ｅ（Ｘ_１Ｘ_２ ^３）＝８．０５１１７（８．０３２２１）である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。

また、準乱数（ソボル列）では表７２（ソボル列）にあるように、それぞれ、ξ_１（歪度０．７４９６１、尖度４．９９６９７）、ξ_２（歪度０．９９９８９、尖度５．４９７１６）、ν（歪度０．４９９８４、尖度３．９９６６１）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．４６５１３（０．４５７０８）、Ｅ（Ｘ_２ ^３）＝０．３４７２４（０．３４６６０）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．１５１９５（０．１４４９９）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．１３１４７（０．１２２８９）、Ｅ（Ｘ_１ ^４）＝２０．９６１４０（２０．７３８６５）、Ｅ（Ｘ_２ ^４）＝１８．９６４９３（１９．００６７８）、Ｅ（Ｘ_１ ^３Ｘ_２）＝８．６２４４４（８．４６５１５）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝７．２３７７９（７．１２３９０）、Ｅ（Ｘ_１Ｘ_２ ^３）＝８．１５８６６（８．０３２２１）である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。

ピアソンＶＩ型の場合、２次元においてξ_１、ξ_２、νすべてに、ピアソンＶＩ型を用いた例を表５７（合同法）、表６６（メルセンヌツイスター）、表７５（準乱数 ソボル列）に示す。まず、目標とする分布の特性はξ_１（歪度１．０、尖度４．８）、ξ_２（歪度１．２５、尖度５．５）、ν（歪度０．７５、尖度４．０）で、それから生成される多変量非正規分布のモーメント特性の理論値は、Ｅ（ｘ_１ ^３）＝０．８５７８３、Ｅ（ｘ_２ ^３）＝０．６９２１１、Ｅ（ｘ_１ ^２ｘ_２）＝０．２７４８８、Ｅ（ｘ_１ｘ_２ ^２）＝０．２４１７４、Ｅ（ｘ_１ ^４）＝２０．７３５０６、Ｅ（ｘ_２ ^４）＝１８．３１０４５、Ｅ（ｘ_１ ^３ｘ_２）＝８．４５１７６、Ｅ（ｘ_１ ^２ｘ_２ ^２）＝７．０７３９１、Ｅ（ｘ_１ｘ_２ ^３）＝７．８４５６３である。

例えば、１０，０００，０００個のシミュレーションの場合、表５７（合同法）では、それぞれ、ξ_１（歪度０．９９１２６、尖度４．６４４９５）、ξ_２（歪度１．２２６１５、尖度５．２２３８６）、ν（歪度０．７５３０４、尖度３．９２１２３）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．８４３０２（０．８５７８３）、Ｅ（Ｘ_２ ^３）＝０．６７３４２（０．６９２１１）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．２６８６５（０．２７４８８）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．２４０２５（０．２４１７４）、Ｅ（Ｘ_１ ^４）＝１９．２００２１（２０．７３５０６）、Ｅ（Ｘ_２ ^４）＝１７．１４７０５（１８．３１０４５）、Ｅ（Ｘ_１ ^３Ｘ_２）＝７．９６７９４（８．４５１７６）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝６．７６２１５（７．０７３９１）、Ｅ（Ｘ_１Ｘ_２ ^３）＝７．４１６７９（７．８４５６３）である。

括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ_１、ξ_２、νの実際のシミュレーション時の歪度、尖度（例えば、上の例ではξ_１（歪度０．９９１２６、尖度４．６４４９５）、ξ_２（歪度１．２２６１５、尖度５．２２３８６）、ν（歪度０．７５３０４、尖度３．９２１２３））から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との間に開きがある。

ところが、これをメルセンヌツイスターでシミュレーションすると、表６６（メルセンヌツイスター）にあるように、それぞれ、ξ_１（歪度０．９９９６８、尖度４．８００６２）、ξ_２（歪度１．２５０４２、尖度５．４９７０３）、ν（歪度０．７５０２２、尖度３．９９８７７）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．８６１８５（０．８５７８３）、Ｅ（Ｘ_２ ^３）＝０．６７７４６（０．６９２１１）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．２７７１３（０．２７４８８）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．２３８６８（０．２４１７４）、Ｅ（Ｘ_１ ^４）＝２０．８７６６８（２０．７３５０６）、Ｅ（Ｘ_２ ^４）＝１８．０１４６４（１８．３１０４５）、Ｅ（Ｘ_１ ^３Ｘ_２）＝８．５９７６０（８．４５１７６）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝７．１０８３７（７．０７３９１）、Ｅ（Ｘ_１Ｘ_２ ^３）＝７．７５２９８（７．８４５６３）である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。特に４次のモーメントがより近づいている。

また、準乱数（ソボル列）では表７５（ソボル列）にあるように、それぞれ、それぞれ、ξ_１（歪度１．０００２５、尖度４．８０８３２）、ξ_２（歪度１．２４９９０、尖度５．４９８６２）、ν（歪度０．７５０５２、尖度４．００５３３）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．８５８３４（０．８５７８３）、Ｅ（Ｘ_２ ^３）＝０．６７５８６（０．６９２１１）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．２６４８６（０．２７４８８）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．２２９５７（０．２４１７４）、Ｅ（Ｘ_１ ^４）＝２０．７５１１１（２０．７３５０６）、Ｅ（Ｘ_２ ^４）＝１８．３５８９８（１８．３１０４５）、Ｅ（Ｘ_１ ^３Ｘ_２）＝８．３４４００（８．４５１７６）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝６．９４１４３（７．０７３９１）、Ｅ（Ｘ_１Ｘ_２ ^３）＝７．７０２９１（７．８４５６３）である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。特に４次のモーメントがより近づいている。

ピアソンＩ型の場合、２次元においてξ_１、ξ_２、νすべてに、ピアソンＩ型を用いた例を表６０（合同法）、表６９（メルセンヌツイスター）、表７８（準乱数 ソボル列）に示す。まず、目標とする分布の特性はξ_１（歪度１．０、尖度４．０）、ξ_２（歪度１．２５、尖度５．０）、ν（歪度０．７５、尖度３．０）で、それから生成される多変量非正規分布のモーメント特性の理論値は、Ｅ（ｘ_１ ^３）＝０．８５７８３、Ｅ（ｘ_２ ^３）＝０．６９２１１、Ｅ（ｘ_１ ^２ｘ_２）＝０．２７４８８、Ｅ（ｘ_１ｘ_２ ^２）＝０．２４１７４、Ｅ（ｘ_１ ^４）＝１４．２３４９０、Ｅ（ｘ_２ ^４）＝１１．６３７１０、Ｅ（ｘ_１ ^３ｘ_２）＝５．９４８８１、Ｅ（ｘ_１ ^２ｘ_２ ^２）＝５．０６１７３、Ｅ（ｘ_１ｘ_２ ^３）＝５．２９９３９である。

例えば、１０，０００，０００個のシミュレーションの場合、表６０（合同法）では、それぞれ、ξ_１（歪度１．０２８０２、尖度４．０４７６０）、ξ_２（歪度１．２９３２３、尖度５．１５２６４）、ν（歪度０．７６２０８、尖度２．９９８７１）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．９０４９１（０．８５７８３）、Ｅ（Ｘ_２ ^３）＝０．７２６３４（０．６９２１１）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．２９２４４（０．２７４８８）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．２５６４２（０．２４１７４）、Ｅ（Ｘ_１ ^４）＝１４．５７５８７（１４．２３４９０）、Ｅ（Ｘ_２ ^４）＝１１．８３１０５（１１．６３７１０）、Ｅ（Ｘ_１ ^３Ｘ_２）＝６．０７５６４（５．９４８８１）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝５．１２７２８（５．０６１７３）、Ｅ（Ｘ_１Ｘ_２ ^３）＝５．３７７８４（５．２９９３９）である。

括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ_１、ξ_２、νの実際のシミュレーション時の歪度、尖度（例えば、上の例ではξ_１（歪度１．０２８０２、尖度４．０４７６０）、ξ_２（歪度１．２９３２３、尖度５．１５２６４）、ν（歪度０．７６２０８、尖度２．９９８７１））から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との間に開きがある。

ところが、これをメルセンヌツイスターでシミュレーションすると、表６９（メルセンヌツイスター）にあるように、ξ_１（歪度１．０００３９、尖度４．００４０２）、ξ_２（歪度１．２４９１４、尖度４．９９５６７）、ν（歪度０．７５００９、尖度３．０００６１）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．８６２１５（０．８５７８３）、Ｅ（Ｘ_２ ^３）＝０．６９７９１（０．６９２１１）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．２７６９６（０．２７４８８）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．２４３８１（０．２４１７４）、Ｅ（Ｘ_１ ^４）＝１４．２８１２８（１４．２３４９０）、Ｅ（Ｘ_２ ^４）＝１１．６６９１０（１１．６３７１０）、Ｅ（Ｘ_１ ^３Ｘ_２）＝５．９３１９９（５．９４８８１）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝５．０５６０２（５．０６１７３）、Ｅ（Ｘ_１Ｘ_２ ^３）＝５．２９９１６（５．２９９３９）である。この結果からすべての目標のモーメント特性の理論値と非常に近くなっていることがわかる。

また、準乱数（ソボル列）では表７８（ソボル列）にあるように、ξ_１（歪度０．９９９７３、尖度３．９９８７８）、ξ_２（歪度１．２５０８２、尖度５．００６４０）、ν（歪度０．７４９５４、尖度２．９９７６６）で、それから生成される多変量非正規分布のモーメント特性の計測値は、Ｅ（Ｘ_１ ^３）＝０．８５１５１ （０．８５７８３）、Ｅ（Ｘ_２ ^３）＝０．６９４２５（０．６９２１１）、Ｅ（Ｘ_１ ^２Ｘ_２）＝０．２７２７１（０．２７４８８）、Ｅ（Ｘ_１Ｘ_２ ^２）＝０．２４３５８（０．２４１７４）、Ｅ（Ｘ_１ ^４）＝１４．１２４５４（１４．２３４９０）、Ｅ（Ｘ_２ ^４）＝１１．６３９４０（１１．６３７１０）、Ｅ（Ｘ_１ ^３Ｘ_２）＝５．８７２８１（５．９４８８１）、Ｅ（Ｘ_１ ^２Ｘ_２ ^２）＝５．０２３８１（５．０６１７３）、Ｅ（Ｘ_１Ｘ_２ ^３）＝５．２８３３５（５．２９９３９）である。この結果からすべての目標のモーメント特性の理論値と非常に近くなっていることがわかる。

２．３ 応用のためのあてはめプロシージャ
統計モデリングの見地から、経験分布の１つに３次と４次のモーメント行列をあてはめるプロシージャが必要である。各要素と標本のモーメントの式の間の偏りの２乗の和を最小にすることによって、Ｔ，ζ，γ，κ及びβを計算するプロシージャを提案する。すなわち、ＴＴ´＝Σという条件の下、次の式（９３）によって定義される関数の値を最小化することによって、適切なＴ，ζ，γ，κ及びβを得ることができる。

ここで、ｍ＝ｐのときＴ＝Σ^１／２（行列の平方根）であり、固定されている。行列の平方根の計算は、実対称行列の固有値問題に帰着され、おもな方法として、ヤコビ法、ハウスホルダー法、二分法、逆反復法、二分法・同時逆反復法、キューエル（ＱＬ）法、ディバイド・アンド・コンクウェア（Divide and Conquer）法がある。

ただし、ｆ_iは、３次と４次のモーメントの各要素の式であり、例えば、ｐ＝２について、ｆ₁＝Ｅ（ｘ₁ ³）＝γｔ₁₁ ³ζ₁＋γｔ₁₂ ³ζ₂，ｋ＝９であり、ｍ_iは、標本分布のモーメントである。ｐ＝ｍ＝２の場合、Ｔは固定されている（ユアン及びベントラー（２０００））。例えば、標本モーメントｍ_iは、表８１に従う。式（９３）の値を最小にすることによって、ζ₁＝０．６５９６７，ζ₂＝１．０３６２４，γ＝０．６３５５３，κ₁＝５．２７４７５，κ₂＝６．６９１７６，β＝４．０２６８０というパラメータを得る。これらのパラメータによって得られた行列要素は、表８１に示されている。ここで、ｗ_iは、重みである。

３ 最尤法
前述のモーメント法にかわる方法として最尤法がある。モーメント法では、３次と４次モーメントのような高次モーメントについては不安定になることがあるが、最尤法ではこの点が改善される可能性がある。

最尤法では、生成した乱数によって密度関数を数値的に近似する。例として、特定の多変量非正規分布を１０，０００，０００点の乱数により非常に正確にシミュレーションする。最尤法のプロシージャは、以下の通りである。

（１）最尤法に用いるパラメータの初期値をモーメント法から採用する。

（２）多変量非正規分布から例えば１０，０００，０００点の乱数を得る。

（３）各次元をステップΔｈによって分割に分ける。超立方体となる各分割内にある点の個数を数え、全ての点の個数で各分割内にある点の個数を除する。

（４）次のように対数尤度を求める。

｛Ｘ₁，・・・，Ｘ_N｝は、ｎ次元データの集合（Ｘ_ｉは、ｎ次元データ・ベクトルである。）を示す。特定の多変量非正規分布のパラメータθを有するＸ_iの分割の確率Ｐ_i（θ）は、次の式（９４）によって定義される。ここで、パラメータθは、ユアンとベントラーの方法Ｉで用いられるパラメータＴ，ζ＝（ζ₁，・・・，ζ_ｍ），γ，κ＝（κ₁，・・・，κ_ｍ）、βの少なくとも１つを表すものとする。

したがって、パラメータθを有する近似した尤度ｆ_i（θ）は、次の式（９５）によって定義される。

対数尤度の和は、次の式（９６）によって計算される。

そして、この対数尤度の和を最大にすることにより、最尤推定量

を得る。このようにして、最尤推定パラメータが得られた。なお、尤度の積を最大にしても同様である。

多変量非正規分布推定量

を、対数尤度の和を最大化することにより推定する。この方法を予想損失額に適用するため、１年の取引日に対応する２５０点を採用する。２５０点のモーメントを表８２に示し、最尤法による推定分布のモーメントを表８３〜８５に示す。

２５０点は、合同法を用いて得られた、目的とする分布の特性が、ξ_１（歪度０．７５、尖度５．０）、ξ_２（歪度１．００、尖度５．５）、ν（歪度０．５、尖度４．０）を持つ（先のピアソンＩＶ型の多変量シミュレーションで使われた目標値）乱数である。２５０点からなる分布の特性の計測値は表８２に示されており、目的とする分布の特性とはかなり違ったものになっている。

次に、多数の乱数からなる（ここでは、１０，０００，０００点）分布による最尤法の推定結果を示す。なお、分散共分散のみを経験分散共分散０．４１８１９と同じに設定するため、ｔ₁₁＝０．２１４８，ｔ₁₂＝０．９７６７に設定する。まず、表８３は合同法による推定結果であるが、２５０点の計測値よりも、目標とされた分布の特性値（表５４）に近づいていることがわかる。このことから、少数のデータの分布の推定には、モーメント法より、最尤法がすぐれていることがわかる。なお、最大対数尤度は、−５６２．７３になった。特に、金融のリスク管理などでは、一年のデータは２５０点程度であり，これから、多変量非正規分布を用い、良好な推定をするためには、最尤法を用いらなければならないといえよう。

さらに、表８４（メルセンヌツイスター）、表８５（準乱数 ソボル列）を示す。前の多変量非正規分布に従う乱数のシミュレーション結果から、これらの方法は合同法に比べ、非常に目標とする分布に近い乱数を発生することがわかっており、これは最尤法を実行するための重要な性質である。結果として、両者の方法も、より目標値に近い値を推定している。なお、メルセンヌツイスターについての最大対数尤度は−５６１．３６、ソボル列についての最大対数尤度は−５６０．３６となった。

本発明では、最尤法について主に連続型確率変数の場合の密度関数を乱数から直接計算している方法を中心にしているが、離散型確率変数の場合は密度関数の代わりに確率関数となる。また、密度関数や確率関数を積分した分布関数として表現する方法も可能である。

本実施の形態は、一般性を失わないように標準化された平均０、分散１の設定（Ｅ［ｘ_ｉ］＝０，Ｖａｒ［ｘ_ｉ］＝σ_ｉｉ＝１（Ｅ［ｘ_ｉ ^２］＝σ_ｉｉ＝１），１≦ｉ≦ｎ）で行われているが、一般の場合でも標準化の逆変換を行えば容易に拡張でき、また、最尤法を使えば、さらに、その標準化のための平均と標準偏差もパラメータとして、直接、生のデータ（標準化などの前処理をしていないという意味）から分布を推定することが容易に可能である。

なお、さらに、分布の１次から４次までの一部または全部のモーメントが解析的に容易にもとまらない場合や、発散などで存在しない場合、あるいは密度関数でなく特性関数で定義される分布などの場合でも、その分布の乱数を発生させてそれからモーメントを計算する方法（必ず計算可能）や、最尤法を行うにあたって分布を特徴付けるパラメータやその他のパラメータ（両パラメータを含めてθと表す。）を適当な間隔で適当に与え、その分布の乱数を用いて幅広く尤度を計算する方法を使えば、そのような分布もこの発明のνやξの分布として採用でき
る。そのような分布として例えば、一般の場合、２次のモーメントが発散し、密度関数が陽表的に書けないが特性関数によって表現できる安定分布などがある。ゆえに、本発明の最尤法は、νやξの分布としては、乱数が発生できれば１次から４次までの一部又は全部のモーメントが求まらない分布でも採用可能となり、ユアン・ベントラーの方法からさらに拡張された方法となっている。

４ 実施例
図１は、２次元を例とした予想損失額計算を示すフローチャートである。ここでは、日本株式リターンと外国株式リターンを２次元の変数に取る。この予想損失額計算の一連の手順は、コンピュータにおいて所定のコンピュータ・プログラムを実行することによって実現される。

ステップＳ１においては、日本株式リターンＸ₁´と外国株式リターンＸ₂´を標準化する。

標準化は、次のように経験分布の平均を０、分散を１とすることにより行う。日本株式リターンと外国株式リターンのある実現値を持つ確率変数Ｘ₁´及びＸ₂´について、それぞれの平均をＥ（Ｘ₁´）及びＥ（Ｘ₂´）、分散をＶ（Ｘ₁´）及びＶ（Ｘ₂´）とする。このとき、標準化日本株式リターンＸ₁と標準化外国株式リターンＸ₂を、次の式（９７）と（９８）により与える。

ステップＳ２においては、標準化日本株式リターンＸ₁と標準化外国株式リターンＸ₂から分散共分散行列Σを計算する。

ステップＳ３においては、標準化日本株式リターンＸ₁と標準化外国株式リターンＸ₂について３次モーメントを計算する。３次モーメントは、次の式（９９）に示すように、ｍ₁〜ｍ₄までの４種類がある。具体的には、Ｅ（・）は、Ｘ_１，Ｘ_２のデータ（実現値）から計算される。

ステップＳ４においては、標準化日本株式リターンＸ₁と標準化外国株式リターンＸ₂について４次モーメントを計算する。４次モーメントは、次の式（１００）に示すように、ｍ₅〜ｍ₉までの５種類がある。

ステップＳ５においては、ステップＳ２において計算した分散共分散行列Σについて、Σ＝ＴＴ´を満たす行列Ｔ＝（ｔ_ij）を計算する。いまξの次元とＸ_ｉの次元がともに２次元なので、行列Ｔは、Ｔ＝Σ^１／２（行列の平方根）と固定される。

具体的には、次の式（１０１）の値を最小にするパラメータζ＝（ζ₁，ζ₂），γ，κ＝（κ₁，κ₂），βを決定する。

ここで、重みｗ_iによって、モーメントに軽重をつけることができる。本実施例では、簡単のためにｗ_i＝１とする。

式（１０１）においては、３次及び４次モーメントの表現として、次の式（１０２）と（１０３）の表現を用いる。行列Ｔは、Ｔ＝Σ^１／２（行列の平方根）と計算されているので、式（１０１）に最小値を与えるパラメータζ＝（ζ₁，ζ₂），γ，κ＝（κ₁，κ₂），βを求める。

［３次モーメント］

ただし、ｆ₁，ｆ₂，ｆ₃，ｆ₄は、Ｘ＝νＴξで発生される確率ベクトルＸ＝（ｘ_１，ｘ_２）´のそれぞれ３次モーメントＥ（ｘ₁ ³），Ｅ（ｘ₂ ³），Ｅ（ｘ₁ ²ｘ₂），Ｅ（ｘ₁ｘ₂ ²）の表現である。

［４次モーメント］

ただし、ｆ₅，ｆ₆，ｆ₇，ｆ₈，ｆ₉は、Ｘ＝νＴξで発生される確率ベクトルＸ＝（ｘ_１，ｘ_２）´のそれぞれ４次モーメントＥ（ｘ₁ ⁴），Ｅ（ｘ₂ ⁴），Ｅ（ｘ₁ ³ｘ₂），Ｅ（ｘ₁ ²ｘ₂ ²），Ｅ（ｘ₁ｘ₂ ³）の表現である。ここで、σ_１１＝σ_２２＝１であり、また、ｔ_１１ ^２＝ｔ_２２ ^２＝１−ｔ_１２ ^２＝１−ｔ_２１ ^２である。

ステップＳ６においては、確率変数ξ₁, ξ₂，νがピアソン分布のどの型に属するか決定する。

確率変数ξ₁は、Ｅ（ξ₁）＝０，Ｅ（ξ₁ ²）＝１，Ｅ（ξ₁ ³）＝ζ₁，Ｅ（ξ₁ ⁴）＝κ₁を満たす。同様に、確率変数ξ₂は、Ｅ（ξ₂）＝０，Ｅ（ξ₂ ²）＝１，Ｅ（ξ₂ ³）＝ζ₂，Ｅ（ξ₂ ⁴）＝κ₂を満たす。確率変数νは、Ｅ（ν）＝０，Ｅ（ν^２）＝１，Ｅ（ν³）＝γ，Ｅ（ν⁴）＝βを満たす。規格化された確率変数ξ₁, ξ₂，νの歪度はそれぞれζ₁，ζ₂，γであり、尖度はそれぞれκ₁，κ₂，βである。

ここで、次の式（１０４）と表８６を用い、確率変数ξ₁，ξ₂，νがピアソン分布のどの型になるか決定する。式（１０４）におけるβ₁とβ₂は、それぞれ歪度の２乗と尖度である。

ステップＳ７においては、確率変数ξ₁，ξ₂，νのそれぞれについて、ステップＳ６で求めた型のピアソン分布に従う乱数を発生する。

具体的に、ピアソンＩＶ型以外の場合には、次の表８７と表８８に示すプロシージャによって乱数を発生する。

ピアソンＩＶ型の場合は、次のアルゴリズムによって乱数を発生する。

このアルゴリズムは、ピアソンＩＶ型分布の正規化定数を解析解又はこの展開（式（１０５））に基づいて計算するステップと、棄却法により乱数を発生するステップを含む。

ステップＳ８においては、ステップＳ７において発生した乱数ξ₁，ξ₂，νを次の式（１０６）によって、標準化日本株式リターン従う乱数ｘ₁と標準化外国株式リターンに従う乱数ｘ₂に変換する。

ステップＳ９においては、ステップＳ８で得られた乱数ｘ₁，ｘ₂を次の式（１０７）と（１０８）によって、標準化前の日本株式リターンに従う乱数ｘ₁´と標準化前の外国株式リターンに従う乱数ｘ₂´に変換する。ここで、分散Ｖ（・）と平均Ｅ（・）には、ステップＳ１で標準化した際の値を用いる。

ステップＳ１０においては、ステップＳ９で得られた乱数ｘ₁´と乱数ｘ₂´を用いて、予想損失額（ＶＡＲ）を計算する。

日本株式と外国株式の資産の残高をＡ₁，Ａ₂とすると、これらの重みを用いて日本株式と外国株式を併せたリターンに従う乱数を得る。次の式（１０９）と（１１０）にそれぞれ日本株式と外国株式の資産の重みを示し、日本株式と外国株式を併せたリターンに従う乱数Ｗを式（１１１）に示す。

乱数Ｗを用いたモンテカルロ法によるシミュレーションによって、日本株式と外国株式を併せたリターン分布をマイナス無限大から積分した領域が、全体の領域の１％になったところのリターンを用い、このリターン（通常、マイナスである。）の絶対値に日本株式と外国株式を併せた資産の残高を乗じたものが標準的な予想損失額である。

具体的には、日本株式と外国株式を併せたリターン分布に従う乱数Ｗを複数発生させ、低いものから並べて下から１％目（通常、銀行規制などで使われている基準）のところのリターンを用いる。なお、予想損失額に用いるリターンは１％に限らず、所望の値を用いることができる。

本実施例においては、日本株式と外国株式の２種類の資産を有する２次元の場合について説明したが、ｎ種類の資産を有する一般のｎ次元の場合には、容易に拡張することができる。以上の方法は、損失額自体の分布によるリスク管理にも適用できる。

この場合、ステップＳ５においては、式（１１２）において、ｆ_ijk（Ｔ，ζ，γ，κ，β）とｆ_ijkl（Ｔ，ζ，γ，κ，β）には、３次と４次モーメントの表現として式（１１３）と（１１４）の関係を用いる。また、ｍ_ijkとｍ_ijklは、式（１１５）と（１１６）によって定義される。具体的には、Ｅ（・）は、Ｘ_ｉのデータ（実現値）から計算される。

そして、式（１１２）を最小にするようなパラメータζ＝（ζ₁，・・・，ζ_ｍ），γ，κ＝（κ₁，・・・，κ_ｍ），βを求める。

ステップＳ６においては、これらのパラメータから、確率変数ν，ξ₁〜ξ_ｍがどの型のピアソン分布か決定し、ステップＳ７においてその型の乱数を発生する。ステップＳ８においては、ステップＳ７において発生した乱数を乱数ｘ₁，ｘ₂・・・，ｘ_nに変換し、さらにステップＳ９においては、乱数ｘ₁´，ｘ₂´・・・，ｘ_n´に変換する。

ステップＳ１０においては、ｎ種類の資産の残高をＡ₁，Ａ₂，・・・，Ａ_nとすると、これらの重みを乗じたｎ種類の資産を併せたリターンに従う乱数を得る。次の式（１１７）にｎ種類の資産の重みを示し、式（１１８）にｎ種類の資産を併せたリターン分布に従う乱数Ｗを示す。

このような乱数Ｗを用いて、ｎ種類の資産を有する場合にも２次元の場合と同様にモンテカルロ法によるシミュレーションによって予想損失額を計算することができる。すなわち、複数の乱数Ｗを発生させ、低いものから並べて下から１％目のところのリターンを取り出し、このリターンの絶対値にｎ種類の資産の残高を乗じたものを予想損失額とする。

さらに、本実施例においては、モーメント法に代わって最尤法によって推定したパラメータを用いて予想損失額を計算することができる。前述したように、株式リターンを用いる場合、年間取引日は２５０日と少数である。最尤法によると、例えば２５０点のような少数の標本点であっても、パラメータを正確に算出することができる。

図２は、最尤法を用いて予想損失額を計算する一連の手順を示すフローチャートである。この一連の手順においては、ステップＳ５とＳ６の間にステップＳ１０１が加わった他は、図１に示した一連の手順と同様であるので、ステップＳ１０１を除いて説明を省略するものとする。

ステップＳ１０１において、ステップＳ５においてモーメント法により計算したパラメータを初期値として、最尤法によりパラメータζ₁，ζ₂，γ，κ₁，κ₂，βの推定値を計算する。

具体的に、ステップＳ１０１においては、次の手順によりパラメータの推定値を計算する。なお、以下このステップＳ１０１においては、パラメータＴ，ζ₁，ζ₂，γ，κ₁，κ₂，βの少なくとも１つをθで表記することにする。

（１）ステップＳ５において、モーメント法から計算したパラメータを初期値として採用する。

（２）多変量非正規分布から例えば１０，０００，０００点の乱数を得る。乱数の発生は、ステップＳ６とＳ７と同様にして行う。

（３）各次元をステップΔｈによって分割に分ける。超立方体となる各分割内にある点の個数を数え、全ての点の個数で各分割内にある点の個数を除する。

（４）次のように対数尤度を求める。

｛Ｘ₁，・・・，Ｘ_N｝は、ｎ次元データの集合（Ｘ_ｉは、ｎ次元データ・ベクトルである。）を示す。特定の多変量非正規分布のパラメータθを有するＸ_iの分割の確率Ｐ_i（θ）は、次の式（１１９）によって定義される。ここで、パラメータθは、ユアンとベントラーの方法Ｉで用いられるパラメータＴ，ζ＝（ζ₁，・・・，ζ_ｍ），γ，κ＝（κ₁，・・・，κ_ｍ）、βの少なくとも１つを表すものとする。

したがって、パラメータθを有する近似した尤度ｆ_i（θ）は、次の式（１２０）によって定義される。

対数尤度の和は、次の式（１２１）によって計算される。

そして、この対数尤度の和を最大にすることにより、最尤推定量

を得る。この最尤推定量をパラメータとして用いる。なお、尤度の積を最大にしても同様である。

本実施の形態によると、従来のリスクをより正確に測定出来るだけでなく、不動産やオルタネティブ投資などの歪んだリターン分布を含んだポートフォリオなどのリスク管理も可能になる。また、従来、リスク管理で使われていた、リターン（収益率）ベースの方法だけでなく、損失額自体の分布を直接評価できるようになり、現物だけでなく、オプションなど様々なデリバティブを含むポートフォリオのリスク管理や、保険・天候・地震などのその他の様々なリスク（金融リスク、オペレーショナル・リスク、保険リスク、災害リスク、事業リスクなど）も統合したリスク管理が可能となる。ここでは、リスク計測対象の損失額自体の分布や収益率の分布など、リスクを計測する際の評価のための分布を総称して損益分布と呼ぶ。

本実施の形態においては、ステップＳ２、Ｓ３及びＳ４で求めた３次及び４次モーメントｍ_１〜ｍ_９を用いてステップＳ５でパラメータζ，γ，κ，βを求め、そのパラメータを初期値として、ステップＳ１０１において前記初期値の近傍に属する最大対数尤度を与えるパラメータを求めた。しかしながら、パラメータの初期値を定めるステップＳ２〜Ｓ５の一部または全部を省略することもできる。この場合、パラメータの初期値は未定として、ステップＳ１０１において、パラメータ空間における離散的な格子点について幅広い領域で対数尤度を計算し、最大対数尤度（又は最大尤度）を与えるパラメータを求める。このような方法は、格子探索法（グリッドサーチ）と称される。

この方法によれば、１次から４次までのモーメントの一部又は全部が存在しない確率分布であっても、乱数さえ発生できればξやνの分布として採用可能となる。そのような場合、最尤法のための、分布を特徴付けるようなパラメータやその他あらゆる必要なパラメータを広くθと記す。

５ 変形例
５．１ 金融商品のシミュレーション方法
次に、本実施の形態の第１の変形例として、本実施の形態の乱数を使った多変量非正規分布の特性を用いた金融商品のシミュレーション方法を説明する。

この方法では、次のアルゴリズムによって発生されるピアソンＩＶ型分布に従う乱数を使ったモンテカルロ法によるシミュレーションを用いることができる。このアルゴリズムは、ピアソンＩＶ型分布の正規化定数の解析解又はこの展開（式（１２２））を計算するステップと、棄却法により乱数を発生するステップとを含む。

この方法は、多資産の高次の相関を考慮した乱数を発生させ、実際のマーケットの特性を正確に表すリターンを発生することができる。例えば、金融商品のパフォーマンス測定などに用いられる。

この方法は、前記多変量非正規分布（ξやνに、ピアソン分布系を使うことに限らず一般的な確率分布を使うことを含む。）を用いることにより、金融の分野における分布に特徴的ないわゆるファット・テールを正確に再現することができる。また、様々なデリバティブの対象データを本実施の形態の方法により、正確な確率分布（ξやνに、ピアソン分布系を使うことに限らず一般的な確率分布を使うことを含む。）をあてはめ（単変量分布の場合も、多変量の特殊な場合として同様に扱える。）、その分布に従うモンテカルロ・シミュレーションをすることで価格が計算できる。デリバティブの設計には、多地点での気象データや様々な指数の相関が必要で、その分析にも本実施の形態が応用される。

例えば、天候・保険デリバティブやリアル・オプションのプライシングや設計などに応用できる。特に、プライシングに多変量非正規分布を使う例としては、コリレーティング・トリガー・オプション（複数の気象要素を選んで、それぞれに補償金の支払基準を設定するもので、オプション購入者は、契約終了時にその中で一番有利なものを選択・行使できるオプション。）のような複合オプションや複数地点オプションなどがある。

５．２ 半導体へのイオン注入
次に、本実施の形態の第２の変形例として、半導体へのイオン注入への応用を説明する。イオン注入への応用では、半導体に注入されるイオンの分布をモンテカルロ法によってシミュレーションする。

イオンの空間分布や電荷特性は、一般に多変量非正規分布となるが、ピアソン分布を含む前記多変量非正規分布に従う乱数を用いた方法により正確に再現することができる。

このように本実施の形態の多変量非正規分布を適用することにより、従来の解析的なモデルを、本発明の多変量非正規分布の乱数表現によって実際の分布により近い正確なモデルに改良できる。

また、現在、イオン注入後のイオン分布を知るための物理モデルのモンテカルロ・シミュレーションは、非常に長時間かかっているが、本実施の形態の多変量非正規分布の乱数とパラメータ推定法によって、より少数の物理シミュレーションの結果から短時間にその分布の全体を正確に推定可能になり、大幅に現在のモンテカルロ・シミュレーションの短時間化を行える。

５．３ アセット・アロケーション
本実施の形態の方法によって、多変量非正規分布に従う資産間のアセット・アロケーションを、従来の平均・分散による最適化だけでなく歪度も考慮できる最適化手法と組み合わせて実行する。この歪度まで考慮できる最適化手法については、参考文献に今野浩「理財工学ＩＩ」、日科技連、１９９８年、コンノ，エイチ．，スズキ，ティー.及びコバヤシ，ディー．（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）巻、第２部。

異なった能力検査や試験の結果に、本発明の実施によって、正確な多変量非正規分布が推定され、能力評価を正確に実施することができる。このことから、入試判定や人事アセスメント（選抜基準の設定や選抜有効性の検証）などをより正確化することができる。また、パーセンタイルなど様々な基準により、総合的な能力値を定義することができる。

さらに、大学入学試験の結果と入学者の卒業成績、就職試験の結果と入社後の評価など、一部の周辺分布が大量のデータであるとき、本実施の形態の方法により、残りの少数のデータの分布を最尤法により推定することもできる。

産業上の利用可能性
本発明によると、ピアソン分布系や一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布の推定方法、及びこれらの手順を格納した媒体を提供することができる。

具体的には、多変量非正規分布に従う乱数発生方法及び多変量非正規分布のパラメータ推定方法は、様々な分野における統計モデリングに有効である。具体的には、統計ソフトウェアのパッケージ、リスク・マネジメントにおける予想損失額の計算への応用、金融商品のシミュレーション、アセット・アロケーションへの応用、天候・保険デリバティブやリアル・オプションのプライシングと設計、半導体へのイオン注入のシミュレーション、教育や人事アセスメントにおける能力評価に有効である。

図１は、予想損失額を計算する一連の手順を示すフローチャートである。 図２は、最尤法を用いて予想損失額を計算する一連の手順を示すフローチャートである。

Claims (14)

1. コンピュータ上で多変量非正規分布に従う乱数を発生する乱数発生方法であって、
   経験データ｛Ｘ _１ ´，・・・，Ｘ _Ｎ ´｝、すなわち、
   

   

   をコンピュータに入力するステップと、
   コンピュータを用いて、前記経験データ｛Ｘ _１ ´，・・・，Ｘ _Ｎ ´｝を標準化して｛Ｘ _１ ，・・・，Ｘ _Ｎ ｝とするステップと、
   コンピュータを用いて、前記標準化されたデータ｛Ｘ _１ ，・・・，Ｘ _Ｎ ｝に基づいて、分散共分散行列Σを計算するステップと、
   コンピュータを用いて、前記標準化されたデータ｛Ｘ _１ ，・・・，Ｘ _Ｎ ｝に基づいて、３次モーメントｍ _ｉｊｋ （１≦ｉ≦ｊ≦ｋ≦ｎ）を計算するステップと、
   コンピュータを用いて、前記標準化されたデータ｛Ｘ _１ ，・・・，Ｘ _Ｎ ｝に基づいて、４次モーメントｍ _ｉｊｋｌ （１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）を計算するステップと、
   コンピュータを用いて、前記分散共分散行列ΣからＴ＝Σ ^１／２ （行列の平方根）を満たす行列Ｔを計算するステップと、
   前記標準化されたデータ｛Ｘ _１ ，・・・，Ｘ _Ｎ ｝が従う分布が、Ｘ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ＝ｎである。）個の確率変数ξ _１ ，・・・，ξ _ｍ からなる確率ベクトル（ξ _１ ，・・・，ξ _ｍ ）´である。独立な確率変数ξ _１ ，・・・，ξ _ｍ は、パラメータζ _ｊ 及びκ _ｊ について、Ｅ（ξ _ｊ ）＝０，Ｅ（ξ _ｊ ^２ ）＝１，Ｅ（ξ _ｊ ^３ ）＝ζｊ，Ｅ（ξ _ｊ ^４ ）＝κ _ｊ （１≦ｊ≦ｍ）を満たす。また、νはν＝１で、パラメータγ及びβについて、Ｅ（ν ^２ ）＝１，Ｅ（ν ^３ ）＝γ＝１，Ｅ（ν ^４ ）＝β＝１を満たす。）を満たす非正規分布を有すると仮定して、コンピュータを用いて、３次モーメントｍ _ｉｊｋ （１≦ｉ≦ｊ≦ｋ≦ｎ）及び４次モーメントｍ _ｉｊｋｌ （１≦ｉ≦ｊ≦ｋ≦ｌ≦ｎ）とｆ _ｉｊｋ （Ｔ，ζ，γ，κ，β），ｆ _ｉｊｋｌ （Ｔ，ζ，γ，κ，β）（ここでζ＝（ζ _１ ，・・・，ζ _ｍ ），κ＝（κ _１ ，・・・，κ _ｍ ）である。）との差を損失とする損失関数を導入し、全体として評価するリスク関数を最小化するようにパラメータＴ，（ζ _１ ，・・・，ζ _ｍ ），γ，（κ _１ ，・・・，κ _ｍ ），βを決定するステップと、
   コンピュータを用いて、一様乱数を発生するとともに、この一様乱数から、前記決定されたパラメータＴ，（ζ _１ ，・・・，ζ _ｍ ），γ，（κ _１ ，・・・，κ _ｍ ），βに基づく前記確率変数ξ _１ ，・・・，ξ _ｍ に従ったｍ個の系列の乱数を計算するステップと、
   コンピュータを用いて、前記確率変数ξ _１ ，・・・，ξ _ｍ に従ったｍ個の系列の乱数から、Ｘ＝νＴξに基づいて、標準化された多変量非正規分布に従った乱数を計算するステップと、
   コンピュータを用いて、前記標準化された多変量非正規分布に従った乱数から、前記標準化の逆変換により、標準化前の経験データ｛Ｘ _１ ´，・・・，Ｘ _Ｎ ´｝に対応する多変量非正規分布に従った乱数に変換するステップと、
   を有する
   多変量非正規分布に従う乱数発生方法。
   ただし、前記損失関数（Ｌ）及びリスク関数（Ｒ）は、式（１）及び（２）、式（３）及び（４）、又は式（５）及び（６）のいずれかの組で与えられるものとする。
   

   

   

   

   

   

   

   ただし、ｆ _ｉｊｋ （Ｔ，ζ，γ，κ，β）及びｆ _ｉｊｋｌ （Ｔ，ζ，γ，κ，β）は、３次モーメントＥ（ｘ _ｉ ｘ _ｊ ｘ _ｋ ）及び４次モーメントＥ（ｘ _ｉ ｘ _ｊ ｘ _ｋ ｘ _ｌ ）にそれぞれ対応する関係式（７）及び（８）による表現とする。また、ｗ _ｉｊｋ 及びｗ _ｉｊｋｌ は、所定の重みである。
   

   

   

   ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄ _ｎ は、ｎ次のデュプリケーション行列であり、Ｄ _ｎ ＋は、Ｄ _ｎ のムーア・ペンローズの一般化逆行列であり、
   

   

   は、クロネッカー積であり、Ｅ _ｉｉ は、ｅ _ｉ を第ｉ列単位としてｅ _ｉ ｅ _ｉ ′である。
2. 前記確率変数ξ _１ ，・・・，ξ _ｍ には、少なくとも２以上の型のピアソン分布系に従う乱数を使い、
   前記リスク関数を最小化するようにパラメータＴ，（ζ _１ ，・・・，ζ _ｍ ），γ，（κ _１ ，・・・，κ _ｍ ），βを決定するステップは、式（９）及び（１０）でそれぞれ与えられるｎ次元の経験分布（確率ベクトル（Ｘ _１ ，・・・，Ｘ _ｎ ）´とする。）の３次及び４次モーメントに関して、式（１１）の値を最小にするパラメータＴ，ζ＝（ζ _１ ，・・・，ζ _ｍ ），κ＝（κ _１ ，・・・，κ _ｍ ）の少なくとも１つを求める請求項１に記載の多変量非正規分布に従う乱数発生方法。
   

   

   

   

   ただし、ｆ _ｉｊｋ （Ｔ，ζ，γ，κ，β）及びｆ _ｉｊｋｌ （Ｔ，ζ，γ，κ，β）は、３次モーメントＥ（ｘ _ｉ ｘ _ｊ ｘ _ｋ ）及び４次モーメントＥ（ｘ _ｉ ｘ _ｊ ｘ _ｋ ｘ _ｌ ）にそれぞれ対応する関係式（１２）及び（１３）による表現とする。また、ｗ _ｉｊｋ 及びｗ _ｉｊｋｌ は、所定の重みである。
   

   

   
3. ｗ_ｉｊｋ＝１，ｗ_ｉｊｋｌ＝１である請求項２に記載の多変量非正規分布に従う乱数発生方法。
4. ｎ＝２であって、前記あてはめるステップは、式（１４）及び（１５）でそれぞれ与えられる２次元の経験分布（確率ベクトル（Ｘ_１，Ｘ_２）´とする。）の３次及び４次モーメントと、式（１６）及び（１７）で与えられる３次モーメントＥ（ｘ_ｉｘ_ｊｘ_ｋ）及び４次モーメントＥ（ｘ_ｉｘ_ｊｘ_ｋｘ_ｌ）にそれぞれ対応する表現とについて、式（１８）で与えられる値を最小にするようにパラメータＴ，ζ＝（ζ_１，ζ_２），κ＝（κ_１，κ_２）の少なくとも１つを求める請求項２又は３に記載の多変量非正規分布に従う乱数発生方法。
   

   

   

   

   

   
5. ｗ_ｉ＝１である請求項４に記載の多変量非正規分布に従う乱数発生方法。
6. パラメータＴ，ζ＝（ζ_１，・・・，ζ_ｍ），κ＝（κ_１，・・・，κ_ｍ）の少なくとも１つを最尤法により推定する請求項１乃至５のいずれか１項に記載の多変量非正規分布に従う乱数発生方法。
7. 式（１９）で与えられるκの値を用いて、表１に示すκの値と分布の型との対応関係の少なくとも１つを用いて前記確率変数の属する型を決定する請求項１乃至６のいずれか１項に記載の多変量非正規分布に従う乱数発生方法。
   

   

   

   ただし、β_１とβ_２は、それぞれ歪度の２乗と尖度であり、式（１９）及び表１におけるκは、前記パラメータκ＝（κ_１，・・・，κ_ｍ）とは異なる。
8. 経験データ｛Ｘ _１ ´，・・・，Ｘ _Ｎ ´｝、すなわち、
   

   

   をコンピュータに入力するステップと、
   コンピュータを用いて、前記経験データ｛Ｘ _１ ´，・・・，Ｘ _Ｎ ´｝を標準化して｛Ｘ _１ ，・・・，Ｘ _Ｎ ｝とするステップと、
   コンピュータを用いて、任意の分散共分散行列ΣからＴ＝Σ ^１／２ （行列の平方根）を満たす行列Ｔを計算するステップと、
   前記標準化されたデータ｛Ｘ _１ ，・・・，Ｘ _Ｎ ｝が従う分布が、Ｘ＝νＴξ（ここで、ξは、ｍ（ここで、ｍ＝ｎである。）個の確率変数ξ_１ ，・・・，ξ _ｍ からなる確率ベクトル（ξ_１，・・・，ξ_ｍ）´である。独立な確率変数ξ_１，・・・，ξ_ｍは、パラメータζ_ｊ及びκ_ｊについて、Ｅ（ξ_ｊ）＝０，Ｅ（ξ_ｊ ^２）＝１，Ｅ（ξ_ｊ ^３）＝ζ_ｊ，Ｅ（ξ_ｊ ^４）＝κ_ｊ（１≦ｊ≦ｍ）を満たす。また、νはν＝１で、パラメータγ及びβについて、Ｅ（ν^２）＝１，Ｅ（ν^３）＝γ＝１，Ｅ（ν^４）＝β＝１を満たす。）を満たす非正規分布を有すると仮定して、コンピュータを用いて、任意のパラメータ（ζ _１ ，・・・，ζ _ｍ ），γ，（κ _１ ，・・・，κ _ｍ ），βに基づいて、式（２０）及び表２により（β _１ とβ _２ は、それぞれ歪度の２乗と尖度である。）、前記確率変数ξ _１ ，・・・，ξ _ｍ が各々ピアソン分布のどの型に属するかを決定するステップと、
   

   

   

   コンピュータを用いて、一様乱数を発生するとともに、この一様乱数から、表３、表４、表５、表６により、ピアソン分布の型に応じて、前記確率変数ξ _１ ，・・・，ξ _ｍ に従ったｍ個の系列の乱数を計算するステップと、
   

   

   

   

   

   コンピュータを用いて、前記確率変数ξ _１ ，・・・，ξ _ｍ に従ったｍ個の系列の乱数から、Ｘ＝νＴξに基づいて、標準化された多変量非正規分布に従った乱数を計算するステップと、
   コンピュータを用いて、前記標準化された多変量非正規分布に従う確率ベクトルＸについての空間を、所定間隔Δｈで分割し、このｎ次元空間を（Δｈ）^ｎの超立方体へ分割するステップと、
   コンピュータを用いて、各データ・ベクトルの属する区間（Δｈ）^ｎ内に存在する乱数の個数を乱数の総数で除して、当該区間（Δｈ）^ｎの確率Ｐ_ｉ（θ）を計算するステップと、
   コンピュータを用いて、前記確率Ｐ_ｉ（θ）をｎ次元体積（Δｈ）^ｎで除して、前記データ・ベクトルの属する区間（Δｈ）^ｎの尤度ｆ_ｉ（θ）を計算するステップと、
   コンピュータを用いて、前記尤度ｆ_ｉ（θ）の積Π_ｉ＝１ ^Ｎｆ_ｉ（θ）又は対数尤度の和
   Σ_ｉ＝１ ^Ｎｌｏｇ ｆ_ｉ（θ）が最大になるようなパラメータθを計算するステップと、
   を有する多変量非正規分布のパラメータ推定方法。
   ただし、パラメータθは、パラメータＴ，ζ＝（ζ _１ ，・・・，ζ _ｍ ），γ，κ＝（κ _１ ，・・・，κ _ｍ ）、βの少なくとも１つを表す。
   ただし、式（２０）及び表２において、β _１ とβ _２ は、それぞれ歪度の２乗と尖度であり、κは、前記パラメータκ＝（κ _１ ，・・・κ _ｍ ）とは異なる。
   ただし、Ｅ（・）は、期待値であり（以下同じ）、ｖｅｃｈ（・）は、対称行列から重複しない行列要素を取り出したベクトルであり、Ｄ_ｎは、ｎ次のデュプリケーション行列であり、Ｄ_ｎ ^＋は、Ｄ_ｎのムーア・ペンローズの一般化逆行列であり、
   

   

   は、クロネッカー積であり、Ｅ_ｉｉは、ｅ_ｉを第ｉ列単位としてｅ_ｉｅ_ｉ ^’である。
   ここで、独立な確率変数ξ_１，・・・，ξ_ｍは、パラメータζ_ｊ及びκ_ｊについて、Ｅ（ξ_ｊ）＝０，Ｅ（ξ_ｊ ^２）＝１，Ｅ（ξ_ｊ ^３）＝ζ_ｊ，Ｅ（ξ_ｊ ^４）＝κ_ｊ（１≦ｊ≦ｍ）を満たす。
   νはν＝１で、パラメータγ及びβについて、Ｅ（ν^２）＝１，Ｅ（ν^３）＝γ＝１，Ｅ（ν^４）＝β＝１を満たす。
   階数ｎのｎ行ｍ列（ｍ≧ｎ）の非確率行列Ｔ＝（ｔ_ｉｊ）は、行列Σ＝（σ_ｉｊ）について、ＴＴ´＝Σを満たす。ただし、行列Ｔ´はＴの転置行列である。
   このとき、次の式（２１）で与えられる確率ベクトルＸ＝（ｘ_１，・・・，ｘ_ｎ）´は、Ｃｏｖ（Ｘ）＝Σを満たす。ただし、Ｃｏｖ（・）は、ベクトルの分散共分散行列であり、ξ＝（ξ_１，・・・，ξ_ｍ）´である。
   

   
9. 前記確率変数ξ_１，・・・，ξ_ｍには、ピアソン分布系に従う乱数を使う請求項８に記載の多変量非正規分布のパラメータ推定方法。
10. 前記確率変数ξ_１，・・・，ξ_ｍには、少なくとも２以上の型のピアソン分布を用いる請求項９に記載の多変量非正規分布のパラメータ推定方法。
11. 請求項１乃至請求項７のいずれか一に記載の乱数発生方法を実現する、コンピュータにより実行可能なコンピュータ・プログラム。
12. 請求項１１記載のコンピュータ・プログラムを記録した記録媒体。
13. 請求項８乃至請求項１０のいずれか一に記載の多変量非正規分布のパラメータ推定方法を実現する、コンピュータにより実行可能なコンピュータ・プログラム。
14. 請求項１３記載のコンピュータ・プログラムを記録した記録媒体。

Images