JP6219528B2 - シミュレーションシステム、及びシミュレーション方法 - Google Patents

Description

本発明は、シミュレーションシステム、及びシミュレーション方法に関する。

興味の対象であるシステム（系）が与えられたとき、当該システムの振る舞いをモデル化し、計算機が当該モデルに基づいて当該システムの振る舞いをシミュレーションし、将来の振る舞いを予測することが行われる。例えば、ある種の物理現象や、多数の要素が絡み合う社会現象などは、単純な決定論的なモデルでは記述されることができず、複雑かつ確率的なモデルとして記述される。

決定論的なシステムは、与えられた初期条件のもとで、システムの将来の振る舞いが一つに定まる。これに対して、確率的なシステムは、システムの将来の振る舞いはある確率分布にしたがって分布しており、計算機は、例えば、注目する指標の期待値等の統計量を求める場合、シミュレーションを行う必要がある。

確率的なシステムの統計的な振る舞いを計算機が予測、シミュレーションする方法として、モンテカルロ・シミュレーション法が知られている。モンテカルロ・シミュレーション法は、対象システムに存在する確率的な振る舞いを乱数によって模擬したシミュレーションを多数回繰り返し行うことで、注目する指標の期待値を推定する方法である。このとき、シミュレーションの回数を増やすことで、大数の法則により、注目する指標の期待値が真の値に近づくと考えられる。

本技術分野の背景技術として、特開２００４−１１７２２８号公報（特許文献１）がある。特許文献１には、「気象物理量を対象とした気象庁発表の長期予報を取得するデータ取得ステップ０２０１、気象物理量の従う時系列モデルを作成する気象時系列モデル作成ステップ０２０２、推定期間における推定地点の気象物理量をシミュレーション回数だけシミュレートする気象物理量シミュレーションステップ０２０３、気象物理量のシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力ステップ０２０４を有する。さらに、気象時系列モデル作成ステップ０２０２を、気象時系列モデルのパラメータ推定ステップ０２０２１および気象時系列モデルのパラメータ補正ステップ０２０２２から構成する。」と記載されている（要約参照）。

特開２００４−１１７２２８号公報

例えば、特許文献１に記載の技術のシミュレーション対象である気象物理量のように、モンテカルロ・シミュレーションの対象である、確率的なモデルには、一般的に多数のパラメータが存在することが多い。このとき、モンテカルロ・シミュレーション実行にあたっては、あるパラメータ値の下での注目する指標の期待値だけではなく、注目指標の期待値のパラメータ値に対する感度を算出したい場合がある。注目指標の期待値のパラメータ値に対する感度とは、与えられたパラメータ値が微小に変化したときに、注目指標の期待値がどれくらい変化するかを表す微分係数である。

例えば、システムのパラメータの最適化を行う、すなわち、許容されるパラメータ値の範囲の中で、ある注目指標の期待値を最大化するパラメータ値を探す場合、与えられたパラメータ値の下での感度を得られれば、勾配法やニュートン法といった公知の最適化手法によって最適化を行うことができる。

ところが、モンテカルロ・シミュレーションにおいて、注目指標の期待値のパラメータ値に対する感度を求めることは、とくにパラメータが多数ある場合には、容易ではない。注目指標の期待値のパラメータに対する感度は、単純には、与えられたパラメータ値の下での注目指標の期待値と、パラメータをそこから微小に変化させたときの注目指標の期待値との差を、パラメータの微小変化量で割ることで算出できる。

特許文献１に記載の技術を、多数のパラメータについて注目指標の期待値の感度の推定に適用する場合、例えば特許文献１の数５、数１２のように、多数のパラメータそれぞれを一つずつ微小変化させた注目指標の期待値それぞれを算出する必要がある。すなわち、特許文献１に記載の技術を適用した感度推定手法は、当該多数パラメータの個数と同数セットのモンテカルロ・シミュレーションを行う必要がある。

そもそも、１セットのモンテカルロ・シミュレーションを行って注目指標の期待値を求めること自体が、大数の法則を満足するような多数回のシミュレーションをする必要がある。したがって、当該多数パラメータの個数と同数セットのモンテカルロ・シミュレーションを行うことは、計算機リソースや計算時間の制限により困難である。

上記課題を解決するために、本発明は、例えば、以下のような構成を採用する。複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムであって、プロセッサと、メモリと、を含み、前記メモリは、１以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、前記プロセッサは、複数の乱数ベクトルと、パラメータベクトルと、を取得し、前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、シミュレーションシステム。

本発明の一態様によれば、パラメータが存在する確率的なシステムにおいて、任意に決めた注目指標の期待値のパラメータに対する感度を、１セットのモンテカルロ・シミュレーションによって得られた結果から求めることができる。これにより、注目指標の期待値のパラメータに対する感度を求めるために必要な、計算機リソースや計算時間を削減することができる。

実施例１におけるモンテカルロ・シミュレーション・システムの構成例を表すブロック図である。 実施例１における乱数重要度算出部の構成例を表すブロック図である。 実施例１のモンテカルロ・シミュレーションの全体処理の例を表すフローチャートである。 実施例１におけるシミュレーション実行処理の一例を表すフローチャートである。 実施例１におけるシミュレーション実行処理の別の一例を表すフローチャートである。 実施例１におけるモンテカルロ・シミュレーション・システムのユーザインタフェースの一例である。 実施例２における６個のトランジスタで構成されるＳＲＡＭセルの回路の例である。 実施例３における、銀行の収支計画表の例である。 実施例４において、乱数ベクトルが２次元一様分布に従う場合における、境界残余を算出する概念の例を示す図である。

以下、本発明の実施形態について述べる。本実施形態は本発明を実現するための一例に過ぎず、本発明の技術的範囲を限定するものではないことに注意すべきである。本実施形態においては、便宜上その必要があるときは、複数のセクションまたは実施形態に分割して説明するが、特に明示した場合を除き、それらは互いに無関係なものではなく、一方は他方の一部または全部の変形例、詳細、補足説明等の関係にある。また、以下の実施形態において、要素の数等（個数、数値、量、範囲等を含む）に言及する場合、特に明示した場合および原理的に明らかに特定の数に限定される場合等を除き、その特定の数に限定されるものではなく、特定の数以上でも以下でもよい。

以下、本発明の実施形態を、添付図面を参照して説明する。なお、各図において、共通の構成には原則として同一の参照符号を付し、繰り返しの説明は省略する。

本実施例はモンテカルロ・シミュレーション・システムについて説明する。まず、本実施例のモンテカルロ・シミュレーション・システムがシミュレーションを実行し、感度を算出するための理論的側面について説明する。

確率変数Ｘは、本実施例のモンテカルロ・シミュレーション・システムによる、シミュレーション対象である確率的なシステムの挙動を表す。つまり、当該シミュレーション対象である確率的なシステムの挙動は確率変数Ｘでモデル化される。以下、シミュレーション対象である確率的なシステムを、単にシステムとも表記する。また、ｘはシステムの挙動Ｘの実現値を表す。

ここで、システムの挙動Ｘは、１次元又は多次元のベクトルのいずれであってもよく、また、「時刻」と呼ばれる変数ｔで区別される無限個の確率変数の組（離散時間確率過程と呼ばれる）であってもよい。したがって、システムの挙動Ｘが離散時間確率過程である場合、Ｘは本来Ｘ（ｔ）と表記されるべきであるが、特に必要がある場合を除いて、単にＸとも表記する。また、ある時刻ｔ_１におけるシステムの挙動を示す値ｘ（ｔ_１）を、時刻ｔ_１におけるシステムの状態と呼ぶ。

システムの挙動Ｘの分布は、システムのパラメータｚに依存し、さらにＸの確率密度関数が、ｆ（ｘ，ｚ）という形で表記されるとする。以下、システムのパラメータｚを単にパラメータｚと表記する。パラメータｚは、Ｎ次元のベクトル（すなわち、システムの挙動はＮ個のパラメータに依存している）であり、ｚの成分は、ｚ＝（ｚ_ｉ）（ただし、ｉ＝１，…，Ｎ）であるとする。また、１つのシステムの挙動の実現値Ｘ＝ｘに対して１つの実数値を定める、関数ａ（ｘ）が与えられているとする。ａ（ｘ）をシステムの挙動の評価関数と呼ぶ。システムの挙動の評価関数は、注目指標の一例である。さらに、ａ（ｘ）の期待値Ａ（ｘ）＝Ｅ［ａ（ｘ）］をシステム評価関数と呼ぶ。

実際に興味あるシステムについて検討する場合、しばしば、Ｘの確率密度関数ｆ（ｘ，ｚ）が未知であることが問題となる。しかしながら、多くの場合、Ｘの確率密度関数ｆ（ｘ，ｚ）自体は未知であっても、Ｘの分布を生み出す挙動自体は既知（又はモデル化可能）である。

この場合、モンテカルロ・シミュレーション・システムは、モンテカルロ・シミュレーションによって、システムの挙動の実現値ｘを多数生成して、Ｘの分布、すなわちｆ（ｘ，ｚ）を模擬することで、システム評価関数Ａ（ｘ）の値を近似的に算出できる。なお、生成されたシステムの挙動の実現値ｘそれぞれをパスとも呼ぶ。モンテカルロ・シミュレーションは、分布が既知である（複数の）乱数から、システムの挙動の実現値ｘを算出する計算を含む。以下、これを定式化する。

Ｗ＝（Ｗ_ｊ）をＭ次元の確率変数とし（ただし、ｊ＝１，…，Ｍ）、確率変数Ｗの実現値をｗ＝（ｗ_ｊ）と書く。また、Ｗの（同時）確率密度関数をｇ（ｗ）と書く。ここで、ｇ（ｗ）は既知であり、また計算機上でこの分布に従う乱数を容易に生成できるものとする。ｗは、例えば、［０，１］上の一様分布、又は、標準正規分布Ｎ（０，１）等に従うＭ個の独立な乱数からなるベクトルである。このとき、パラメータｚ＝（ｚ_ｊ）のもとで、乱数ベクトルｗから、対応するＸの実現値ｘを計算する方法、すなわち、関数ｘ＝ｘ（ｗ,ｚ）が与えられているとする。このとき、以下の数１が成立する。

ただし、Ωは、乱数ベクトルｗの分布範囲、つまりｇ（ｗ）のサポート（台）を表す。また、最右辺の記号ｄｗはｄｗ_１・ｄｗ_２・…・ｄｗ_Ｍの略記であり、積分記号はＭ次元の多重積分を表す。モンテカルロ・シミュレーションにおいて１つのパスをシミュレーションすることは、関数ｘ＝ｘ（ｗ,ｚ）を計算することに相当する。

さらに、Ｌ個のパスによりモンテカルロ・シミュレーションを行う場合、生成されたＬ個の乱数ベクトルをそれぞれｗ^ｋ＝（ｗ_ｊ ^ｋ）と書けば（ただし、ｊ＝１，…，Ｍ、ｋ＝１，…，Ｌ、Ｌは大数の法則を満足する数）、大数の法則から、システム評価関数は、例えば、以下の数２のように、期待値を多数（Ｌ個）のサンプルｘ（ｗ^ｋ，ｚ）の平均値で近似することにより算出される。

本実施例の目的の１つは、以下の数３で表される、システム評価関数をＮ個のパラメータ（ｚ_ｉ）それぞれで微分した微分係数を効率的に算出する方法を提供することである。当該微分係数を、システム評価関数のパラメータに対する感度と呼ぶ。

まず、本実施例の基本的なアイディアを説明する。本実施例の基本的なアイディアは、与えられたパラメータのセットｚのもとで生成したシステムの挙動（実現値）ｘを、パラメータｚ_ｉを微小変化させた後にも利用しつづけることである。多くのシステムにおいてｘ（ｗ，ｚ）は連続であり、このとき、パラメータの所定の閾値以下である微小変化Δｚ_ｉに対応して、それを打ち消すような乱数ベクトルの所定の閾値以下である微小変化Δｗが存在する場合を考える。つまり、与えられたΔｚ_ｉに対して以下の数４を満たすようなΔｗが定まる場合を考える。

このとき、ｗ＋Δｗの分布は、確率密度関数ｇ（ｗ）には従っていないため、例えば数２のような単純な平均によってＡ（ｚ＋Δｚ_ｉ）を計算することはできない。しかしながら、ｗ＋Δｗの分布がｇ（ｗ）からどれくらい異なっているのかを勘案して、ａ（ｘ（ｗ，ｚ））にそれぞれ適切な重みをつけたうえで平均をとることで、Ａ（ｚ＋Δｚ_ｉ）を算出できるというのが、本実施例の基本的なアイディアである。本実施例の手法は、重点サンプリング（Ｉｍｐｏｒｔａｎｃｅ Ｓａｍｐｌｉｎｇ）手法の一例と考えることもでき、当該重みは重要度に相当する。

具体的には、以下の数５で表される、各パスｘに対してｘを一定に保つようなΔｚ_ｉとΔｗとの関係を算出し、当該関係に基づいて重みを算出する。

数５は、ベクトル（∂ｗ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}の各要素が右辺で表されることを意味している。ただし、前述したように、ｘは無限次元の変数になりえるため、この表記は形式的なものである。この点については後述する。

実際にモンテカルロ・シミュレーションの対象であるシステムにおいては、多くの場合、ｗ、ｚからｘを計算する構成的なアルゴリズム、即ち関数ｘ＝ｘ（ｗ，ｚ）の具体的な形が分かっている。このとき、ユーザ等は数５の右辺が示示す（ｗとｚの）関数の具体的な形を解析的に算出することができる。したがって、（∂ｗ／∂ｚ_ｉ）｜_{（ｘ＝ｃｏｎｓｔ）}は、例えば、当該解析的に算出された関数に乱数ベクトルとパラメータを代入することにより算出される。また、関数ｘ（ｗ,ｚ）の具体的な形がわからないため、数５の右辺が示す関数の具体的な形を解析的に算出することはできない場合でも、パラメータの微小変化Δｚ_ｉに対して、ｘを一定に保つような、Δｗを求めることにより、（∂ｗ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}＝Δｗ／Δｚ_ｉと近似的に計算することができる。

以下、上記したアイディアを具体的に定式化する。Δｚ_ｉをｉ番目のパラメータの微小変化として、パラメータがｚ＋Δｚ_ｉ＝（ｚ_１,ｚ_２,…,ｚ_ｉ＋Δｚ_ｉ,…,ｚ_Ｎ）であるときの、システム評価関数Ａ（ｚ＋Δｚ_ｉ）は、以下の数６のように表される。

ただし、（∂ｗ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}は、前述したようにｘ＝ｘ（ｗ,ｚ）＝ｘ（ｗ＋Δｗ,ｚ＋Δｚ_ｉ）を満たす、パラメータの微小変化Δｚ_ｉと、乱数ベクトルの微小変化Δｗとの関係を表す。また、Ｉは単位行列を、ｔｒは行列の跡（ｔｒａｃｅ）を、中点・は内積を表す。｜ｄｗ’／ｄｗ｜はｗからｗ’への変数変換によるヤコビ行列式を、Ω’は変数変換によるΩの像を表す。したがって、以下の数７が成立する。

ただし、数７の積分が存在するものとする。ここで、以下の数８で定義されるＮ変数関数ｈ（ｗ，ｚ）＝（ｈ_ｉ（ｗ,ｚ））を（ただし、ｉ＝１，…，Ｎ）、等価微分重要度と呼ぶ。

数７より、システム評価関数のパラメータに対する微分係数、すなわち感度は、等価微分重要度によって重みづけされたシステムの挙動の評価値、の期待値に基づいて算出される。また、数７中のＲ_ｉは、数７の積分範囲をΩ’ではなくΩとしたことにともなう差し引きを表す項であり，以下の数９のように表される。

Ｒ_ｉを境界残余と呼ぶ。ただし∂ΩはΩの境界（表面）を、ｎは境界∂Ωの外向き単位法線ベクトルを表す。数９より、Ωの境界（表面）∂Ωで（∂ｗ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}＝０（ベクトル）であれば、Ｒ_ｉ＝０である。この場合、数７により、以下の数１０が成立する。

ここで、代表的な乱数ベクトルｗの分布における等価微分重要度ｈ_ｉ（ｗ,ｚ）について述べる。ｗが一様分布の乱数ベクトルであるとき、ｈ_ｉ（ｗ，ｚ）は以下の数１１で表される。

また、ｗがＭ次元独立標準正規分布の乱数ベクトルであるとき，ｈ_ｉ（ｗ，ｚ）は以下の数１２で表される。

また、モンテカルロ・シミュレーションにより、数７が示す、システム評価関数Ａ（ｚ）のパラメータｚに対する微分係数は、以下の数１３のように近似的に算出される。

とくに，Ωの境界（表面）∂Ωで（∂ｗ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}＝０の場合は、数１０より、システム評価関数Ａ（ｚ）のパラメータｚに対する感度は、以下の数１４で表される。

数２、数１３、数１４より、初期パラメータｚについてのみ、モンテカルロ・シミュレーションを１セット行って、Ｌ個のパス｛ｘ（ｗ^ｋ，ｚ）｝と、Ｌ個のパスそれぞれにおける等価微分重要度｛ｈ（ｗ^ｋ,ｚ）｝を計算すれば（ただし、ｋ＝１，…，Ｌ）、ｚにおけるシステム評価関数Ａ（ｚ）と同時に、システム評価関数Ａ（ｚ）の全てのパラメータに対する微分係数を計算できる。すなわち、本実施例のモンテカルロ・シミュレーション・システムは、微小変化後の各パラメータｚ＋Δｚ_ｉについてのモンテカルロ・シミュレーションを行う必要がない。

さらに、数２および数１４より、システムの挙動の評価関数ａ（ｘ）はモンテカルロ・シミュレーション終了後に自由に設定できることがわかる。例えば、複雑なシステムをシミュレーション対象とする場合、実際にシステムの挙動をシミュレーションしてみてはじめて何について評価するべきか、どのような現象に注目するべきかが分かる、ということがしばしばあるため、この特徴は望ましい。

また，数１３、数１４は、通常の期待値計算の形をしているため、Ｌ回のシミュレーションによる微分係数の推定値の分散は、Ｏ（Ｌ^−１）で減少する。一方、パラメータの数だけ単純にモンテカルロ・シミュレーションをやり直す方法による微分係数の推定値の分散は、シミュレーション回数に対してＯ（Ｌ^−１/４）、又はＯ（Ｌ^−１/３）でしか減少しない。つまり、本実施例のシミュレーションによる微分係数の推定値の分散の収束のスピードが速い。

なお、数９で定義される境界残余Ｒ_ｉは、ｘ（ｗ,ｚ）＝ｘ（ｗ＋Δｗ,ｚ＋Δｚ_ｉ）を満たすようなｗ＋Δｗの分布範囲が、ｗの分布範囲（すなわち、ｗが従う確率密度関数の定義域）Ωとずれてしまうことで生じる。つまり、出現する可能性がある全ての乱数ベクトルの実現値ｗ∈Ωに対して、パラメータの微小変化Δｚ_ｉを打ち消すような、乱数ベクトルの実現値ｗ＋Δｗを考えたときに、ｗ＋Δｗの分布範囲Ω’＝｛ｗ＋Δｗ；ｗ∈Ω｝がΩと一致していれば境界残余Ｒ_ｉは生じない。

モンテカルロ・シミュレーションの対象となる実際のシステムを想定すると、システムのどのような挙動ｘについても、パラメータｚ_ｉが微小変化したときに、それを打ち消すようなノイズ（乱数ベクトルの実現値）が存在する場合が多い。この場合には、境界残余Ｒ_ｉは生じない。例えば、正規分布に従う乱数ベクトルのように、ｗの分布範囲Ωが実数全体にわたっている場合には、Ωに境界が存在しないため境界残余Ｒ_ｉはゼロである。

このとき、モンテカルロ・シミュレーション・システムは、境界残余Ｒ_ｉを算出しなくてもよいため、計算時間及び計算リソースを削減することができる。以下、特に断らない限り乱数ベクトルｗの分布範囲ΩはΩ＝Ω’を満たす、即ち境界残余Ｒ_ｉがゼロであるものとする。

ここで、本実施例の感度計算手法が適用できる条件の例を説明する。前述の基本的なアイディアの説明において、パラメータの微小変化Δｚ_ｉに対してｘを一定に保つようなΔｗを考える、と書いたが、実際には、与えられたΔｚ_ｉに対して数４を満たすようなΔｗが存在するとは限らない。

前述したように、例えば、ｘは無限次元のベクトルにもなりえる、つまりｘの次元はｗの次元Ｍよりも大きくなりえる。この場合、数４において、制約条件の数が変数Δｗの数（次元数）よりも多くなるため、数４を満たすΔｗは一般には存在しない。このように、与えられたパラメータの微小変化Δｚ_ｉに対して、数４を満たす乱数ベクトルの微小変化Δｗが存在しない場合には本実施例の手法は適用できない。

以下、本実施例の手法が適用可能、すなわち、パラメータの微小変化Δｚ_ｉに対して数４を満たすような乱数ベクトルの微小変化Δｗが存在する例を説明する。

システムの挙動を表す時系列ｘ（ｔ）が、例えば、以下の数１５のように表せる場合、本実施例の手法が適用可能である。

数１５は、乱数ベクトルｗとパラメータｚで定まる値ｙ（ｗ,ｚ）が定まると、時刻ｔでのシステムの状態ｘ（ｔ）が、時刻ｔとｙから定まることを示す。なお、この場合、明らかに、以下の数１６が成立する。

また、システムの挙動を表す時系列ｘ（ｔ）が、例えば、以下の数１７のような漸化式で記述される場合、本実施例の手法が適用可能である。

数１７は、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）が、時刻ｔと、過去（時刻ｔ以前）のシステムの挙動と、時刻ｔで新たに生成され、時刻毎に異なる乱数ｗ_ｊ（ｗ_ｊは複数個の乱数を含む乱数ベクトルでもよい）と、パラメータｚと、から定まることを示す。なお、この場合、以下の数１８が成立する。

数１８は、パラメータの微小変化Δｚ_ｉを打ち消す乱数の微小変化Δｗ_ｊが、初期時刻ｔ＝０からシミュレーションの終了時刻ｔ＝Ｔまで順に定まることを示す。

以下、本実施例のモンテカルロ・シミュレーション・システムがシミュレーションを実行し、感度を算出する実践的側面について説明する。図１は、本実施例のモンテカルロ・シミュレーション・システム１００の構成例を示す。モンテカルロ・シミュレーション・システム１００は、例えば、ＣＰＵ１１０、メモリ１２０、入出力インタフェース１３０、二次記憶装置１４０を含む計算機上に構成される。また、図１において、モンテカルロ・シミュレーション・システムは１つの計算機から構成されているが、複数の計算機から構成されてもよい。

ＣＰＵ１１０は、プログラムに従って動作するプロセッサ及び／又は論理回路を含み、データの入力／出力、読み込み／書き込みを行い、さらに、後述する各プログラムを実行する。メモリ１２０は、ＣＰＵ１１０が実行するプログラム及びデータを一時的にロードして記憶し、さらに各プログラム及び各データを保持する。入出力インタフェース１３０は、外部装置等からデータ等の入力、及び外部装置等に対してデータ等の出力を行うインタフェースである。二次記憶装置１４０は、不揮発性の記憶装置であり、例えばＨＤＤ等である。プログラムやデータ等は二次記憶装置１４０に格納されてもよい。

プログラムはＣＰＵ１１０によって実行されることで、定められた処理をメモリ１２０及び入出力インタフェース１３０を用いながら行う。従って、本実施例及び他の実施例においてプログラムを主語とする説明は、ＣＰＵ１１０を主語とした説明でもよい。若しくは、プログラムが実行する処理は、そのプログラムが動作する計算機及び計算機システムが行う処理である。

ＣＰＵ１１０は、プログラムに従って動作することによって、所定の機能を実現する機能部として動作する。例えば、ＣＰＵ１１０は、後述する乱数発生部２１０に従って動作することで乱数発生部として機能し、後述するモデル挙動更新部２２０に従って動作することでモデル挙動更新部として機能する。さらに、ＣＰＵ１１０は、各プログラムが実行する複数の処理のそれぞれを実現する機能部としても動作する。計算機及び計算機システムは、これらの機能部を含む装置及びシステムである。プログラムは二次記憶装置１４０に格納されていてもよい。

メモリ１２０は、１又は複数のシナリオ・シミュレータ２００と、記憶部３００と、後処理計算部４００と、モンテカルロ・シミュレーション全体制御部５００と、を含む。シナリオ・シミュレータ２００は、異なる乱数ベクトルを用いてＬ回のモンテカルロ・シミュレーションを行う。記憶部３００は、シナリオ・シミュレータ２００から出力された数値等を記憶する。後処理計算部４００は、シミュレーション終了後に期待値および感度の計算を行う。

モンテカルロ・シミュレーション全体制御部５００は、モンテカルロ・シミュレーション・システム１００全体の制御を行う。また、モンテカルロ・シミュレーション全体制御部５００は、ユーザ等が指定したシミュレーション回数によって算出されるシステム評価関数及び感度の推定値の計算精度を算出し、出力してもよい。また、モンテカルロ・シミュレーション全体制御部５００は、ユーザ等が指定した、システム評価関数及び感度の推定値の計算精度を保証するために必要なシミュレーション回数等を算出、出力してもよい。モンテカルロ・シミュレーション全体制御部５００は、例えば、指定された又は算出したシミュレーション回数、例えばＬ回、シナリオ・シミュレータ２００を起動し、Ｌ回のシミュレーションが終了すると、後処理計算部４００を起動する。

シナリオ・シミュレータ２００は、ユーザ等から与えられるシミュレーションモデルの１次元又は多次元のパラメータのベクトルｚを入力として、シミュレーションに必要な乱数ベクトルを内部で生成しながらモデルのシミュレーションを実行する。また、シナリオ・シミュレータ２００は、システムの挙動の評価関数ａ^（ｋ）と、感度計算に必要な等価微分重要度ｈ_ｉ ^（ｋ）を出力する。

シナリオ・シミュレータ２００は、それぞれプログラムである、乱数発生部２１０、モデル挙動更新部２２０、乱数重要度算出部２３０、指標算出部２６０、等価微分重要度算出部２７０、及びシミュレーション制御部２８０を含む。また、シナリオ・シミュレータ２００は、それぞれデータを格納する領域である、モデル挙動記憶部２４０、及び乱数重要度記憶部２５０を含む。

乱数発生部２１０は、予め定められた、又はユーザ等によって指定されたシミュレーションモデルに対応する確率密度関数に従う乱数ベクトルを生成する。モデル挙動更新部２２０は、ユーザ等によって与えられたパラメータと乱数発生部２１０が生成した乱数ベクトルとを入力として、入力されたパラメータと乱数ベクトルとを予め定められた、又はユーザ等によって指定されたシミュレーションモデル（関数）に適用し、システムの挙動の実現値をシミュレートする。

乱数重要度算出部２３０は、パラメータと乱数発生部２１０が生成した乱数ベクトルとを入力として、入力された各パラメータについて、ｗ_ｊ ^（ｋ）の各パラメータに関する乱数重要度、すなわち数８の総和記号の中の項を算出する。モデル挙動記憶部２４０は、モデル挙動更新部２２０が算出したシステムの挙動の実現値を格納する。乱数重要度記憶部２５０は、乱数重要度算出部２３０が算出した乱数重要度を格納する。指標算出部２６０は、例えば、モデル挙動記憶部２４０に格納されたシステムの挙動の実現値の平均を算出することにより、システムの挙動の評価関数を算出する。

等価微分重要度算出部２７０は、例えば、乱数重要度記憶部２５０に格納された乱数重要度の和を算出することにより等価微分重要度を算出する。シミュレーション制御部２８０は、シナリオ・シミュレータ２００全体の動作を制御する。

記憶部３００は、指標算出部２６０が算出したシステムの挙動の評価関数ａ^（ｋ）を格納する指標値記憶部３１０と、等価微分重要度算出部２７０が算出した等価微分重要度を格納する等価微分重要度記憶部３２０と、を含む。

後処理計算部４００は、それぞれプログラムである、期待値計算部４１０、及び感度計算部４２０を含む。期待値計算部４１０は、指標値記憶部３１０に格納されたＬ個のシステムの挙動の評価関数ａ^（ｋ）の平均値を計算し、当該平均値を期待値、すなわちシステム評価関数Ａ（ｚ）として出力する。このとき、期待値計算部４１０は、さらにＬ個のシステムの挙動の評価関数ａ^（ｋ）の分散を計算し、Ａ（ｚ）の推定値の精度（信頼区間）を同時に出力してもよい。

感度計算部４２０は、指標値記憶部３１０に格納されたシステムの挙動の評価関数ａ^（ｋ）と、等価微分重要度記憶部３２０に格納されたＮ個のパラメータそれぞれに関する等価微分重要度ｈ_ｉ ^（ｋ）の組について、それぞれ対応する値同士の積を計算する。さらに、感度計算部４２０は、算出した積の平均値を計算して、Ｎ個のパラメータそれぞれに関するシステム評価関数Ａ（ｚ）の感度∂Ａ／∂ｚ_ｉとして出力する。このとき、感度計算部４２０は、さらにａ^（ｋ）とｈ_ｉ ^（ｋ）の積の分散を計算し、感度∂Ａ／∂ｚ_ｉの推定値の精度（信頼区間）を同時に出力してもよい。

図２は、乱数重要度算出部２３０の構成例を示す。乱数重要度算出部２３０は、いずれもプログラムである、サンプルパス固定微分計算部２３１、サンプルパス固定２階微分計算部２３２、及び乱数重要度計算部２３３を含む。

サンプルパス固定微分計算部２３１は、パラメータと、乱数発生部２１０が生成した乱数ベクトルと、を入力として（∂ｗ_ｊ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}を算出する。サンプルパス固定２階微分計算部２３２は、パラメータと、乱数発生部２１０が生成した乱数ベクトルと、を入力として∂Ｄ_ｉ，ｊ／∂ｗ_ｊを算出する（Ｄ_ｉ，ｊ ^（ｋ）＝（∂ｗ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}）。サンプルパス固定２階微分計算部２３２は、サンプルパス固定微分計算部２３１が算出した（∂ｗ_ｊ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}と、パラメータと、を入力として、∂Ｄ_ｉ，ｊ／∂ｗ_ｊを算出してもよい。また、サンプルパス固定２階微分計算部２３２が（∂ｗ_ｊ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}を算出してもよく、このとき、サンプルパス固定微分計算部２３１は（∂ｗ_ｊ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}を算出しなくてもよい。

乱数重要度計算部２３３は、乱数発生部２１０から入力された乱数ベクトルと、サンプルパス固定微分計算部２３１又はサンプルパス固定２階微分計算部２３２が算出した（∂ｗ_ｊ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}と、サンプルパス固定２階微分計算部２３２が算出した∂Ｄ_ｉ，ｊ／∂ｗ_ｊと、を入力として、各パラメータについて、ｗ_ｊ ^（ｋ）の各パラメータに関する乱数重要度、すなわち数８の総和記号の中の項を算出する。

図３は、本実施例におけるモンテカルロ・シミュレーションの全体処理の一例を示す。モンテカルロ・シミュレーションの全体処理は、シナリオ・シミュレータ２００が異なる乱数ベクトルを用いてＬ回のシミュレーションを行うシミュレーション実行処理と、シミュレーション終了後に後処理計算部４００が行う期待値・感度計算処理と、を含む。

シナリオ・シミュレータ２００は、シミュレーション実行処理を行う（Ｓ３０１）。シミュレーション実行処理は、１セットのモンテカルロ・シミュレーションを構成するＬ回のシミュレーションを実行する処理である。シナリオ・シミュレータ２００は、異なる乱数ベクトルを用いて個々のシミュレーションを実行する。

なお、本実施例のモンテカルロ・シミュレーション・システムは、ｋ回目のシミュレーションにより得られたシステムの挙動の実現値ｘ^（ｋ）から、システムの挙動の評価関数ａ^（ｋ）に加えて、Ｎ個のパラメータそれぞれについて等価微分重要度ｈ_ｉ ^（ｋ）を計算する。なお、シミュレーション実行処理におけるＬ回のシミュレーションそれぞれは、互いに独立であるため、Ｌ回のシミュレーションを逐次的に行うのではなく、例えば、多数のシナリオ・シミュレータ２００又は多数の計算機を用いて同時に並列実行してもよい。等価微分重要度ｈ_ｉ ^（ｋ）を算出する処理を定める一連の関数を、重み関数と呼ぶ。

続いて、後処理計算部４００は、期待値・感度計算処理を行う（Ｓ３０２）。期待値・感度計算処理は、１セットのモンテカルロ・シミュレーションを構成するＬ回のシミュレーションの実行後に、後処理計算部４００が、注目指標の期待値の推定値と、注目指標の期待値のパラメータに対する感度の推定値と、を計算する処理である。

注目指標の期待値、すなわちシステム評価関数Ａ（ｚ）の推定値は、期待値計算部４１０が、指標値記憶部に格納されたＬ個のシステムの挙動の評価関数ａ^（ｋ）の平均を計算することで得られる。本実施例におけるモンテカルロ・シミュレーション・システムは、これに加えて、注目指標の期待値のＮ個のパラメータそれぞれに対する感度の推定値を計算する。

具体的には、感度計算部４２０は、システム評価関数Ａ（ｚ）のｉ番目のパラメータそれぞれに対する感度の推定値を、指標値記憶部３１０に格納された対応するシステムの挙動の評価関数ａ^（ｋ）に、等価微分重要度記憶部３２０に格納されたｉ番目のパラメータについての等価微分重要度ｈ_ｉ ^（ｋ）を掛けることで重みづけをしたものの平均を計算することで得る。なお、モンテカルロ・シミュレーション・システムは、例えば、ユーザ等が指定したパラメータについてのみ、ステップＳ３０１における等価微分重要度の算出処理、及びステップＳ３０２における感度の推定値の算出処理を行ってもよい。

なお、前述のように、例えば、システムの挙動を表す時系列ｘ（ｔ）が、数１５又は数１７の形で表される場合、本実施例における注目指標の期待値のパラメータに対する感度の推定手法が適用できる。図３、図４は、それぞれ、システムの挙動を表す時系列ｘ（ｔ）が、数１５、数１７の形で表されたモデルで算出される場合における、ｋ回目のシミュレーション実行処理の詳細な処理内容を示す。

図４は、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、すなわち、乱数ベクトルｗとパラメータｚで定まる値ｙ（ｗ,ｚ）を定めると、時刻ｔでのシステムの状態ｘ（ｔ）が，時刻ｔとｙのみで決定論的に定まる場合について、シミュレーション実行処理の詳細な処理内容の一例を示す。

まず、シナリオ・シミュレータ２００は、乱数ベクトルｗ^（ｋ）を生成し、ｗ^（ｋ）とパラメータｚとからｙ^（ｋ）を計算する（Ｓ４０１）。具体的には、乱数発生部２１０は、予め定められた確率密度関数ｇ（ｗ）にしたがうＭ次元の乱数ベクトルｗ^（ｋ）を生成し、生成した乱数をモデル挙動更新部２２０及び乱数重要度算出部２３０に送信する。なお、乱数発生部２１０は、例えば、既知の疑似乱数発生アルゴリズム、又は物理乱数発生装置といった公知の手段を用いて乱数ベクトルを生成すればよい。モデル挙動更新部２２０は、受信したｗ^（ｋ）と、与えられたＮ次元のパラメータ値ｚを、予め与えられた関係式ｙ（ｗ，ｚ）に代入してｙ^（ｋ）を計算する。

続いて、シナリオ・シミュレータ２００は、Ｎ個のパラメータそれぞれについて等価微分重要度ｈ_ｉ ^（ｋ）を計算する（Ｓ４０２）。これは、注目指標の期待値のパラメータに対する感度の推定のために必要な処理であり、重み関数が定める処理は当該処理を含む。なお、シナリオ・シミュレータ２００は、ユーザ等が指定したパラメータについてのみ、ステップＳ４０２における処理を行ってもよい。ｉ番目のパラメータについての等価微分重要度ｈ_ｉ ^（ｋ）は、具体的には以下の４つのステップによって算出される、すなわちステップＳ４０２は以下の４つのステップを含む。

第１のステップは、サンプルパス固定微分計算部２３１が、システムの挙動ｘを一定に保つような、乱数ベクトルｗ^（ｋ）のｊ番目の要素とｉ番目のパラメータの微小変化の関係Ｄ_ｉ，ｊ ^（ｋ）＝（∂ｗ_ｊ ^（ｋ）／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}を算出するステップである。ここで、添え字ｊは、Ｍ次元の乱数ベクトルｗ^（ｋ）のｊ番目の要素についての計算を表す。

システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、数１６が成り立つから、サンプルパス固定微分計算部２３１は、（∂ｗ_ｊ／∂ｚ_ｉ）｜_{（ｙ＾（ｋ）＝ｃｏｎｓｔ）}を求めればよい。つまり、サンプルパス固定微分計算部２３１は、例えば、関数−（∂ｙ^（ｋ）／∂ｚ_ｉ）／（∂ｙ^（ｋ）／∂ｗ_ｊ ^（ｋ））がユーザ等によって与えられた場合、入力されたパラメータｚ_ｉ及び入力された乱数ベクトルｗ_ｊ ^（ｋ）を当該関数に代入することにより、Ｄ_ｉ，ｊ ^（ｋ）を算出すればよい。サンプルパス固定微分計算部２３１は、当該関数から高速にＤ_ｉ，ｊ ^（ｋ）を算出することができる。

また、サンプルパス固定微分計算部２３１は、例えば、パラメータの微小変化Δｚ_ｉに対して、ｙを一定に保つようなΔｗを数値的に探索してΔｗ／Δｚ_ｉを計算することにより、（∂ｗ_ｊ／∂ｚ_ｉ）｜_{（ｙ＾（ｋ）＝ｃｏｎｓｔ）}、すなわちＤ_ｉ，ｊ ^（ｋ）を近似的に算出してもよい。例えば、ユーザ等が関数ｘ及びｙの具体的な形がわからない場合であっても、サンプルパス固定微分計算部２３１は、上述の方法により、Ｄ_ｉ，ｊ ^（ｋ）を算出することができる。サンプルパス固定微分計算部２３１は算出したＤ_ｉ，ｊ ^（ｋ）を、サンプルパス固定２階微分計算部２３２及び乱数重要度計算部２３３に送信する。

第２のステップは、サンプルパス固定２階微分計算部２３２が、受信したＤ_ｉ，ｊ ^（ｋ）を、ｗ_ｊ ^（ｋ）で微分するステップである。Ｄ_ｉ，ｊ ^（ｋ）を、ｗ_ｊ ^（ｋ）で微分した結果をＤＤ_ｉ，ｊ ^（ｋ）と表記する。サンプルパス固定２階微分計算部２３２は、例えば、受信したＤ_ｉ，ｊ ^（ｋ）を数値的に微分することによって、ＤＤ_ｉ，ｊ ^（ｋ）を算出すればよい。また、サンプルパス固定２階微分計算部２３２は、例えば、ＤＤ_ｉ，ｊ ^（ｋ）を算出する関数が与えられている場合、入力されたパラメータｚ_ｉ及び入力された乱数ベクトルｗ_ｊ ^（ｋ）当該関数に代入することによって、ＤＤ_ｉ，ｊ ^（ｋ）を算出してもよい。サンプルパス固定２階微分計算部２３２は、算出したＤＤ_ｉ，ｊ ^（ｋ）を、乱数重要度計算部２３３に送信する。

第３のステップは、乱数重要度計算部２３３が、乱数ベクトルｗ^（ｋ）のｊ番目の要素ｗ_ｊ ^（ｋ）に関する乱数重要度Ｈ_ｉ，ｊ ^（ｋ）を計算するステップである。第３のステップは、数８における総和記号の中の項を計算することに相当する。具体的には、例えば、乱数重要度計算部２３３は、受信したＤ_ｉ，ｊ ^（ｋ）と、受信したＤＤ_ｉ，ｊ ^（ｋ）と、予め定められた確率密度関数ｇ（ｗ）から一意に定まる関数ＬＧに乱数ベクトルｗ_ｊ ^（ｋ）を代入した値ＬＧ_j（ｗ_ｊ ^（ｋ））＝（１／ｇ（ｗ_ｊ ^（ｋ）））×（∂ｇ／∂ｗ_ｊ ^（ｋ））と、を図中に示した関係式に代入して乱数重要度Ｈ_ｉ，ｊ ^（ｋ）を計算すればよい。乱数重要度計算部２３３は、算出したＨ_ｉ，ｊ ^（ｋ）を乱数重要度記憶部２５０に格納する。

なお、前述のように、第３のステップにおける計算は、Ｍ次元のｗが一様分布の乱数ベクトルである場合、又は、ｗがＭ次元独立標準正規分布の乱数ベクトルである場合には、それぞれ数１１又は数１２のように表される。従って、例えば、乱数ｗ_ｊ ^（ｋ）が一様分布に従う場合、数１１より、乱数重要度はＤＤ_ｉ，ｊ ^（ｋ）で表されるため、乱数重要度算出部２３０は、第１のステップのみを行うことにより乱数重要度を算出することができる。シナリオ・シミュレータ２００は、第１のステップから第３のステップにおける処理を、全てのｊ（ただし、ｊ＝１…Ｍ）について行う。

第４のステップは、等価微分重要度算出部２７０がｉ番目のパラメータに関する等価微分重要度ｈ_ｉ ^（ｋ）を計算するステップである。第４のステップは、前述の数８における総和計算に相当する。等価微分重要度算出部２７０は、乱数重要度記憶部２５０に格納された乱数重要度Ｈ_ｉ，ｊ ^（ｋ）を取得し、ｊについての総和を計算する。なお、各パラメータｚ_ｉについての等価微分重要度ｈ_ｉ ^（ｋ）の計算は互いに独立であるため、１つのシナリオ・シミュレータ２００によってＮ個のパラメータｚ_１，…，ｚ_Nについて等価微分重要度の計算を逐次的に行うのではなく、例えば、多数の計算機を用いて同時に並列に当該計算を実行してもよい。

続いて、モデル挙動更新部２２０は、は、実際にシステムの挙動ｘ^（ｋ）（ｔ）をシミュレートし計算する（Ｓ４０３）。モデル挙動更新部２２０は、時刻ｔを初期時刻から終了時刻に向かって増加させながら、時刻ｔでのシステムの状態を逐次的に計算することで当該処理を行うことができるが、これによらずその他の公知の効率的なシミュレーションの計算手法を用いてもよい。また、モデル挙動更新部２２０は算出したｘ^（ｋ）（ｔ）をモデル挙動記憶部２４０に格納する。なお、ステップＳ４０３における処理は、例えば、ステップＳ４０２における処理と並行して行われてもよいし、ステップＳ４０２における処理の前に行われてもよい。

続いて、指標算出部２６０は、モデル挙動記憶部２４０から、システムの挙動ｘ^（ｋ）（ｔ）を取得し、注目指標の評価値、すなわち、システムの挙動の評価関数ａ^（ｋ）を計算する（Ｓ４０４）。指標算出部２６０は、算出したａ^（ｋ）を指標値記憶部３１０に格納する。

以上のように、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、図４に示した処理により、シナリオ・シミュレータ２００はシステム挙動の評価関数ａ^（ｋ）、およびＮ個のパラメータそれぞれについての等価微分重要度ｈ_ｉ ^（ｋ）を計算することができる。

図５は、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される場合、すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）が、時刻ｔと、過去（時刻ｔ以前）のシステムの状態と、時刻ｔで新たに生成する（時刻毎に異なる）乱数ｗ_ｊ（複数個の乱数を含む乱数ベクトルもよい）と、パラメータｚと、によって定まる場合の、シミュレーション実行処理の詳細な処理内容を示す。

シナリオ・シミュレータ２００は、時刻ｔを、シミュレーション開始時刻ｔ＿ｓｔａｒｔから、シミュレーション終了時刻ｔ＿ｅｎｄまで、シミュレーション時間ステップΔｔずつ増加させながら、時刻ｔにおけるシステムの状態ｘ^（ｋ）（ｔ）と、Ｎ個のパラメータ全てについてｊ番目の乱数ｗ_ｊ ^（ｋ）に関する乱数重要度Ｈ_ｉ，ｊ ^（ｋ）と、を計算する処理を繰り返し行うことでシミュレーションを行う。さらに、シナリオ・シミュレータ２００は、シミュレーションの結果得られたシステムの挙動ｘ^（ｋ）（ｔ）の時系列を用いて、システムの挙動の評価関数ａ^（ｋ）、およびＮ個のパラメータそれぞれについての等価微分重要度ｈ_ｉ ^（ｋ）を計算する。

まず、シミュレーション制御部２８０は、時刻ｔをシミュレーション開始時刻ｔ＿ｓｔａｒｔにセットし、時刻ｔに関するループ処理の開始をシナリオ・シミュレータ２００の各部に通知する（Ｓ５０１）。時刻ｔに関するループ処理は、モデル挙動更新部２２０が時刻ｔでのシステムの状態ｘ^（ｋ）（ｔ）を計算する処理（Ｓ５０２）と、乱数重要度算出部２３０が時刻ｔで生成した乱数のＮ個のパラメータそれぞれに関して乱数重要度Ｈ_ｉ，ｊ ^（ｋ）を計算する処理（Ｓ５０３）と、を含む。

ステップＳ５０２における具体的な処理について説明する。まず、乱数発生部２１０は、時刻ｔでのシステムの状態ｘ^（ｋ）（ｔ）の計算に必要な乱数ｗ_ｊ ^（ｋ）を、予め定められた確率密度関数ｇ（ｗ）にしたがって生成し、モデル挙動更新部２２０、及び乱数重要度算出部２３０に生成した乱数を送信する。モデル挙動更新部２２０は、受信した乱数ｗ_ｊ ^（ｋ）と過去（時刻ｔ以前）のシステムの挙動とパラメータｚと時刻ｔとを与えられた漸化式（数１７）に代入して、時刻ｔでのシステムの状態ｘ^（ｋ）（ｔ）を計算する。

続いて、乱数重要度算出部２３０は、Ｎ個のパラメータｚ_１，…，ｚ_Nそれぞれについて、時刻ｔで生成した乱数ｗ_ｊ ^（ｋ）のパラメータに関する乱数重要度Ｈ_ｉ，ｊ ^（ｋ）を計算する（Ｓ５０３）。なお、シナリオ・シミュレータ２００は、ユーザ等が指定したパラメータについてのみ、ステップＳ５０３における処理を行ってもよい。

ｉ番目のパラメータについての乱数重要度Ｈ_ｉ、ｊ ^（ｋ）は、具体的には以下の３つのステップにより算出される、すなわちステップＳ５０３は以下の３つのステップを含む。なお、各パラメータについての等価微分重要度ｈ_ｉ ^（ｋ）の計算は互いに独立であるため、１つのシナリオ・シミュレータ２００によってＮ個のパラメータｚ_１，…，ｚ_Nについて等価微分重要度の計算を逐次的に行うのではなく、例えば多数の計算機を用いて当該計算を同時に並列実行してもよい。

第１のステップは、サンプルパス固定微分計算部２３１が、システムの挙動ｘを一定に保つような、乱数ｗ_ｊ ^（ｋ）とｉ番目のパラメータの微小変化の関係Ｄ_ｉ，ｊ ^（ｋ）＝（∂ｗ_ｊ ^（ｋ）／∂ｚ_ｉ）｜_{ｘ＿ｊ＋１＾（ｋ）＝ｃｏｎｓｔ}を求めるステップである。

サンプルパス固定微分計算部２３１は、例えば、関数−（∂ｆ／∂ｚ_ｉ）／（∂ｆ／∂ｗ_ｊ ^（ｋ））が与えられた場合、入力されたパラメータｚ_ｉ及び入力された乱数ｗ_ｊ ^（ｋ）を当該関数に代入することにより、Ｄ_ｉ，ｊ ^（ｋ）を算出すればよい。また、サンプルパス固定微分計算部２３１は、例えば、パラメータの微小変化Δｚ_ｉに対して、ｆを一定に保つようなΔｗを数値的に探索してΔｗ／Δｚ_ｉを計算することにより、（∂ｗ_ｊ／∂ｚ_ｉ）｜_{（x(t)＝ｃｏｎｓｔ）}、すなわちＤ_ｉ，ｊ ^（ｋ）を近似的に算出してもよい。サンプルパス固定微分計算部２３１は、算出したＤ_ｉ，ｊ ^（ｋ）を、サンプルパス固定２階微分計算部２３２、及び乱数重要度計算部２３３に送信する。

第２のステップは、サンプルパス固定２階微分計算部２３２が、受信したＤ_ｉ，ｊ ^（ｋ）を、ｗ_ｊ ^（ｋ）で微分し、ＤＤ_ｉ，ｊ ^（ｋ）を算出するステップである。第２のステップは、ステップＳ４０２における第２のステップと同様である。

第３のステップは、乱数重要度計算部２３３が、時刻ｔにおける乱数ｗ_ｊ ^（ｋ）に関する乱数重要度Ｈ_ｉ，ｊ ^（ｋ）を計算するステップである。第３のステップは、ステップＳ４０２における第３のステップと同様である。

続いて、シミュレーション制御部２８０は、時刻をｔ＋Δｔにセットする（Ｓ５０４）。時刻ｔ＋Δｔがシミュレーション終了時刻ｔ＿ｅｎｄ以前である場合（Ｓ５０５：Ｙｅｓ）、モンテカルロ・シミュレーション全体制御部５００は、シナリオ・シミュレータ２００を再度起動し、シナリオ・シミュレータ２００は、ステップＳ５０２〜ステップＳ５０４の処理を行う。

時刻ｔ＋Δｔがシミュレーション終了時刻ｔ＿ｅｎｄより後である場合（Ｓ５０５：Ｎｏ）、指標算出部２６０は、モデル挙動記憶部２４０から、システムの挙動ｘ^（ｋ）（ｔ）の時系列を取得し、注目指標の評価値、すなわちシステムの挙動の評価関数ａ^（ｋ）を計算する（Ｓ５０６）。指標算出部２６０は、算出したａ^（ｋ）を指標値記憶部３１０に格納する。

続いて、等価微分重要度算出部２７０は、乱数重要度記憶部２５０から、全てのｊについてＨ_ｉ，ｊ ^（ｋ）を取得し、ｊについての総和を計算することで、ｉ番目のパラメータに関する等価微分重要度ｈ_ｉ ^（ｋ）を算出する（Ｓ５０７）。なお、等価微分重要度算出部２７０は、は、ユーザ等が指定したパラメータについてのみ、ステップＳ５０７における処理を行ってもよい。なお、ステップＳ５０６における処理とステップＳ５０７における処理は並行して行われてもよいし、順序を入れ替えて行われてもよい。重み関数が定める処理は、Ｓ５０３及びＳ５０７における処理を含む。

以上のように、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される場合に、図５に示した処理により、シナリオ・シミュレータ２００は、システム挙動の評価関数ａ^（ｋ）およびＮ個のパラメータそれぞれについて等価微分重要度ｈ_ｉ ^（ｋ）を計算することができる。

図６は、実施例１におけるモンテカルロ・シミュレーション・システムのユーザインタフェースの一例である。ユーザインタフェース６００は、シミュレーション条件設定セクション６１０、期待値計算セクション６２０、及び感度計算セクション６３０を含む。

シミュレーション条件設定セクション６１０は、入力受付セクション６１１、６１２、及び表示セクション６１３を含む。入力受付セクション６１１は、シミュレーションモデル、すなわちシステムの挙動の実現値ｘを定める関数ｘ（ｗ，ｔ）等の入力を受け付ける。入力受付セクション６１２は、当該シミュレーションモデルによるシミュレーション回数の入力を受け付ける。表示セクション６１３は、入力受付セクション６１１〜６１２が入力を受け付けた内容のシミュレーションを実行するために必要な計算時間を表示する。

期待値計算セクション６２０は、入力受付セクション６２１、及び表示セクション６２２を含む。入力受付セクション６２１は、注目指標、すなわちシステムの挙動の評価関数ａ（ｘ）の入力を受け付ける。表示セクション６２２は、システム評価関数Ａ（ｘ）の計算精度を表示する。

感度計算セクション６３０は、入力受付セクション６３２、６３４、チェックボックス６３１、６３３、及び表示セクション６３５を含む。入力受付セクション６３２は、感度計算の対象パラメータ、すなわち感度∂Ａ／∂ｚ_ｉを算出する対象となるパラメータｚ_ｉの入力を受け付ける。チェックボックス６３１は、全てのパラメータｚ_ｉを感度計算の対象パラメータとして選択するためのチェックボックスである。入力受付セクション６３４は、Ｄ_ｉ，ｊを計算する方法の入力を受け付ける。チェックボックス６３３は、Ｄ_ｉ，ｊを計算する方法を自動で選択するためのチェックボックスである。表示セクション６３５は、感度∂Ａ／∂ｚ_ｉの計算精度を表示する。

なお、ユーザインタフェース６００における、入力セクション、表示セクションは、それぞれ表示セクション、入力セクションであってもよいし、それぞれ表示セクション、入力セクションを兼ねてもよい。例えば、表示セクション６２２がシステム評価関数Ａ（ｘ）の計算精度の入力を受け付け、表示セクション６３５が感度∂Ａ／∂ｚ_ｉの計算精度の入力を受け付けてもよい。このとき、入力受付セクション６１１は、例えば、入力されたシステム評価関数Ａ（ｘ）の計算精度と、入力された感度∂Ａ／∂ｚ_ｉの計算精度と、を保証するために必要なシミュレーション回数を表示してもよい。

以上、本実施例のモンテカルロ・シミュレーション・システムは、１セットのモンテカルロ・シミュレーションの結果を用いて、システム評価関数の各パラメータに関する感度を算出することができ、ひいては、計算時間及び計算リソースを削減することができる。

また、本実施例のモンテカルロ・シミュレーション・システムが、算出した感度を用いて、ユーザ等は、パラメータの最適化等、すなわちパラメータ値の許容される範囲で、システム評価関数を最大化するパラメータ値の探索等を行うことができる。

本実施例は、多数のパラメータが存在する確率的なモデルにおいて、注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって同時に求める手段を、半導体集積回路の電子回路における製造ばらつきを考慮した設計に応用する例を説明する。

本実施例は、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される、すなわち、乱数ベクトルｗとパラメータｚで定まる値ｙ（ｗ,ｚ）が定まると、時刻ｔでのシステムの状態ｘ（ｔ）が、時刻ｔとｙのみで決定論的に定まる（図３に相当）一例である。

近年、半導体集積回路の微細化が進むにしたがって、集積回路を構成する個々のトランジスタの特性のばらつきが大きくなっている。集積回路を構成するトランジスタの特性ばらつきは、例えば、集積回路の製造時において、形成される回路要素の物理的なサイズ、基板に注入される不純物イオンの数、又は格子欠陥の数等の製造時に不可避的に生じる確率的な現象に起因する。つまり、当該回路は、個々のトランジスタの特性が確率的にばらつくことを想定して設計される必要がある。

したがって、トランジスタの特性のばらつきを確率分布でモデル化した上で、例えば、回路全体の特性の期待値やばらつき（分散）の大きさ、又は集積回路の歩留まり、すなわち集積回路を構成する回路全体の特性が望ましい範囲におさまる確率、等をモンテカルロ・シミュレーションによって推定することが行われる。回路設計者は、例えば、設計した回路のばらつきを考慮した特性が与えられた仕様を満たすまで、設計パラメータをチューニングしては、モンテカルロ・シミュレーションを行ってばらつきを考慮した特性を確認する、というプロセスを繰り返し行う。

個々のトランジスタにおける設計パラメータは、例えば、ゲート幅Ｗ、ゲート長Ｌ、拡散層幅ＬＯＤ、トランジスタ間距離ＰＤＸといった、ばらつき・電力・レイアウト面積等に関係するものを１０程度含む。したがって、一般に、回路ブロックの設計パラメータの数は、回路を構成するトランジスタ数の１０倍程度になる。このため、多数のトランジスタで構成される複雑な回路ブロックにおける設計パラメータは膨大な数となりえるため、上記のように個々の設計パラメータのチューニングを人手で行うことは困難である。

これに対して、本実施例における手法では、任意の注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって求めることができる。本実施例の手法を用いることで、回路設計者は、単に、現在のパラメータ設定のもとでの回路ブロック全体の特性（性能期待値や歩留まり）や電力の値のみならず、全てのトランジスタの全てのパラメータについて、それらが回路ブロック全体の特性（性能期待値や歩留まり）や電力にどのように影響するのか、その感度を同時に得ることができる。

したがって、回路設計者は、例えば、どのトランジスタのどのパラメータを変えることで、消費電力を増やさずに特性を向上させることができるか等を、やみくもにパラメータをチューニングすることなく知ることができ、設計の効率を大幅に向上できる。

図７は、６個のトランジスタで構成されるＳＲＡＭセルの回路の例を示す。ＳＲＡＭセルは、通常時にはＳＲＡＭセルに記憶された値が失われないことと、記憶の書き換え時にはＳＲＡＭセルの記憶値が意図した値にきちんと書き変わること、という互いに相反する特性を満たす必要がある。

ＳＲＡＭに含まれる各トランジスタのパラメータは、集積回路の製造時のばらつきによって、製造された集積回路における個々のトランジスタの実際のパラメータ値は設計時に指定した値にはならず、例えば、ある確率分布にしたがって分布する。その結果、例えば、通常にＳＲＡＭセルに記憶していた値が失われる、又は、記憶の書き換え時にはＳＲＡＭセルの記憶値が意図した値にきちんと書き変わらない、等の不良が起こりえる。

そこで、回路設計者は、設計パラメータ値をいろいろな値に変えてモンテカルロ・シミュレーションを行い、実際に製造されるＳＲＡＭセル回路の不良率が与えられた閾値以下になるように設計パラメータを最適化する必要がある。この例に、実施例１のモンテカルロ・シミュレーション・システムのシミュレーション手法を適用する。

本実施例におけるパラメータベクトルｚは、トランジスタの設計パラメータからなる。例えば、個々のトランジスタに１０の設計パラメータがある場合、６個のトランジスタで構成されるＳＲＡＭセル回路におけるパラメータベクトルｚは、６０次元のベクトルである。

ところで、前述のように、ＳＲＡＭ回路の特性ばらつきは、トランジスタの実際のパラメータ値が製造時にばらつくことに起因する。すなわち、トランジスタの特性を決める６０次元の設計パラメータｚに対して、実際には製造時のばらつきにより、例えば、６０次元のパラメータｙが生成される。

乱数発生部２１０は、シミュレーション毎に確率密度関数ｇ（ｗ）に従う多次元の乱数ベクトルｗを生成する。モデル挙動更新部２２０は、乱数ベクトルｗと６０次元の設計パラメータベクトルｚと、から製造後の実際のパラメータ値ｙを関数ｙ＝ｙ（ｗ，ｚ）によって定める。このとき、製造後の実際のパラメータ値ｙの確率分布が、設計パラメータベクトルがｚであるときの、製造ばらつきを表すように、乱数ベクトルｗが従う確率密度関数ｇ（ｗ）と、関数ｙ（ｗ，ｚ）が定められているものとする。

モデル挙動更新部２２０は、製造後の実際のパラメータｙを定めれば、ＳＲＡＭセル回路の挙動ｘ（ｔ）を、通常の回路シミュレーションによって決定論的に定めることができる。したがって、本実施例におけるシステムの挙動を表す時系列ｘ（ｔ）は、数１５の形で表される、すなわち、乱数ベクトルｗとパラメータベクトルｚで定まる値ｙ（ｗ,ｚ）を定めると、時刻ｔでのシステムの状態ｘ（ｔ）が，時刻ｔとｙのみで決定論的に定まる。また、システムの挙動の評価関数ａ（ｘ）を、そのセルの挙動が正常動作の範囲内であれば０、不良であれば１と定義すれば、システム評価関数Ａ（ｘ）、すなわち、ａ（ｘ）の期待値Ｅ［ａ（ｘ）］は、ＳＲＡＭセルの不良率を示す。

したがって、モンテカルロ・シミュレーション・システムは、図４で説明した処理を行うことにより、ＳＲＡＭセルの不良率の６０個の設計パラメータ全てに対する感度を算出することができる。また、モンテカルロ・シミュレーション・システムは、は、算出した感度から、勾配法やニュートン法といった公知の最適化手法を用いてＳＲＡＭセル不良率を最小にするような６０個の設計パラメータ値の組み合わせを、算出してもよい。

以上のように、本実施例のモンテカルロ・シミュレーション・システムは、多数のトランジスタで構成される半導体集積回路の電子回路における、製造ばらつきを考慮した設計パラメータ値の最適値を効率的に求めることができる。

本実施例は、多数のパラメータが存在する確率的なモデルにおいて、注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって同時に求める手段を、金融分野、とくに、銀行の収益予測シミュレーションシステムに適用した例を説明する。本実施例における、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される、すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）が、時刻ｔと、過去（時刻ｔ以前）のシステムの挙動と、時刻ｔで新たに生成する（時刻毎に異なる）乱数ｗ_ｊと、パラメータｚと、で定まる。

金融機関においては、将来の事業計画の立案時に、取り扱っている金融商品の将来の取扱高やその収益について計画を立てて、当該事業計画の妥当性を経営判断する。

図８は、銀行における収支計画表の例を示す。収支計画表の最左列は、当該銀行における取扱金融商品の科目を示す。取扱金融商品は、例えば、企業向けの事業融資、個人向けの住宅ローン貸付、カードローン貸付、および銀行が自ら保有している有価証券配当について、それぞれＡからＣまでの３区分の科目を含む。

実際の銀行においては、数千程度の金融商品の科目を保有することがある。表の各行は、最上行に示される期間における、金融商品の科目それぞれについての取扱高を示す。例えば、現時点が２０１４年２Ｑ終了時とすれば、‘１４／１Ｑおよび‘１４／２Ｑの取扱高は実績値を示すが、‘１４／３Ｑ以降の取扱高は現時点（２０１４年２Ｑ終了時点）における計画値を示す。図８の例では四半期単位で計６期間についての計画を示しているが、実際の銀行においては、例えば、月単位又は週単位で５年程度の計画を立てることもあり、このとき、収支計画表は数百程度の期間を含む。

各科目における取扱高と収益との関係は、例えば、主に銀行の資金調達金利、すなわちスポットレートに依存する。スポットレートと収益の関係については、科目毎にある程度の確度をもったモデル化が可能であるのに対して、将来のスポットレートは不確定であり確率的な予測しかできない。

そこで、銀行の経営においては、モンテカルロ・シミュレーションによって収益の期待値やばらつきをシミュレーションし、事業計画の妥当性を判断する。具体的には、例えば、乱数を用いて将来のスポットレートの実現パスを生成し、各科目の計画取扱高をもとに収益値を計算し、全科目の収益を合計することで銀行全体の収益値を計算する、という処理を多数回行う。当該多数回の処理によって、与えられた事業計画の下での銀行全体の収益の期待値やばらつきを推定する。

ところで銀行において、事業計画の妥当性を判断するにあたっては、将来の各科目の取扱高の計画値が完全に達成された場合の銀行全体の期待収益やばらつきを求めるだけでは不足である。各科目の将来の実際の取扱高は、計画値からずれることが当然に想定されるためである。このとき、計画値と将来の実績がずれたときの影響について、事業計画立案時に予め考慮しておく必要がある。

具体的には、例えば、数千程度の科目それぞれについての数百程度の未来の各期間における取扱高の計画値、すなわち、合計十万個程度の計画値それぞれについて、将来の実績値が計画値からずれたときの銀行全体の収益に対する影響を計算する必要がある。もし実績値が計画値からほんの少しだけずれるだけで、銀行全体の収益に大きな影響を及ぼすような計画となっていることが分かれば、適切な引当金を積んでおく、又は、事業計画そのものを立て直す、といった対策が必要になるであろう。

以下、実施例１のモンテカルロ・シミュレーション・システムによるシミュレーション手法を銀行の収益予測シミュレーションに適用する例を説明する。

本実施例におけるパラメータベクトルｚは、例えば、全ての科目それぞれについて未来の各期間における取扱高の計画値である。例えば、科目数は数千程度、期間は数百期間である場合、パラメータベクトルｚは十万次元程度のベクトルである。システムの挙動の時系列ｘ（ｔ）は、例えば、将来の各期間における各科目の収益値の集合（例えば、数千次元のベクトルになりえる）と、金利のスポットレートと、を合わせた多次元のベクトルである。

ｔ番目の期間における各科目の収益値は、例えば、スポットレートと、その時刻における各科目の取扱高の計画値であるパラメータｚと、に依存する。将来のスポットレートは、例えば、Ｈｕｌｌ−Ｗｈｉｔｅモデル、ＣＩＲモデル、又はＢＤＴモデルなどといった確率的なモデルで表現される。これらのモデルは、時刻ｔ＋１でのスポットレートを、時刻ｔでのスポットレートと、（時刻ｔで新たに生成する）乱数と、に基づいて定める。

したがって、本実施例におけるシステムの挙動ｘ（ｔ）は、全体として、数１７の形で表される。すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）は、時刻ｔと、過去（時刻ｔ以前）のシステムの挙動と、時刻ｔで新たに生成される（時刻毎に異なる）乱数ｗ_ｊと、パラメータｚと、からで定まる。

システムの挙動の評価関数ａ（ｘ）は、例えば、全科目の全期間の収益の合計値である。このとき、システム評価関数Ａ（ｘ）、すなわちａ（ｘ）の期待値Ｅ［ａ（ｘ）］は、与えられた計画のもとでの銀行全体の将来の収益の期待値を示す。

本実施例のモンテカルロ・シミュレーション・システムは、上述した前提条件のもとで、図５に示した処理を行うことで、例えば、数千程度ある全ての科目それぞれについて数百程度ある未来の各期間における取扱高の計画値、すなわち、合計十万個程度の計画値それぞれについて、銀行全体の収益の期待値に対する感度、を同時に推定することが可能になる。

銀行は、推定した感度に基づいて、例えば、計画値と実績値のずれの影響を予め考慮して適切な引当金を積んでおく、又は事業計画そのものを立て直す、といった対策を立てることが可能となる。

以上のように、本実施例のモンテカルロ・シミュレーション・システムは、金融分野、とくに銀行の収益予測シミュレーションシステムにおいて計画値の収益に対する感度を推定することが可能となるため、銀行は事業計画立案のプロセスを効率化することができる。

本実施例のモンテカルロ・シミュレーション・システムは、実施例１の数９における境界残余Ｒ_ｉを算出する。前述のように、境界残余Ｒ_ｉは、乱数ベクトルｗの分布範囲Ωの境界面（表面）∂Ωで（∂ｗ／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}＝０（ベクトル）であれば、Ｒ_ｉ＝０となる。現実的にモンテカルロ・シミュレーションの考察対象となるシステムでは、この条件が成立して、Ｒ_ｉ＝０となる場合が多いと考えられる。

しかしながら、例えば、乱数ベクトルｗが一様分布乱数である場合などに、この条件が成立せず、具体的に境界残余Ｒ_ｉを算出する必要にせまられる場合がある。境界残余Ｒ_ｉの算出には、数９の１行目に基づく方法と、数９の３行目に基づく方法と、の２通りの方法が考えられる。

図９は、乱数ベクトルｗが、区間［ａ，ｂ］上の一様分布乱数ベクトルｗ_１と、区間［ｃ，ｄ］上の一様分布乱数ｗ_２と、の２つの要素をもつ２次元のベクトルである場合の境界残余Ｒ_ｉを算出する概念の例を示す。図中に矩形で示された範囲がΩである。確率密度関数ｇ（ｗ）は、以下の数１９で表される。

なお、本実施例における記憶部３００は、サンプルパス固定微分記憶部３３０（不図示）をさらに含む。サンプルパス固定微分記憶部３３０は、サンプルパス固定微分計算部２３１が算出したＤ_ｉ，ｊを保持する。また、本実施例における後処理計算部４００は、プログラムである境界残余算出部４３０（不図示）をさらに含む。境界残余算出部４３０は、乱数ベクトルｗと、パラメータベクトルｚと、サンプルパス固定微分記憶部３３０から受信するＤ_ｉ，ｊと、を入力として境界残余Ｒ_ｉを算出し、算出した境界残余Ｒ_ｉを期待値計算部４１０に送信する。

以下、乱数ベクトルｗが従う確率密度関数ｇが数１９で表される例を用いて、境界残余Ｒ_ｉの２通りの算出方法を説明する。境界残余算出部４３０が、境界残余Ｒ_ｉを算出する第１の方法は、数９の１行目に基づくものである。図９の条件において、数９の１行目は、以下の数２０のように展開される。

したがって、境界残余算出部４３０は、ｗ_１をΩの境界上の値であるｂに固定したうえでｗ_２を動かしてそれぞれ（∂ｗ_２／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}を計算した積分値から、ｗ_１を境界値ａに固定したうえでｗ_２を動かしてそれぞれ（∂ｗ_２／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}を計算した積分値を減算した値を算出する。また、境界残余算出部４３０は、ｗ_２を分布の境界上の値であるｄに固定したうえでｗ_１を動かしてそれぞれ（∂ｗ_１／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}を計算した積分値から、ｗ_２を境界値ｃに固定したうえでｗ_１を動かしてそれぞれ（∂ｗ_１／∂ｚ_ｉ）｜_{ｘ＝ｃｏｎｓｔ}を計算した積分値を減算したものと、を加算した値を算出する。

境界残余算出部４３０は、算出した２つの値を加算することで、境界残余Ｒ_ｉを算出する。なお、境界残余算出部４３０は、例えば、既知の数値積分の手法を用いて、上述した積分値それぞれを算出すればよい。

境界残余算出部４３０が、境界残余Ｒ_ｉを算出する第２の方法は、数９の３行目に基づくものである。数９の３行目は、乱数ベクトルの分布範囲Ω全域における積分であるから、境界残余算出部４３０は、以下の数２１より、モンテカルロ・シミュレーションを用いた期待値として境界残余Ｒ_ｉを算出することができる。

すなわち、乱数発生部２１０が、与えられた確率密度関数ｇ（ｗ）に従う乱数を多数発生させる。境界残余算出部４３０は、発生した乱数ベクトルそれぞれに対して、数２１の２行目の角括弧内の値を計算し、計算した値の平均値をＲ_ｉの推定値とすればよい。

なお、境界残余算出部４３０は、数値的な微分により、数２１中のｄｉｖを算出すればよい。また、ユーザ等によって関数としてｄｉｖが与えられた場合、境界残余算出部４３０は、当該関数にパラメータベクトルｚと乱数ベクトルｗとを代入することにより、ｄｉｖを算出してもよい。なお、境界残余算出部４３０は、確率密度関数ｇ（ｗ）が、２次元一様分布以外である場合も同様に境界残余Ｒ_ｉを算出することができる。感度計算部４２０は、上述の方法で算出された境界残余Ｒ_ｉを、算出した平均値から引いた差を算出する。当該算出した差が当該感度の推定値である。

以上のように、本実施例のモンテカルロ・シミュレーション・システムは、境界残余Ｒ_ｉがゼロでな場合に、その値を算出することが可能となる。また、モンテカルロ・シミュレーション・システムは、境界残余Ｒ_ｉを算出することにより、システムの挙動の実現値ｘを算出するモデルに適用する乱数ベクトルがどのような分布であっても、即ちｘを算出する全てのモデルに対して、本実施例の感度計算手法を適用することができる。

なお、本発明は上記した実施例に限定されるものではなく、様々な変形例が含まれる。例えば、上記した実施例は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、ある実施例の構成の一部を他の実施例の構成に置き換えることも可能であり、また、ある実施例の構成に他の実施例の構成を加えることも可能である。また、各実施例の構成の一部について、他の構成の追加・削除・置換をすることが可能である。

また、上記の各構成、機能、処理部、処理手段等は、それらの一部又は全部を、例えば集積回路で設計する等によりハードウェアで実現してもよい。また、上記の各構成、機能等は、プロセッサがそれぞれの機能を実現するプログラムを解釈し、実行することによりソフトウェアで実現してもよい。各機能を実現するプログラム、テーブル、ファイル等の情報は、メモリや、ハードディスク、ＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｒｉｖｅ）等の記録装置、または、ＩＣカード、ＳＤカード、ＤＶＤ等の記録媒体に置くことができる。

また、制御線や情報線は説明上必要と考えられるものを示しており、製品上必ずしも全ての制御線や情報線を示しているとは限らない。実際には殆ど全ての構成が相互に接続されていると考えてもよい。

Claims (15)

1. 複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムであって、
   プロセッサと、
   メモリと、を含み、
   前記メモリは、
   １以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、
   前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、
   乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、
   前記プロセッサは、
   複数の乱数ベクトルと、パラメータベクトルと、を取得し、
   前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、
   前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、
   前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、
   前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、シミュレーションシステム。
2. 請求項１に記載のシミュレーションシステムであって、
   前記重み関数は、
   前記第１の微分係数それぞれの、前記第１の微分係数それぞれに対応する乱数による、第２の微分係数を決定し、
   前記第２の微分係数に基づき前記パラメータに対する前記乱数ベクトルの重みを決定する、シミュレーションシステム。
3. 請求項２に記載のシミュレーションシステムであって、
   前記重み関数における重みは、下記数式で表わされ、
   

   

   上記数式における、ｗは前記乱数ベクトル、ｚは前記パラメータベクトル、Ｍは前記乱数ベクトルの次元、ｗ_ｊは前記乱数ベクトルｗのｊ番目の要素_、ｚ_ｉは前記パラメータベクトルのｉ番目の要素、ｘは前記乱数ベクトルｗに対応する前記確率的システムの挙動値、ｇは前記乱数ベクトルｗが従う確率密度関数である、シミュレーションシステム。
4. 請求項２に記載のシミュレーションシステムであって、
   前記取得した複数の乱数ベクトルは、一様分布に従い、
   前記重み関数において、前記重みは、前記乱数ベクトルの全ての乱数の前記第２の微分係数の和である、シミュレーションシステム。
5. 請求項２に記載のシミュレーションシステムであって、
   前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度は、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、の積の平均値で表わされる、シミュレーションシステム。
6. 請求項１に記載のシミュレーションシステムであって、
   前記挙動関数の独立変数は、時刻及び時刻非依存の内部関数であり、
   前記内部関数は、確率分布に従う乱数ベクトルと、パラメータベクトルと、に基づいて、第１ベクトルを決定し、
   前記挙動関数は、一つの乱数ベクトルから、前記確率的システムの挙動の実現値の一つの時系列を決定し、
   前記挙動評価関数は、前記一つの時系列に基づいて、前記確率的システムの挙動の一つの評価値を決定し、
   前記重み関数は、前記内部関数による第１ベクトルを一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定し、
   前記プロセッサは、
   前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値の時系列を算出し、
   前記算出した実現値の時系列それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出する、シミュレーションシステム。
7. 請求項１に記載のシミュレーションシステムであって、
   前記挙動関数は、前記確率的システムの挙動の実現値の時系列を決定し、
   前記挙動関数は、前記時系列における各時刻の実現値を、各時刻に対応付けられている乱数と各時刻より前の時刻の実現値とに基づき決定し、
   前記挙動評価関数は、乱数ベクトルに対応する前記確率的システムの挙動の実現値の時系列に基づいて、前記乱数ベクトルに対応する前記確率的システムの挙動の評価値を決定し、
   前記重み関数は、前記乱数ベクトルの乱数それぞれに対応する時刻における挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定し、
   前記プロセッサは、
   前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値の時系列を算出し、
   前記算出した実現値の時系列それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出する、シミュレーションシステム。
8. 請求項５に記載のシミュレーションシステムであって、
   前記プロセッサは、
   前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値、前記取得した複数の乱数ベクトルが発生する確率、及び前記取得した複数の乱数ベクトルの乱数それぞれの第１の微分係数からなるベクトルと前記取得した複数の乱数ベクトルが従う確率密度関数の分布範囲の境界の外向き単位法線ベクトルとの内積、の積の前記境界上における積分値を算出し、
   前記平均値と前記積分値とに基づいて、前記確率的システムの挙動の評価値の期待値の感度を算出する、シミュレーションシステム。
9. 請求項５に記載のシミュレーションシステムであって、
   前記プロセッサは、
   前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値、前記取得した複数の乱数ベクトルが発生する確率、及び前記取得した複数の乱数ベクトルの乱数それぞれ第１の微分係数からなるベクトル、の積の発散と、前記取得した複数の乱数ベクトルが発生する確率の逆数と、の積の第１の期待値を算出し、
   前記平均値と前記第１の期待値とに基づいて、前記確率的システムの挙動の評価値の期待値の感度を算出する、シミュレーションシステム。
10. 請求項１に記載のシミュレーションシステムであって、
    前記取得した複数の乱数ベクトルが従う確率密度関数は、前記取得したパラメータベクトルのパラメータそれぞれが微小変化する条件において、前記乱数ベクトルに対応する挙動関数の実現値を一定に保つ、乱数ベクトルの微小変化に対して、定義域が一定である、シミュレーションシステム。
11. 複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムにおけるシミュレーション方法であって、
    前記シミュレーションシステムは、プロセッサと、メモリと、を含み、
    前記メモリは、
    １以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、
    前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、
    乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、
    前記方法は、
    前記プロセッサが、
    複数の乱数ベクトルと、パラメータベクトルと、を取得し、
    前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、
    前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、
    前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、
    前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、ことを含む方法。
12. 請求項１１に記載のシミュレーション方法であって、
    前記重み関数は、
    前記第１の微分係数それぞれの、前記第１の微分係数それぞれに対応する乱数による、第２の微分係数を決定し、
    前記第２の微分係数に基づき前記パラメータに対する前記乱数ベクトルの重みを決定する、方法。
13. 請求項１１に記載のシミュレーション方法であって、
    前記重み関数における重みは、下記数式で表わされ、
    

    

    上記数式における、ｗは前記乱数ベクトル、ｚは前記パラメータベクトル、Ｍは前記乱数ベクトルの次元、ｗ_ｊは前記乱数ベクトルｗのｊ番目の要素_、ｚ_ｉは前記パラメータベクトルのｉ番目の要素、ｘは前記乱数ベクトルｗに対応する前記確率的システムの挙動値、ｇは前記乱数ベクトルｗが従う確率密度関数である、方法。
14. 請求項１２に記載のシミュレーション方法であって、
    前記取得した複数の乱数ベクトルは、一様分布に従い、
    前記重み関数において、前記重みは、前記乱数ベクトルの全ての乱数の前記第２の微分係数の和である、方法。
15. 請求項１２に記載のシミュレーション方法であって、
    前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度は、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、の積の平均値で表わされる、方法。

Images