(57)【要約】  シミュレーションシステムは、複数の乱数ベクトルと、パラメータベクトルと、を取得し、取得した複数の乱数ベクトルそれぞれに対応する確率的システムの挙動の実現値を算出し、取得した複数の乱数ベクトルそれぞれと、取得したパラメータベクトルと、重み関数と、に基づき、取得したパラメータベクトルのパラメータそれぞれに対する、取得した乱数ベクトルそれぞれの重みを算出し、取得した複数の乱数ベクトルそれぞれに対応する確率的システムの挙動の評価値を算出し、取得した複数の乱数ベクトルそれぞれに対応する算出した挙動の評価値と、取得したパラメータベクトルから選択したパラメータに対する取得した乱数ベクトルそれぞれの重みと、に基づいて、選択したパラメータに対する、確率的システムの挙動の評価値の期待値の感度、を算出する。
2016030938.tif 000002
【特許請求の範囲】
【請求項１】
複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムであって、
プロセッサと、
メモリと、を含み、
前記メモリは、
１以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、
前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、
乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、
前記プロセッサは、
複数の乱数ベクトルと、パラメータベクトルと、を取得し、
前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、
前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、
前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、
前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、シミュレーションシステム。
【請求項２】
請求項１に記載のシミュレーションシステムであって、
前記重み関数は、
前記第１の微分係数それぞれの、前記第１の微分係数それぞれに対応する乱数による、第２の微分係数を決定し、
前記第２の微分係数に基づき前記パラメータに対する前記乱数ベクトルの重みを決定する、シミュレーションシステム。
【請求項３】
請求項２に記載のシミュレーションシステムであって、
前記重み関数における重みは、下記数式で表わされ、
【数１】
2016030938.tif 000024
上記数式における、ｗは前記乱数ベクトル、ｚは前記パラメータベクトル、Ｍは前記乱数ベクトルの次元、ｗｊは前記乱数ベクトルｗのｊ番目の要素、ｚｉは前記パラメータベクトルのｉ番目の要素、ｘは前記乱数ベクトルｗに対応する前記確率的システムの挙動値、ｇは前記乱数ベクトルｗが従う確率密度関数である、シミュレーションシステム。
【請求項４】
請求項２に記載のシミュレーションシステムであって、
前記取得した複数の乱数ベクトルは、一様分布に従い、
前記重み関数において、前記重みは、前記乱数ベクトルの全ての乱数の前記第２の微分係数の和である、シミュレーションシステム。
【請求項５】
請求項２に記載のシミュレーションシステムであって、
前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度は、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、の積の平均値で表わされる、シミュレーションシステム。
【請求項６】
請求項１に記載のシミュレーションシステムであって、
前記挙動関数の独立変数は、時刻及び時刻非依存の内部関数であり、
前記内部関数は、確率分布に従う乱数ベクトルと、パラメータベクトルと、に基づいて、第１ベクトルを決定し、
前記挙動関数は、一つの乱数ベクトルから、前記確率的システムの挙動の実現値の一つの時系列を決定し、
前記挙動評価関数は、前記一つの時系列に基づいて、前記確率的システムの挙動の一つの評価値を決定し、
前記重み関数は、前記内部関数による第１ベクトルを一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定し、
前記プロセッサは、
前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値の時系列を算出し、
前記算出した実現値の時系列それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出する、シミュレーションシステム。
【請求項７】
請求項１に記載のシミュレーションシステムであって、
前記挙動関数は、前記確率的システムの挙動の実現値の時系列を決定し、
前記挙動関数は、前記時系列における各時刻の実現値を、各時刻に対応付けられている乱数と各時刻より前の時刻の実現値とに基づき決定し、
前記挙動評価関数は、乱数ベクトルに対応する前記確率的システムの挙動の実現値の時系列に基づいて、前記乱数ベクトルに対応する前記確率的システムの挙動の評価値を決定し、
前記重み関数は、前記乱数ベクトルの乱数それぞれに対応する時刻における挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定し、
前記プロセッサは、
前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値の時系列を算出し、
前記算出した実現値の時系列それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出する、シミュレーションシステム。
【請求項８】
請求項５に記載のシミュレーションシステムであって、
前記プロセッサは、
前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値、前記取得した複数の乱数ベクトルが発生する確率、及び前記取得した複数の乱数ベクトルの乱数それぞれの第１の微分係数からなるベクトルと前記取得した複数の乱数ベクトルが従う確率密度関数の分布範囲の境界の外向き単位法線ベクトルとの内積、の積の前記境界上における積分値を算出し、
前記平均値と前記積分値とに基づいて、前記確率的システムの挙動の評価値の期待値の感度を算出する、シミュレーションシステム。
【請求項９】
請求項５に記載のシミュレーションシステムであって、
前記プロセッサは、
前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値、前記取得した複数の乱数ベクトルが発生する確率、及び前記取得した複数の乱数ベクトルの乱数それぞれ第１の微分係数からなるベクトル、の積の発散と、前記取得した複数の乱数ベクトルが発生する確率の逆数と、の積の第１の期待値を算出し、
前記平均値と前記第１の期待値とに基づいて、前記確率的システムの挙動の評価値の期待値の感度を算出する、シミュレーションシステム。
【請求項１０】
請求項１に記載のシミュレーションシステムであって、
前記取得した複数の乱数ベクトルが従う確率密度関数は、前記取得したパラメータベクトルのパラメータそれぞれが微小変化する条件において、前記乱数ベクトルに対応する挙動関数の実現値を一定に保つ、乱数ベクトルの微小変化に対して、定義域が一定である、シミュレーションシステム。
【請求項１１】
複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムにおけるシミュレーション方法であって、
前記シミュレーションシステムは、プロセッサと、メモリと、を含み、
前記メモリは、
１以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、
前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、
乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、
前記方法は、
前記プロセッサが、
複数の乱数ベクトルと、パラメータベクトルと、を取得し、
前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、
前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、
前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、
前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、ことを含む方法。
【請求項１２】
請求項１１に記載のシミュレーション方法であって、
前記重み関数は、
前記第１の微分係数それぞれの、前記第１の微分係数それぞれに対応する乱数による、第２の微分係数を決定し、
前記第２の微分係数に基づき前記パラメータに対する前記乱数ベクトルの重みを決定する、方法。
【請求項１３】
請求項１１に記載のシミュレーション方法であって、
前記重み関数における重みは、下記数式で表わされ、
【数２】
2016030938.tif 000025
上記数式における、ｗは前記乱数ベクトル、ｚは前記パラメータベクトル、Ｍは前記乱数ベクトルの次元、ｗｊは前記乱数ベクトルｗのｊ番目の要素、ｚｉは前記パラメータベクトルのｉ番目の要素、ｘは前記乱数ベクトルｗに対応する前記確率的システムの挙動値、ｇは前記乱数ベクトルｗが従う確率密度関数である、方法。
【請求項１４】
請求項１２に記載のシミュレーション方法であって、
前記取得した複数の乱数ベクトルは、一様分布に従い、
前記重み関数において、前記重みは、前記乱数ベクトルの全ての乱数の前記第２の微分係数の和である、方法。
【請求項１５】
請求項１２に記載のシミュレーション方法であって、
前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度は、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、の積の平均値で表わされる、方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーションシステム、及びシミュレーション方法に関する。
【背景技術】
【０００２】
興味の対象であるシステム（系）が与えられたとき、当該システムの振る舞いをモデル化し、計算機が当該モデルに基づいて当該システムの振る舞いをシミュレーションし、将来の振る舞いを予測することが行われる。例えば、ある種の物理現象や、多数の要素が絡み合う社会現象などは、単純な決定論的なモデルでは記述されることができず、複雑かつ確率的なモデルとして記述される。【０００３】
決定論的なシステムは、与えられた初期条件のもとで、システムの将来の振る舞いが一つに定まる。これに対して、確率的なシステムは、システムの将来の振る舞いはある確率分布にしたがって分布しており、計算機は、例えば、注目する指標の期待値等の統計量を求める場合、シミュレーションを行う必要がある。【０００４】
確率的なシステムの統計的な振る舞いを計算機が予測、シミュレーションする方法として、モンテカルロ・シミュレーション法が知られている。モンテカルロ・シミュレーション法は、対象システムに存在する確率的な振る舞いを乱数によって模擬したシミュレーションを多数回繰り返し行うことで、注目する指標の期待値を推定する方法である。このとき、シミュレーションの回数を増やすことで、大数の法則により、注目する指標の期待値が真の値に近づくと考えられる。【０００５】
本技術分野の背景技術として、特開２００４−１１７２２８号公報（特許文献１）がある。特許文献１には、「気象物理量を対象とした気象庁発表の長期予報を取得するデータ取得ステップ０２０１、気象物理量の従う時系列モデルを作成する気象時系列モデル作成ステップ０２０２、推定期間における推定地点の気象物理量をシミュレーション回数だけシミュレートする気象物理量シミュレーションステップ０２０３、気象物理量のシミュレーション結果を気象物理量推定結果として出力する気象物理量推定結果出力ステップ０２０４を有する。さらに、気象時系列モデル作成ステップ０２０２を、気象時系列モデルのパラメータ推定ステップ０２０２１および気象時系列モデルのパラメータ補正ステップ０２０２２から構成する。」と記載されている（要約参照）。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００４−１１７２２８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
例えば、特許文献１に記載の技術のシミュレーション対象である気象物理量のように、モンテカルロ・シミュレーションの対象である、確率的なモデルには、一般的に多数のパラメータが存在することが多い。このとき、モンテカルロ・シミュレーション実行にあたっては、あるパラメータ値の下での注目する指標の期待値だけではなく、注目指標の期待値のパラメータ値に対する感度を算出したい場合がある。注目指標の期待値のパラメータ値に対する感度とは、与えられたパラメータ値が微小に変化したときに、注目指標の期待値がどれくらい変化するかを表す微分係数である。【０００８】
例えば、システムのパラメータの最適化を行う、すなわち、許容されるパラメータ値の範囲の中で、ある注目指標の期待値を最大化するパラメータ値を探す場合、与えられたパラメータ値の下での感度を得られれば、勾配法やニュートン法といった公知の最適化手法によって最適化を行うことができる。【０００９】
ところが、モンテカルロ・シミュレーションにおいて、注目指標の期待値のパラメータ値に対する感度を求めることは、とくにパラメータが多数ある場合には、容易ではない。注目指標の期待値のパラメータに対する感度は、単純には、与えられたパラメータ値の下での注目指標の期待値と、パラメータをそこから微小に変化させたときの注目指標の期待値との差を、パラメータの微小変化量で割ることで算出できる。【００１０】
特許文献１に記載の技術を、多数のパラメータについて注目指標の期待値の感度の推定に適用する場合、例えば特許文献１の数５、数１２のように、多数のパラメータそれぞれを一つずつ微小変化させた注目指標の期待値それぞれを算出する必要がある。すなわち、特許文献１に記載の技術を適用した感度推定手法は、当該多数パラメータの個数と同数セットのモンテカルロ・シミュレーションを行う必要がある。【００１１】
そもそも、１セットのモンテカルロ・シミュレーションを行って注目指標の期待値を求めること自体が、大数の法則を満足するような多数回のシミュレーションをする必要がある。したがって、当該多数パラメータの個数と同数セットのモンテカルロ・シミュレーションを行うことは、計算機リソースや計算時間の制限により困難である。
【課題を解決するための手段】
【００１２】
上記課題を解決するために、本発明は、例えば、以下のような構成を採用する。複数の乱数ベクトルそれぞれを入力とする確率的システムのシミュレーションを繰り返して、前記確率的システムの挙動の複数の実現値を算出する、シミュレーションシステムであって、プロセッサと、メモリと、を含み、前記メモリは、１以上の乱数からなる乱数ベクトルと、１以上のパラメータからなるパラメータベクトルと、に基づいて、確率的システムの挙動の実現値を決定する、挙動関数と、前記確率的システムの挙動の実現値に基づいて、前記確率的システムの挙動の評価値を決定する、挙動評価関数と、乱数ベクトルに対応する挙動関数の実現値を一定に保つ条件において、前記乱数ベクトルの乱数それぞれのパラメータによる第１の微分係数を決定し、当該第１の微分係数に基づき当該パラメータに対する当該乱数ベクトルの重みを決定する、重み関数と、を保持し、前記プロセッサは、複数の乱数ベクトルと、パラメータベクトルと、を取得し、前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記挙動関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の実現値を算出し、前記取得した複数の乱数ベクトルそれぞれと、前記取得したパラメータベクトルと、前記重み関数と、に基づき、前記取得したパラメータベクトルのパラメータそれぞれに対する、前記取得した乱数ベクトルそれぞれの重みを算出し、前記算出した実現値それぞれと、前記挙動評価関数と、に基づいて、前記取得した複数の乱数ベクトルそれぞれに対応する前記確率的システムの挙動の評価値を算出し、前記取得した複数の乱数ベクトルそれぞれに対応する前記算出した挙動の評価値と、前記取得したパラメータベクトルから選択したパラメータに対する前記取得した乱数ベクトルそれぞれの重みと、に基づいて、前記選択したパラメータに対する、前記確率的システムの挙動の評価値の期待値の感度、を算出する、シミュレーションシステム。
【発明の効果】
【００１３】
本発明の一態様によれば、パラメータが存在する確率的なシステムにおいて、任意に決めた注目指標の期待値のパラメータに対する感度を、１セットのモンテカルロ・シミュレーションによって得られた結果から求めることができる。これにより、注目指標の期待値のパラメータに対する感度を求めるために必要な、計算機リソースや計算時間を削減することができる。
【図面の簡単な説明】
【００１４】
【図１】実施例１におけるモンテカルロ・シミュレーション・システムの構成例を表すブロック図である。
【図２】実施例１における乱数重要度算出部の構成例を表すブロック図である。
【図３】実施例１のモンテカルロ・シミュレーションの全体処理の例を表すフローチャートである。
【図４】実施例１におけるシミュレーション実行処理の一例を表すフローチャートである。
【図５】実施例１におけるシミュレーション実行処理の別の一例を表すフローチャートである。
【図６】実施例１におけるモンテカルロ・シミュレーション・システムのユーザインタフェースの一例である。
【図７】実施例２における６個のトランジスタで構成されるＳＲＡＭセルの回路の例である。
【図８】実施例３における、銀行の収支計画表の例である。
【図９】実施例４において、乱数ベクトルが２次元一様分布に従う場合における、境界残余を算出する概念の例を示す図である。
【発明を実施するための形態】
【００１５】
以下、本発明の実施形態について述べる。本実施形態は本発明を実現するための一例に過ぎず、本発明の技術的範囲を限定するものではないことに注意すべきである。本実施形態においては、便宜上その必要があるときは、複数のセクションまたは実施形態に分割して説明するが、特に明示した場合を除き、それらは互いに無関係なものではなく、一方は他方の一部または全部の変形例、詳細、補足説明等の関係にある。また、以下の実施形態において、要素の数等（個数、数値、量、範囲等を含む）に言及する場合、特に明示した場合および原理的に明らかに特定の数に限定される場合等を除き、その特定の数に限定されるものではなく、特定の数以上でも以下でもよい。
【００１６】
以下、本発明の実施形態を、添付図面を参照して説明する。なお、各図において、共通の構成には原則として同一の参照符号を付し、繰り返しの説明は省略する。
【実施例１】
【００１７】
本実施例はモンテカルロ・シミュレーション・システムについて説明する。まず、本実施例のモンテカルロ・シミュレーション・システムがシミュレーションを実行し、感度を算出するための理論的側面について説明する。【００１８】
確率変数Ｘは、本実施例のモンテカルロ・シミュレーション・システムによる、シミュレーション対象である確率的なシステムの挙動を表す。つまり、当該シミュレーション対象である確率的なシステムの挙動は確率変数Ｘでモデル化される。以下、シミュレーション対象である確率的なシステムを、単にシステムとも表記する。また、ｘはシステムの挙動Ｘの実現値を表す。【００１９】
ここで、システムの挙動Ｘは、１次元又は多次元のベクトルのいずれであってもよく、また、「時刻」と呼ばれる変数ｔで区別される無限個の確率変数の組（離散時間確率過程と呼ばれる）であってもよい。したがって、システムの挙動Ｘが離散時間確率過程である場合、Ｘは本来Ｘ（ｔ）と表記されるべきであるが、特に必要がある場合を除いて、単にＸとも表記する。また、ある時刻ｔ１におけるシステムの挙動を示す値ｘ（ｔ１）を、時刻ｔ１におけるシステムの状態と呼ぶ。【００２０】
システムの挙動Ｘの分布は、システムのパラメータｚに依存し、さらにＸの確率密度関数が、ｆ（ｘ，ｚ）という形で表記されるとする。以下、システムのパラメータｚを単にパラメータｚと表記する。パラメータｚは、Ｎ次元のベクトル（すなわち、システムの挙動はＮ個のパラメータに依存している）であり、ｚの成分は、ｚ＝（ｚｉ）（ただし、ｉ＝１，…，Ｎ）であるとする。また、１つのシステムの挙動の実現値Ｘ＝ｘに対して１つの実数値を定める、関数ａ（ｘ）が与えられているとする。ａ（ｘ）をシステムの挙動の評価関数と呼ぶ。システムの挙動の評価関数は、注目指標の一例である。さらに、ａ（ｘ）の期待値Ａ（ｘ）＝Ｅ［ａ（ｘ）］をシステム評価関数と呼ぶ。【００２１】
実際に興味あるシステムについて検討する場合、しばしば、Ｘの確率密度関数ｆ（ｘ，ｚ）が未知であることが問題となる。しかしながら、多くの場合、Ｘの確率密度関数ｆ（ｘ，ｚ）自体は未知であっても、Ｘの分布を生み出す挙動自体は既知（又はモデル化可能）である。【００２２】
この場合、モンテカルロ・シミュレーション・システムは、モンテカルロ・シミュレーションによって、システムの挙動の実現値ｘを多数生成して、Ｘの分布、すなわちｆ（ｘ，ｚ）を模擬することで、システム評価関数Ａ（ｘ）の値を近似的に算出できる。なお、生成されたシステムの挙動の実現値ｘそれぞれをパスとも呼ぶ。モンテカルロ・シミュレーションは、分布が既知である（複数の）乱数から、システムの挙動の実現値ｘを算出する計算を含む。以下、これを定式化する。【００２３】
Ｗ＝（Ｗｊ）をＭ次元の確率変数とし（ただし、ｊ＝１，…，Ｍ）、確率変数Ｗの実現値をｗ＝（ｗｊ）と書く。また、Ｗの（同時）確率密度関数をｇ（ｗ）と書く。ここで、ｇ（ｗ）は既知であり、また計算機上でこの分布に従う乱数を容易に生成できるものとする。ｗは、例えば、［０，１］上の一様分布、又は、標準正規分布Ｎ（０，１）等に従うＭ個の独立な乱数からなるベクトルである。このとき、パラメータｚ＝（ｚｊ）のもとで、乱数ベクトルｗから、対応するＸの実現値ｘを計算する方法、すなわち、関数ｘ＝ｘ（ｗ,ｚ）が与えられているとする。このとき、以下の数１が成立する。【００２４】
【数１】
2016030938.tif 000003
【００２５】
ただし、Ωは、乱数ベクトルｗの分布範囲、つまりｇ（ｗ）のサポート（台）を表す。また、最右辺の記号ｄｗはｄｗ１・ｄｗ２・…・ｄｗＭの略記であり、積分記号はＭ次元の多重積分を表す。モンテカルロ・シミュレーションにおいて１つのパスをシミュレーションすることは、関数ｘ＝ｘ（ｗ,ｚ）を計算することに相当する。【００２６】
さらに、Ｌ個のパスによりモンテカルロ・シミュレーションを行う場合、生成されたＬ個の乱数ベクトルをそれぞれｗｋ＝（ｗｊｋ）と書けば（ただし、ｊ＝１，…，Ｍ、ｋ＝１，…，Ｌ、Ｌは大数の法則を満足する数）、大数の法則から、システム評価関数は、例えば、以下の数２のように、期待値を多数（Ｌ個）のサンプルｘ（ｗｋ，ｚ）の平均値で近似することにより算出される。【００２７】
【数２】
2016030938.tif 000004
【００２８】
本実施例の目的の１つは、以下の数３で表される、システム評価関数をＮ個のパラメータ（ｚｉ）それぞれで微分した微分係数を効率的に算出する方法を提供することである。当該微分係数を、システム評価関数のパラメータに対する感度と呼ぶ。【００２９】
【数３】
2016030938.tif 000005
【００３０】
まず、本実施例の基本的なアイディアを説明する。本実施例の基本的なアイディアは、与えられたパラメータのセットｚのもとで生成したシステムの挙動（実現値）ｘを、パラメータｚｉを微小変化させた後にも利用しつづけることである。多くのシステムにおいてｘ（ｗ，ｚ）は連続であり、このとき、パラメータの所定の閾値以下である微小変化Δｚｉに対応して、それを打ち消すような乱数ベクトルの所定の閾値以下である微小変化Δｗが存在する場合を考える。つまり、与えられたΔｚｉに対して以下の数４を満たすようなΔｗが定まる場合を考える。【００３１】
【数４】
2016030938.tif 000006
【００３２】
このとき、ｗ＋Δｗの分布は、確率密度関数ｇ（ｗ）には従っていないため、例えば数２のような単純な平均によってＡ（ｚ＋Δｚｉ）を計算することはできない。しかしながら、ｗ＋Δｗの分布がｇ（ｗ）からどれくらい異なっているのかを勘案して、ａ（ｘ（ｗ，ｚ））にそれぞれ適切な重みをつけたうえで平均をとることで、Ａ（ｚ＋Δｚｉ）を算出できるというのが、本実施例の基本的なアイディアである。本実施例の手法は、重点サンプリング（Ｉｍｐｏｒｔａｎｃｅ  Ｓａｍｐｌｉｎｇ）手法の一例と考えることもでき、当該重みは重要度に相当する。【００３３】
具体的には、以下の数５で表される、各パスｘに対してｘを一定に保つようなΔｚｉとΔｗとの関係を算出し、当該関係に基づいて重みを算出する。【００３４】
【数５】
2016030938.tif 000007
【００３５】
数５は、ベクトル（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔの各要素が右辺で表されることを意味している。ただし、前述したように、ｘは無限次元の変数になりえるため、この表記は形式的なものである。この点については後述する。【００３６】
実際にモンテカルロ・シミュレーションの対象であるシステムにおいては、多くの場合、ｗ、ｚからｘを計算する構成的なアルゴリズム、即ち関数ｘ＝ｘ（ｗ，ｚ）の具体的な形が分かっている。このとき、ユーザ等は数５の右辺が示示す（ｗとｚの）関数の具体的な形を解析的に算出することができる。したがって、（∂ｗ／∂ｚｉ）｜（ｘ＝ｃｏｎｓｔ）は、例えば、当該解析的に算出された関数に乱数ベクトルとパラメータを代入することにより算出される。また、関数ｘ（ｗ,ｚ）の具体的な形がわからないため、数５の右辺が示す関数の具体的な形を解析的に算出することはできない場合でも、パラメータの微小変化Δｚｉに対して、ｘを一定に保つような、Δｗを求めることにより、（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔ＝Δｗ／Δｚｉと近似的に計算することができる。【００３７】
以下、上記したアイディアを具体的に定式化する。Δｚｉをｉ番目のパラメータの微小変化として、パラメータがｚ＋Δｚｉ＝（ｚ１,ｚ２,…,ｚｉ＋Δｚｉ,…,ｚＮ）であるときの、システム評価関数Ａ（ｚ＋Δｚｉ）は、以下の数６のように表される。【００３８】
【数６】
2016030938.tif 000008
【００３９】
ただし、（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔは、前述したようにｘ＝ｘ（ｗ,ｚ）＝ｘ（ｗ＋Δｗ,ｚ＋Δｚｉ）を満たす、パラメータの微小変化Δｚｉと、乱数ベクトルの微小変化Δｗとの関係を表す。また、Ｉは単位行列を、ｔｒは行列の跡（ｔｒａｃｅ）を、中点・は内積を表す。｜ｄｗ’／ｄｗ｜はｗからｗ’への変数変換によるヤコビ行列式を、Ω’は変数変換によるΩの像を表す。したがって、以下の数７が成立する。【００４０】
【数７】
2016030938.tif 000009
【００４１】
ただし、数７の積分が存在するものとする。ここで、以下の数８で定義されるＮ変数関数ｈ（ｗ，ｚ）＝（ｈｉ（ｗ,ｚ））を（ただし、ｉ＝１，…，Ｎ）、等価微分重要度と呼ぶ。【００４２】
【数８】
2016030938.tif 000010
【００４３】
数７より、システム評価関数のパラメータに対する微分係数、すなわち感度は、等価微分重要度によって重みづけされたシステムの挙動の評価値、の期待値に基づいて算出される。また、数７中のＲｉは、数７の積分範囲をΩ’ではなくΩとしたことにともなう差し引きを表す項であり，以下の数９のように表される。【００４４】
【数９】
2016030938.tif 000011
【００４５】
Ｒｉを境界残余と呼ぶ。ただし∂ΩはΩの境界（表面）を、ｎは境界∂Ωの外向き単位法線ベクトルを表す。数９より、Ωの境界（表面）∂Ωで（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔ＝０（ベクトル）であれば、Ｒｉ＝０である。この場合、数７により、以下の数１０が成立する。【００４６】
【数１０】
2016030938.tif 000012
【００４７】
ここで、代表的な乱数ベクトルｗの分布における等価微分重要度ｈｉ（ｗ,ｚ）について述べる。ｗが一様分布の乱数ベクトルであるとき、ｈｉ（ｗ，ｚ）は以下の数１１で表される。【００４８】
【数１１】
2016030938.tif 000013
【００４９】
また、ｗがＭ次元独立標準正規分布の乱数ベクトルであるとき，ｈｉ（ｗ，ｚ）は以下の数１２で表される。【００５０】
【数１２】
2016030938.tif 000014
【００５１】
また、モンテカルロ・シミュレーションにより、数７が示す、システム評価関数Ａ（ｚ）のパラメータｚに対する微分係数は、以下の数１３のように近似的に算出される。【００５２】
【数１３】
2016030938.tif 000015
【００５３】
とくに，Ωの境界（表面）∂Ωで（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔ＝０の場合は、数１０より、システム評価関数Ａ（ｚ）のパラメータｚに対する感度は、以下の数１４で表される。【００５４】
【数１４】
2016030938.tif 000016
【００５５】
数２、数１３、数１４より、初期パラメータｚについてのみ、モンテカルロ・シミュレーションを１セット行って、Ｌ個のパス｛ｘ（ｗｋ，ｚ）｝と、Ｌ個のパスそれぞれにおける等価微分重要度｛ｈ（ｗｋ,ｚ）｝を計算すれば（ただし、ｋ＝１，…，Ｌ）、ｚにおけるシステム評価関数Ａ（ｚ）と同時に、システム評価関数Ａ（ｚ）の全てのパラメータに対する微分係数を計算できる。すなわち、本実施例のモンテカルロ・シミュレーション・システムは、微小変化後の各パラメータｚ＋Δｚｉについてのモンテカルロ・シミュレーションを行う必要がない。【００５６】
さらに、数２および数１４より、システムの挙動の評価関数ａ（ｘ）はモンテカルロ・シミュレーション終了後に自由に設定できることがわかる。例えば、複雑なシステムをシミュレーション対象とする場合、実際にシステムの挙動をシミュレーションしてみてはじめて何について評価するべきか、どのような現象に注目するべきかが分かる、ということがしばしばあるため、この特徴は望ましい。【００５７】
また，数１３、数１４は、通常の期待値計算の形をしているため、Ｌ回のシミュレーションによる微分係数の推定値の分散は、Ｏ（Ｌ−１）で減少する。一方、パラメータの数だけ単純にモンテカルロ・シミュレーションをやり直す方法による微分係数の推定値の分散は、シミュレーション回数に対してＯ（Ｌ−１/４）、又はＯ（Ｌ−１/３）でしか減少しない。つまり、本実施例のシミュレーションによる微分係数の推定値の分散の収束のスピードが速い。【００５８】
なお、数９で定義される境界残余Ｒｉは、ｘ（ｗ,ｚ）＝ｘ（ｗ＋Δｗ,ｚ＋Δｚｉ）を満たすようなｗ＋Δｗの分布範囲が、ｗの分布範囲（すなわち、ｗが従う確率密度関数の定義域）Ωとずれてしまうことで生じる。つまり、出現する可能性がある全ての乱数ベクトルの実現値ｗ∈Ωに対して、パラメータの微小変化Δｚｉを打ち消すような、乱数ベクトルの実現値ｗ＋Δｗを考えたときに、ｗ＋Δｗの分布範囲Ω’＝｛ｗ＋Δｗ；ｗ∈Ω｝がΩと一致していれば境界残余Ｒｉは生じない。【００５９】
モンテカルロ・シミュレーションの対象となる実際のシステムを想定すると、システムのどのような挙動ｘについても、パラメータｚｉが微小変化したときに、それを打ち消すようなノイズ（乱数ベクトルの実現値）が存在する場合が多い。この場合には、境界残余Ｒｉは生じない。例えば、正規分布に従う乱数ベクトルのように、ｗの分布範囲Ωが実数全体にわたっている場合には、Ωに境界が存在しないため境界残余Ｒｉはゼロである。【００６０】
このとき、モンテカルロ・シミュレーション・システムは、境界残余Ｒｉを算出しなくてもよいため、計算時間及び計算リソースを削減することができる。以下、特に断らない限り乱数ベクトルｗの分布範囲ΩはΩ＝Ω’を満たす、即ち境界残余Ｒｉがゼロであるものとする。【００６１】
ここで、本実施例の感度計算手法が適用できる条件の例を説明する。前述の基本的なアイディアの説明において、パラメータの微小変化Δｚｉに対してｘを一定に保つようなΔｗを考える、と書いたが、実際には、与えられたΔｚｉに対して数４を満たすようなΔｗが存在するとは限らない。【００６２】
前述したように、例えば、ｘは無限次元のベクトルにもなりえる、つまりｘの次元はｗの次元Ｍよりも大きくなりえる。この場合、数４において、制約条件の数が変数Δｗの数（次元数）よりも多くなるため、数４を満たすΔｗは一般には存在しない。このように、与えられたパラメータの微小変化Δｚｉに対して、数４を満たす乱数ベクトルの微小変化Δｗが存在しない場合には本実施例の手法は適用できない。【００６３】
以下、本実施例の手法が適用可能、すなわち、パラメータの微小変化Δｚｉに対して数４を満たすような乱数ベクトルの微小変化Δｗが存在する例を説明する。【００６４】
システムの挙動を表す時系列ｘ（ｔ）が、例えば、以下の数１５のように表せる場合、本実施例の手法が適用可能である。【００６５】
【数１５】
2016030938.tif 000017
【００６６】
数１５は、乱数ベクトルｗとパラメータｚで定まる値ｙ（ｗ,ｚ）が定まると、時刻ｔでのシステムの状態ｘ（ｔ）が、時刻ｔとｙから定まることを示す。なお、この場合、明らかに、以下の数１６が成立する。【００６７】
【数１６】
2016030938.tif 000018
【００６８】
また、システムの挙動を表す時系列ｘ（ｔ）が、例えば、以下の数１７のような漸化式で記述される場合、本実施例の手法が適用可能である。【００６９】
【数１７】
2016030938.tif 000019
【００７０】
数１７は、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）が、時刻ｔと、過去（時刻ｔ以前）のシステムの挙動と、時刻ｔで新たに生成され、時刻毎に異なる乱数ｗｊ（ｗｊは複数個の乱数を含む乱数ベクトルでもよい）と、パラメータｚと、から定まることを示す。なお、この場合、以下の数１８が成立する。【００７１】
【数１８】
2016030938.tif 000020
【００７２】
数１８は、パラメータの微小変化Δｚｉを打ち消す乱数の微小変化Δｗｊが、初期時刻ｔ＝０からシミュレーションの終了時刻ｔ＝Ｔまで順に定まることを示す。【００７３】
以下、本実施例のモンテカルロ・シミュレーション・システムがシミュレーションを実行し、感度を算出する実践的側面について説明する。図１は、本実施例のモンテカルロ・シミュレーション・システム１００の構成例を示す。モンテカルロ・シミュレーション・システム１００は、例えば、ＣＰＵ１１０、メモリ１２０、入出力インタフェース１３０、二次記憶装置１４０を含む計算機上に構成される。また、図１において、モンテカルロ・シミュレーション・システムは１つの計算機から構成されているが、複数の計算機から構成されてもよい。【００７４】
ＣＰＵ１１０は、プログラムに従って動作するプロセッサ及び／又は論理回路を含み、データの入力／出力、読み込み／書き込みを行い、さらに、後述する各プログラムを実行する。メモリ１２０は、ＣＰＵ１１０が実行するプログラム及びデータを一時的にロードして記憶し、さらに各プログラム及び各データを保持する。入出力インタフェース１３０は、外部装置等からデータ等の入力、及び外部装置等に対してデータ等の出力を行うインタフェースである。二次記憶装置１４０は、不揮発性の記憶装置であり、例えばＨＤＤ等である。プログラムやデータ等は二次記憶装置１４０に格納されてもよい。【００７５】
プログラムはＣＰＵ１１０によって実行されることで、定められた処理をメモリ１２０及び入出力インタフェース１３０を用いながら行う。従って、本実施例及び他の実施例においてプログラムを主語とする説明は、ＣＰＵ１１０を主語とした説明でもよい。若しくは、プログラムが実行する処理は、そのプログラムが動作する計算機及び計算機システムが行う処理である。【００７６】
ＣＰＵ１１０は、プログラムに従って動作することによって、所定の機能を実現する機能部として動作する。例えば、ＣＰＵ１１０は、後述する乱数発生部２１０に従って動作することで乱数発生部として機能し、後述するモデル挙動更新部２２０に従って動作することでモデル挙動更新部として機能する。さらに、ＣＰＵ１１０は、各プログラムが実行する複数の処理のそれぞれを実現する機能部としても動作する。計算機及び計算機システムは、これらの機能部を含む装置及びシステムである。プログラムは二次記憶装置１４０に格納されていてもよい。【００７７】
メモリ１２０は、１又は複数のシナリオ・シミュレータ２００と、記憶部３００と、後処理計算部４００と、モンテカルロ・シミュレーション全体制御部５００と、を含む。シナリオ・シミュレータ２００は、異なる乱数ベクトルを用いてＬ回のモンテカルロ・シミュレーションを行う。記憶部３００は、シナリオ・シミュレータ２００から出力された数値等を記憶する。後処理計算部４００は、シミュレーション終了後に期待値および感度の計算を行う。【００７８】
モンテカルロ・シミュレーション全体制御部５００は、モンテカルロ・シミュレーション・システム１００全体の制御を行う。また、モンテカルロ・シミュレーション全体制御部５００は、ユーザ等が指定したシミュレーション回数によって算出されるシステム評価関数及び感度の推定値の計算精度を算出し、出力してもよい。また、モンテカルロ・シミュレーション全体制御部５００は、ユーザ等が指定した、システム評価関数及び感度の推定値の計算精度を保証するために必要なシミュレーション回数等を算出、出力してもよい。モンテカルロ・シミュレーション全体制御部５００は、例えば、指定された又は算出したシミュレーション回数、例えばＬ回、シナリオ・シミュレータ２００を起動し、Ｌ回のシミュレーションが終了すると、後処理計算部４００を起動する。【００７９】
シナリオ・シミュレータ２００は、ユーザ等から与えられるシミュレーションモデルの１次元又は多次元のパラメータのベクトルｚを入力として、シミュレーションに必要な乱数ベクトルを内部で生成しながらモデルのシミュレーションを実行する。また、シナリオ・シミュレータ２００は、システムの挙動の評価関数ａ（ｋ）と、感度計算に必要な等価微分重要度ｈｉ（ｋ）を出力する。【００８０】
シナリオ・シミュレータ２００は、それぞれプログラムである、乱数発生部２１０、モデル挙動更新部２２０、乱数重要度算出部２３０、指標算出部２６０、等価微分重要度算出部２７０、及びシミュレーション制御部２８０を含む。また、シナリオ・シミュレータ２００は、それぞれデータを格納する領域である、モデル挙動記憶部２４０、及び乱数重要度記憶部２５０を含む。【００８１】
乱数発生部２１０は、予め定められた、又はユーザ等によって指定されたシミュレーションモデルに対応する確率密度関数に従う乱数ベクトルを生成する。モデル挙動更新部２２０は、ユーザ等によって与えられたパラメータと乱数発生部２１０が生成した乱数ベクトルとを入力として、入力されたパラメータと乱数ベクトルとを予め定められた、又はユーザ等によって指定されたシミュレーションモデル（関数）に適用し、システムの挙動の実現値をシミュレートする。【００８２】
乱数重要度算出部２３０は、パラメータと乱数発生部２１０が生成した乱数ベクトルとを入力として、入力された各パラメータについて、ｗｊ（ｋ）の各パラメータに関する乱数重要度、すなわち数８の総和記号の中の項を算出する。モデル挙動記憶部２４０は、モデル挙動更新部２２０が算出したシステムの挙動の実現値を格納する。乱数重要度記憶部２５０は、乱数重要度算出部２３０が算出した乱数重要度を格納する。指標算出部２６０は、例えば、モデル挙動記憶部２４０に格納されたシステムの挙動の実現値の平均を算出することにより、システムの挙動の評価関数を算出する。【００８３】
等価微分重要度算出部２７０は、例えば、乱数重要度記憶部２５０に格納された乱数重要度の和を算出することにより等価微分重要度を算出する。シミュレーション制御部２８０は、シナリオ・シミュレータ２００全体の動作を制御する。【００８４】
記憶部３００は、指標算出部２６０が算出したシステムの挙動の評価関数ａ（ｋ）を格納する指標値記憶部３１０と、等価微分重要度算出部２７０が算出した等価微分重要度を格納する等価微分重要度記憶部３２０と、を含む。【００８５】
後処理計算部４００は、それぞれプログラムである、期待値計算部４１０、及び感度計算部４２０を含む。期待値計算部４１０は、指標値記憶部３１０に格納されたＬ個のシステムの挙動の評価関数ａ（ｋ）の平均値を計算し、当該平均値を期待値、すなわちシステム評価関数Ａ（ｚ）として出力する。このとき、期待値計算部４１０は、さらにＬ個のシステムの挙動の評価関数ａ（ｋ）の分散を計算し、Ａ（ｚ）の推定値の精度（信頼区間）を同時に出力してもよい。【００８６】
感度計算部４２０は、指標値記憶部３１０に格納されたシステムの挙動の評価関数ａ（ｋ）と、等価微分重要度記憶部３２０に格納されたＮ個のパラメータそれぞれに関する等価微分重要度ｈｉ（ｋ）の組について、それぞれ対応する値同士の積を計算する。さらに、感度計算部４２０は、算出した積の平均値を計算して、Ｎ個のパラメータそれぞれに関するシステム評価関数Ａ（ｚ）の感度∂Ａ／∂ｚｉとして出力する。このとき、感度計算部４２０は、さらにａ（ｋ）とｈｉ（ｋ）の積の分散を計算し、感度∂Ａ／∂ｚｉの推定値の精度（信頼区間）を同時に出力してもよい。【００８７】
図２は、乱数重要度算出部２３０の構成例を示す。乱数重要度算出部２３０は、いずれもプログラムである、サンプルパス固定微分計算部２３１、サンプルパス固定２階微分計算部２３２、及び乱数重要度計算部２３３を含む。【００８８】
サンプルパス固定微分計算部２３１は、パラメータと、乱数発生部２１０が生成した乱数ベクトルと、を入力として（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出する。サンプルパス固定２階微分計算部２３２は、パラメータと、乱数発生部２１０が生成した乱数ベクトルと、を入力として∂Ｄｉ，ｊ／∂ｗｊを算出する（Ｄｉ，ｊ（ｋ）＝（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔ）。サンプルパス固定２階微分計算部２３１は、サンプルパス固定微分計算部２３１が算出した（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔと、パラメータと、を入力として、∂Ｄｉ，ｊ／∂ｗｊを算出してもよい。また、サンプルパス固定２階微分計算部２３２が（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出してもよく、このとき、サンプルパス固定微分計算部２３１は（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出しなくてもよい。【００８９】
乱数重要度計算部２３３は、乱数発生部２１０から入力された乱数ベクトルと、サンプルパス固定微分計算部２３１又はサンプルパス固定２階微分計算部２３２が算出した（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔと、サンプルパス固定２階微分計算部２３２が算出した∂Ｄｉ，ｊ／∂ｗｊと、を入力として、各パラメータについて、ｗｊ（ｋ）の各パラメータに関する乱数重要度、すなわち数８の総和記号の中の項を算出する。【００９０】
図３は、本実施例におけるモンテカルロ・シミュレーションの全体処理の一例を示す。モンテカルロ・シミュレーションの全体処理は、シナリオ・シミュレータ２００が異なる乱数ベクトルを用いてＬ回のシミュレーションを行うシミュレーション実行処理と、シミュレーション終了後に後処理計算部４００が行う期待値・感度計算処理と、を含む。【００９１】
シナリオ・シミュレータ２００は、シミュレーション実行処理を行う（Ｓ３０１）。シミュレーション実行処理は、１セットのモンテカルロ・シミュレーションを構成するＬ回のシミュレーションを実行する処理である。シナリオ・シミュレータ２００は、異なる乱数ベクトルを用いて個々のシミュレーションを実行する。【００９２】
なお、本実施例のモンテカルロ・シミュレーション・システムは、ｋ回目のシミュレーションにより得られたシステムの挙動の実現値ｘ（ｋ）から、システムの挙動の評価関数ａ（ｋ）に加えて、Ｎ個のパラメータそれぞれについて等価微分重要度ｈｉ（ｋ）を計算する。なお、シミュレーション実行処理におけるＬ回のシミュレーションそれぞれは、互いに独立であるため、Ｌ回のシミュレーションを逐次的に行うのではなく、例えば、多数のシナリオ・シミュレータ２００又は多数の計算機を用いて同時に並列実行してもよい。等価微分重要度ｈｉ（ｋ）を算出する処理を定める一連の関数を、重み関数と呼ぶ。【００９３】
続いて、後処理計算部４００は、期待値・感度計算処理を行う（Ｓ３０２）。期待値・感度計算処理は、１セットのモンテカルロ・シミュレーションを構成するＬ回のシミュレーションの実行後に、後処理計算部４００が、注目指標の期待値の推定値と、注目指標の期待値のパラメータに対する感度の推定値と、を計算する処理である。【００９４】
注目指標の期待値、すなわちシステム評価関数Ａ（ｚ）の推定値は、期待値計算部４１０が、指標値記憶部に格納されたＬ個のシステムの挙動の評価関数ａ（ｋ）の平均を計算することで得られる。本実施例におけるモンテカルロ・シミュレーション・システムは、これに加えて、注目指標の期待値のＮ個のパラメータそれぞれに対する感度の推定値を計算する。【００９５】
具体的には、感度計算部４２０は、システム評価関数Ａ（ｚ）のｉ番目のパラメータそれぞれに対する感度の推定値を、指標値記憶部３１０に格納された対応するシステムの挙動の評価関数ａ（ｋ）に、等価微分重要度記憶部３２０に格納されたｉ番目のパラメータについての等価微分重要度ｈｉ（ｋ）を掛けることで重みづけをしたものの平均を計算することで得る。なお、モンテカルロ・シミュレーション・システムは、例えば、ユーザ等が指定したパラメータについてのみ、ステップＳ３０１における等価微分重要度の算出処理、及びステップＳ３０２における感度の推定値の算出処理を行ってもよい。【００９６】
なお、前述のように、例えば、システムの挙動を表す時系列ｘ（ｔ）が、数１５又は数１７の形で表される場合、本実施例における注目指標の期待値のパラメータに対する感度の推定手法が適用できる。図３、図４は、それぞれ、システムの挙動を表す時系列ｘ（ｔ）が、数１５、数１７の形で表されたモデルで算出される場合における、ｋ回目のシミュレーション実行処理の詳細な処理内容を示す。【００９７】
図４は、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、すなわち、乱数ベクトルｗとパラメータｚで定まる値ｙ（ｗ,ｚ）を定めると、時刻ｔでのシステムの状態ｘ（ｔ）が，時刻ｔとｙのみで決定論的に定まる場合について、シミュレーション実行処理の詳細な処理内容の一例を示す。【００９８】
まず、シナリオ・シミュレータ２００は、乱数ベクトルｗ（ｋ）を生成し、ｗ（ｋ）とパラメータｚとからｙ（ｋ）を計算する（Ｓ４０１）。具体的には、乱数発生部２１０は、予め定められた確率密度関数ｇ（ｗ）にしたがうＭ次元の乱数ベクトルｗ（ｋ）を生成し、生成した乱数をモデル挙動更新部２２０及び乱数重要度算出部２３０に送信する。なお、乱数発生部２１０は、例えば、既知の疑似乱数発生アルゴリズム、又は物理乱数発生装置といった公知の手段を用いて乱数ベクトルを生成すればよい。モデル挙動更新部２２０は、受信したｗ（ｋ）と、与えられたＮ次元のパラメータ値ｚを、予め与えられた関係式ｙ（ｗ，ｚ）に代入してｙ（ｋ）を計算する。【００９９】
続いて、シナリオ・シミュレータ２００は、Ｎ個のパラメータそれぞれについて等価微分重要度ｈｉ（ｋ）を計算する（Ｓ４０２）。これは、注目指標の期待値のパラメータに対する感度の推定のために必要な処理であり、重み関数が定める処理は当該処理を含む。なお、シナリオ・シミュレータ２００は、ユーザ等が指定したパラメータについてのみ、ステップＳ４０２における処理を行ってもよい。ｉ番目のパラメータについての等価微分重要度ｈｉ（ｋ）は、具体的には以下の４つのステップによって算出される、すなわちステップＳ４０２は以下の４つのステップを含む。【０１００】
第１のステップは、サンプルパス固定微分計算部２３１が、システムの挙動ｘを一定に保つような、乱数ベクトルｗ（ｋ）のｊ番目の要素とｉ番目のパラメータの微小変化の関係Ｄｉ，ｊ（ｋ）＝（∂ｗｊ（ｋ）／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出するステップである。ここで、添え字ｊは、Ｍ次元の乱数ベクトルｗ（ｋ）のｊ番目の要素についての計算を表す。【０１０１】
システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、数１６が成り立つから、サンプルパス固定微分計算部２３１は、（∂ｗｊ／∂ｚｉ）｜（ｙ＾（ｋ）＝ｃｏｎｓｔ）を求めればよい。つまり、サンプルパス固定微分計算部２３１は、例えば、関数−（∂ｙ（ｋ）／∂ｚｉ）／（∂ｙ（ｋ）／∂ｗｊ（ｋ））がユーザ等によって与えられた場合、入力されたパラメータｚｉ及び入力された乱数ベクトルｗｊ（ｋ）を当該関数に代入することにより、Ｄｉ，ｊ（ｋ）を算出すればよい。サンプルパス固定微分計算部２３１は、当該関数から高速にＤｉ，ｊ（ｋ）を算出することができる。【０１０２】
また、サンプルパス固定微分計算部２３１は、例えば、パラメータの微小変化Δｚｉに対して、ｙを一定に保つようなΔｗを数値的に探索してΔｗ／Δｚｉを計算することにより、（∂ｗｊ／∂ｚｉ）｜（ｙ＾（ｋ）＝ｃｏｎｓｔ）、すなわちＤｉ，ｊ（ｋ）を近似的に算出してもよい。例えば、ユーザ等が関数ｘ及びｙの具体的な形がわからない場合であっても、サンプルパス固定微分計算部２３１は、上述の方法により、Ｄｉ，ｊ（ｋ）を算出することができる。サンプルパス固定微分計算部２３１は算出したＤｉ，ｊ（ｋ）を、サンプルパス固定２階微分計算部２３２及び乱数重要度計算部２３３に送信する。【０１０３】
第２のステップは、サンプルパス固定２階微分計算部２３２が、受信したＤｉ，ｊ（ｋ）を、ｗｊ（ｋ）で微分するステップである。Ｄｉ，ｊ（ｋ）を、ｗｊ（ｋ）で微分した結果をＤＤｉ，ｊ（ｋ）と表記する。サンプルパス固定２階微分計算部２３２は、例えば、受信したＤｉ，ｊ（ｋ）を数値的に微分することによって、ＤＤｉ，ｊ（ｋ）を算出すればよい。また、サンプルパス固定２階微分計算部２３２は、例えば、ＤＤｉ，ｊ（ｋ）を算出する関数が与えられている場合、入力されたパラメータｚｉ及び入力された乱数ベクトルｗｊ（ｋ）当該関数に代入することによって、ＤＤｉ，ｊ（ｋ）を算出してもよい。サンプルパス固定２階微分計算部２３２は、算出したＤＤｉ，ｊ（ｋ）を、乱数重要度計算部２３３に送信する。【０１０４】
第３のステップは、乱数重要度計算部２３３が、乱数ベクトルｗ（ｋ）のｊ番目の要素ｗｊ（ｋ）に関する等価微分重要度Ｈｉ，ｊ（ｋ）を計算するステップである。第３のステップは、数８における総和記号の中の項を計算することに相当する。具体的には、例えば、乱数重要度計算部２３３は、受信したＤｉ，ｊ（ｋ）と、受信したＤＤｉ，ｊ（ｋ）と、予め定められた確率密度関数ｇ（ｗ）から一意に定まる関数ＬＧに乱数ベクトルｗｊ（ｋ）を代入した値ＬＧj（ｗｊ（ｋ））＝（１／ｇ（ｗｊ（ｋ）））×（∂ｇ／∂ｗｊ（ｋ））と、を図中に示した関係式に代入して乱数重要度Ｈｉ，ｊ（ｋ）を計算すればよい。乱数重要度計算部２３３は、算出したＨｉ，ｊ（ｋ）を乱数重要度記憶部２５０に格納する。【０１０５】
なお、前述のように、第３のステップにおける計算は、Ｍ次元のｗが一様分布の乱数ベクトルである場合、又は、ｗがＭ次元独立標準正規分布の乱数ベクトルである場合には、それぞれ数１１又は数１２のように表される。従って、例えば、乱数ｗｊ（ｋ）が一様分布に従う場合、数１１より、乱数重要度はＤＤｉ，ｊ（ｋ）で表されるため、乱数重要度算出部２３０は、第１のステップのみを行うことにより乱数重要度を算出することができる。シナリオ・シミュレータ２００は、第１のステップから第３のステップにおける処理を、全てのｊ（ただし、ｊ＝１…Ｍ）について行う。【０１０６】
第４のステップは、等価微分重要度算出部２７０がｉ番目のパラメータに関する等価微分重要度ｈｉ（ｋ）を計算するステップである。第４のステップは、前述の数８における総和計算に相当する。等価微分重要度算出部２７０は、乱数重要度記憶部２５０に格納された乱数重要度Ｈｉ，ｊ（ｋ）を取得し、ｊについての総和を計算する。なお、各パラメータｚｉについての等価微分重要度ｈｉ（ｋ）の計算は互いに独立であるため、１つのシナリオ・シミュレータ２００によってＮ個のパラメータｚ１，…，ｚNについて等価微分重要度の計算を逐次的に行うのではなく、例えば、多数の計算機を用いて同時に並列に当該計算を実行してもよい。【０１０７】
続いて、モデル挙動更新部２２０は、は、実際にシステムの挙動ｘ（ｋ）（ｔ）をシミュレートし計算する（Ｓ４０３）。モデル挙動更新部２２０は、時刻ｔを初期時刻から終了時刻に向かって増加させながら、時刻ｔでのシステムの状態を逐次的に計算することで当該処理を行うことができるが、これによらずその他の公知の効率的なシミュレーションの計算手法を用いてもよい。また、モデル挙動更新部２２０は算出したｘ（ｋ）（ｔ）をモデル挙動記憶部２４０に格納する。なお、ステップＳ４０３における処理は、例えば、ステップＳ４０２における処理と並行して行われてもよいし、ステップＳ４０２における処理の前に行われてもよい。【０１０８】
続いて、指標算出部２６０は、モデル挙動記憶部２４０から、システムの挙動ｘ（ｋ）（ｔ）を取得し、注目指標の評価値、すなわち、システムの挙動の評価関数ａ（ｋ）を計算する（Ｓ４０４）。指標算出部２６０は、算出したａ（ｋ）を指標値記憶部３１０に格納する。【０１０９】
以上のように、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される場合、図４に示した処理により、シナリオ・シミュレータ２００はシステム挙動の評価関数ａ（ｋ）、およびＮ個のパラメータそれぞれについての等価微分重要度ｈｉ（ｋ）を計算することができる。【０１１０】
図５は、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される場合、すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）が、時刻ｔと、過去（時刻ｔ以前）のシステムの状態と、時刻ｔで新たに生成する（時刻毎に異なる）乱数ｗｊ（複数個の乱数を含む乱数ベクトルもよい）と、パラメータｚと、によって定まる場合の、シミュレーション実行処理の詳細な処理内容を示す。【０１１１】
シナリオ・シミュレータ２００は、時刻ｔを、シミュレーション開始時刻ｔ＿ｓｔａｒｔから、シミュレーション終了時刻ｔ＿ｅｎｄまで、シミュレーション時間ステップΔｔずつ増加させながら、時刻ｔにおけるシステムの状態ｘ（ｋ）（ｔ）と、Ｎ個のパラメータ全てについてｊ番目の乱数ｗｊ（ｋ）に関する乱数重要度Ｈｉ，ｊ（ｋ）と、を計算する処理を繰り返し行うことでシミュレーションを行う。さらに、シナリオ・シミュレータ２００は、シミュレーションの結果得られたシステムの挙動ｘ（ｋ）（ｔ）の時系列を用いて、システムの挙動の評価関数ａ（ｋ）、およびＮ個のパラメータそれぞれについての等価微分重要度ｈｉ（ｋ）を計算する。【０１１２】
まず、シミュレーション制御部２８０は、時刻ｔをシミュレーション開始時刻ｔ＿ｓｔａｒｔにセットし、時刻ｔに関するループ処理の開始をシナリオ・シミュレータ２００の各部に通知する（Ｓ５０１）。時刻ｔに関するループ処理は、モデル挙動更新部２２０が時刻ｔでのシステムの状態ｘ（ｋ）（ｔ）を計算する処理（Ｓ５０２）と、乱数重要度算出部２３０が時刻ｔで生成した乱数のＮ個のパラメータそれぞれに関して乱数重要度Ｈｉ，ｊ（ｋ）を計算する処理（Ｓ５０３）と、を含む。【０１１３】
ステップＳ５０２における具体的な処理について説明する。まず、乱数発生部２１０は、時刻ｔでのシステムの状態ｘ（ｋ）（ｔ）の計算に必要な乱数ｗｊ（ｋ）を、予め定められた確率密度関数ｇ（ｗ）にしたがって生成し、モデル挙動更新部２２０、及び乱数重要度算出部２３０に生成した乱数を送信する。モデル挙動更新部２２０は、受信した乱数ｗｊ（ｋ）と過去（時刻ｔ以前）のシステムの挙動とパラメータｚと時刻ｔとを与えられた漸化式（数１７）に代入して、時刻ｔでのシステムの状態ｘ（ｋ）（ｔ）を計算する。【０１１４】
続いて、乱数重要度算出部２３０は、Ｎ個のパラメータｚ１，…，ｚNそれぞれについて、時刻ｔで生成した乱数ｗｊ（ｋ）のパラメータに関する乱数重要度Ｈｉ，ｊ（ｋ）を計算する（Ｓ５０３）。なお、シナリオ・シミュレータ２００は、ユーザ等が指定したパラメータについてのみ、ステップＳ５０３における処理を行ってもよい。【０１１５】
ｉ番目のパラメータについての乱数重要度Ｈｉ、ｊ（ｋ）は、具体的には以下の３つのステップにより算出される、すなわちステップＳ５０３は以下の３つのステップを含む。なお、各パラメータについての等価微分重要度ｈｉ（ｋ）の計算は互いに独立であるため、１つのシナリオ・シミュレータ２００によってＮ個のパラメータｚ１，…，ｚNについて等価微分重要度の計算を逐次的に行うのではなく、例えば多数の計算機を用いて当該計算を同時に並列実行してもよい。【０１１６】
第１のステップは、サンプルパス固定微分計算部２３１が、システムの挙動ｘを一定に保つような、乱数ｗｊ（ｋ）とｉ番目のパラメータの微小変化の関係Ｄｉ，ｊ（ｋ）＝（∂ｗｊ（ｋ）／∂ｚｉ）｜ｘ＿ｊ＋１＾（ｋ）＝ｃｏｎｓｔを求めるステップである。【０１１７】
サンプルパス固定微分計算部２３１は、例えば、関数−（∂ｆ／∂ｚｉ）／（∂ｆ／∂ｗｊ（ｋ））が与えられた場合、入力されたパラメータｚｉ及び入力された乱数ｗｊ（ｋ）を当該関数に代入することにより、Ｄｉ，ｊ（ｋ）を算出すればよい。また、サンプルパス固定微分計算部２３１は、例えば、パラメータの微小変化Δｚｉに対して、ｆを一定に保つようなΔｗを数値的に探索してΔｗ／Δｚｉを計算することにより、（∂ｗｊ／∂ｚｉ）｜（x(t)＝ｃｏｎｓｔ）、すなわちＤｉ，ｊ（ｋ）を近似的に算出してもよい。サンプルパス固定微分計算部２３１は、算出したＤｉ，ｊ（ｋ）を、サンプルパス固定２階微分計算部２３２、及び乱数重要度計算部２３３に送信する。【０１１８】
第２のステップは、サンプルパス固定２階微分計算部２３２が、受信したＤｉ，ｊ（ｋ）を、ｗｊ（ｋ）で微分し、ＤＤｉ，ｊ（ｋ）を算出するステップである。第２のステップは、ステップＳ４０２における第２のステップと同様である。【０１１９】
第３のステップは、乱数重要度計算部２３３が、時刻ｔにおける乱数ｗｊ（ｋ）に関する乱数重要度Ｈｉ，ｊ（ｋ）を計算するステップである。第３のステップは、ステップＳ４０２における第３のステップと同様である。【０１２０】
続いて、シミュレーション制御部２８０は、時刻をｔ＋Δｔにセットする（Ｓ５０４）。時刻ｔ＋Δｔがシミュレーション終了時刻ｔ＿ｅｎｄ以前である場合（Ｓ５０５：Ｙｅｓ）、モンテカルロ・シミュレーション全体制御部５００は、シナリオ・シミュレータ２００を再度起動し、シナリオ・シミュレータ２００は、ステップＳ５０２〜ステップＳ５０４の処理を行う。【０１２１】
時刻ｔ＋Δｔがシミュレーション終了時刻ｔ＿ｅｎｄより後である場合（Ｓ５０５：Ｎｏ）、指標算出部２６０は、モデル挙動記憶部２４０から、システムの挙動ｘ（ｋ）（ｔ）の時系列を取得し、注目指標の評価値、すなわちシステムの挙動の評価関数ａ（ｋ）を計算する（Ｓ５０６）。指標算出部２６０は、算出したａ（ｋ）を指標値記憶部３１０に格納する。【０１２２】
続いて、等価微分重要度算出部２７０は、乱数重要度記憶部２５０から、全てのｊについてＨｉ，ｊ（ｋ）を取得し、ｊについての総和を計算することで、ｉ番目のパラメータに関する等価微分重要度ｈｉ（ｋ）を算出する（Ｓ５０７）。なお、等価微分重要度算出部２７０は、は、ユーザ等が指定したパラメータについてのみ、ステップＳ５０７における処理を行ってもよい。なお、ステップＳ５０６における処理とステップＳ５０７における処理は並行して行われてもよいし、順序を入れ替えて行われてもよい。重み関数が定める処理は、Ｓ５０３及びＳ５０７における処理を含む。【０１２３】
以上のように、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される場合に、図５に示した処理により、シナリオ・シミュレータ２００は、システム挙動の評価関数ａ（ｋ）およびＮ個のパラメータそれぞれについて等価微分重要度ｈｉ（ｋ）を計算することができる。【０１２４】
図６は、実施例１におけるモンテカルロ・シミュレーション・システムのユーザインタフェースの一例である。ユーザインタフェース６００は、シミュレーション条件設定セクション６１０、期待値計算セクション６２０、及び感度計算セクション６３０を含む。【０１２５】
シミュレーション条件設定セクション６１０は、入力受付セクション６１１、６１２、及び表示セクション６１３を含む。入力受付セクション６１１は、シミュレーションモデル、すなわちシステムの挙動の実現値ｘを定める関数ｘ（ｗ，ｔ）等の入力を受け付ける。入力受付セクション６１２は、当該シミュレーションモデルによるシミュレーション回数の入力を受け付ける。表示セクション６１３は、入力受付セクション６１１〜６１２が入力を受け付けた内容のシミュレーションを実行するために必要な計算時間を表示する。【０１２６】
期待値計算セクション６２０は、入力受付セクション６２１、及び表示セクション６２２を含む。入力受付セクション６２１は、注目指標、すなわちシステムの挙動の評価関数ａ（ｘ）の入力を受け付ける。表示セクション６２２は、システム評価関数Ａ（ｘ）の計算精度を表示する。【０１２７】
感度計算セクション６３０は、入力受付セクション６３２、６３４、チェックボックス６３１、６３３、及び表示セクション６３５を含む。入力受付セクション６３２は、感度計算の対象パラメータ、すなわち感度∂Ａ／∂ｚｉを算出する対象となるパラメータｚｉの入力を受け付ける。チェックボックス６３１は、全てのパラメータｚｉを感度計算の対象パラメータとして選択するためのチェックボックスである。入力受付セクション６３４は、Ｄｉ，ｊを計算する方法の入力を受け付ける。チェックボックス６３３は、Ｄｉ，ｊを計算する方法を自動で選択するためのチェックボックスである。表示セクション６３５は、感度∂Ａ／∂ｚｉの計算精度を表示する。【０１２８】
なお、ユーザインタフェース６００における、入力セクション、表示セクションは、それぞれ表示セクション、入力セクションであってもよいし、それぞれ表示セクション、入力セクションを兼ねてもよい。例えば、表示セクション６２２がシステム評価関数Ａ（ｘ）の計算精度の入力を受け付け、表示セクション６３５が感度∂Ａ／∂ｚｉの計算精度の入力を受け付けてもよい。このとき、入力受付セクション６１１は、例えば、入力されたシステム評価関数Ａ（ｘ）の計算精度と、入力された感度∂Ａ／∂ｚｉの計算精度と、を保証するために必要なシミュレーション回数を表示してもよい。【０１２９】
以上、本実施例のモンテカルロ・シミュレーション・システムは、１セットのモンテカルロ・シミュレーションの結果を用いて、システム評価関数の各パラメータに関する感度を算出することができ、ひいては、計算時間及び計算リソースを削減することができる。【０１３０】
また、本実施例のモンテカルロ・シミュレーション・システムが、算出した感度を用いて、ユーザ等は、パラメータの最適化等、すなわちパラメータ値の許容される範囲で、システム評価関数を最大化するパラメータ値の探索等を行うことができる。
【実施例２】
【０１３１】
本実施例は、多数のパラメータが存在する確率的なモデルにおいて、注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって同時に求める手段を、半導体集積回路の電子回路における製造ばらつきを考慮した設計に応用する例を説明する。【０１３２】
本実施例は、システムの挙動を表す時系列ｘ（ｔ）が数１５の形で表される、すなわち、乱数ベクトルｗとパラメータｚで定まる値ｙ（ｗ,ｚ）が定まると、時刻ｔでのシステムの状態ｘ（ｔ）が、時刻ｔとｙのみで決定論的に定まる（図３に相当）一例である。【０１３３】
近年、半導体集積回路の微細化が進むにしたがって、集積回路を構成する個々のトランジスタの特性のばらつきが大きくなっている。集積回路を構成するトランジスタの特性ばらつきは、例えば、集積回路の製造時において、形成される回路要素の物理的なサイズ、基板に注入される不純物イオンの数、又は格子欠陥の数等の製造時に不可避的に生じる確率的な現象に起因する。つまり、当該回路は、個々のトランジスタの特性が確率的にばらつくことを想定して設計される必要がある。【０１３４】
したがって、トランジスタの特性のばらつきを確率分布でモデル化した上で、例えば、回路全体の特性の期待値やばらつき（分散）の大きさ、又は集積回路の歩留まり、すなわち集積回路を構成する回路全体の特性が望ましい範囲におさまる確率、等をモンテカルロ・シミュレーションによって推定することが行われる。回路設計者は、例えば、設計した回路のばらつきを考慮した特性が与えられた仕様を満たすまで、設計パラメータをチューニングしては、モンテカルロ・シミュレーションを行ってばらつきを考慮した特性を確認する、というプロセスを繰り返し行う。【０１３５】
個々のトランジスタにおける設計パラメータは、例えば、ゲート幅Ｗ、ゲート長Ｌ、拡散層幅ＬＯＤ、トランジスタ間距離ＰＤＸといった、ばらつき・電力・レイアウト面積等に関係するものを１０程度含む。したがって、一般に、回路ブロックの設計パラメータの数は、回路を構成するトランジスタ数の１０倍程度になる。このため、多数のトランジスタで構成される複雑な回路ブロックにおける設計パラメータは膨大な数となりえるため、上記のように個々の設計パラメータのチューニングを人手で行うことは困難である。【０１３６】
これに対して、本実施例における手法では、任意の注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって求めることができる。本実施例の手法を用いることで、回路設計者は、単に、現在のパラメータ設定のもとでの回路ブロック全体の特性（性能期待値や歩留まり）や電力の値のみならず、全てのトランジスタの全てのパラメータについて、それらが回路ブロック全体の特性（性能期待値や歩留まり）や電力にどのように影響するのか、その感度を同時に得ることができる。【０１３７】
したがって、回路設計者は、例えば、どのトランジスタのどのパラメータを変えることで、消費電力を増やさずに特性を向上させることができるか等を、やみくもにパラメータをチューニングすることなく知ることができ、設計の効率を大幅に向上できる。【０１３８】
図７は、６個のトランジスタで構成されるＳＲＡＭセルの回路の例を示す。ＳＲＡＭセルは、通常時にはＳＲＡＭセルに記憶された値が失われないことと、記憶の書き換え時にはＳＲＡＭセルの記憶値が意図した値にきちんと書き変わること、という互いに相反する特性を満たす必要がある。【０１３９】
ＳＲＡＭに含まれる各トランジスタのパラメータは、集積回路の製造時のばらつきによって、製造された集積回路における個々のトランジスタの実際のパラメータ値は設計時に指定した値にはならず、例えば、ある確率分布にしたがって分布する。その結果、例えば、通常にＳＲＡＭセルに記憶していた値が失われる、又は、記憶の書き換え時にはＳＲＡＭセルの記憶値が意図した値にきちんと書き変わらない、等の不良が起こりえる。【０１４０】
そこで、回路設計者は、設計パラメータ値をいろいろな値に変えてモンテカルロ・シミュレーションを行い、実際に製造されるＳＲＡＭセル回路の不良率が与えられた閾値以下になるように設計パラメータを最適化する必要がある。この例に、実施例１のモンテカルロ・シミュレーション・システムのシミュレーション手法を適用する。【０１４１】
本実施例におけるパラメータベクトルｚは、トランジスタの設計パラメータからなる。例えば、個々のトランジスタに１０の設計パラメータがある場合、６個のトランジスタで構成されるＳＲＡＭセル回路におけるパラメータベクトルｚは、６０次元のベクトルである。【０１４２】
ところで、前述のように、ＳＲＡＭ回路の特性ばらつきは、トランジスタの実際のパラメータ値が製造時にばらつくことに起因する。すなわち、トランジスタの特性を決める６０次元の設計パラメータｚに対して、実際には製造時のばらつきにより、例えば、６０次元のパラメータｙが生成される。【０１４３】
乱数発生部２１０は、シミュレーション毎に確率密度関数ｇ（ｗ）に従う多次元の乱数ベクトルｗを生成する。モデル挙動更新部２２０は、乱数ベクトルｗと６０次元の設計パラメータベクトルｚと、から製造後の実際のパラメータ値ｙを関数ｙ＝ｙ（ｗ，ｚ）によって定める。このとき、製造後の実際のパラメータ値ｙの確率分布が、設計パラメータベクトルがｚであるときの、製造ばらつきを表すように、乱数ベクトルｗが従う確率密度関数ｇ（ｗ）と、関数ｙ（ｗ，ｚ）が定められているものとする。【０１４４】
モデル挙動更新部２２０は、製造後の実際のパラメータｙを定めれば、ＳＲＡＭセル回路の挙動ｘ（ｔ）を、通常の回路シミュレーションによって決定論的に定めることができる。したがって、本実施例におけるシステムの挙動を表す時系列ｘ（ｔ）は、数１５の形で表される、すなわち、乱数ベクトルｗとパラメータベクトルｚで定まる値ｙ（ｗ,ｚ）を定めると、時刻ｔでのシステムの状態ｘ（ｔ）が，時刻ｔとｙのみで決定論的に定まる。また、システムの挙動の評価関数ａ（ｘ）を、そのセルの挙動が正常動作の範囲内であれば０、不良であれば１と定義すれば、システム評価関数Ａ（ｘ）、すなわち、ａ（ｘ）の期待値Ｅ［ａ（ｘ）］は、ＳＲＡＭセルの不良率を示す。【０１４５】
したがって、モンテカルロ・シミュレーション・システムは、図４で説明した処理を行うことにより、ＳＲＡＭセルの不良率の６０個の設計パラメータ全てに対する感度を算出することができる。また、モンテカルロ・シミュレーション・システムは、は、算出した感度から、勾配法やニュートン法といった公知の最適化手法を用いて最急降下法ＳＲＡＭセル不良率を最小にするような６０個の設計パラメータ値の組み合わせを、算出してもよい。【０１４６】
以上のように、本実施例のモンテカルロ・シミュレーション・システムは、多数のトランジスタで構成される半導体集積回路の電子回路における、製造ばらつきを考慮した設計パラメータ値の最適値を効率的に求めることができる。
【実施例３】
【０１４７】
本実施例は、多数のパラメータが存在する確率的なモデルにおいて、注目指標の期待値の多数のパラメータ全てに対する感度（微分係数）を、１セットのモンテカルロ・シミュレーションによって同時に求める手段を、金融分野、とくに、銀行の収益予測シミュレーションシステムに適用した例を説明する。本実施例における、システムの挙動を表す時系列ｘ（ｔ）が数１７の形で表される、すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）が、時刻ｔと、過去（時刻ｔ以前）のシステムの挙動と、時刻ｔで新たに生成する（時刻毎に異なる）乱数ｗｊと、パラメータｚと、で定まる。【０１４８】
金融機関においては、将来の事業計画の立案時に、取り扱っている金融商品の将来の取扱高やその収益について計画を立てて、当該事業計画の妥当性を経営判断する。【０１４９】
図８は、銀行における収支計画表の例を示す。収支計画表の最左列は、当該銀行における取扱金融商品の科目を示す。取扱金融商品は、例えば、企業向けの事業融資、個人向けの住宅ローン貸付、カードローン貸付、および銀行が自ら保有している有価証券配当について、それぞれＡからＣまでの３区分の科目を含む。【０１５０】
実際の銀行においては、数千程度の金融商品の科目を保有することがある。表の各行は、最上行に示される期間における、金融商品の科目それぞれについての取扱高を示す。例えば、現時点が２０１４年２Ｑ終了時とすれば、‘１４／１Ｑおよび‘１４／２Ｑの取扱高は実績値を示すが、‘１４／３Ｑ以降の取扱高は現時点（２０１４年２Ｑ終了時点）における計画値を示す。図８の例では四半期単位で計６期間についての計画を示しているが、実際の銀行においては、例えば、月単位又は週単位で５年程度の計画を立てることもあり、このとき、収支計画表は数百程度の期間を含む。【０１５１】
各科目における取扱高と収益との関係は、例えば、主に銀行の資金調達金利、すなわちスポットレートに依存する。スポットレートと収益の関係については、科目毎にある程度の確度をもったモデル化が可能であるのに対して、将来のスポットレートは不確定であり確率的な予測しかできない。【０１５２】
そこで、銀行の経営においては、モンテカルロ・シミュレーションによって収益の期待値やばらつきをシミュレーションし、事業計画の妥当性を判断する。具体的には、例えば、乱数を用いて将来のスポットレートの実現パスを生成し、各科目の計画取扱高をもとに収益値を計算し、全科目の収益を合計することで銀行全体の収益値を計算する、という処理を多数回行う。当該多数回の処理によって、与えられた事業計画の下での銀行全体の収益の期待値やばらつきを推定する。【０１５３】
ところで銀行において、事業計画の妥当性を判断するにあたっては、将来の各科目の取扱高の計画値が完全に達成された場合の銀行全体の期待収益やばらつきを求めるだけでは不足である。各科目の将来の実際の取扱高は、計画値からずれることが当然に想定されるためである。このとき、計画値と将来の実績がずれたときの影響について、事業計画立案時に予め考慮しておく必要がある。【０１５４】
具体的には、例えば、数千程度の科目それぞれについての数百程度の未来の各期間における取扱高の計画値、すなわち、合計十万個程度の計画値それぞれについて、将来の実績値が計画値からずれたときの銀行全体の収益に対する影響を計算する必要がある。もし実績値が計画値からほんの少しだけずれるだけで、銀行全体の収益に大きな影響を及ぼすような計画となっていることが分かれば、適切な引当金を積んでおく、又は、事業計画そのものを立て直す、といった対策が必要になるであろう。【０１５５】
以下、実施例１のモンテカルロ・シミュレーション・システムによるシミュレーション手法を銀行の収益予測シミュレーションに適用する例を説明する。【０１５６】
本実施例におけるパラメータベクトルｚは、例えば、全ての科目それぞれについて未来の各期間における取扱高の計画値である。例えば、科目数は数千程度、期間は数百期間である場合、パラメータベクトルｚは十万次元程度のベクトルである。システムの挙動の時系列ｘ（ｔ）は、例えば、将来の各期間における各科目の収益値の集合（例えば、数千次元のベクトルになりえる）と、金利のスポットレートと、を合わせた多次元のベクトルである。【０１５７】
ｔ番目の期間における各科目の収益値は、例えば、スポットレートと、その時刻における各科目の取扱高の計画値であるパラメータｚと、に依存する。将来のスポットレートは、例えば、Ｈｕｌｌ−Ｗｈｉｔｅモデル、ＣＩＲモデル、又はＢＤＴモデルなどといった確率的なモデルで表現される。これらのモデルは、時刻ｔ＋１でのスポットレートを、時刻ｔでのスポットレートと、（時刻ｔで新たに生成する）乱数と、に基づいて定める。【０１５８】
したがって、本実施例におけるシステムの挙動ｘ（ｔ）は、全体として、数１７の形で表される。すなわち、時刻ｔ＋１でのシステムの状態ｘ（ｔ＋１）は、時刻ｔと、過去（時刻ｔ以前）のシステムの挙動と、時刻ｔで新たに生成される（時刻毎に異なる）乱数ｗｊと、パラメータｚと、からで定まる。【０１５９】
システムの挙動の評価関数ａ（ｘ）は、例えば、全科目の全期間の収益の合計値である。このとき、システム評価関数Ａ（ｘ）、すなわちａ（ｘ）の期待値Ｅ［ａ（ｘ）］は、与えられた計画のもとでの銀行全体の将来の収益の期待値を示す。【０１６０】
本実施例のモンテカルロ・シミュレーション・システムは、上述した前提条件のもとで、図５に示した処理を行うことで、例えば、数千程度ある全ての科目それぞれについて数百程度ある未来の各期間における取扱高の計画値、すなわち、合計十万個程度の計画値それぞれについて、銀行全体の収益の期待値に対する感度、を同時に推定することが可能になる。【０１６１】
銀行は、推定した感度に基づいて、例えば、計画値と実績値のずれの影響を予め考慮して適切な引当金を積んでおく、又は事業計画そのものを立て直す、といった対策を立てることが可能となる。【０１６２】
以上のように、本実施例のモンテカルロ・シミュレーション・システムは、金融分野、とくに銀行の収益予測シミュレーションシステムにおいて計画値の収益に対する感度を推定することが可能となるため、銀行は事業計画立案のプロセスを効率化することができる。
【実施例４】
【０１６３】
本実施例のモンテカルロ・シミュレーション・システムは、実施例１の数９における境界残余Ｒｉを算出する。前述のように、境界残余Ｒｉは、乱数ベクトルｗの分布範囲Ωの境界面（表面）∂Ωで（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔ＝０（ベクトル）であれば、Ｒｉ＝０となる。現実的にモンテカルロ・シミュレーションの考察対象となるシステムでは、この条件が成立して、Ｒｉ＝０となる場合が多いと考えられる。【０１６４】
しかしながら、例えば、乱数ベクトルｗが一様分布乱数である場合などに、この条件が成立せず、具体的に境界残余Ｒｉを算出する必要にせまられる場合がある。境界残余Ｒｉの算出には、数９の１行目に基づく方法と、数９の３行目に基づく方法と、の２通りの方法が考えられる。【０１６５】
図９は、乱数ベクトルｗが、区間［ａ，ｂ］上の一様分布乱数ベクトルｗ１と、区間［ｃ，ｄ］上の一様分布乱数ｗ２と、の２つの要素をもつ２次元のベクトルである場合の境界残余Ｒｉを算出する概念の例を示す。図中に矩形で示された範囲がΩである。確率密度関数ｇ（ｗ）は、以下の数１９で表される。【０１６６】
【数１９】
2016030938.tif 000021
【０１６７】
なお、本実施例における記憶部３００は、サンプルパス固定微分記憶部３３０（不図示）をさらに含む。サンプルパス固定微分記憶部３３０は、サンプルパス固定微分計算部２３１が算出したＤｉ，ｊを保持する。また、本実施例における後処理計算部４００は、プログラムである境界残余算出部４３０（不図示）をさらに含む。境界残余算出部４３０は、乱数ベクトルｗと、パラメータベクトルｚと、サンプルパス固定微分記憶部３３０から受信するＤｉ，ｊと、を入力として境界残余Ｒｉを算出し、算出した境界残余Ｒｉを期待値計算部４１０に送信する。【０１６８】
以下、乱数ベクトルｗが従う確率密度関数ｇが数１９で表される例を用いて、境界残余Ｒｉの２通りの算出方法を説明する。境界残余算出部４３０が、境界残余Ｒｉを算出する第１の方法は、数９の１行目に基づくものである。図９の条件において、数９の１行目は、以下の数２０のように展開される。【０１６９】
【数２０】
2016030938.tif 000022
【０１７０】
したがって、境界残余算出部４３０は、ｗ１をΩの境界上の値であるｂに固定したうえでｗ２を動かしてそれぞれ（∂ｗ２／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを計算した積分値から、ｗ１を境界値ａに固定したうえでｗ２を動かしてそれぞれ（∂ｗ２／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを計算した積分値を減算した値を算出する。また、境界残余算出部４３０は、ｗ２を分布の境界上の値であるｄに固定したうえでｗ１を動かしてそれぞれ（∂ｗ１／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを計算した積分値から、ｗ２を境界値ｃに固定したうえでｗ１を動かしてそれぞれ（∂ｗ１／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを計算した積分値を減算したものと、を加算した値を算出する。【０１７１】
境界残余算出部４３０は、算出した２つの値を加算することで、境界残余Ｒｉを算出する。なお、境界残余算出部４３０は、例えば、既知の数値積分の手法を用いて、上述した積分値それぞれを算出すればよい。【０１７２】
境界残余算出部４３０が、境界残余Ｒｉを算出する第２の方法は、数９の３行目に基づくものである。数９の３行目は、乱数ベクトルの分布範囲Ω全域における積分であるから、境界残余算出部４３０は、以下の数２１より、モンテカルロ・シミュレーションを用いた期待値として境界残余Ｒｉを算出することができる。【０１７３】
【数２１】
2016030938.tif 000023
【０１７４】
すなわち、乱数発生部２１０が、与えられた確率密度関数ｇ（ｗ）に従う乱数を多数発生させる。境界残余算出部４３０は、発生した乱数ベクトルそれぞれに対して、数２１の２行目の角括弧内の値を計算し、計算した値の平均値をＲｉの推定値とすればよい。【０１７５】
なお、境界残余算出部４３０は、数値的な微分により、数２１中のｄｉｖを算出すればよい。また、ユーザ等によって関数としてｄｉｖが与えられた場合、境界残余算出部４３０は、当該関数にパラメータベクトルｚと乱数ベクトルｗとを代入することにより、ｄｉｖを算出してもよい。なお、境界残余算出部４３０は、確率密度関数ｇ（ｗ）が、２次元一様分布以外である場合も同様に境界残余Ｒｉを算出することができる。感度計算部４２０は、上述の方法で算出された境界残余Ｒｉを、算出した平均値から引いた差を算出する。当該算出した差が当該感度の推定値である。【０１７６】
以上のように、本実施例のモンテカルロ・シミュレーション・システムは、境界残余Ｒｉがゼロでな場合に、その値を算出することが可能となる。また、モンテカルロ・シミュレーション・システムは、境界残余Ｒｉを算出することにより、システムの挙動の実現値ｘを算出するモデルに適用する乱数ベクトルがどのような分布であっても、即ちｘを算出する全てのモデルに対して、本実施例の感度計算手法を適用することができる。【０１７７】
なお、本発明は上記した実施例に限定されるものではなく、様々な変形例が含まれる。例えば、上記した実施例は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、ある実施例の構成の一部を他の実施例の構成に置き換えることも可能であり、また、ある実施例の構成に他の実施例の構成を加えることも可能である。また、各実施例の構成の一部について、他の構成の追加・削除・置換をすることが可能である。【０１７８】
また、上記の各構成、機能、処理部、処理手段等は、それらの一部又は全部を、例えば集積回路で設計する等によりハードウェアで実現してもよい。また、上記の各構成、機能等は、プロセッサがそれぞれの機能を実現するプログラムを解釈し、実行することによりソフトウェアで実現してもよい。各機能を実現するプログラム、テーブル、ファイル等の情報は、メモリや、ハードディスク、ＳＳＤ（Ｓｏｌｉｄ  Ｓｔａｔｅ  Ｄｒｉｖｅ）等の記録装置、または、ＩＣカード、ＳＤカード、ＤＶＤ等の記録媒体に置くことができる。【０１７９】
また、制御線や情報線は説明上必要と考えられるものを示しており、製品上必ずしも全ての制御線や情報線を示しているとは限らない。実際には殆ど全ての構成が相互に接続されていると考えてもよい。
【図１】
2016030938.tif 000026
【図２】
2016030938.tif 000027
【図３】
2016030938.tif 000028
【図４】
2016030938.tif 000029
【図５】
2016030938.tif 000030
【図６】
2016030938.tif 000031
【図７】
2016030938.tif 000032
【図８】
2016030938.tif 000033
【図９】
2016030938.tif 000034
【手続補正書】
【提出日】平成29年1月24日(2017.1.24)
【手続補正１】
【補正対象書類名】明細書
【補正対象項目名】００８８
【補正方法】変更
【補正の内容】
【００８８】
サンプルパス固定微分計算部２３１は、パラメータと、乱数発生部２１０が生成した乱数ベクトルと、を入力として（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出する。サンプルパス固定２階微分計算部２３２は、パラメータと、乱数発生部２１０が生成した乱数ベクトルと、を入力として∂Ｄｉ，ｊ／∂ｗｊを算出する（Ｄｉ，ｊ（ｋ）＝（∂ｗ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔ）。サンプルパス固定２階微分計算部２３２は、サンプルパス固定微分計算部２３１が算出した（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔと、パラメータと、を入力として、∂Ｄｉ，ｊ／∂ｗｊを算出してもよい。また、サンプルパス固定２階微分計算部２３２が（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出してもよく、このとき、サンプルパス固定微分計算部２３１は（∂ｗｊ／∂ｚｉ）｜ｘ＝ｃｏｎｓｔを算出しなくてもよい。
【手続補正２】
【補正対象書類名】明細書
【補正対象項目名】０１０４
【補正方法】変更
【補正の内容】
【０１０４】
第３のステップは、乱数重要度計算部２３３が、乱数ベクトルｗ（ｋ）のｊ番目の要素ｗｊ（ｋ）に関する乱数重要度Ｈｉ，ｊ（ｋ）を計算するステップである。第３のステップは、数８における総和記号の中の項を計算することに相当する。具体的には、例えば、乱数重要度計算部２３３は、受信したＤｉ，ｊ（ｋ）と、受信したＤＤｉ，ｊ（ｋ）と、予め定められた確率密度関数ｇ（ｗ）から一意に定まる関数ＬＧに乱数ベクトルｗｊ（ｋ）を代入した値ＬＧj（ｗｊ（ｋ））＝（１／ｇ（ｗｊ（ｋ）））×（∂ｇ／∂ｗｊ（ｋ））と、を図中に示した関係式に代入して乱数重要度Ｈｉ，ｊ（ｋ）を計算すればよい。乱数重要度計算部２３３は、算出したＨｉ，ｊ（ｋ）を乱数重要度記憶部２５０に格納する。
【手続補正３】
【補正対象書類名】明細書
【補正対象項目名】０１４５
【補正方法】変更
【補正の内容】
【０１４５】
したがって、モンテカルロ・シミュレーション・システムは、図４で説明した処理を行うことにより、ＳＲＡＭセルの不良率の６０個の設計パラメータ全てに対する感度を算出することができる。また、モンテカルロ・シミュレーション・システムは、は、算出した感度から、勾配法やニュートン法といった公知の最適化手法を用いてＳＲＡＭセル不良率を最小にするような６０個の設計パラメータ値の組み合わせを、算出してもよい。
【国際調査報告】
2016030938.tif 000035
2016030938.tif 000036