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

明細書 :非正規分布に従う乱数発生方法及びそのパラメータの推定方法

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4350062号 (P4350062)
公開番号 特開2005-259170 (P2005-259170A)
登録日 平成21年7月31日(2009.7.31)
発行日 平成21年10月21日(2009.10.21)
公開日 平成17年9月22日(2005.9.22)
発明の名称または考案の名称 非正規分布に従う乱数発生方法及びそのパラメータの推定方法
国際特許分類 G06F   7/58        (2006.01)
FI G06F 7/58 B
請求項の数または発明の数 14
全頁数 100
出願番号 特願2005-130854 (P2005-130854)
分割の表示 特願2003-082816 (P2003-082816)の分割、【原出願日】平成15年3月25日(2003.3.25)
出願日 平成17年4月28日(2005.4.28)
優先権出願番号 2002093355
2002093356
優先日 平成14年3月28日(2002.3.28)
平成14年3月28日(2002.3.28)
優先権主張国 日本国(JP)
日本国(JP)
審査請求日 平成17年4月28日(2005.4.28)
特許権者または実用新案権者 【識別番号】801000027
【氏名又は名称】学校法人明治大学
発明者または考案者 【氏名】永原 裕一
個別代理人の代理人 【識別番号】100083806、【弁理士】、【氏名又は名称】三好 秀和
【識別番号】100100712、【弁理士】、【氏名又は名称】岩▲崎▼ 幸邦
【識別番号】100100929、【弁理士】、【氏名又は名称】川又 澄雄
【識別番号】100095500、【弁理士】、【氏名又は名称】伊藤 正和
【識別番号】100101247、【弁理士】、【氏名又は名称】高橋 俊一
【識別番号】100098327、【弁理士】、【氏名又は名称】高松 俊雄
審査官 【審査官】高瀬 勤
参考文献・文献 特開2002-015265(JP,A)
特開平8-139049(JP,A)
特開2002-288437(JP,A)
国際公開第96/018144(WO,A1)
特開平10-270376(JP,A)
特開平10-125602(JP,A)
特開平11-110192(JP,A)
実開平05-036449(JP,U)
特開平08-139049(JP,A)
特開2001-051974(JP,A)
調査した分野 G06F 7/58
G06F 17/18
JSTPlus(JDreamII)
特許請求の範囲 【請求項1】
コンピュータ上で多変量非正規分布に従う乱数を発生する乱数発生方法であって
経験データ{X´,・・・,X´}、すなわち、
【数1】
JP0004350062B2_000195t.gif
をコンピュータに入力するステップと、
コンピュータを用いて、前記経験データ{X´,・・・,X´}を標準化して{X,・・・,X}とするステップと、
コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、分散共分散行列Σを計算するステップと、
コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、3次モーメントmijk(1≦i≦j≦k≦n)を計算するステップと、
コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、4次モーメントmijkl(1≦i≦j≦k≦l≦n)を計算するステップと、
コンピュータを用いて、前記分散共分散行列ΣからT=Σ1/2(行列の平方根)を満たす行列Tを計算するステップと、
前記標準化されたデータ{X,・・・,Xが従う分布が、X=νTξ(ここで、ξは、m(ここで、m=nである。)個の確率変数ξ,・・・,ξからなる確率ベクトル(ξ,・・・,ξ)´である。独立な確率変数ξ,・・・,ξは、パラメータζ及びκについて、E(ξ)=0,E(ξ)=1,E(ξ)=ζj,E(ξ)=κ(1≦j≦m)を満たす。また、νはν=1で、パラメータγ及びβについて、E(ν)=1,E(ν)=γ=1,E(ν)=β=1を満たす。)を満たす非正規分布を有すると仮定して、コンピュータを用いて、3次モーメントmijk(1≦i≦j≦k≦n)及び4次モーメントmijkl(1≦i≦j≦k≦l≦n)とfijk(T,ζ,γ,κ,β),fijkl(T,ζ,γ,κ,β)(ここでζ=(ζ,・・・,ζ),κ=(κ,・・・,κ)である。)との差を損失とする損失関数を導入し、全体として評価するリスク関数を最小化するようにパラメータT,(ζ,・・・,ζ),γ,(κ,・・・,κ),βを決定するステップと、
コンピュータを用いて、一様乱数を発生するとともに、この一様乱数から、前記決定されたパラメータT,(ζ,・・・,ζ),γ,(κ,・・・,κ),βに基づく前記確率変数ξ,・・・,ξに従ったm個の系列の乱数を計算するステップと、
コンピュータを用いて、前記確率変数ξ,・・・,ξに従ったm個の系列の乱数から、X=νTξに基づいて、標準化された多変量非正規分布に従った乱数を計算するステップと、
コンピュータを用いて、前記標準化された多変量非正規分布に従った乱数から、前記標準化の逆変換により、標準化前の経験データ{X´,・・・,X´}に対応する多変量非正規分布に従った乱数に変換するステップと、
を有する
多変量非正規分布に従う乱数発生方法。
ただし、前記損失関数(L)及びリスク関数(R)は、式(1)及び(2)、式(3)及び(4)、又は式(5)及び(6)のいずれかの組で与えられるものとする。
【数2】
JP0004350062B2_000196t.gif
【数3】
JP0004350062B2_000197t.gif
【数4】
JP0004350062B2_000198t.gif
【数5】
JP0004350062B2_000199t.gif
【数6】
JP0004350062B2_000200t.gif
【数7】
JP0004350062B2_000201t.gif
ただし、fijk(T,ζ,γ,κ,β)及びfijkl(T,ζ,γ,κ,β)は、3次モーメントE(x)及び4次モーメントE(x)にそれぞれ対応する関係式(7)及び(8)による表現とする。また、wijk及びwijklは、所定の重みである。
【数8】
JP0004350062B2_000202t.gif
【数9】
JP0004350062B2_000203t.gif
ただし、E(・)は、期待値であり(以下同じ)、vech(・)は、対称行列から重複しない行列要素を取り出したベクトルであり、Dは、n次のデュプリケーション行列であり、D+は、Dのムーア・ペンローズの一般化逆行列であり、
【数10】
JP0004350062B2_000204t.gif
は、クロネッカー積であり、Eiiは、eを第i列単位としてe′である。
【請求項2】
前記確率変数ξ,・・・,ξには、少なくとも2以上の型のピアソン分布系に従う乱数を使い、
前記リスク関数を最小化するようにパラメータT,(ζ,・・・,ζ),γ,(κ,・・・,κ),βを決定するステップは、式(9)及び(10)でそれぞれ与えられるn次元の経験分布(確率ベクトル(X,・・・,X)´とする。)の3次及び4次モーメントに関して、式(11)の値を最小にするパラメータT,ζ=(ζ,・・・,ζ),κ=(κ,・・・,κ)の少なくとも1つを求める請求項1に記載の多変量非正規分布に従う乱数発生方法。
【数11】
JP0004350062B2_000205t.gif
【数12】
JP0004350062B2_000206t.gif
【数13】
JP0004350062B2_000207t.gif
ただし、fijk(T,ζ,γ,κ,β)及びfijkl(T,ζ,γ,κ,β)は、3次モーメントE(x)及び4次モーメントE(x)にそれぞれ対応する関係式(12)及び(13)による表現とする。また、wijk及びwijklは、所定の重みである。
【数14】
JP0004350062B2_000208t.gif
【数15】
JP0004350062B2_000209t.gif

【請求項3】
ijk=1,wijkl=1である請求項に記載の多変量非正規分布に従う乱数発生方法。
【請求項4】
n=2であって、前記あてはめるステップは、式(14)及び(15)でそれぞれ与えられる2次元の経験分布(確率ベクトル(X,X)´とする。)の3次及び4次モーメントと、式(16)及び(17)で与えられる3次モーメントE(x)及び4次モーメントE(x)にそれぞれ対応する表現とについて、式(18)で与えられる値を最小にするようにパラメータT,ζ=(ζ,ζ),κ=(κ,κ)の少なくとも1つを求める請求項又はに記載の多変量非正規分布に従う乱数発生方法。
【数16】
JP0004350062B2_000210t.gif
【数17】
JP0004350062B2_000211t.gif
【数18】
JP0004350062B2_000212t.gif
【数19】
JP0004350062B2_000213t.gif
【数20】
JP0004350062B2_000214t.gif

【請求項5】
=1である請求項に記載の多変量非正規分布に従う乱数発生方法。
【請求項6】
パラメータT,ζ=(ζ,・・・,ζ),κ=(κ,・・・,κ)の少なくとも1つを最尤法により推定する請求項1乃至のいずれか1項に記載の多変量非正規分布に従う乱数発生方法。
【請求項7】
式(19)で与えられるκの値を用いて、表1に示すκの値と分布の型との対応関係の少なくとも1つを用いて前記確率変数の属する型を決定する請求項1乃至のいずれか1項に記載の多変量非正規分布に従う乱数発生方法。
【数21】
JP0004350062B2_000215t.gif
【表1】
JP0004350062B2_000216t.gif
ただし、βとβは、それぞれ歪度の2乗と尖度であり、式(19)及び表1におけるκは、前記パラメータκ=(κ,・・・,κ)とは異なる。
【請求項8】
経験データ{X´,・・・,X´}、すなわち、
【数1】
JP0004350062B2_000217t.gif
をコンピュータに入力するステップと、
コンピュータを用いて、前記経験データ{X´,・・・,X´}を標準化して{X,・・・,X}とするステップと、
コンピュータを用いて、任意の分散共分散行列ΣからT=Σ1/2(行列の平方根)を満たす行列Tを計算するステップと、
前記標準化されたデータ{X,・・・,X}が従う分布が、X=νTξ(ここで、ξは、m(ここで、m=nである。)個の確率変数ξ,・・・,ξからなる確率ベクトル(ξ,・・・,ξ)´である。独立な確率変数ξ,・・・,ξは、パラメータζ及びκについて、E(ξ)=0,E(ξ)=1,E(ξ)=ζ,E(ξ)=κ(1≦j≦m)を満たす。また、νはν=1で、パラメータγ及びβについて、E(ν)=1,E(ν)=γ=1,E(ν)=β=1を満たす。)を満たす非正規分布を有すると仮定して、コンピュータを用いて、任意のパラメータ(ζ,・・・,ζ),γ,(κ,・・・,κ),βに基づいて、式(20)及び表2により(βとβは、それぞれ歪度の2乗と尖度である。)、前記確率変数ξ,・・・,ξが各々ピアソン分布のどの型に属するかを決定するステップと、
【数23】
JP0004350062B2_000218t.gif
【表2】
JP0004350062B2_000219t.gif
コンピュータを用いて、一様乱数を発生するとともに、この一様乱数から、表3、表4、表5、表6により、ピアソン分布の型に応じて、前記確率変数ξ,・・・,ξに従ったm個の系列の乱数を計算するステップと、
【表3】
JP0004350062B2_000220t.gif
【表4】
JP0004350062B2_000221t.gif
【表5】
JP0004350062B2_000222t.gif
【表6】
JP0004350062B2_000223t.gif
コンピュータを用いて、前記確率変数ξ,・・・,ξに従ったm個の系列の乱数から、X=νTξに基づいて、標準化された多変量非正規分布に従った乱数を計算するステップと、
コンピュータを用いて、前記標準化された多変量非正規分布に従う確率ベクトルXについての空間を、所定間隔Δhで分割し、このn次元空間を(Δh)の超立方体へ分割するステップと、
コンピュータを用いて、各データ・ベクトルの属する区間(Δh)内に存在する乱数の個数を乱数の総数で除して、当該区間(Δh)の確率P(θ)を計算するステップと、
コンピュータを用いて、前記確率P(θ)をn次元体積(Δh)で除して、前記データ・ベクトルの属する区間(Δh)の尤度f(θ)を計算するステップと、
コンピュータを用いて、前記尤度f(θ)の積Πi=1(θ)又は対数尤度の和
Σi=1log f(θ)が最大になるようなパラメータθを計算するステップと、
を有する多変量非正規分布のパラメータ推定方法。
ただし、パラメータθは、パラメータT,ζ=(ζ,・・・,ζ),γ,κ=(κ,・・・,κ)、βの少なくとも1つを表す。
ただし、式(20)及び表2において、βとβは、それぞれ歪度の2乗と尖度であり、κは、前記パラメータκ=(κ,・・・κ)とは異なる。
ただし、E(・)は、期待値であり(以下同じ)、vech(・)は、対称行列から重複しない行列要素を取り出したベクトルであり、Dは、n次のデュプリケーション行列であり、Dは、Dのムーア・ペンローズの一般化逆行列であり、
【数24】
JP0004350062B2_000224t.gif
は、クロネッカー積であり、Eiiは、eを第i列単位としてeである。
ここで、独立な確率変数ξ,・・・,ξは、パラメータζ及びκについて、E(ξ)=0,E(ξ)=1,E(ξ)=ζ,E(ξ)=κ(1≦j≦m)を満たす。
νはν=1で、パラメータγ及びβについて、E(ν)=1,E(ν)=γ=1,E(ν)=β=1を満たす。
階数nのn行m列(m≧n)の非確率行列T=(tij)は、行列Σ=(σij)について、TT´=Σを満たす。ただし、行列T´はTの転置行列である。
このとき、次の式(21)で与えられる確率ベクトルX=(x,・・・,x)´は、Cov(X)=Σを満たす。ただし、Cov(・)は、ベクトルの分散共分散行列であり、ξ=(ξ,・・・,ξ)´である。
【数25】
JP0004350062B2_000225t.gif

【請求項9】
前記確率変数ξ,・・・,ξには、ピアソン分布系に従う乱数を使う請求項に記載の多変量非正規分布のパラメータ推定方法。
【請求項10】
前記確率変数ξ,・・・,ξには、少なくとも2以上の型のピアソン分布を用いる請求項に記載の多変量非正規分布のパラメータ推定方法。
【請求項11】
請求項1乃至請求項7のいずれか一に記載の乱数発生方法を実現する、コンピュータにより実行可能なコンピュータ・プログラム。
【請求項12】
請求項1記載のコンピュータ・プログラムを記録した記録媒体。
【請求項13】
請求項8乃至請求項10のいずれか一に記載の多変量非正規分布のパラメータ推定方法を実現する、コンピュータにより実行可能なコンピュータ・プログラム。
【請求項14】
請求項13記載のコンピュータ・プログラムを記録した記録媒体。
発明の詳細な説明 【技術分野】
【0001】
本発明は、多変量非正規分布に従う乱数発生方法及び多変量非正規分布のパラメータ推定方法、ピアソンIV型分布に従う乱数発生方法及びコンピュータ・プログラムに関する。
【背景技術】
【0002】
従来、多変量非正規分布を生成する方法は、数えるしか存在していない。このような分布は、多変量解析と統計モデリングの分野において、切望されている。
【0003】
フライシュマン(Fleishmann)(1978)、タジカマラ(Tadikamalla)(1980)、ベイルとマブレリ(Vale and Mavrelli)(1983)は、多変量非正規分布を生成する方法を提案した。これらの方法は、シミュレーションした母集団の確率ベクトルの4次のモーメント行列を計算することが難しいという欠点がある。ユアン及びベントラー(1997)は、4次のモーメント行列を考慮し、特定の周辺の歪度と尖度を有する多変量非正規分布を生成する手法を提案した。
【0004】
本明細書では、ユアン等が方法Iと称したものは、次の通りである。ξ1,・・・,ξmは、独立な確率変数であり、E(・)を期待値として、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj及びξ=(ξ1,・・・,ξ1)´である。νは、ξに独立な確率変数であり、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。また、T=(tij)は、階数pのp行m列の非確率行列であり、TT´=Σとm≧pを満たす。すると、次の式(35)で与えられる確率ベクトルは、一般に、非楕円分布に従い、Cov(X)=Σを満たす。ここで、Cov(・)は、確率ベクトルX=(x1,・・・,xn)´の分散共分散行列である。
【0005】
【数40】
JP0004350062B2_000002t.gif
iの周辺の歪度(skew)と尖度(kurt)は、それぞれ次の式(36)及び(37)によって与えられる。
【0006】
【数41】
JP0004350062B2_000003t.gif
【数42】
JP0004350062B2_000004t.gif
Xは、ξ~Nm(0,I)及びPr(ν≧0)=1(ただし、Nmはm次元正規分布、Prは確率を表す。)が満たされるとき、楕円分布に従う(ファン、コッツ及びング(Fang, Kotz and Ng)(1990))。ν=1とし、z1,・・・,zpを独立な標準正規変数とし、ξ=|z1|-E(|z1|),ξj=zj,j=2,・・・,pとすると、式(35)におけるXは、アザリニ及びバレ(Azzalini and Valle)(1996)によって定義された歪んだ正規分布に従う。
【0007】
この方法には、母集団の4次モーメント行列を容易に計算できるという利点がある。Sを式(35)のXからの大きさnの標本に基づいた標本共分散とする。vech(・)を対称行列から重複しない要素を取り出してベクトルに変換する演算子とし、s=vech(S)とする。すると、sの漸近的分散共分散行列はΓ/nによって与えられる。ここで、Γ=var(vech(XX´))であり、次の式(38)によって与えられる。varは、分散である。
【0008】
【数43】
JP0004350062B2_000005t.gif
ただし、tjは、T=(tij)の第j列ベクトルであり、
【数44】
JP0004350062B2_000006t.gif
は、クロネッカー積である。Dpは、p次のデュプリケーション行列であり、Dp+は、Dpのムーア・ペンローズの一般化逆行列である。
【0009】
この後、ユアン等は、検定統計の理論を研究した(ユアン及びベントラー(1999a,1999b,2000))。ユアン等は、確率ベクトルを生成するにはフライシュマン(1978)の方法を推奨している。しかし、この分布は、あまり広い歪度及び尖度を有していない。これに対して、ピアソン分布系は、次のような特徴を有する。
【0010】
(1)ピアソン分布系は、様々な歪度と尖度を有する広いクラスの分布を表現することができる。ピアソン分布系は、ガンマ分布、ベータ分布、t分布等のいくつかの良く知られた分布を含む。
【0011】
(2)ピアソン分布系は、4つのモーメントによって特徴付けられる。この特徴は、ユアン等の方法を適用するためには十分である。
【0012】
(3)IV型を除いたピアソン分布系に従う乱数は、正規分布とガンマ分布に従う乱数を用いて生成される。
【0013】
(4)最近、実現することが難しかったピアソンIV型分布を含む殆ど全ての種類を用いる実用的な手法が開発された(永原(2002))。
【0014】
ここで、ピアソン分布系を導入する。ケイ・ピアソン(K. Pearson)(1895)は、確率密度関数pに関する次の式(39)の微分方程式に従い、ピアソン分布系を定義した。
【0015】
【数45】
JP0004350062B2_000007t.gif
ピアソン第IV型分布は、正規化定数を計算することが難しいので、実現することが難しかったが、論文(永原(1999))において、計算に用いる再帰的公式を導出し、様々な問題を解決した。この後、非ガウスフィルタに用いるピアソンIV型分布を含む分布系の殆ど全ての種類を用いる実用的な手法を開発した(永原(2002))。これらの研究により、ピアソン分布系が実用的に利用できるようになった。
【0016】
ここで、ピアソン分布系に関して参照される研究者と文献を掲げる(ピアソン(1914)、ピアソン及びハートリー(1954)、ピアソン(1963)、ホードリ-(1968)、エルダートン及びジョンソン(1969)、オード(1972)、パリシュ(1983)、ジョンソン、コッツ及びバラカリシュナン(1994)、スチュアート及びオード(1994))。
【0017】
参考文献としては、次のものが挙げられる。
【0018】
[1]アザリニ,エイ.及びバレ,エイ.ディー.(Azzalini, A. and Valle, A. D.)(1996)「多変量の歪んだ正規分布(The multivariate skew-normal distribution)」、バイオメトリカ(Biometrika)、83巻、715~726頁
[2]ドブロイ,エル.(Devroye. L.)(1984)「対数凹密度を有する確率変数を生成する簡単なアルゴリズム(A simple algorithm for generating random variates with a log-concave density)」、コンピューティング(Computing)、33巻、247~257頁
[3]ドブロイ,エル.(Devroye. L.)(1986)「非一様な確率変数の生成(Non-Uniform Random Variate Generation)」、スプリンガー(Springer)、ニューヨーク(New York)
[4]エルダートン,ダブリュ.ピー.及びジョンソン,エヌ.エル.(Elderton, W. P. and Johnson, N. P.)(1969)「頻度曲線の系(Systems of Frequency Curves)」、ケンブリッジ大学出版(Cambridge University Press)
[5]ファング ケイ.ティー,コッツ,エス及びング,ケイ.ダブリュ.(Fang K. -T, Kotz, S. and Ng, K. W.)(1990)「対称多変量及び関連する分布(Symmetric Multivariate and Related Distributions)」、チャップマン及びホール(Chapman & Hall)、ロンドン(London)
[6]フライシュマン,エイ.(Fleishmann, A.)(1978)「非正規分布のシミュレーション方法(A method for simulating non-normal distributions)」、サイコメトリカ(Psychometrika)、43巻、521~532頁
[7]ホードリー,エイ.ビー.(Hoadley, A. B.)(1968)「左端と最初の3つのモーメントが知られた歪んだ密度を近似するためのピアソン分布の使用(Use of the Pearson densities for approximating a skew density whose left terminal and first three moments are known)」、バイオメトリカ(Biometrika)、55巻、3号、559~563頁
[8]ジョンソン,エヌ.エル,コッツ,エス.及びバラクリシュナン,エヌ.(Johnson, N. L, Kotz, S. and Balakrishnan, N.)(1994)「連続単変量分布1(Continuous univariate distribution)」、第2版、ジョン・ワイリー(John Wiley)、
[9]マグナス,ジェイ.アール.及びニューデッカー,エイチ.(Magnus, J. R. and Neudecker, H.)(1988)「統計及び計量経済学における応用を有する行列微分法(Matrix Differential Calculus with Application in Statistics and Econometrics)」、ワイリー(Wiley)、ニューヨーク(New York)
[10]マクグラス,イー.ジェイ.及びアービング,ディー.シー(McGrath, E. J. and Irving, D. C.)(1973)「効果的なモンテカルロ・シミュレーションの技術 第2巻:選択した確率分布についての乱数の生成(Techniques for efficient Monte Carlo simulation: volume 2: random number generation for selected probability distributions)」、テクニカル・リポート(Technical Report SAI-72-590-LJ)、サイエンス・アプリケーションズ社(Science Applications, Inc.),ラヨラ(La Jolla)、カリフォルニア(CA)
[11]ムーアヘッド,アール.ジェー.(Muirhead, R. J.)(1982)「多変量統計理論の諸問題(Aspscts of Multivariate Statistical Theory)」、ワイリー(Wiley)、ニューヨーク(New York)
[12]永原,ワイ.(Nagahara, Y.)(1999)「ピアソンIV型分布の確率密度関数及び分布関数とパラメータの最尤推定(The PDF and CF of Pearson type IV distributions and the ML estimation of the parameters)」、スタティスティクス・アンド・プロバビリティー・レターズ(Statistics & Probability Letters)、43巻、251~264頁
[13]永原,ワイ(Nagahara, Y)(2002)「ピアソン分布系に基づく非ガウスフィルタ及びスムーザー(Non-Gaussian filter and smoother based on the Pearson distribution system)」、ジャーナル・オブ・タイム・シリーズ・アナリシス(Journal of Time Series Analysis)、近刊
[14]オード,ジェイ.ケイ.(Ord, J. K.)(1972)「頻度分布の族(Family of Frequency Distributions)」、グリフィン(Griffin)出版
[15]パリシュ,アール.(Parrish, R.)(1983)「ピアソン分布についてのメンバー選択とパラメータ推定への積分アプローチ(On an integrated approach to member selection and parameter estimation for Pearson distributions)」、コンピュテイショナル・スタティスティクス及びデータ解析(Computational Statistics & Data Analysis)、1巻、239~255頁
[16]ピアソン,イー.エス.(Pearson, E. S.)(1963)「モーメントを用いた確率分布近似で生じるいくつかの問題(Some problems arising in approximating to probability distributions, using moments)」、バイオメトリカ(Biometrika)、50巻、95~111頁
[17]ピアソン,ケイ(Pearson, K.)(1895)「均質な物質における歪んだ分散についての研究(Memoir on skew variation in homogeneous material)」、フィロソフィカル・トランザクションズ・オブ・ザ・ローヤル・ササイアティー(Phil. Trans. Roy. Soc.)、エイ(A)、186巻、343~414頁
[18]ピアソン,ケイ(Pearson, K.)(1914)「統計及び計量生物学のための数表(Tables for Statistics and Biometrics)」、1巻、ケンブリッジ大学出版(Cambridge University Press)
[19]ピアソン,イー.エス.及びハートリー,エイチ.オー.(Pearson, E. S. and Hartrey, H. O. )(1954)「統計のための計量生物学数表(Biometrika Tables for Statistics)」、1巻、ケンブリッジ大学出版(Cambridge University Press)
[20]ショラック,ジー.アール.(Shorack, G. R.)(2000)「統計のための確率(Probability for Statistics)」、スプリンガー(Springer)、ニューヨーク(New York)
[21]スチュアート,エイ及びオード,ジェイ.ケイ.(Stuart, A. and Ord, J. K.)(1994)「ケンダルによる統計の進んだ理論(Kendall's Acvanced Theory of Statistics)」、1巻、分布理論(Distribution Theory)、6版、エドワード・アーノルド(Edward Arnold)
[22]タジカマラ,ピー.アール.(Tadikamalla, P. R.)(1980)「非正規分布のシミュレーションについて(On simulating non-normal distributions)」、サイコメトリカ(Psychometrika)、45巻、273~279頁
[23]ユアン,ケイ.エイチ.及びベントラー,ピー.エム.(Yuan K, -H. and Bentler P. M.)(1997)「特定の周辺の歪度及び尖度を有する多変量分布の生成(Generating multivariate distributions with specified marginal skewness and kurtosis)」、ソフトスタット‘97-統計ソフトウェアにおける発展(ダブリュ.バンディラ及びエフ.ファウルバウム編集)内(in SoftStat'97-Advances in Statistical Software 6, (W. Bandilla and F. Faulbaum, Eds.))、385~391頁、ルシアス・アンド・ルシアス(Lucius & Lucius)、シュトゥットガルト(Stuttgart)、ドイツ
[24]ユアン,ケイ.エイチ.及びベントラー,ピー.エム.(Yuan K, -H. and Bentler P. M.)(1999a)「非正規分布における2つのクラスの下での共分散構造解析における正規理論及び付随する検定統計(On normal theory and associated test statistics in covariance structure analysis under two classes of nonnormal distributions)」、スタティスティカ・シニカ(Statist. Sinica)、9巻、831~853頁
[25]ユアン,ケイ.エイチ.及びベントラー,ピー.エム.(Yuan K, -H. and Bentler P. M.)(1999b)「いくつかの非正規分布の下での共分散構造解析における正規理論MLEの漸近的分布について(On asymptotic distributions of normal theory MLE in covariance structure analysis under some nonnormal distributions)」、スタティスティクス・アンド・プロバビリティー・レターズ(Statistics & Probability Letters)、42巻、107~113頁
[26]ユアン,ケイ.エイチ.及びベントラー,ピー.エム.(Yuan K, -H. and Bentler P. M.)(2000)「非正規分布のいくつかのクラスにおける相関係数の推定Inferences on correlation coefficients in some classes of nonnormal distributions)」、ジャーナル・オブ・マルチヴァリエイト・アナリシス(Journal of Multivariate Analysis)、72巻、230~248頁
[27]ベイル,ディー及びモーレイ,ブイ.エイ.(Vale, D. and Maurelli, V. A.)(1983)「多変量非正規分布シミュレーション(Simulating Multivariate Nonormal Distributions)」、サイコメトリカ(Psycometrika)、48巻、465~471頁
【発明の開示】
【発明が解決しようとする課題】
【0019】
ところで、検定統計のみならず統計モデリングの見地からも、多変量非正規分布を導出し、経験分布を近似する分布を構築することがさらに求められている。また、少数の標本からも多変量非正規分布を安定して推定することが求められている。
【0020】
本発明は、前述の実情に鑑みて提案されるものであって、ピアソン分布系及び一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンIV型分布に従う乱数発生方法及びコンピュータ・プログラムを提供することを目的とする。
【課題を解決するための手段】
【0021】
前述の課題を解決するために、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータ上でユアン及びベントラーの方法Iに基づいて多変量非正規分布に従う乱数を発生する乱数発生方法において、コンピュータを用いてn次元経験分布にn次元多変量非正規分布をあてはめるステップと、コンピュータを用いて乱数を発生するステップとを有し、前記あてはめるステップは、前記経験分布の3次及び4次モーメントに関するあてはめについて、次の関係式(40)及び(41)を用いる。
【0022】
【数46】
JP0004350062B2_000008t.gif
【数47】
JP0004350062B2_000009t.gif
ただし、E(・)は、期待値であり(以下同じ)、vech(・)は、対称行列から重複しない行列要素を取り出したベクトルであり、Dは、n次のデュプリケーション行列であり、D+は、Dのムーア・ペンローズの一般化逆行列であり、
【数48】
JP0004350062B2_000010t.gif
は、クロネッカー積であり、Eiiは、eiを第i列単位としてeii'である。
【0023】
ここで、前記ユアン及びベントラーの方法Iは、次の通りである。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。
【0024】
ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。
【0025】
階数nのn行m列(m≧n)の非確率行列T=(tij)は、行列Σ=(σij)について、TT´=Σを満たす。ただし、行列T´はTの転置行列である。
【0026】
このとき、次の式(42)で与えられる確率ベクトルX=(x1,・・・,xn)´は、Cov(X)=Σを満たす。ただし、Cov(・)は、ベクトルの分散共分散行列であり、ξ=(ξ1,・・・,ξm)´である。
【0027】
【数49】
JP0004350062B2_000011t.gif
好ましくは、確率変数ξ1,・・・,ξm及びνには、ピアソン分布系に従う乱数を使う。
【0028】
好ましくは、確率変数ξ1,・・・,ξm及びνには、少なくとも2以上の型のピアソン分布を用いる。
【0029】
好ましくは、前記あてはめるステップは、式(43)及び(44)でそれぞれ与えられるn次元の経験分布(確率ベクトル(X,・・・,X)´とする。)の3次及び4次モーメントに関して、式(45)の値を最小にするパラメータT,ζ=(ζ1,・・・,ζ),γ,κ=(κ1,・・・,κ),βの少なくとも1つを求める。
【0030】
【数50】
JP0004350062B2_000012t.gif
【数51】
JP0004350062B2_000013t.gif
【数52】
JP0004350062B2_000014t.gif
ただし、fijk(T,ζ,γ,κ,β)及びfijkl(T,ζ,γ,κ,β)は、3次モーメントE(xijk)及び4次モーメントE(xijk)にそれぞれ対応する関係式(40)及び(41)による表現とする。また、wijk及びwijklは、所定の重みである。
ここで、記号
【数53】
JP0004350062B2_000015t.gif
は、
【数54】
JP0004350062B2_000016t.gif
(ただし、i≦j≦k)の意味で、1≦i≦n,1≦j≦n,1≦k≦n,i≦j≦kの各条件を満たす全てのi,j,kの組み合わせについて式に代入し、それらの総和を取ったものである。同様に、記号
【数55】
JP0004350062B2_000017t.gif
は、
【数56】
JP0004350062B2_000018t.gif
(ただし、i≦j≦k≦l)の意味で、1≦i≦n,1≦j≦n,1≦k≦n,1≦l≦n,i≦j≦k≦lの各条件を満たす全てのi,j,k,lの組み合わせについて式に代入し、それらの総和を取ったものである。
また、前記mijk=E(Xijk)及びmijkl=E(Xijkl)の意味について説明する。データの分布(X,・・・,X)´は、次の式(A1)のようにN個のデータ・ベクトル
【数57】
JP0004350062B2_000019t.gif
(ただし、
【数58】
JP0004350062B2_000020t.gif
は各々のデータ(実現値)を示す。1≦q≦N)から構成されている。
【数59】
JP0004350062B2_000021t.gif
前記mijk及びmijklは、このようなデータ・ベクトルの要素となるデータを用いて次の式(A2)及び(A3)によって計算される。
【数60】
JP0004350062B2_000022t.gif
【数61】
JP0004350062B2_000023t.gif

【0031】
好ましくは、wijk=1,wijkl=1である。
【0032】
好ましくは、n=2であって、前記あてはめるステップは、式(46)及び(47)でそれぞれ与えられる2次元経験分布(確率ベクトル(X,X)´とする。)の3次及び4次モーメントと、式(48)及び(49)で与えられる3次モーメントE(xijk)及び4次モーメントE(xijk)にそれぞれ対応する表現とについて、式(50)で与えられる値を最小にするようにパラメータT,ζ=(ζ1,ζ),γ,κ=(κ1,κ),βの少なくとも1つを求める。
【0033】
【数62】
JP0004350062B2_000024t.gif
【数63】
JP0004350062B2_000025t.gif
【数64】
JP0004350062B2_000026t.gif
【数65】
JP0004350062B2_000027t.gif
【数66】
JP0004350062B2_000028t.gif
好ましくは、wi=1である。
【0034】
好ましくは、パラメータT,ζ=(ζ,・・・,ζ),γ,κ=(κ,・・・,κ),βの少なくとも1つを最尤法により推定する。
【0035】
好ましくは、式(51)で与えられるκの値を用いて、表8に示すκの値と分布の型との対応関係の少なくとも1つを用いて前記確率変数の属する型を決定する。
【0036】
【数67】
JP0004350062B2_000029t.gif
【表8】
JP0004350062B2_000030t.gif
ただし、βとβは、それぞれ歪度の2乗と尖度であり、式(51)及び表8におけるκは、前記パラメータκ=(κ1,・・・,κ)とは異なる。
【0037】
好ましくは、前記決定した型に応じて、次の表9及び表10に示す生成法の少なくとも1つを用いて乱数Zを発生する。
【0038】
【表9】
JP0004350062B2_000031t.gif
【表10】
JP0004350062B2_000032t.gif
好ましくは、コンピュータを用いて、式(52)に示すピアソンIV型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンIV型分布に従う乱数を発生させる。
【0039】
【数68】
JP0004350062B2_000033t.gif
好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、式(53)に示すピアソンIV型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンIV型分布に従う乱数を発生させる。
【0040】
【数69】
JP0004350062B2_000034t.gif
好ましくは、前記乱数を発生するステップは、コンピュータを用いて発生した、加算生成法、M系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数に基づく。
【0041】
好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、n個のデータ
【数70】
JP0004350062B2_000035t.gif
からなるデータ・ベクトルX´について、経験分布のデータ{X´,・・・,X´}を取得するステップと、コンピュータを用いて、前記データ{X´,・・・,X´}を標準化して{X,・・・,X}とするステップと、コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、分散共分散行列Σを計算するステップと、コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、3次モーメントmijk(1≦i≦j≦k≦n)を計算するステップと、コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、4次モーメントmijkl(1≦i≦j≦k≦l≦n)を計算するステップと、コンピュータを用いて、前記分散共分散行列Σから行列Tを計算するステップと、前記データ・ベクトルが従う分布が、X=νTξ(ここで、ξは、m(ここで、m>nである。)個の確率変数ξ~ξからなる確率ベクトル(ξ,・・・,ξ)´である。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。また、ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。)を満たす非正規分布を有すると仮定して、コンピュータを用いて、3次モーメントmijk(1≦i≦j≦k≦n)及び4次モーメントmijkl(1≦i≦j≦k≦l≦n)とfijk(T,ζ,γ,κ,β),fijkl(T,ζ,γ,κ,β)(ここでζ=(ζ,・・・,ζ),κ=(κ,・・・,κ)である。)との差を損失とする損失関数を導入し、TT´=Σ(T´はTの転置行列である。)という条件の下で、全体として評価するリスク関数を最小化するようにパラメータT,(ζ,・・・,ζ),γ,(κ,・・・,κ),βを決定するステップと、コンピュータを用いて、前記決定されたパラメータ(ζ,・・・,ζ),γ,(κ,・・・,κ),βに基づいて、式(54)及び表11により(βとβは、それぞれ歪度の2乗と尖度である。)、前記確率ベクトル(ξ,・・・,ξ)´及び確率変数νがピアソン分布のどの型に属するかを決定するステップと、
【数71】
JP0004350062B2_000036t.gif
【表11】
JP0004350062B2_000037t.gif
コンピュータを用いて乱数を発生し、この乱数に基づいて前記確率ベクトル(ξ,・・・,ξ)´及び確率変数νについて乱数を計算するステップと、
コンピュータを用いて、X=νTξに基づいて、標準化されたXの乱数を計算するステップと、コンピュータを用いて、標準化されたXの乱数を標準化前の乱数に変換するステップと、を有する。
【0042】
ただし、式(54)及び表11において、βとβは、それぞれ歪度の2乗と尖度であり、κは、前記パラメータκ=(κ1,・・・,κ)とは異なる。
【0043】
好ましくは、本発明に係る多変量非正規分布に従う乱数発生方法は、コンピュータを用いて、n個のデータ
【数72】
JP0004350062B2_000038t.gif
からなるデータ・ベクトルX´について、経験分布のデータ{X´,・・・,X´}を取得するステップと、コンピュータを用いて、前記データ{X´,・・・,X´}を標準化して{X,・・・,X}とするステップと、コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、分散共分散行列Σを計算するステップと、
コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、3次モーメントmijk(1≦i≦j≦k≦n)を計算するステップと、コンピュータを用いて、前記標準化されたデータ{X,・・・,X}に基づいて、4次モーメントmijkl(1≦i≦j≦k≦l≦n)を計算するステップと、コンピュータを用いて、前記分散共分散行列ΣからT=Σ1/2(行列の平方根)を満たす行列Tを計算するステップと、前記データ・ベクトルが従う分布が、X=νTξ(ここで、ξは、m(ここで、m=nである。)個の確率変数ξ~ξからなる確率ベクトル(ξ,・・・,ξ)´である。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。また、ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。)を満たす非正規分布を有すると仮定して、コンピュータを用いて、3次モーメントmijk(1≦i≦j≦k≦n)及び4次モーメントmijkl(1≦i≦j≦k≦l≦n)とfijk(ζ,γ,κ,β),fijkl(ζ,γ,κ,β)(ここでζ=(ζ,・・・,ζ),κ=(κ,・・・,κ)である。)との差を損失とする損失関数を導入し、全体として評価するリスク関数を最小化するようにパラメータ(ζ,・・・,ζ),γ,(κ,・・・,κ),βを決定するステップと、コンピュータを用いて、前記決定されたパラメータ(ζ,・・・,ζ),γ,(κ,・・・,κ),βに基づいて、式(55)及び表12により(βとβは、それぞれ歪度の2乗と尖度である。)、前記確率ベクトル(ξ,・・・,ξ)´及び確率変数νがピアソン分布のどの型に属するかを決定するステップと、
【数73】
JP0004350062B2_000039t.gif
【表12】
JP0004350062B2_000040t.gif
コンピュータを用いて乱数を発生し、この乱数に基づいて前記確率ベクトル(ξ,・・・,ξ)´及び確率変数νについて乱数を計算するステップと、コンピュータを用いて、X=νTξに基づいて、標準化されたXの乱数を計算するステップと、コンピュータを用いて、標準化されたXを標準化前の乱数に変換するステップと、を有する。
ここで、経験分布のデータ{X´,・・・,X´}とは、次の意味である

【数74】
JP0004350062B2_000041t.gif
ここでは、X´(1≦r≦N)はデータ・ベクトルを表し、先のデータの分布(X,・・・,X)´での表記Xとは異なる。標準化されたデータ{X,・・・,X}も同様の意味である。
【0044】
ただし、式(55)及び表12において、βとβは、それぞれ歪度の2乗と尖度であり、κは、前記パラメータκ=(κ1,・・・κ)とは異なる。
【0045】
好ましくは、前記損失関数(L)及びリスク関数(R)は、式(56)及び(57)、式(58)及び(59)、又は式(60)及び(61)のいずれかの組で与えられる請求項13又は14に記載の多変量非正規分布のパラメータ推定方法。
【0046】
【数75】
JP0004350062B2_000042t.gif
【数76】
JP0004350062B2_000043t.gif
【数77】
JP0004350062B2_000044t.gif
【数78】
JP0004350062B2_000045t.gif
【数79】
JP0004350062B2_000046t.gif
【数80】
JP0004350062B2_000047t.gif
好ましくは、前記決定した型に応じて、次の表13及び表14に示す生成法の少なくとも1つを用いて乱数Zを発生する。
【0047】
【表13】
JP0004350062B2_000048t.gif
【表14】
JP0004350062B2_000049t.gif
好ましくは、コンピュータを用いて、式(62)に示すピアソンIV型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンIV型分布に従う乱数を発生させる。
【0048】
【数81】
JP0004350062B2_000050t.gif
好ましくは、前記乱数を発生するステップは、加算生成法、M系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。
【0049】
本発明に係る多変量非正規分布のパラメータ推定方法は、コンピュータを用いて、式(63)に示すように、所定の確率分布に従う確率変数ν、階数nのn行m列(m≧n)の非確率行列T、所定の確率分布に従う確率ベクトルξ=(ξ1,・・・,ξm)´の積として与えられる確率ベクトルX=(x1,・・・,xn)´で与えられる多次元非正規分布を計算するステップと、
【数82】
JP0004350062B2_000051t.gif
最尤法によりパラメータを推定するステップと、を有する。
【0050】
好ましくは、前記最尤法によりパラメータを推定するステップは、コンピュータを用いて、n次元空間を分割するステップと、コンピュータを用いて、経験分布のデータ{X1,・・・,XN}について、Xiが属する各分割に属する乱数の数を乱数の総数で除して、次の式(64)で与えられる当該分割の確率を求めるステップと、
【数83】
JP0004350062B2_000052t.gif
コンピュータを用いて、前記分割の確率を当該分割のN次元体積で除して前記Xiが属する、次の式(65)で与えられる分割の尤度fi(θ)を求めるステップと、
【数84】
JP0004350062B2_000053t.gif
コンピュータを用いて、尤度の積Πi=1Ni(θ)又は対数尤度の和Σi=1Nlog fi(θ)が最大になるようなパラメータθを推定するステップと、を有する。
【0051】
好ましくは、確率変数ξ1,・・・,ξm及びνには、ピアソン分布系に従う乱数を使う。
【0052】
好ましくは、確率変数ξ1,・・・,ξm及びνには、少なくとも2以上の型のピアソン分布を用いる。
【0053】
好ましくは、コンピュータを用いて、式(66)に示すピアソンIV型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンIV型分布に従う乱数を発生させる。
【0054】
【数85】
JP0004350062B2_000054t.gif
好ましくは、コンピュータを用いて、ユアン及びベントラーの方法Iに基づいて、次の関係式(67)及び(68)を用いてn次元経験分布にn次元多変量非正規分布をあてはめるステップを有し、このあてはめるステップで求めたパラメータの近傍について前記パラメータθを推定する。
【0055】
【数86】
JP0004350062B2_000055t.gif
【数87】
JP0004350062B2_000056t.gif
ただし、E(・)は、期待値であり(以下同じ)、vech(・)は、対称行列から重複しない行列要素を取り出したベクトルであり、Dは、n次のデュプリケーション行列であり、D+は、Dのムーア・ペンローズの一般化逆行列であり、
【数88】
JP0004350062B2_000057t.gif
は、クロネッカー積であり、Eiiは、eiを第i列単位としてeii'である。
【0056】
ここで、前記ユアン及びベントラーの方法Iは、次の通りである。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。
【0057】
ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。
【0058】
階数nのn行m列(m≧n)の非確率行列T=(tij)は、行列Σ=(σij)について、TT´=Σを満たす。ただし、行列T´はTの転置行列である。
【0059】
このとき、次の式(69)で与えられる確率ベクトルX=(x1,・・・,xn)´は、Cov(X)=Σを満たす。ただし、Cov(・)は、ベクトルの分散共分散行列であり、ξ=(ξ1,・・・,ξm)´である。
【0060】
【数89】
JP0004350062B2_000058t.gif
好ましくは、前記パラメータθは、パラメータT,ζ=(ζ1,・・・,ζ),γ,κ=(κ1,・・・,κ),βの少なくとも1つを表す。
【0061】
好ましくは、コンピュータを用いて加算生成法、M系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。
【0062】
好ましくは、本発明に係る多変量非正規分布のパラメータ推定方法は、コンピュータを用いて、ユアン及びベントラーの方法Iに基づいて、次の関係式(70)及び(71)を用いてn次元経験分布にn次元多変量非正規分布をあてはめるステップと、コンピュータを用いて、n次元経験分布がX=νTξ(ここで、ξは、m(ここで、m≧nである。)個の確率変数ξ~ξからなる確率ベクトル(ξ,・・・,ξ)´である。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。また、ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。)を満たす非正規分布を有すると仮定して、前記Xについて乱数を計算するステップと、コンピュータを用いて、各確率変数x~xについての空間を、所定間隔Δhで分割し、このn次元空間を(Δh)の超立方体へ分割するステップと、コンピュータを用いて、各データ・ベクトルの属する区間(Δh)内に存在する乱数の個数を乱数の総数で除して、当該区間(Δh)の確率P(θ)を計算するステップと、コンピュータを用いて、前記確率P(θ)をn次元体積(Δh)で除して、前記データ・ベクトルの属する区間(Δh)の尤度f(θ)を計算するステップと、コンピュータを用いて、前記尤度f(θ)の積Πi=1(θ)又は対数尤度の和Σi=1Nlog fi(θ)が最大になるようなパラメータθを計算するステップと、を有する。
【0063】
【数90】
JP0004350062B2_000059t.gif
【数91】
JP0004350062B2_000060t.gif
ただし、E(・)は、期待値であり(以下同じ)、vech(・)は、対称行列から重複しない行列要素を取り出したベクトルであり、Dは、n次のデュプリケーション行列であり、D+は、Dのムーア・ペンローズの一般化逆行列であり、
【数92】
JP0004350062B2_000061t.gif
は、クロネッカー積であり、Eiiは、eiを第i列単位としてeii'である。
【0064】
ここで、前記ユアン及びベントラーの方法Iは、次の通りである。独立な確率変数ξ1,・・・,ξmは、パラメータζj及びκjについて、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj(1≦j≦m)を満たす。
【0065】
ξjに独立な確率変数νは、パラメータγ及びβについて、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。
【0066】
階数nのn行m列(m≧n)の非確率行列T=(tij)は、行列Σ=(σij)について、TT´=Σを満たす。ただし、行列T´はTの転置行列である。
【0067】
このとき、次の式(72)で与えられる確率ベクトルX=(x1,・・・,xn)´は、Cov(X)=Σを満たす。ただし、Cov(・)は、ベクトルの分散共分散行列であり、ξ=(ξ1,・・・,ξm)´である。
【0068】
【数93】
JP0004350062B2_000062t.gif
好ましくは、確率変数ξ1,・・・,ξm及びνには、ピアソン分布系に従う乱数を使う。
【0069】
好ましくは、確率変数ξ1,・・・,ξm及びνには、少なくとも2以上の型のピアソン分布を用いる。
【0070】
好ましくは、コンピュータを用いて、式(73)に示すピアソンIV型分布の正規化定数の解析解又はこの展開に基づいて正規化定数を計算するステップと、コンピュータを用いて、棄却法により乱数を発生するステップと、によってピアソンIV型分布に従う乱数を発生させる。
【0071】
【数94】
JP0004350062B2_000063t.gif
好ましくは、コンピュータを用いて加算生成法、M系列、一般化フィードバック・シフトレジスタ法及びメルセンヌツイスターを含む、合同法を除く擬似乱数、準乱数、低ディスクレパンシー列、物理乱数を含む乱数を発生する。
【0072】
本発明に係る予想損失額の計算方法は、前記手順に従い、ポートフォリオを構成するn種類の資産のn次元損益分布にn次元多変量非正規分布をあてはめ、シミュレーションによって予想損失額を計算する。
【0073】
本発明に係る金融商品のシミュレーション方法は前記手順を用い、前記多変量非正規分布に従う乱数に基づいて金融商品の特性をコンピュータを用いてシミュレーションする。
【0074】
本発明に係る天候・保険デリバティブの設計又は価格付け方法は、前記手順を用いて天候・保険デリバティブを設計又は価格付けする。
【0075】
本発明に係るアセット・アロケーションの方法は、前記手順を用いてアセット・アロケーションを行う。
【0076】
本発明に係るイオン注入のシミュレーション方法は、前記手順を用い、前記多変量非正規分布に従う乱数に基づいてイオン注入によるイオンの分布をコンピュータを用いてシミュレーションする。
【0077】
本発明に係る能力評価の方法は、前記手順を用いて能力評価を行う。
【0078】
本発明に係るコンピュータにより実行可能なコンピュータ・プログラムは、前記手順を実現するものである。好ましくは、前記コンピュータ・プログラムは、光ディスク、磁気ディスク、光磁気ディスク又は磁気テープの少なくとも1つによって提供される。
【発明の効果】
【0079】
本発明によると、ピアソン分布系及び一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンIV型分布に従う乱数発生方法及びコンピュータ・プログラムを提供することができる。
【発明を実施するための最良の形態】
【0080】
以下、本発明に係る多変量非正規分布に従う乱数発生方法、多変量非正規分布のパラメータ推定方法、ピアソンIV型分布に従う乱数発生方法及びコンピュータ・プログラムの実施の形態並びにリスク・マネジメントにおける予想損失額の計算方法について、図面を参照して詳細に説明する。
【0081】
以下の本実施の形態は、大きくは次の構成からなっている。
【0082】
1.ピアソン分布系及び乱数生成法
1.1 ピアソン分布系
1.2 ピアソン分布IV型に従う乱数の生成
1.3 ピアソン分布系に従う乱数の生成
2. ピアソン分布系を用いた多変量非正規分布の生成
2.1 多変量正規分布の確率ベクトル
2.2 シミュレーション
2.3 応用のためのあてはめプロシージャ
3. 最尤法
4. 実施例
5. 変形例
5.1 金融商品のシミュレーション方法
5.2 半導体へのイオン注入
5.3 アセット・アロケーション
5.4 能力評価
本実施の形態は、様々な歪度と尖度を有する広いクラスの分布を表すことができるピアソン分布系を用い、多変量非正規分布をシミュレーションする方法を提案する。第1に、実現が困難であったピアソンIV型分布を含むピアソン分布系に従う乱数を生成するプロシージャを導く。第2に、ピアソン分布系を適用し、ユアンとベントラー(1997)によって開発された方法を用い、多変量非正規分布に従う乱数を生成する。第3に、本願発明者によって求められた3次のモーメント行列を用い、種々の標本分布について3次と4次のモーメント行列をあてはめ、近似的な多変量非正規分布を構築することができる。第4に、最尤法によるパラメータ推定方法を提案する。このパラメータ推定方法によって、少数のデータからモーメント法によっては実現できない正確さで、多変量非正規分布を推定することができることを示す。そして、これらの方法が、様々なビジネスに応用できることを、リスク・マネジメントにおける予想損失額の計算について具体的に説明する。最後に、変形例として、金融商品のシミュレーション方法、半導体へのイオン注入、アセット・アロケーション及び能力評価について簡単に触れる。
【0083】
1 ピアソン分布系及び乱数生成法
1.1 ピアソン分布系
ピアソン分布系の型は、次のパラメータκによって分類される。パラメータκは、次の式(74)によって定義される。
【0084】
【数95】
JP0004350062B2_000064t.gif
ただし、β1とβ2は、それぞれ歪度の2乗と尖度である。
【0085】
ピアソン分布は、κ<0,0<κ<1及び1<κについて、それぞれI型、IV型及びVI型と呼ばれる。これらの3つの場合は、「主要型」と呼ばれる(エルダートン及びジョンソン(Elderton and Johnson)(1969))。III型(κ=±∞)は、I型とVI型の境界にある。V型(κ=1)は、第IV型と第VI型の境界にある。κ=0の場合については、VII型、II型及び正規分布が含まれる。
【0086】
ピアソン分布系においては、型が指定されると、この分布の3つ又は4つのパラメータは、平均、分散、歪度及び尖度から決定される。歪度の2乗と尖度の間の典型的な関係は、多数の教科書に示されている(ジョンソン、コッツ及びバラカリシュナン(Johnson, Kotz and Balakrishnan)(1994)、スチュアート及びオード(Stuart and Ord)(1994)等)。
【0087】
まず、表15~18に示すように、ピアソン分布系は、最も適切な表現によると、8つの型に分類され、密度関数、及び平均、分散、歪度及び尖度による陽な形が与えられる。
【0088】
次に、モーメントから各型のパラメータを得るための変換式とプロシージャを示す。各型の変換式を表15~18に示す。変換式によって、全てのパラメータを順に容易に計算することができる。
【0089】
【表15】
JP0004350062B2_000065t.gif
【表16】
JP0004350062B2_000066t.gif
表16において、ピアソンI型分布について、p,q>1のとき単峰型、0<q≦1<pのときJ型、0<p≦1<qのときJ型、0<p,q≦1のときU型になる。なお、p=q=1のときには一様分布になる。また、ピアソンII型分布について、p>1のとき単峰型、0<p<1のときU型になる。なお、p=1のときには一様分布になる。
【0090】
【表17】
JP0004350062B2_000067t.gif
【表18】
JP0004350062B2_000068t.gif
1.2 ピアソン分布IV型に従う乱数の生成
本実施の形態においては、ピアソン分布IV型の確率密度関数は、次の式(75)によって定義される。ただし、Γ(・)はガンマ関数、B(・)はベータ関数である。
【0091】
【数96】
JP0004350062B2_000069t.gif
ドブロイ(Devrye)(1986)によると、確率変数XがピアソンIV型分布であってb>1のとき、次の式(76)により分布するとすると、arctan(x)は、式(77)によって定義される対数凹密度を有する。
【0092】
【数97】
JP0004350062B2_000070t.gif
【数98】
JP0004350062B2_000071t.gif
ドブロイ(1984)は、対数凹密度についての棄却アルゴリズム(指数版)が、この分布に従う乱数を生成するのに適用できることを示唆した。しかし、その正規化定数は、著書において問題として残されていた(ドブロイ(1986)の練習2.7)。論文(永原(1999))では、この問題を解析的に解き、その正規化定数が容易に計算できることを示した(式(78))。ピアソンIV型分布に従う乱数の生成のアルゴリズムを次に示す。Cは、ピアソンIV型分布の正規化定数を示す。
【0093】
【数99】
JP0004350062B2_000072t.gif
JP0004350062B2_000073t.gif このアルゴリズムは、正規化定数の解析解又はこの展開を計算するステップと、棄却法により乱数を発生するステップを含んでいる。
【0094】
ここで、乱数のシーズについて説明する。ある確率分布に従う乱数を発生させるにあたっては、一様乱数がその基本になっている。
【0095】
原理的には、一様乱数が与えられれば、逆関数法によって、様々な確率分布に従う乱数が発生できる。本発明においても、この一様乱数をもとにして、ピアソン分布系の確率分布に従う乱数の発生のために使う指数分布、ガンマ分布、正規分布に従う乱数を発生している。この一様乱数を発生させる方法には、大別すると物理乱数、擬似乱数及び準乱数及び低ディスクレパンシー列が含まれる。準乱数及び低ディスクレパンシー列とは、擬似乱数とは違って、ランダムな性質はもっていないが、一様に分布する点列である。
【0096】
擬似乱数には、合同法、加算生成法、M系列、一般化フィードバック・シフトレジスタ法(Generalized Feedback Shift-register; GFSR)法及びメルセンヌツイスターが含まれる。準乱数には、リチトメヤー(Ritchtmyer)列、ヘイセルグローブ(Haselgrove)列、ハマーズレイ(Hammersley)列、低ディスクレパンシー(Low Discrepancy)列である、ホルトン(Halton)列、一般化ホルトン(Halton)列、ソボル(Sobol)列、改良ソボル(Sobol)列、フォーレ(Faure)列、一般化フォーレ(Faure)列、一般化ニーダーライター(Niederreiter)列が含まれる。近年、低ディスクレパンシー列は、多く研究されてきて、一つの定義として定着してきている。
【0097】
なお、メルセンヌツイスターの参考文献には、マツモト,エム.及びニシムラ,ティー.(Matsumoto, M. and Nishimura, T)「メルセンヌツイスター: 623次元等分布一様擬似乱数発生(Mersenne twister: A 623-dimensionally equidistributed uniform pseudorandom number generator)」、1998年、エイシーエム・トランザクションズ・オン・モデリング・アンド・コンピューター・シミュレーション(ACM Transaction on Modeling and Computer Simulation)、8巻、3~30頁がある。
【0098】
また、ソボル列の参考文献には、プレス,ダブリュ.,ツコルスキー,エス.,ヴェタリング,ダブリュ.,フラネリー,ビー.(Press, W., Teukolsky, S., Vetterling, W., Flannery, B)「フォートラン77の数値処理法(Numerical Recipes in Fortran 77)」、1992年、第2版、ケンブリッジ大学出版(Cambridge Univ. Press.)がある。
【0099】
この中で、一般的な方法は、擬似乱数の合同法である。ところが、最近では、合同法よりはるかに長い周期を持つ、GFSR法やメルセンヌツイスターが開発されてきている。さらに、準乱数や低ディスクレパンシー列が研究されていて、実用にも供されてきている。さて、これらの周期のより長い擬似乱数(例えばメルセンヌツイスター)や準乱数(低ディスクレパンシー列)を本実施の形態に適用したところ、本実施の形態の実施に重要な乱数の性質を大幅に改善することができた。従来の合同法では実現できなかった正確性で、目標とする確率分布の平均、分散、歪度、尖度を再現した乱数を得ることができた。
【0100】
本実施の形態の一番の重要な点は、非常に多数の乱数で、もはや解析的には表現できないような複雑な性質を持った確率分布を表現することである。この目的のためには、従来の合同法よりさらに正確な確率分布を再現できる乱数が望まれ、特に本実施の形態の特徴である、歪度、尖度が正規分布と違う分布(ピアソンIV型等)については、擬似乱数では合同法より長い周期のもの(例えばメルセンヌツイスター)、又は準乱数を使うことによって大幅に改善された、より実用的な正確性が達成された。
【0101】
以下に1次元の例を述べる。まず、ピアソンIV型に関する例では、乱数の個数が多くなると、合同法に比べ、メルセンヌツイスター、準乱数(ここでは、ソボル列を使用する。ソボル列は、低ディスクレパンシー列でもある。)を使った方が、平均、分散、歪度、尖度すべてにわたって目標とする値(表の一段目の数)に大幅に近くなっている。
【0102】
例えば、表21(合同法)、表25(メルセンヌツイスター)、表29(ソボル列)において、100000個の乱数で1000回シミュレーションした結果、平均について、合同法では0.01023、メルセンヌツイスターでは0.00018、準乱数では、-0.00001となり、合同法での小数点以下2桁目の正確さが、メルセンヌツイスターでは4桁、準乱数では、5桁までと大幅に改善している。同様に、分散について、合同法では1.03433、メルセンヌツイスターでは1.00009、準乱数では、0.99999となり、合同法での小数点以下2桁目の正確さが、メルセンヌツイスターでは5桁、準乱数では、6桁までと大幅に改善している。
【0103】
また、本実施の形態で重要な歪度と尖度についてはさらに改善している。同じ表で、歪度について目標値は0.75で、合同法では0.79651、メルセンヌツイスターでは0.75041、準乱数では、0.75007となり、合同法での小数点以下2桁目の正確さが、メルセンヌツイスターでは4桁、準乱数では、5桁までと大幅に改善している。尖度について目標値は5.0で、合同法では5.36334、メルセンヌツイスターでは5.00057、準乱数では、5.00257となり、合同法での小数点以下1桁目の正確さが、メルセンヌツイスターでは4桁、準乱数でも、3桁までと大幅に改善している。また、メルセンヌツイスター及び準乱数では、乱数の数が増えれば、それに従って目標の値にそれぞれ収束している点が合同法より優れており、この性質は、本実施の形態において乱数が持つべき最も重要な性質である。
【0104】
ここで、モンテカルロ・シミュレーション技法として対照変量(負の相関)法、制御変量法、モーメント適合法、回帰分析法、非線形回帰分析法、加重抽出法、層化抽出法、ラテン・ハイパーキューブ法、マルチンゲール分散減少法、条件付モンテカルロ法のような分散減少法が容易に適用できる。これら技法は、乱数の発生法に依らず、計算を改良する方法である。
【0105】
表19~30において、数は、乱数の個数を表し、値は、1000回の平均である。各括弧内の値は各値の標準偏差である。この結果、アルゴリズムを用いて妥当な乱数を得た。
【0106】
【表19】
JP0004350062B2_000074t.gif
【表20】
JP0004350062B2_000075t.gif
【表21】
JP0004350062B2_000076t.gif
【表22】
JP0004350062B2_000077t.gif
【表23】
JP0004350062B2_000078t.gif
【表24】
JP0004350062B2_000079t.gif
【表25】
JP0004350062B2_000080t.gif
【表26】
JP0004350062B2_000081t.gif
【表27】
JP0004350062B2_000082t.gif
【表28】
JP0004350062B2_000083t.gif
【表29】
JP0004350062B2_000084t.gif
【表30】
JP0004350062B2_000085t.gif
1.3 ピアソン分布系に従う乱数の生成
IV型を除くピアソン分布系については、正規分布とガンマ分布に従う乱数を用いることで、容易に乱数を生成できることが知られている(マクグラス及びアービング(McGrath and Irving)(1973))。本実施の形態では、表15~18のような定義を実現する実用的なプロシージャを開発する。これらの結果は、表31と表32に示されている。典型的な例は、表33~53に示されている。数は、乱数の個数を示し、値は1000回の平均である。各括弧内の値は、各値の標準偏差である。この結果、アルゴリズムを用いて妥当な乱数を得る。
【0107】
ピアソンIV型以外の分布系については、合同法に比べ、メルセンヌツイスター、準乱数(ここでは、ソボル列を使用。ソボル列は低ディスクレパンシー列でもある。)を使った方が、平均、分散では、ほとんどすべてで大幅に改善した。重要な歪度、尖度については、幅の広い領域を受け持つピアソンVI型及びピアソンI型では、100000個の乱数で1000回シミュレーションの結果で、まず、ピアソンVI型で、表33(合同法)、表40(メルセンヌツイスター)及び表44(ソボル列)において、歪度について目標値は1.0で、合同法では0.99323、メルセンヌツイスターでは0.99975、準乱数では0.99968となり、合同法での小数点以下3桁目の正確さが、メルセンヌツイスターでは4桁、準乱数でも、同様に4桁までと大幅に改善している。同じ表で、尖度について目標値は4.8で、合同法では4.73386、メルセンヌツイスターでは4.79863、準乱数では4.79836となり、合同法での小数点以下2桁目の正確さが、メルセンヌツイスターでは3桁、準乱数でも、同様に3桁までと改善している。
【0108】
次に、ピアソンI型について、表34(合同法)、表41(メルセンヌツイスター)及び表48(ソボル列)において、歪度について目標値は1.0で、合同法では1.00610、メルセンヌツイスターでは0.99945、準乱数では1.00004となり、合同法での小数点以下3桁目の正確さが、メルセンヌツイスターでは4桁、準乱数は5桁までと大幅に改善している。同じ表で、尖度について目標値は4.0で、合同法では4.00580、メルセンヌツイスターでは3.99737、準乱数では4.00037となり、合同法での小数点以下3桁目の正確さが、メルセンヌツイスターでは同様の3桁、準乱数では、4桁までと改善している。
【0109】
ピアソンII型では、表35(合同法)、表42(メルセンヌツイスター)及び表49(ソボル列)において、歪度について目標値は0.0で、合同法では0.00012、メルセンヌツイスターでは-0.00005、準乱数では0.00005となり、合同法での小数点以下4桁目の正確さが、メルセンヌツイスターでは5桁、準乱数でも、同様に5桁までと大幅に改善している。尖度について目標値は2.5で、同じ表で、合同法では2.49490、メルセンヌツイスターでは2.49992、準乱数では2.50014となり、合同法での小数点以下3桁目の正確さが、メルセンヌツイスターでは5桁、準乱数でも、4桁までと改善している。
【0110】
正規分布では、表36(合同法)、表43(メルセンヌツイスター)及び表50(ソボル列)において、歪度について目標値は0.0で、合同法では0.00439、メルセンヌツイスターでは0.00023、準乱数では0.00000となり、合同法での小数点以下3桁目の正確さが、メルセンヌツイスターでは4桁、準乱数では6桁までと大幅に改善している。尖度について目標値は3.0で、合同法では2.98130、メルセンヌツイスターでは3.00018、準乱数では3.00001となり、合同法での小数点以下2桁目の正確さが、メルセンヌツイスターでは4桁、準乱数でも、5桁までと改善している。
【0111】
ピアソンVII型では、表37(合同法)、表44(メルセンヌツイスター)、表51(ソボル列)において、歪度について目標値は0.0で、合同法では0.00282、メルセンヌツイスターでは0.00062、準乱数では-0.00047となり、合同法での小数点以下3桁目の正確さが、メルセンヌツイスターでは4桁、準乱数でも4桁までと大幅に改善している。尖度について目標値は4.0で、合同法では3.96639、メルセンヌツイスターでは3.99690、準乱数では3.99894となり、合同法での小数点以下2桁目の正確さが、メルセンヌツイスターでは3桁、準乱数でも、3桁までと改善している。
【0112】
ピアソンIII型では、表38(合同法)、表45(メルセンヌツイスター)、表52(ソボル列)において、歪度について目標値は1.0で、合同法では1.00128、メルセンヌツイスターでは0.99951、準乱数では0.99996となり、合同法での小数点以下3桁目の正確さが、メルセンヌツイスターでは4桁、準乱数では5桁までと大幅に改善している。尖度について目標値は4.5で、合同法では4.45711、メルセンヌツイスターでは4.49864、準乱数では4.49978となり、合同法での小数点以下2桁目の正確さが、メルセンヌツイスターでは3桁、準乱数では、4桁までと改善している。
【0113】
ピアソンV型では、表39(合同法)、表46(メルセンヌツイスター)、表53(ソボル列)において、歪度について目標値は1.0で、合同法では1.00456、メルセンヌツイスターでは0.99842、準乱数では0.99992となり、合同法での小数点以下3桁目の正確さが、メルセンヌツイスターでは同様の3桁、準乱数では5桁までと大幅に改善している。尖度について目標値は4.970で、合同法では5.15505、メルセンヌツイスターでは4.95829、準乱数では4.96925となり、合同法での小数点以下1桁目の正確さが、メルセンヌツイスターでは2桁、準乱数では、4桁までと大幅に改善している。また、以上分布についても、ピアソンIV型のときと同様に、メルセンヌツイスター及び準乱数では、乱数の数が増えれば、それに従って目標の値にそれぞれ収束している。この点が合同法より優れており、この性質は本発明において乱数が持つべき最も重要な性質である。
【0114】
【表31】
JP0004350062B2_000086t.gif
【表32】
JP0004350062B2_000087t.gif
【表33】
JP0004350062B2_000088t.gif
【表34】
JP0004350062B2_000089t.gif
【表35】
JP0004350062B2_000090t.gif
【表36】
JP0004350062B2_000091t.gif
【表37】
JP0004350062B2_000092t.gif
【表38】
JP0004350062B2_000093t.gif
【表39】
JP0004350062B2_000094t.gif
【表40】
JP0004350062B2_000095t.gif
【表41】
JP0004350062B2_000096t.gif
【表42】
JP0004350062B2_000097t.gif
【表43】
JP0004350062B2_000098t.gif
【表44】
JP0004350062B2_000099t.gif
【表45】
JP0004350062B2_000100t.gif
【表46】
JP0004350062B2_000101t.gif
【表47】
JP0004350062B2_000102t.gif
【表48】
JP0004350062B2_000103t.gif
【表49】
JP0004350062B2_000104t.gif
【表50】
JP0004350062B2_000105t.gif
【表51】
JP0004350062B2_000106t.gif
【表52】
JP0004350062B2_000107t.gif
【表53】
JP0004350062B2_000108t.gif
2 ピアソン分布系を用いた多変量非正規分布の生成
2.1 多変量非正規分布の確率ベクトル
ユアン等(1997)によると、ユアン等が方法Iと称したものは、次の通りである。ξ1,・・・,ξmは、独立な確率変数であり、E(ξj)=0,E(ξj2)=1,E(ξj3)=ζj,E(ξj4)=κj及びξ=(ξ1,・・・,ξm)´である。νは、ξに独立な確率変数であり、E(ν)=0,E(ν2)=1,E(ν3)=γ,E(ν4)=βを満たす。ここで、この分布は、νが1のとき、ピアソン分布系の確率密度関数によって表されることが特記される。また、T=(tij)は、階数pの非確率p行m列の行列であり、TT´=Σ及びm≧pを満たす。すると、次の式(79)で与えられる確率ベクトルは、一般に、Cov(X)=Σを満たす非楕円分布に従うようになる。ここで、Cov(・)は、ベクトルX=(x1,・・・,xn)´の分散共分散行列である。
【0115】
【数100】
JP0004350062B2_000109t.gif
iの周辺の歪度及び尖度は、それぞれ次の式(80)及び(81)によって与えられる。
【0116】
【数101】
JP0004350062B2_000110t.gif
【数102】
JP0004350062B2_000111t.gif
そして、4次のモーメント行列Γ=var(vech(XX´))は、次の式(82)によって与えられる。
【0117】
【数103】
JP0004350062B2_000112t.gif
ただし、tjは、T=(tij)の第j列ベクトルであり、
【数104】
JP0004350062B2_000113t.gif
は、クロネッカー積である。Dpは、p次のデュプリケーション行列であり、Dp+は、Dpのムーア・ペンローズの一般化逆行列である。
【0118】
統計モデリングの見地から、3次のモーメント行列が必要である。本願発明者は、3次のモーメント行列について、次の定理を得ることができた。
【0119】
定理3.1(3次のモーメント行列)3次のモーメント行列Λは、次の式(83)によって与えられる。
【0120】
【数105】
JP0004350062B2_000114t.gif
ここで、eiは、第i列単位であり、Eii=eii'である。
【0121】
証明 3次のモーメント行列Λは、Y=Tξとして、次の式(84)のようになる。
【0122】
【数106】
JP0004350062B2_000115t.gif
YY´=Tξξ´T´及びvech(YY´)=Dp+vec(YY´)である
ので、
【数107】
JP0004350062B2_000116t.gif
となる。ここで、記号
【数108】
JP0004350062B2_000117t.gif
は、行列のクロネッカー積を表す。以下でも同様である。
【0123】
したがって、次の式(86)を得る。
【0124】
【数109】
JP0004350062B2_000118t.gif
次の式(87)は、明らかである。
【0125】
【数110】
JP0004350062B2_000119t.gif
この定理は、式(84),(86)及び(87)から得られる。
【0126】
証明終わり
ユアン及びベントラー(1999a,1999b及び2000)は、検定統計理論を研究した。ユアン等は、次の式(88)を示唆した。
【0127】
【数111】
JP0004350062B2_000120t.gif
ただし、Π=∂f/∂vech´(Σ)の階数sの行列として、Ω=ΠΓΠ´であり、Sは、標本分散共分散行列である。また、ユアン等は、以前、Πvech(tjj´)=0,j=1,・・・,mを仮定すると、次の式(89)が成立することを示唆している。
【0128】
【数112】
JP0004350062B2_000121t.gif
ただし、
【数113】
JP0004350062B2_000122t.gif
である。
【0129】
本実施の形態は、様々な3次と4次のモーメント行列を用い、多変量非正規分布に従う乱数を得ることを目的としている。T,ξ,νを変化させることによって、異なる3次と4次のモーメント行列を有する様々な多変量非正規分布を得ることができる。まず、X1,X2,ξ1及びξ2を有するp=m=2の典型的な場合についてシミュレーションを行う。次に、3次と4次のモーメントを経験データの1つにあてはめるプロシージャを提案する。
【0130】
なお、ピアソン分布系のほかには、4次までのモーメントで表せるジョンソン分布系やラムダ分布などがνやξの分布として採用できる。
【0131】
2.2 シミュレーション
本実施の形態では、x1,x2,ξ1及びξ2を有するp=m=2の場合についてシミュレーションする。このシミュレーションにおいて、3次と4次のモーメント行列の要素は、次の通りである。
【0132】
[3次モーメント]
【数114】
JP0004350062B2_000123t.gif
ただし、f1,f2,f3,f4は、X=νTξで発生される確率ベクトルX=(x,x)´のそれぞれ3次モーメントE(x13),E(x23),E(x122),E(x122)の表現である。
【0133】
[4次モーメント]
【数115】
JP0004350062B2_000124t.gif
ただし、f5,f6,f7,f8,f9は、X=νTξで発生される確率ベクトルX=(x,x)´のそれぞれ4次モーメントE(x14),E(x24),E(x132),E(x1222),E(x123)の表現である。ここで、σ11=σ22=1であり、また、t11=t22=1-t12=1-t21である。
【0134】
ユアンとベントラー(2000)の例を採用し、ρ=0.5及びt11=t22=0.2588,t12=t21=0.9659とする。結果は、表54~80に示されている。特定の歪度及び尖度を有する、各分布のξ1,ξ2,νから乱数を得た後、この乱数を平均と分散がそれぞれ0と1になるように標準化する。続いて、式X=νTξによって、x1とx2のための乱数X,Xを得る。各表において、E(・)は、これらの乱数を用いて計算される。括弧内は、3次と4次のモーメント行列の式を用いて計算された理論値である。これらの結果によると、10,000,000点と1,000,000点の場合、この方法による乱数は、理論的なモーメント、すなわち分散共分散、3次と4次のモーメント行列を有する目的の分布を正確に表すことができる。
【0135】
ここで、異なるピアソン分布の型の組み合わせも可能であることは特記されるべきであることである。例えば、金融リスク、オペレーショナル・リスク、保険リスク、災害リスク、事業リスクは、一般に異なったピアソン分布の型で記述される。本実施の形態によると、異なったピアソン分布の組み合わせも可能なので、このような異なった種類のリスクの組み合わせも統合して管理することができる。例えば、各種の資産を含んだポートフォリオの予想損失額を統合して計算することができる。
【0136】
なお、金融での最大損失額の計算では、極値分布(ガンベル分布、フレッシュ分布、ワイブル分布)や一般化パレート分布・対数正規分布が使われるが、これらの分布もこの発明の方法で、νやξの分布として採用できる。このことによって、多変量化や、よりデータを正確に表す分布を推定できる。
【0137】
ピアソンIV型の場合、2次元においてξ、ξ、νすべてに、ピアソンIV型を用いた例を表54(合同法)、表63(メルセンヌツイスター)、表72(準乱数 ソボル列)に示す。以下、x、xは確率変数、X,Xは乱数を表す。
【0138】
まず、目標とする分布の特性はξ(歪度0.75、尖度5.0)、ξ(歪度1.0、尖度5.5)、ν(歪度0.50、尖度4.0)で、それから生成される多変量非正規分布のモーメント特性の理論値は、E(x)=0.45708、E(x)=0.34660、E(x)=0.14499、E(x)=0.12289、E(x)=20.73865、E(x)=19.00678、E(x)=8.46515、E(x)=7.12390、E(x)=8.03221である。
【0139】
例えば、10,000,000個のシミュレーションの場合、表54(合同法)では、それぞれ、ξ(歪度0.83082、尖度5.28928)、ξ(歪度1.05906、尖度6.19495)、ν(歪度0.58807、尖度4.17774)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.57507(0.45708)、E(X)=0.45499(0.34660)、E(X)=0.18715(0.14499)、E(X)=0.16383(0.12289)、E(X)=24.14493(20.73865)、E(X)=20.94897(19.00678)、E(X)=9.53054(8.46515)、E(X)=7.71168(7.12390)、E(X)=8.75500(8.03221)である。
【0140】
括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ、ξ、νの実際のシミュレーション時の歪度、尖度(例えば、上の例ではξ(歪度0.83082、尖度5.28928)、ξ(歪度1.05906、尖度6.19495)、ν(歪度0.58807、尖度4.17774))から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との開きは大きい。
【0141】
ところが、これをメルセンヌツイスターでシミュレーションすると、表63(メルセンヌツイスター)にあるように、それぞれ、ξ(歪度0.75279、尖度5.02123)、ξ(歪度1.00004、尖度5.50831)、ν(歪度0.50414、尖度4.00780)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.45185(0.45708)、E(X)=0.35019(0.34660)、E(X)=0.14202(0.14499)、E(X)=0.12404(0.12289)、E(X)=20.03210(20.73865)、E(X)=19.16341(19.00678)、E(X)=8.16249(8.46515)、E(X)=7.02127(7.12390)、E(X)=8.05117(8.03221)である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。
【0142】
また、準乱数(ソボル列)では表72(ソボル列)にあるように、それぞれ、ξ(歪度0.74961、尖度4.99697)、ξ(歪度0.99989、尖度5.49716)、ν(歪度0.49984、尖度3.99661)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.46513(0.45708)、E(X)=0.34724(0.34660)、E(X)=0.15195(0.14499)、E(X)=0.13147(0.12289)、E(X)=20.96140(20.73865)、E(X)=18.96493(19.00678)、E(X)=8.62444(8.46515)、E(X)=7.23779(7.12390)、E(X)=8.15866(8.03221)である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。
【0143】
ピアソンVI型の場合、2次元においてξ、ξ、νすべてに、ピアソンVI型を用いた例を表57(合同法)、表66(メルセンヌツイスター)、表75(準乱数 ソボル列)に示す。まず、目標とする分布の特性はξ(歪度1.0、尖度4.8)、ξ(歪度1.25、尖度5.5)、ν(歪度0.75、尖度4.0)で、それから生成される多変量非正規分布のモーメント特性の理論値は、E(x)=0.85783、E(x)=0.69211、E(x)=0.27488、E(x)=0.24174、E(x)=20.73506、E(x)=18.31045、E(x)=8.45176、E(x)=7.07391、E(x)=7.84563である。
【0144】
例えば、10,000,000個のシミュレーションの場合、表57(合同法)では、それぞれ、ξ(歪度0.99126、尖度4.64495)、ξ(歪度1.22615、尖度5.22386)、ν(歪度0.75304、尖度3.92123)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.84302(0.85783)、E(X)=0.67342(0.69211)、E(X)=0.26865(0.27488)、E(X)=0.24025(0.24174)、E(X)=19.20021(20.73506)、E(X)=17.14705(18.31045)、E(X)=7.96794(8.45176)、E(X)=6.76215(7.07391)、E(X)=7.41679(7.84563)である。
【0145】
括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ、ξ、νの実際のシミュレーション時の歪度、尖度(例えば、上の例ではξ(歪度0.99126、尖度4.64495)、ξ(歪度1.22615、尖度5.22386)、ν(歪度0.75304、尖度3.92123))から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との間に開きがある。
【0146】
ところが、これをメルセンヌツイスターでシミュレーションすると、表66(メルセンヌツイスター)にあるように、それぞれ、ξ(歪度0.99968、尖度4.80062)、ξ(歪度1.25042、尖度5.49703)、ν(歪度0.75022、尖度3.99877)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.86185(0.85783)、E(X)=0.67746(0.69211)、E(X)=0.27713(0.27488)、E(X)=0.23868(0.24174)、E(X)=20.87668(20.73506)、E(X)=18.01464(18.31045)、E(X)=8.59760(8.45176)、E(X)=7.10837(7.07391)、E(X)=7.75298(7.84563)である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。特に4次のモーメントがより近づいている。
【0147】
また、準乱数(ソボル列)では表75(ソボル列)にあるように、それぞれ、それぞれ、ξ(歪度1.00025、尖度4.80832)、ξ(歪度1.24990、尖度5.49862)、ν(歪度0.75052、尖度4.00533)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.85834(0.85783)、E(X)=0.67586(0.69211)、E(X)=0.26486(0.27488)、E(X)=0.22957(0.24174)、E(X)=20.75111(20.73506)、E(X)=18.35898(18.31045)、E(X)=8.34400(8.45176)、E(X)=6.94143(7.07391)、E(X)=7.70291(7.84563)である。この結果からすべての目標のモーメント特性の理論値に非常に近くなっていることがわかる。特に4次のモーメントがより近づいている。
【0148】
ピアソンI型の場合、2次元においてξ、ξ、νすべてに、ピアソンI型を用いた例を表60(合同法)、表69(メルセンヌツイスター)、表78(準乱数 ソボル列)に示す。まず、目標とする分布の特性はξ(歪度1.0、尖度4.0)、ξ(歪度1.25、尖度5.0)、ν(歪度0.75、尖度3.0)で、それから生成される多変量非正規分布のモーメント特性の理論値は、E(x)=0.85783、E(x)=0.69211、E(x)=0.27488、E(x)=0.24174、E(x)=14.23490、E(x)=11.63710、E(x)=5.94881、E(x)=5.06173、E(x)=5.29939である。
【0149】
例えば、10,000,000個のシミュレーションの場合、表60(合同法)では、それぞれ、ξ(歪度1.02802、尖度4.04760)、ξ(歪度1.29323、尖度5.15264)、ν(歪度0.76208、尖度2.99871)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.90491(0.85783)、E(X)=0.72634(0.69211)、E(X)=0.29244(0.27488)、E(X)=0.25642(0.24174)、E(X)=14.57587(14.23490)、E(X)=11.83105(11.63710)、E(X)=6.07564(5.94881)、E(X)=5.12728(5.06173)、E(X)=5.37784(5.29939)である。
【0150】
括弧内は、理論値を再掲した。また、各表中の括弧内の数字は、各ξ、ξ、νの実際のシミュレーション時の歪度、尖度(例えば、上の例ではξ(歪度1.02802、尖度4.04760)、ξ(歪度1.29323、尖度5.15264)、ν(歪度0.76208、尖度2.99871))から計算した理論値である。乱数の数を増やすと、表の中の括弧内の値と計測値は非常に近くなるが、実用上の目標値は、上で計算されたモーメント特性の理論値であり、合同法では、目標の理論値とシミュレーションの計測値との間に開きがある。
【0151】
ところが、これをメルセンヌツイスターでシミュレーションすると、表69(メルセンヌツイスター)にあるように、ξ(歪度1.00039、尖度4.00402)、ξ(歪度1.24914、尖度4.99567)、ν(歪度0.75009、尖度3.00061)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.86215(0.85783)、E(X)=0.69791(0.69211)、E(X)=0.27696(0.27488)、E(X)=0.24381(0.24174)、E(X)=14.28128(14.23490)、E(X)=11.66910(11.63710)、E(X)=5.93199(5.94881)、E(X)=5.05602(5.06173)、E(X)=5.29916(5.29939)である。この結果からすべての目標のモーメント特性の理論値と非常に近くなっていることがわかる。
【0152】
また、準乱数(ソボル列)では表78(ソボル列)にあるように、ξ(歪度0.99973、尖度3.99878)、ξ(歪度1.25082、尖度5.00640)、ν(歪度0.74954、尖度2.99766)で、それから生成される多変量非正規分布のモーメント特性の計測値は、E(X)=0.85151 (0.85783)、E(X)=0.69425(0.69211)、E(X)=0.27271(0.27488)、E(X)=0.24358(0.24174)、E(X)=14.12454(14.23490)、E(X)=11.63940(11.63710)、E(X)=5.87281(5.94881)、E(X)=5.02381(5.06173)、E(X)=5.28335(5.29939)である。この結果からすべての目標のモーメント特性の理論値と非常に近くなっていることがわかる。
【0153】
【表54】
JP0004350062B2_000125t.gif
【表55】
JP0004350062B2_000126t.gif
【表56】
JP0004350062B2_000127t.gif
【表57】
JP0004350062B2_000128t.gif
【表58】
JP0004350062B2_000129t.gif
【表59】
JP0004350062B2_000130t.gif
【表60】
JP0004350062B2_000131t.gif
【表61】
JP0004350062B2_000132t.gif
【表62】
JP0004350062B2_000133t.gif
【表63】
JP0004350062B2_000134t.gif
【表64】
JP0004350062B2_000135t.gif
【表65】
JP0004350062B2_000136t.gif
【表66】
JP0004350062B2_000137t.gif
【表67】
JP0004350062B2_000138t.gif
【表68】
JP0004350062B2_000139t.gif
【表69】
JP0004350062B2_000140t.gif
【表70】
JP0004350062B2_000141t.gif
【表71】
JP0004350062B2_000142t.gif
【表72】
JP0004350062B2_000143t.gif
【表73】
JP0004350062B2_000144t.gif
【表74】
JP0004350062B2_000145t.gif
【表75】
JP0004350062B2_000146t.gif
【表76】
JP0004350062B2_000147t.gif
【表77】
JP0004350062B2_000148t.gif
【表78】
JP0004350062B2_000149t.gif
【表79】
JP0004350062B2_000150t.gif
【表80】
JP0004350062B2_000151t.gif
2.3 応用のためのあてはめプロシージャ
統計モデリングの見地から、経験分布の1つに3次と4次のモーメント行列をあてはめるプロシージャが必要である。各要素と標本のモーメントの式の間の偏りの2乗の和を最小にすることによって、T,ζ,γ,κ及びβを計算するプロシージャを提案する。すなわち、TT´=Σという条件の下、次の式(93)によって定義される関数の値を最小化することによって、適切なT,ζ,γ,κ及びβを得ることができる。
【0154】
ここで、m=pのときT=Σ1/2(行列の平方根)であり、固定されている。行列の平方根の計算は、実対称行列の固有値問題に帰着され、おもな方法として、ヤコビ法、ハウスホルダー法、二分法、逆反復法、二分法・同時逆反復法、キューエル(QL)法、ディバイド・アンド・コンクウェア(Divide and Conquer)法がある。
【0155】
【数116】
JP0004350062B2_000152t.gif
ただし、fiは、3次と4次のモーメントの各要素の式であり、例えば、p=2について、f1=E(x13)=γt113ζ1+γt123ζ2,k=9であり、miは、標本分布のモーメントである。p=m=2の場合、Tは固定されている(ユアン及びベントラー(2000))。例えば、標本モーメントmiは、表81に従う。式(93)の値を最小にすることによって、ζ1=0.65967,ζ2=1.03624,γ=0.63553,κ1=5.27475,κ2=6.69176,β=4.02680というパラメータを得る。これらのパラメータによって得られた行列要素は、表81に示されている。ここで、wiは、重みである。
【0156】
【表81】
JP0004350062B2_000153t.gif
3 最尤法
前述のモーメント法にかわる方法として最尤法がある。モーメント法では、3次と4次モーメントのような高次モーメントについては不安定になることがあるが、最尤法ではこの点が改善される可能性がある。
【0157】
最尤法では、生成した乱数によって密度関数を数値的に近似する。例として、特定の多変量非正規分布を10,000,000点の乱数により非常に正確にシミュレーションする。最尤法のプロシージャは、以下の通りである。
【0158】
(1)最尤法に用いるパラメータの初期値をモーメント法から採用する。
【0159】
(2)多変量非正規分布から例えば10,000,000点の乱数を得る。
【0160】
(3)各次元をステップΔhによって分割に分ける。超立方体となる各分割内にある点の個数を数え、全ての点の個数で各分割内にある点の個数を除する。
【0161】
(4)次のように対数尤度を求める。
【0162】
{X1,・・・,XN}は、n次元データの集合(Xは、n次元データ・ベクトルである。)を示す。特定の多変量非正規分布のパラメータθを有するXiの分割の確率Pi(θ)は、次の式(94)によって定義される。ここで、パラメータθは、ユアンとベントラーの方法Iで用いられるパラメータT,ζ=(ζ1,・・・,ζ),γ,κ=(κ1,・・・,κ)、βの少なくとも1つを表すものとする。
【0163】
【数117】
JP0004350062B2_000154t.gif
したがって、パラメータθを有する近似した尤度fi(θ)は、次の式(95)によって定義される。
【0164】
【数118】
JP0004350062B2_000155t.gif
対数尤度の和は、次の式(96)によって計算される。
【0165】
【数119】
JP0004350062B2_000156t.gif
そして、この対数尤度の和を最大にすることにより、最尤推定量
【数146】
JP0004350062B2_000157t.gif
を得る。このようにして、最尤推定パラメータが得られた。なお、尤度の積を最大にしても同様である。
【0166】
多変量非正規分布推定量
【数147】
JP0004350062B2_000158t.gif
を、対数尤度の和を最大化することにより推定する。この方法を予想損失額に適用するため、1年の取引日に対応する250点を採用する。250点のモーメントを表82に示し、最尤法による推定分布のモーメントを表83~85に示す。
【0167】
250点は、合同法を用いて得られた、目的とする分布の特性が、ξ(歪度0.75、尖度5.0)、ξ(歪度1.00、尖度5.5)、ν(歪度0.5、尖度4.0)を持つ(先のピアソンIV型の多変量シミュレーションで使われた目標値)乱数である。250点からなる分布の特性の計測値は表82に示されており、目的とする分布の特性とはかなり違ったものになっている。
【0168】
次に、多数の乱数からなる(ここでは、10,000,000点)分布による最尤法の推定結果を示す。なお、分散共分散のみを経験分散共分散0.41819と同じに設定するため、t11=0.2148,t12=0.9767に設定する。まず、表83は合同法による推定結果であるが、250点の計測値よりも、目標とされた分布の特性値(表54)に近づいていることがわかる。このことから、少数のデータの分布の推定には、モーメント法より、最尤法がすぐれていることがわかる。なお、最大対数尤度は、-562.73になった。特に、金融のリスク管理などでは、一年のデータは250点程度であり,これから、多変量非正規分布を用い、良好な推定をするためには、最尤法を用いらなければならないといえよう。
【0169】
さらに、表84(メルセンヌツイスター)、表85(準乱数 ソボル列)を示す。前の多変量非正規分布に従う乱数のシミュレーション結果から、これらの方法は合同法に比べ、非常に目標とする分布に近い乱数を発生することがわかっており、これは最尤法を実行するための重要な性質である。結果として、両者の方法も、より目標値に近い値を推定している。なお、メルセンヌツイスターについての最大対数尤度は-561.36、ソボル列についての最大対数尤度は-560.36となった。
【0170】
本発明では、最尤法について主に連続型確率変数の場合の密度関数を乱数から直接計算している方法を中心にしているが、離散型確率変数の場合は密度関数の代わりに確率関数となる。また、密度関数や確率関数を積分した分布関数として表現する方法も可能である。
【0171】
本実施の形態は、一般性を失わないように標準化された平均0、分散1の設定(E[x]=0,Var[x]=σii=1(E[x]=σii=1),1≦i≦n)で行われているが、一般の場合でも標準化の逆変換を行えば容易に拡張でき、また、最尤法を使えば、さらに、その標準化のための平均と標準偏差もパラメータとして、直接、生のデータ(標準化などの前処理をしていないという意味)から分布を推定することが容易に可能である。
【0172】
【表82】
JP0004350062B2_000159t.gif
【表83】
JP0004350062B2_000160t.gif
【表84】
JP0004350062B2_000161t.gif
【表85】
JP0004350062B2_000162t.gif
なお、さらに、分布の1次から4次までの一部または全部のモーメントが解析的に容易にもとまらない場合や、発散などで存在しない場合、あるいは密度関数でなく特性関数で定義される分布などの場合でも、その分布の乱数を発生させてそれからモーメントを計算する方法(必ず計算可能)や、最尤法を行うにあたって分布を特徴付けるパラメータやその他のパラメータ(両パラメータを含めてθと表す。)を適当な間隔で適当に与え、その分布の乱数を用いて幅広く尤度を計算する方法を使えば、そのような分布もこの発明のνやξの分布として採用でき
る。そのような分布として例えば、一般の場合、2次のモーメントが発散し、密度関数が陽表的に書けないが特性関数によって表現できる安定分布などがある。ゆえに、本発明の最尤法は、νやξの分布としては、乱数が発生できれば1次から4次までの一部又は全部のモーメントが求まらない分布でも採用可能となり、ユアン・ベントラーの方法からさらに拡張された方法となっている。
【0173】
4 実施例
図1は、2次元を例とした予想損失額計算を示すフローチャートである。ここでは、日本株式リターンと外国株式リターンを2次元の変数に取る。この予想損失額計算の一連の手順は、コンピュータにおいて所定のコンピュータ・プログラムを実行することによって実現される。
【0174】
ステップS1においては、日本株式リターンX1´と外国株式リターンX2´を標準化する。
【0175】
標準化は、次のように経験分布の平均を0、分散を1とすることにより行う。日本株式リターンと外国株式リターンのある実現値を持つ確率変数X1´及びX2´について、それぞれの平均をE(X1´)及びE(X2´)、分散をV(X1´)及びV(X2´)とする。このとき、標準化日本株式リターンX1と標準化外国株式リターンX2を、次の式(97)と(98)により与える。
【0176】
【数120】
JP0004350062B2_000163t.gif
【数121】
JP0004350062B2_000164t.gif
ステップS2においては、標準化日本株式リターンX1と標準化外国株式リターンX2から分散共分散行列Σを計算する。
【0177】
ステップS3においては、標準化日本株式リターンX1と標準化外国株式リターンX2について3次モーメントを計算する。3次モーメントは、次の式(99)に示すように、m1~m4までの4種類がある。具体的には、E(・)は、X,Xのデータ(実現値)から計算される。
【0178】
【数122】
JP0004350062B2_000165t.gif
ステップS4においては、標準化日本株式リターンX1と標準化外国株式リターンX2について4次モーメントを計算する。4次モーメントは、次の式(100)に示すように、m5~m9までの5種類がある。
【0179】
【数123】
JP0004350062B2_000166t.gif
ステップS5においては、ステップS2において計算した分散共分散行列Σについて、Σ=TT´を満たす行列T=(tij)を計算する。いまξの次元とXの次元がともに2次元なので、行列Tは、T=Σ1/2(行列の平方根)と固定される。
【0180】
具体的には、次の式(101)の値を最小にするパラメータζ=(ζ1,ζ2),γ,κ=(κ1,κ2),βを決定する。
【0181】
【数124】
JP0004350062B2_000167t.gif
ここで、重みwiによって、モーメントに軽重をつけることができる。本実施例では、簡単のためにwi=1とする。
【0182】
式(101)においては、3次及び4次モーメントの表現として、次の式(102)と(103)の表現を用いる。行列Tは、T=Σ1/2(行列の平方根)と計算されているので、式(101)に最小値を与えるパラメータζ=(ζ1,ζ2),γ,κ=(κ1,κ2),βを求める。
【0183】
[3次モーメント]
【数125】
JP0004350062B2_000168t.gif
ただし、f1,f2,f3,f4は、X=νTξで発生される確率ベクトルX=(x,x)´のそれぞれ3次モーメントE(x13),E(x23),E(x122),E(x122)の表現である。
【0184】
[4次モーメント]
【数126】
JP0004350062B2_000169t.gif
ただし、f5,f6,f7,f8,f9は、X=νTξで発生される確率ベクトルX=(x,x)´のそれぞれ4次モーメントE(x14),E(x24),E(x132),E(x1222),E(x123)の表現である。ここで、σ11=σ22=1であり、また、t11=t22=1-t12=1-t21である。
【0185】
ステップS6においては、確率変数ξ1, ξ2,νがピアソン分布のどの型に属するか決定する。
【0186】
確率変数ξ1は、E(ξ1)=0,E(ξ12)=1,E(ξ13)=ζ1,E(ξ14)=κ1を満たす。同様に、確率変数ξ2は、E(ξ2)=0,E(ξ22)=1,E(ξ23)=ζ2,E(ξ24)=κ2を満たす。確率変数νは、E(ν)=0,E(ν)=1,E(ν3)=γ,E(ν4)=βを満たす。規格化された確率変数ξ1, ξ2,νの歪度はそれぞれζ1,ζ2,γであり、尖度はそれぞれκ1,κ2,βである。
【0187】
ここで、次の式(104)と表86を用い、確率変数ξ1,ξ2,νがピアソン分布のどの型になるか決定する。式(104)におけるβ1とβ2は、それぞれ歪度の2乗と尖度である。
【0188】
【数127】
JP0004350062B2_000170t.gif
【表86】
JP0004350062B2_000171t.gif
ステップS7においては、確率変数ξ1,ξ2,νのそれぞれについて、ステップS6で求めた型のピアソン分布に従う乱数を発生する。
【0189】
具体的に、ピアソンIV型以外の場合には、次の表87と表88に示すプロシージャによって乱数を発生する。
【0190】
【表87】
JP0004350062B2_000172t.gif
【表88】
JP0004350062B2_000173t.gif
ピアソンIV型の場合は、次のアルゴリズムによって乱数を発生する。
【0191】
【数128】
JP0004350062B2_000174t.gif
JP0004350062B2_000175t.gif このアルゴリズムは、ピアソンIV型分布の正規化定数を解析解又はこの展開(式(105))に基づいて計算するステップと、棄却法により乱数を発生するステップを含む。
【0192】
ステップS8においては、ステップS7において発生した乱数ξ1,ξ2,νを次の式(106)によって、標準化日本株式リターン従う乱数x1と標準化外国株式リターンに従う乱数x2に変換する。
【0193】
【数129】
JP0004350062B2_000176t.gif
ステップS9においては、ステップS8で得られた乱数x1,x2を次の式(107)と(108)によって、標準化前の日本株式リターンに従う乱数x1´と標準化前の外国株式リターンに従う乱数x2´に変換する。ここで、分散V(・)と平均E(・)には、ステップS1で標準化した際の値を用いる。
【0194】
【数130】
JP0004350062B2_000177t.gif
【数131】
JP0004350062B2_000178t.gif
ステップS10においては、ステップS9で得られた乱数x1´と乱数x2´を用いて、予想損失額(VAR)を計算する。
【0195】
日本株式と外国株式の資産の残高をA1,A2とすると、これらの重みを用いて日本株式と外国株式を併せたリターンに従う乱数を得る。次の式(109)と(110)にそれぞれ日本株式と外国株式の資産の重みを示し、日本株式と外国株式を併せたリターンに従う乱数Wを式(111)に示す。
【0196】
【数132】
JP0004350062B2_000179t.gif
【数133】
JP0004350062B2_000180t.gif
【数134】
JP0004350062B2_000181t.gif
乱数Wを用いたモンテカルロ法によるシミュレーションによって、日本株式と外国株式を併せたリターン分布をマイナス無限大から積分した領域が、全体の領域の1%になったところのリターンを用い、このリターン(通常、マイナスである。)の絶対値に日本株式と外国株式を併せた資産の残高を乗じたものが標準的な予想損失額である。
【0197】
具体的には、日本株式と外国株式を併せたリターン分布に従う乱数Wを複数発生させ、低いものから並べて下から1%目(通常、銀行規制などで使われている基準)のところのリターンを用いる。なお、予想損失額に用いるリターンは1%に限らず、所望の値を用いることができる。
【0198】
本実施例においては、日本株式と外国株式の2種類の資産を有する2次元の場合について説明したが、n種類の資産を有する一般のn次元の場合には、容易に拡張することができる。以上の方法は、損失額自体の分布によるリスク管理にも適用できる。
【0199】
この場合、ステップS5においては、式(112)において、fijk(T,ζ,γ,κ,β)とfijkl(T,ζ,γ,κ,β)には、3次と4次モーメントの表現として式(113)と(114)の関係を用いる。また、mijkとmijklは、式(115)と(116)によって定義される。具体的には、E(・)は、Xのデータ(実現値)から計算される。
【0200】
【数135】
JP0004350062B2_000182t.gif
【数136】
JP0004350062B2_000183t.gif
【数137】
JP0004350062B2_000184t.gif
【数138】
JP0004350062B2_000185t.gif
【数139】
JP0004350062B2_000186t.gif
そして、式(112)を最小にするようなパラメータζ=(ζ1,・・・,ζ),γ,κ=(κ1,・・・,κ),βを求める。
【0201】
ステップS6においては、これらのパラメータから、確率変数ν,ξ1~ξがどの型のピアソン分布か決定し、ステップS7においてその型の乱数を発生する。ステップS8においては、ステップS7において発生した乱数を乱数x1,x2・・・,xnに変換し、さらにステップS9においては、乱数x1´,x2´・・・,xn´に変換する。
【0202】
ステップS10においては、n種類の資産の残高をA1,A2,・・・,Anとすると、これらの重みを乗じたn種類の資産を併せたリターンに従う乱数を得る。次の式(117)にn種類の資産の重みを示し、式(118)にn種類の資産を併せたリターン分布に従う乱数Wを示す。
【0203】
【数140】
JP0004350062B2_000187t.gif
【数141】
JP0004350062B2_000188t.gif
このような乱数Wを用いて、n種類の資産を有する場合にも2次元の場合と同様にモンテカルロ法によるシミュレーションによって予想損失額を計算することができる。すなわち、複数の乱数Wを発生させ、低いものから並べて下から1%目のところのリターンを取り出し、このリターンの絶対値にn種類の資産の残高を乗じたものを予想損失額とする。
【0204】
さらに、本実施例においては、モーメント法に代わって最尤法によって推定したパラメータを用いて予想損失額を計算することができる。前述したように、株式リターンを用いる場合、年間取引日は250日と少数である。最尤法によると、例えば250点のような少数の標本点であっても、パラメータを正確に算出することができる。
【0205】
図2は、最尤法を用いて予想損失額を計算する一連の手順を示すフローチャートである。この一連の手順においては、ステップS5とS6の間にステップS101が加わった他は、図1に示した一連の手順と同様であるので、ステップS101を除いて説明を省略するものとする。
【0206】
ステップS101において、ステップS5においてモーメント法により計算したパラメータを初期値として、最尤法によりパラメータζ1,ζ2,γ,κ1,κ2,βの推定値を計算する。
【0207】
具体的に、ステップS101においては、次の手順によりパラメータの推定値を計算する。なお、以下このステップS101においては、パラメータT,ζ1,ζ2,γ,κ1,κ2,βの少なくとも1つをθで表記することにする。
【0208】
(1)ステップS5において、モーメント法から計算したパラメータを初期値として採用する。
【0209】
(2)多変量非正規分布から例えば10,000,000点の乱数を得る。乱数の発生は、ステップS6とS7と同様にして行う。
【0210】
(3)各次元をステップΔhによって分割に分ける。超立方体となる各分割内にある点の個数を数え、全ての点の個数で各分割内にある点の個数を除する。
【0211】
(4)次のように対数尤度を求める。
【0212】
{X1,・・・,XN}は、n次元データの集合(Xは、n次元データ・ベクトルである。)を示す。特定の多変量非正規分布のパラメータθを有するXiの分割の確率Pi(θ)は、次の式(119)によって定義される。ここで、パラメータθは、ユアンとベントラーの方法Iで用いられるパラメータT,ζ=(ζ1,・・・,ζ),γ,κ=(κ1,・・・,κ)、βの少なくとも1つを表すものとする。
【0213】
【数142】
JP0004350062B2_000189t.gif
したがって、パラメータθを有する近似した尤度fi(θ)は、次の式(120)によって定義される。
【0214】
【数143】
JP0004350062B2_000190t.gif
対数尤度の和は、次の式(121)によって計算される。
【0215】
【数144】
JP0004350062B2_000191t.gif
そして、この対数尤度の和を最大にすることにより、最尤推定量
【数148】
JP0004350062B2_000192t.gif
を得る。この最尤推定量をパラメータとして用いる。なお、尤度の積を最大にしても同様である。
【0216】
本実施の形態によると、従来のリスクをより正確に測定出来るだけでなく、不動産やオルタネティブ投資などの歪んだリターン分布を含んだポートフォリオなどのリスク管理も可能になる。また、従来、リスク管理で使われていた、リターン(収益率)ベースの方法だけでなく、損失額自体の分布を直接評価できるようになり、現物だけでなく、オプションなど様々なデリバティブを含むポートフォリオのリスク管理や、保険・天候・地震などのその他の様々なリスク(金融リスク、オペレーショナル・リスク、保険リスク、災害リスク、事業リスクなど)も統合したリスク管理が可能となる。ここでは、リスク計測対象の損失額自体の分布や収益率の分布など、リスクを計測する際の評価のための分布を総称して損益分布と呼ぶ。
【0217】
本実施の形態においては、ステップS2、S3及びS4で求めた3次及び4次モーメントm~mを用いてステップS5でパラメータζ,γ,κ,βを求め、そのパラメータを初期値として、ステップS101において前記初期値の近傍に属する最大対数尤度を与えるパラメータを求めた。しかしながら、パラメータの初期値を定めるステップS2~S5の一部または全部を省略することもできる。この場合、パラメータの初期値は未定として、ステップS101において、パラメータ空間における離散的な格子点について幅広い領域で対数尤度を計算し、最大対数尤度(又は最大尤度)を与えるパラメータを求める。このような方法は、格子探索法(グリッドサーチ)と称される。
【0218】
この方法によれば、1次から4次までのモーメントの一部又は全部が存在しない確率分布であっても、乱数さえ発生できればξやνの分布として採用可能となる。そのような場合、最尤法のための、分布を特徴付けるようなパラメータやその他あらゆる必要なパラメータを広くθと記す。
【0219】
5 変形例
5.1 金融商品のシミュレーション方法
次に、本実施の形態の第1の変形例として、本実施の形態の乱数を使った多変量非正規分布の特性を用いた金融商品のシミュレーション方法を説明する。
【0220】
この方法では、次のアルゴリズムによって発生されるピアソンIV型分布に従う乱数を使ったモンテカルロ法によるシミュレーションを用いることができる。このアルゴリズムは、ピアソンIV型分布の正規化定数の解析解又はこの展開(式(122))を計算するステップと、棄却法により乱数を発生するステップとを含む。
【0221】
【数145】
JP0004350062B2_000193t.gif
JP0004350062B2_000194t.gif この方法は、多資産の高次の相関を考慮した乱数を発生させ、実際のマーケットの特性を正確に表すリターンを発生することができる。例えば、金融商品のパフォーマンス測定などに用いられる。
【0222】
この方法は、前記多変量非正規分布(ξやνに、ピアソン分布系を使うことに限らず一般的な確率分布を使うことを含む。)を用いることにより、金融の分野における分布に特徴的ないわゆるファット・テールを正確に再現することができる。また、様々なデリバティブの対象データを本実施の形態の方法により、正確な確率分布(ξやνに、ピアソン分布系を使うことに限らず一般的な確率分布を使うことを含む。)をあてはめ(単変量分布の場合も、多変量の特殊な場合として同様に扱える。)、その分布に従うモンテカルロ・シミュレーションをすることで価格が計算できる。デリバティブの設計には、多地点での気象データや様々な指数の相関が必要で、その分析にも本実施の形態が応用される。
【0223】
例えば、天候・保険デリバティブやリアル・オプションのプライシングや設計などに応用できる。特に、プライシングに多変量非正規分布を使う例としては、コリレーティング・トリガー・オプション(複数の気象要素を選んで、それぞれに補償金の支払基準を設定するもので、オプション購入者は、契約終了時にその中で一番有利なものを選択・行使できるオプション。)のような複合オプションや複数地点オプションなどがある。
【0224】
5.2 半導体へのイオン注入
次に、本実施の形態の第2の変形例として、半導体へのイオン注入への応用を説明する。イオン注入への応用では、半導体に注入されるイオンの分布をモンテカルロ法によってシミュレーションする。
【0225】
イオンの空間分布や電荷特性は、一般に多変量非正規分布となるが、ピアソン分布を含む前記多変量非正規分布に従う乱数を用いた方法により正確に再現することができる。
【0226】
このように本実施の形態の多変量非正規分布を適用することにより、従来の解析的なモデルを、本発明の多変量非正規分布の乱数表現によって実際の分布により近い正確なモデルに改良できる。
【0227】
また、現在、イオン注入後のイオン分布を知るための物理モデルのモンテカルロ・シミュレーションは、非常に長時間かかっているが、本実施の形態の多変量非正規分布の乱数とパラメータ推定法によって、より少数の物理シミュレーションの結果から短時間にその分布の全体を正確に推定可能になり、大幅に現在のモンテカルロ・シミュレーションの短時間化を行える。
【0228】
5.3 アセット・アロケーション
本実施の形態の方法によって、多変量非正規分布に従う資産間のアセット・アロケーションを、従来の平均・分散による最適化だけでなく歪度も考慮できる最適化手法と組み合わせて実行する。この歪度まで考慮できる最適化手法については、参考文献に今野浩「理財工学II」、日科技連、1998年、コンノ,エイチ.,スズキ,ティー.及びコバヤシ,ディー.(Konno,H and Suzuki, T. and Kobayashi, D.)(1998)「平均リスク歪度ポートフォリオ問題を解くための分岐と境界アルゴリズム(A Branch and Bound Algorithm for Solving Mean-Risk-Skewness Portfolio Models)」、最適化方法及びソフトウェア(Optimization Methods and Software)、10巻、297~317頁。コンノ,エイチ.,シラカワ,エイチ.,及びヤマザキ,エイチ.,(Konno, H., Shirakawa, H., and Yamazaki, H.)(1993)「平均絶対偏差ポートフォリオ最適化モデル(A Mean-Absolute Deviation Skewness Portfolio Optimization Model)」、オペレーションズリサーチ年報(Annals of Operations Research)、45巻、205~220頁、コンノ,エイチ.及びスズキ,ケイ.(Konno,H and Suzuki, K.)(1995)「平均分散歪度ポートフォリオ最適化モデル(A Mean-Variance-Skewness Portfolio Optimization Model)」、日本オペレーションズリサーチ学会誌(Journal of the Operations Research Society of Japan)、38巻、173~187頁がある。さらに、尖度まで考慮した最適化を使えば、拡張も容易である。
【0229】
具体的には、本実施の形態の方法により、実際の資産のリターンや価格の分布を推定し、ポートフォリオとしての平均、分散、歪度、尖度を導出する。そして、平均、分散だけでなく、歪度まで考慮できる最適化手法によって、最適な資産配分を計算する。なお、歪度だけでなく尖度まで考慮できる最適化手法によって、最適な資産配分を計算することができる。
【0230】
この方法により、不動産やオルタネティブ投資などの従来の正規分布の近似では無理がある資産も含めて、正確で有効なアセット・アロケーションが可能となる。また、年金コンサルテーションや資産運用への応用としてパフォーマンスが改善するアセット・アロケーションや、多変量非正規分布の資産モデル構築も可能である。
【0231】
5.4 能力評価
例えば知能検査など能力の測定結果は、正規分布を形成しないことが知られている。例えば知能検査の場合には、1次元の場合にはピアソン分布が使われている。能力評価の統計的分析については、次の参考文献がある。バート,シー.(1963)「知能の分布は正規分布に従うか?(Is Intelligence Distributed Normally ?)」、英国統計心理学雑誌(The British Journal of Statistical Psychology)、16(XVI)巻、第2部。
【0232】
異なった能力検査や試験の結果に、本発明の実施によって、正確な多変量非正規分布が推定され、能力評価を正確に実施することができる。このことから、入試判定や人事アセスメント(選抜基準の設定や選抜有効性の検証)などをより正確化することができる。また、パーセンタイルなど様々な基準により、総合的な能力値を定義することができる。
【0233】
さらに、大学入学試験の結果と入学者の卒業成績、就職試験の結果と入社後の評価など、一部の周辺分布が大量のデータであるとき、本実施の形態の方法により、残りの少数のデータの分布を最尤法により推定することもできる。
【0234】
産業上の利用可能性
本発明によると、ピアソン分布系や一般の確率分布に従う乱数を使った多変量非正規分布を用い、経験分布を近似する分布を構築できるような多変量非正規分布に従う乱数発生方法、多変量非正規分布の推定方法、及びこれらの手順を格納した媒体を提供することができる。
【0235】
具体的には、多変量非正規分布に従う乱数発生方法及び多変量非正規分布のパラメータ推定方法は、様々な分野における統計モデリングに有効である。具体的には、統計ソフトウェアのパッケージ、リスク・マネジメントにおける予想損失額の計算への応用、金融商品のシミュレーション、アセット・アロケーションへの応用、天候・保険デリバティブやリアル・オプションのプライシングと設計、半導体へのイオン注入のシミュレーション、教育や人事アセスメントにおける能力評価に有効である。
【図面の簡単な説明】
【0236】
【図1】図1は、予想損失額を計算する一連の手順を示すフローチャートである。
【図2】図2は、最尤法を用いて予想損失額を計算する一連の手順を示すフローチャートである。
図面
【図1】
0
【図2】
1