JP6798487B2 - シミュレーション装置、シミュレーション方法及びプログラムを記憶する記録媒体 - Google Patents

Description

本発明は、実世界で起きている現象や仮説的な状況を数理的にモデル化してコンピュータ上で数値的に計算するシミュレーション技術に関する。

シミュレーションは、実世界で起きている現象や仮説的な状況を数理的にモデル化してコンピュータ上で数値的に計算する技術である。数理的にモデル化することにより、時間や空間を自在に設定して計算することが出来る。このようなシミュレーションにより、現実の結果を得ることが困難な状況（例えば、観測が困難な場所の状況）や、未来に起きうることを予測することが出来る。また、計算の条件を意図的に変化させることにより、現実では見ることが難しい状況での特性や振る舞いを調べることができる。そのようなシミュレーションの結果は、因果関係の理論的解明、設計、及び、計画などの指標として役立てることができる。

特に数理的なモデルに基づく予測では、実際の観測値によってキャリブレーション（校正）された後の予測結果のみならず、未観測及び未知の対象に対する予測能力（以下、汎化能力と呼ぶ）を向上させることが課題である。汎化能力はモデルの出力の汎化誤差を低減することで向上する。ここでの汎化誤差とは分散（ばらつき）とバイアス（偏り）によって表される。未知の対象を扱うにあたって、モデルがその対象の特徴を表す、例えば初期条件や境界条件などの条件、及び、可変なパラメータなどの入力に応じて、その対象の振る舞いを精度良く模擬出来ることが重要である。言い換えると、モデルには外部から設定可能な自由度が残されている。そして、その外部から設定可能な自由度の下で分散とバイアスを低減することが目的となる。しかし、特に物理法則や理論などに基づいた特定の領域（ドメイン）を扱うモデルは、研究グループや企業、機関、または分野、対象、地域ごとにモデル化や改良が行われている場合が多いため、直接の研究対象以外については、十分な精度が期待できない。典型的には、モデルが特定の対象に適合し過ぎている、すなわち自由度が小さい場合は、分散は小さいがバイアスが大きくなる傾向がある。一方で、モデルの汎用性が高い、すなわち自由度が大きい場合は、バイアスは小さいが分散が大きくなる傾向がある。したがって、汎化誤差を低減するためには、モデルの特性で左右される分散とバイアスのトレードオフを解決する必要がある。

そこで、このような汎化誤差を低減させる仕組みとして、アンサンブル（集団）を用いる概念が提案されている。その１つにモデルの変数をアンサンブルとして扱うデータ同化がある。データ同化は、現実で得られた観測データを数値シミュレーションに取り入れる手法として知られており、特に地球科学や海洋学、気象学の分野で発展してきた。データ同化は、シミュレーションによって計算される変数をアンサンブルとして扱うことにより、現実で得られた観測データにもっとも良く適合するシミュレーション結果を、アンサンブルの中から探しだす。同時に、データ同化は、モデルやシミュレーション条件も更新する。従って、データ同化を用いることによって、モデルの自由度が変数アンサンブルによって増大するためにバイアスは低減する。

また、複数のモデルの組み合わせをアンサンブルとして扱うモデルアンサンブルという手法が、例えば非特許文献１において提案されている。この方法では、同種の物理量に相当する出力を有する、少なくとも２つ以上のモデルが用意され、同一の対象に対してシミュレーションが行われ、それぞれの出力結果が平均される。前述したとおり、モデルは対象や地域ごとに開発されている場合が多いことから、汎化誤差を低減させるためには、より多様な対象や地域に適合する必要がある。そのために、同一の対象に対するシミュレーションであっても出力結果の異なる複数のモデル、言い換えると得手不得手が対象や地域に依存するような複数のモデルが用意される。これらの複数の出力結果の組み合わせを平均化することによって、その組み合わせに含まれる出力結果の数、すなわちアンサンブル数の増加によって分散を低減することができる。

Tao Li, Toshihiro Hasegawa, et al, "Uncertainties in predicting rice yield by current crop models under a wide range of climatic conditions", Global Change Biology(2014), DOI: 10.1111/gcb.12758

上述したデータ同化の技術及び非特許文献１に記載されている技術では、未観測及び未知の対象に対する予測能力（汎化能力）を向上させるために、アンサンブル（集団）が用いられる。しかし、その予測能力を左右する汎化誤差の内訳として分散とバイアスが存在し、それらがモデルの特性である自由度に依存したトレードオフの関係になっている。単純にアンサンブルを用いるだけでは、そのトレードオフの関係を打破できない点が課題である。例えば、上述したデータ同化を用いる方法では、誤差のバイアスは低減するが、モデルに内在しているモデル化誤差と観測データに含まれる観測誤差の存在が、モデルの自由度を大きくすることに相当するため、分散の低減が律速される。また、上述したモデルアンサンブルを用いる方法では、誤差の分散は低減するが、モデル１つ１つを見ると、特定の対象に適合している場合、すなわち自由度が小さい場合に相当するため、バイアスは低減しない。以上より、現在提案されている技術では、分散とバイアスのトレードオフ解決には至らない。

また、誤差低減を目的としてアンサンブルを用いる手法は、アンサンブル内の要素（アンサンブルメンバ）ごとに計算が必要となるため、アンサンブル数に比例して必要な計算リソースが増大することが課題である。従って、予測精度と計算リソースがトレードオフとなる。

本発明の目的の１つは、予測能力及び計算リソース効率の双方を向上させることができるシミュレーション装置などを提供することにある。

本発明の一態様に係るシミュレーション装置は、入力としての状態ベクトルを基に、少なくとも２つの異なるモデルについて、当該モデルに従ったシミュレーションによって、前記状態ベクトルの時間発展を算出するモデルアンサンブル手段と、前記状態ベクトルの時間発展と観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とを生成する事後分布生成手段と、前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する事後分布重み決定手段と、を備える。

本発明の一態様に係るシミュレーション方法は、入力としての状態ベクトルを基に、少なくとも２つの異なるモデルについて、当該モデルに従ったシミュレーションによって、前記状態ベクトルの時間発展を算出し、前記状態ベクトルの時間発展と観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とを生成し、前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する。

本発明の一態様に係る記録媒体は、コンピュータを、入力としての状態ベクトルを基に、少なくとも２つの異なるモデルについて、当該モデルに従ったシミュレーションによって、前記状態ベクトルの時間発展を算出するモデルアンサンブル手段と、前記状態ベクトルの時間発展と観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とを生成する事後分布生成手段と、前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する事後分布重み決定手段と、して動作させるプログラムを記憶する。本発明の一態様は、上述のプログラムによっても実現される。

本発明には、予測能力及び計算リソース効率の双方を向上させることができるという効果がある。

本発明の第１の実施の形態にかかるシミュレーション装置１００の構成を模式的に示すブロック図である。 本発明の第１の実施の形態にかかるシミュレーション装置１００における事後分布重み決定部３３の動作例を模式的に示す図である。 本発明の第１の実施の形態にかかるシミュレーション装置１００の動作を説明するフローチャートである。 本発明の第２の実施の形態にかかるシミュレーション装置２００の構成を模式的に示すブロック図である。 図５は、本発明の第３の実施形態のシミュレーション装置３００の構成の例を表すブロック図である。 図６は、シミュレーション装置１００、２００及び３００を実現することができる、コンピュータ１０００のハードウェア構成の一例を表す図である。

以下、図面を参照して本発明の実施の形態について説明する。各図面においては、同一要素には同一の符号が付されており、必要に応じて重複説明は省略される。

（第１の実施の形態）
本発明の第１の実施の形態としてのシミュレーション装置１００について説明する。シミュレーション装置１００は、物理法則に基づいた連続時間・空間の偏微分方程式を解くことによって時間発展を追うシミュレーションに適用可能である。そのような偏微分方程式には、例えば、運動を記述する運動方程式、流体を記述するナビエ−ストークス方程式、熱変化を記述する熱力学方程式、津波を記述する浅水波方程式などがある。また、シミュレーション装置１００は、有限要素法を用いるシミュレーションにも適用可能である。なお、本実施の形態では、シミュレーション対象となる系は、時間変化を追う状態ベクトルが、実際の観測データと何らかの関係式で結ばれる系、すなわち、シミュレーション結果と観測データとが比較可能な系であるものとする。

まず、シミュレーション装置１００の構成を図１に示す。図１において、シミュレーション装置１００は、データ入力部１０と、モデルアンサンブル部２０と、アンサンブル更新部３０と、データ出力部４０とを備える。データ入力部１０は、初期設定部１３と、取得部１４とを含む。モデルアンサンブル部２０は、ｍ個のシステムモデル２１（すなわち、第１のシステムモデル２１−１〜第ｍのシステムモデル２１−ｍ）を含む。アンサンブル更新部３０は、ｍ個の事後分布生成部３１（すなわち、第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）と、事後分布重み決定部３３とを含む。データ出力部４０は、出力部４２を含む。なお、ｍは、２以上の整数である。また、図１において、シミュレーション装置１００は、各機能ブロック間で入出力されるデータの格納領域として、例えば以下に具体的に示す記憶部を含む。具体的には、シミュレーション装置１００のデータ入力部１０は、シミュレーション条件が格納されるシミュレーション条件記憶部１１と、観測データが格納される観測データ記憶部１２をさらに含む。シミュレーション装置１００のモデルアンサンブル部２０は、事前分布が格納される、ｍ個の事前分布記憶部２２（すなわち、第１の事前分布記憶部２２−１〜第ｍの事前分布記憶部２２−ｍ）をさらに含む。シミュレーション装置１００のアンサンブル更新部３０は、事後分布が格納される、ｍ個の事後分布記憶部３２（すなわち、第１の事後分布記憶部３２−１〜第ｍの事後分布記憶部３２−ｍ）を更に含む。シミュレーション装置１００のアンサンブル更新部３０は、入力分布が格納される入力分布記憶部３４と、モデル尤度が格納されるモデル尤度記憶部３５とを更に含む。入力分布記憶部３４は、ｍ個の記憶部（ｍ個の入力分布記憶部３４（すなわち、第１の入力分布記憶部３４−１〜第ｍの入力分布記憶部３４−ｍ）とも表記）によって実現されていてもよい。シミュレーション装置１００のデータ出力部４０は、シミュレーションの結果である時系列状態ベクトルが格納される、時系列状態ベクトル記憶部４１を含む。

次に、シミュレーション装置１００の各機能ブロックの詳細について説明する。

まず、データ入力部１０について説明する。データ入力部１０の取得部１４は、シミュレーションで扱う状態ベクトルの初期状態およびパラメータを含むシミュレーション条件と、観測データとを取得する。取得部１４は、シミュレーション条件と、観測データとを、例えば制御装置（図示されない）及びサーバ装置（図示されない）などから取得すればよい。取得部１４は、取得したシミュレーション条件を、シミュレーション条件の格納領域である、シミュレーション条件記憶部１１に格納する。データ入力部１０は、さらに、取得した観測データを、観測データの格納領域である、観測データ記憶部１２に格納する。ここで、観測データとは、センサ等による観測値である。シミュレーション条件が含むパラメータは、例えば、対象とするシミュレーションの実施に必要な条件である、時間ステップ、終了までの規定時間、及び時間発展が求められる状態変数等である。初期設定部１３は、シミュレーション条件記憶部１１に格納された情報を元に、初期状態のアンサンブルとパラメータとを生成する。初期設定部１３は、生成した初期状態のアンサンブルとパラメータとを、ｍ個のシステムモデル２１のそれぞれに入力する。

次に、モデルアンサンブル部２０について説明する。モデルアンサンブル部２０は、ｍ個のシステムモデル２１（すなわち、第１のシステムモデル２１−１〜第ｍのシステムモデル２１−ｍ）を含む。本発明の各実施形態では、ｍ個のシステムモデル２１は、それぞれ、後述されるように写像によって表されるモデルに従って、入力された状態ベクトルを次のステップの時刻における状態ベクトルに変換する、例えば回路等の専用のハードウェアである。ｍ個のシステムモデル２１は、それぞれ、コンピュータ及びそのコンピュータを制御するプログラムによって実現されていてもよい。ｍ個のシステムモデル２１は、コンピュータ及びそのコンピュータを制御するプログラムと、専用のハードウェアとの組み合わせによって実現されていてもよい。モデルアンサンブル部２０は、さらに、ｍ個のシステムモデル２１にそれぞれ関連付けられているデータ格納領域であるｍ個の事前分布記憶部２２（第１の事前分布記憶部２２−１〜第ｍの事前分布記憶部２２−ｍ）を含む。シミュレーション装置１００では、現実の現象の時間発展を追う状態変数の組は状態ベクトルＸによって表される。ある時刻ｔにおける状態ベクトルＸ_ｔから１ステップ進んだ時刻ｔ＋１における状態ベクトルをＸ_ｔ＋１と表す。例えば、時刻ｔ−１における状態ベクトルＸ_ｔ−１が、第１のシステムモデル２１−１及び第２のシステムモデル２１−２によって、時刻ｔにおける状態ベクトルＸ_１，ｔ及び状態ベクトルＸ_２，ｔにそれぞれ更新されるとする。そして、第１のシステムモデル２１−１及び第２のシステムモデル２１−２による更新を表す写像をそれぞれｆ１、ｆ２とすると、第１のシステムモデル２１−１及び第２のシステムモデル２１−２による状態ベクトルＸ_ｔ−１の更新は、式（１）に示す関係式によって記述される。

式（１）は、複数のシステムモデルを用いることによって、時刻ｔ−１における状態が同一の状態Ｘ_ｔ−１であっても、その同一の状態Ｘ_ｔ−１が時刻ｔではそれぞれ異なる状態Ｘ_１，ｔ、及びＸ_２，ｔへ遷移することを表している。ここで、θ_１及びθ_２はそれぞれのモデルの計算に必要なパラメータベクトルを表し、ｖ_１ｔ及びｖ_２ｔはそれぞれのモデル内の不完全性の効果を数値的に表すために、状態ベクトルに作用を及ぼす確率的な駆動項として導入したシステムノイズを表す。なおここでの写像ｆ_１及びｆ_２は、対象とする現象に応じて、線形であってもよく、または非線形であっても良い。また式（１）から明らかなように、時刻ｔでの状態ベクトルＸ_ｔが、時刻ｔ−１の状態ベクトルＸ_ｔ−１によって陽に記述されている必要はない。すなわち、本実施の形態の第１のシステムモデル２１−１及び第２のシステムモデル２１−２は、時刻ｔ−１における状態ベクトルＸ_ｔ−１を入力として、時刻ｔにおける状態ベクトルＸ_１ｔ及びＸ_２ｔが計算できれば良い。以下の説明では、時刻ｔ−１から時刻ｔへ更新するステップについて示す。またここではシステムモデルの数が２（ｍ＝２）の例を示したが、システムモデルの数は、対象や適用可能なモデルの数に応じて適宜決められていれば良い。

ここで、モデルの変数をアンサンブルとして扱うアンサンブル近似について説明する。以降、時刻ｔの各システムモデル２１−ｊ（ｊは１からｍまでの整数）の状態ベクトルＸ_ｊｔは、システムモデル２１−ｊを表す写像ｆ_ｊの不完全性や入力されるパラメータθ_ｊの不完全性を反映して、確定値Ｘ_ｊｔではなく、確率分布ｐ（Ｘ_ｔ）として扱われる。その確率分布は、式（２）に示すアンサンブルによって、アンサンブル近似を使用して表すことができる。

式（２）に示す例では、アンサンブルはＮ個のアンサンブルメンバ（アンサンブルメンバＸ_ｊ，ｔ ^（ｉ）、以下、アンサンブルメンバｉとも表記）の集合である。各アンサンブルメンバは、状態ベクトルである。式（１）によって表される第１のシステムモデル２１−１及び第２のシステムモデル２１−２は、式（３）に示すアンサンブルメンバｉの時間発展をすべてのアンサンブルメンバｉについて求める。

これによって、それぞれのシステムモデルにて計算された時刻ｔにおける状態ベクトルＸ_１，ｔ、及びＸ_２，ｔの確率分布を計算することが出来る。なお、複数のシステムモデル、及び式（３）によって表されるアンサンブルの計算は、モデル及びアンサンブルメンバごとに独立している点が特徴である。従って、例えばシステムモデルの数がｍ個であり、アンサンブルに含まれるアンサンブルメンバの数がＮ個である場合、ｍ×Ｎ回の計算が繰り返されても良いし、並列計算が行われてもよい。このように、計算リソースに応じて柔軟に計算方法を変えることが出来る。なお以降では、第ｊ（ｊ＝１，．．．，ｍ）のシステムモデルにより出力された状態ベクトルＸ_ｊ，ｔの確率分布を、「ｐ_ｊ（Ｘ_ｔ）」と表し、「事前分布」と表記する。また、このｍ個の事前分布ｐ_ｊ（Ｘ_ｔ）は、ｍ個の事前分布記憶部２２（第１の事前分布記憶部２２−１〜第ｍの事前分布記憶部２２−ｍ）にそれぞれ保存される。アンサンブル近似が使用される場合、事前分布として、その事前分布の近似であるアンサンブルがｍ個の事前分布記憶部２２に保存される。

次に、アンサンブル更新部３０について説明する。アンサンブル更新部３０は、ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）を含む。ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）は、第１の事前分布記憶部２２−１〜第ｍの事前分布記憶部２２−ｍより、時刻ｔにおける事前分布ｐ_１（Ｘ_ｔ）〜ｐ_ｍ（Ｘ_ｔ）をそれぞれ読み出す。ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）のそれぞれは、さらに、観測データ記憶部１２に格納された観測データｙを読み出す。ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）のそれぞれに、第１の事前分布記憶部２２−１〜第ｍの事前分布記憶部２２−ｍより、時刻ｔにおける事前分布ｐ_１（Ｘ_ｔ）〜ｐ_ｍ（Ｘ_ｔ）が入力される。さらに、ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）のそれぞれに、データ入力部１０の観測データ記憶部１２に格納された観測データｙが入力される。ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）のそれぞれは、入力された事前分布と観測データとに基づいて、事後分布を算出する。例えば第１の事前分布記憶部２２−１及び第２の事前分布記憶部２２−２に格納される事前分布ｐ_１（Ｘ_ｔ）及びｐ_２（Ｘ_ｔ）と観測データｙに基づく観測データの分布ｐ（ｙ）とが入力されたときの事後分布を、ｐ_１（Ｘ_ｔ｜ｙ）及びｐ_２（Ｘ_ｔ｜ｙ）と表記する。事後分布ｐ_１（Ｘ_ｔ｜ｙ）及びｐ_２（Ｘ_ｔ｜ｙ）は、ベイズの定理より、それぞれ式（４）によって表される。

ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）によって生成されるｍ個の事後分布も同様の式によって表される。ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）は、生成した事後分布を、事後分布記憶部３２（第１の事後分布記憶部３２−１〜第ｍの事後分布記憶部３２−ｍ）にそれぞれ格納する。これらのｍ個の事後分布は、それぞれ事後分布記憶部３２（第１の事後分布記憶部３２−１〜第ｍの事後分布記憶部３２−ｍ）に保存される。アンサンブル近似が使用される場合、事後分布として、その事後分布の近似であるアンサンブルがｍ個の事後分布記憶部３２に保存される。さらに式（４）の右辺の分子の第１項ｐ_１（ｙ｜Ｘ_ｔ）及びｐ_２（ｙ｜Ｘ_ｔ）は、尤度と呼ばれ、状態ベクトルＸ_ｔの観測データｙへの当てはまり度合の指標である。この尤度は、観測データｙを分散σの正規（ガウス）分布とみなした場合、第ｊ（ｊ＝１，…，ｍ）のシステムモデルから出力されたアンサンブルメンバｉについて、式（５）に従って計算できる。

式（５）に示す尤度λをもとに、時刻ｋ＝１からＴまでの期間における、アンサンブル全体（Ｎ個のアンサンブルメンバ）についての対数尤度Ｌ（ｊ）は、式（６）によって近似的に表すことが出来る。

この値は、第ｊのシステムモデルから出力された時刻ｋ＝１からＴまでの事前分布ｐ_ｊ（Ｘ_１）〜ｐ_ｊ（Ｘ_Ｔ）の観測データへの当てはまり度合を表すことから、モデル尤度と呼ばれる。ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）のそれぞれは、例えば式（６）によって表される値を、モデル尤度として計算する。ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）のそれぞれは、モデル尤度記憶部３５に、計算したモデル尤度を格納する。第１の事前分布記憶部２２−１〜第ｍの事前分布記憶部２２−ｍに格納されたｍ種のシステムモデルから出力された事前分布に対して、この式（６）で表されるモデル尤度をそれぞれ計算した結果Ｌ（１）〜Ｌ（ｍ）が、モデル尤度記憶部３５に保存される。なお、モデル尤度を比較する際、式（６）によって表される、時刻ｋ＝１からＴまでの全ての時系列データの和をモデル尤度として算出するのではなく、式（６−２）に示す、ある時刻ｔごとの逐次的なモデル尤度を算出することもできる。

すなわち、ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）のそれぞれは、式（６−２）によって表される値を、モデル尤度として計算してもよい。この方法は、尤度の計算負荷が高い場合や、観測データの頻度や信頼性などに依存して有効となる場合がある。なお、尤度、及び、モデル尤度は、それぞれ、上述した通り式（５）、及び、式（６）または（６−２）で計算することが出来るが、これらの式はあくまでも例である。特に、観測データが正規（ガウス）分布でない場合など、実用上はこの限りではない。尤度、及び、モデル尤度を算出する方法として、その都度、適した算出方法を使うことができる。

事後分布重み決定部３３は、事後分布記憶部３２（第１の事後分布記憶部３２−１〜第ｍの事後分布記憶部３２−ｍ）に保存された、それぞれのシステムモデルで時間発展演算によって導出された事後分布を読み出す。事後分布重み決定部３３は、さらに、モデル尤度記憶部３５に保存されたそれぞれのシステムモデルに対するモデル尤度を読み出す。事後分布重み決定部３３には、それぞれのシステムモデルで時間発展演算によって導出された事後分布と、モデル尤度記憶部３５に保存されたそれぞれのシステムモデルに対するモデル尤度とが入力される。そして、事後分布重み決定部３３は、モデル尤度Ｌ（１）〜Ｌ（ｍ）を指標として各事後分布ｐ_１（Ｘ_ｔ｜ｙ）〜ｐ_ｍ（Ｘ_ｔ｜ｙ）に重み付けをした新しい事後分布ｐ’_１（Ｘ_ｔ｜ｙ）〜ｐ’_ｍ（Ｘ_ｔ｜ｙ）を生成する。事後分布重み決定部３３は、新しい事後分布ｐ’_１（Ｘ_ｔ｜ｙ）〜ｐ’_ｍ（Ｘ_ｔ｜ｙ）を入力分布記憶部３４（第１の入力分布記憶部３４−１〜第ｍの入力分布記憶部３４−ｍ）にそれぞれ格納する。また、事後分布重み決定部３３は、データ出力部４０へ時間ステップごとの新しい事後分布の重ね合せを出力し、時系列状態ベクトル記憶部４１にその新しい事後分布の重ね合せを格納する。新しい事後分布の重ね合わせについては後述する。

事後分布重み決定部３３において、モデル尤度と新しい事後分布とを生成する方法について、例を示して説明する。図２に動作例の概念図を示す。横軸は時間ステップ、縦軸は状態ベクトルＸの値をそれぞれ表す。以下、モデルアンサンブル部２０が含むシステムモデル２１が第１のシステムモデル２１−１及び第２のシステムモデル２１−２の２つである場合を例として示す。まず、時刻ｔ−１での状態ベクトルがそれぞれのシステムモデル（第１のシステムモデル２１−１及び第２のシステムモデル２１−２）に入力される。それぞれのシステムモデルに入力される状態ベクトルの分布である入力分布を、ｐ_１（Ｘ_ｔ−１）及びｐ_２（Ｘ_ｔ−１）と表す。第１のシステムモデル２１−１は、入力分布ｐ_１（Ｘ_ｔ−１）をもとに、第１のシステムモデル２１−１による変換を表す写像ｆ_１及びｆ_２によって、時間ステップが進められた時刻ｔでのそれぞれの事前分布ｐ_１（Ｘ_ｔ）を導出する。第２のシステムモデル２１−２は、入力分布ｐ_２（Ｘ_ｔ−１）をもとに、第２のシステムモデル２１−２による変換を表す写像ｆ１及びｆ２によって、時間ステップが進められた時刻ｔでのそれぞれの事前分布ｐ_２（Ｘ_ｔ）を導出する。事前分布ｐ_１（Ｘ_ｔ）と、時刻ｔでの観測データｙ_ｔとを第１の事後分布生成部３１−１に入力することにより、事後分布ｐ_１（Ｘ_ｔ｜ｙ_ｔ）と、モデル尤度Ｌ（１）とが得られる。事前分布ｐ_１（Ｘ_ｔ）と、時刻ｔでの観測データｙｔとを第２の事後分布生成部３１−２に入力することにより、事後分布ｐ_２（Ｘ_ｔ｜ｙ_ｔ）と、モデル尤度Ｌ（２）とが得られる。第１の事後分布生成部３１−１は、事前分布ｐ_１（Ｘ_ｔ）と時刻ｔでの観測データｙ_ｔとをもとに、事後分布ｐ_１（Ｘ_ｔ｜ｙ_ｔ）とモデル尤度Ｌ（１）とを導出する。第２の事後分布生成部３１−２は、事前分布ｐ_２（Ｘ_ｔ）と時刻ｔでの観測データｙ_ｔとをもとに、事後分布ｐ_２（Ｘ_ｔ｜ｙ_ｔ）とモデル尤度Ｌ（２）とを導出する。ここで図２に示すように、第２のシステムモデル２１−２を表す写像ｆ_２による事後分布ｐ_２（Ｘ_ｔ｜ｙ_ｔ）の方が、第１のシステムモデル２１−１を表す写像ｆ_１による事後分布ｐ_１（Ｘ_ｔ｜ｙ_ｔ）よりも観測データｙ_ｔに近いと仮定する。すなわち、モデル尤度の大小関係が、Ｌ（１）＜Ｌ（２）の場合を仮定する。すると、元の事後分布と事後分布重み決定部３３にて生成される新しい事後分布ｐ’_１（Ｘ_ｔ｜ｙ）、ｐ’_２（Ｘ_ｔ｜ｙ）とは、モデル尤度の関数である重み係数α_１、α_２によって、式（７）に従って算出される。

特に、ここでは、モデル尤度の大小関係Ｌ（１）＜Ｌ（２）によって、重み係数の大小関係もα_１＜α_２となっている。さらに全ての重み係数の和は、式（８）に示す関係を満たす。

式（８）は、生成された新しい事後分布の積分値が元の事後分布の積分値と同じか、または小さいことを意味する。式（２）に示すアンサンブル近似によれば、生成された新しい事後分布のアンサンブルメンバの、全てのシステムモデルについての総数（以下、合計とも表記）が元の事後分布のアンサンブルメンバの総数と同じか、または小さい。言い換えると、アンサンブルメンバの合計は、保存されるか、または減少が許される。

重み係数α_ｊに関する式（７）及び（８）について、より一般的に説明する。重み係数α_ｊは、ｍ種のモデル尤度Ｌ（ｊ）（ｊ＝１，…，ｍ）に基づいて決定される。重み係数は、例えば前述のとおり、モデル尤度が大きいほど重み係数が大きいように定義されていればよい。重み係数は、例えば、モデル尤度と重み係数との比例関係によって定義されていても良い。モデル尤度と重み係数との関係は、例えば、単調増加する線形な関数又は非線形な関数によって表されていても良い。また、重み係数αｊは、式（８）によって表される重み係数の総和が、例えばモデル尤度の総和が十分大きい場合は減少し、それ以外の場合は維持されるように、事後分布重み決定部３３によって決定されてもよい。重み係数の総和を減少させることは、アンサンブルメンバの総数（以下、アンサンブルメンバ数とも表記）を減少させることに相当する。重み係数の総和を維持することは、アンサンブルメンバの総数を維持することに相当する。このアンサンブルメンバ数は、繰り返し計算の回数、すなわち計算量と関係する。従って、システムモデルのモデル尤度が大きい場合、すなわちそのシステムモデルが観測データと比較して確からしい場合、そのシステムモデルに入力されるアンサンブルメンバの数を減少させることにより計算量を減少させることができる。システムモデルのモデル尤度が小さい場合、すなわち、そのシステムモデルが観測データと比較して明確に確からしいと判断できない場合は、そのシステムモデルに入力されるアンサンブルメンバの数を維持することにより計算の精度を保つことができる。なお、上述した重み係数の決定方法はあくまでも例示である。重み係数の決定方法は、実用上はこの限りではない。

次に、入力分布記憶部３４に格納される入力分布、及び、ｍ個のシステムモデル２１（第１のシステムモデル２１−１〜第ｍのシステムモデル２１−ｍ）へ入力される入力分布について、図２に示す例を使用して説明する。図２に示す例では、事後分布重み決定部３３によって生成される新しい事後分布ｐ’_１（Ｘ_ｔ｜ｙ）及びｐ’_２（Ｘ_ｔ｜ｙ）には、そのままでは状態ベクトルの値に対して大きな偏りが存在する。このように、事後分布重み決定部３３によって生成された新しい事後分布には、大きな偏りが存在する可能性がある。このまま時間ステップを進めた計算を継続した場合、偏った事後分布がそれぞれのシステムモデルに入力され続ける。そのため、計算の継続によって、状態ベクトルの分散が大きくなる。従って、汎化誤差の増大を招く。そこで本実施形態では、事後分布重み決定部３３が、新しい事後分布ｐ’_１（Ｘ_ｔ｜ｙ）及びｐ’_２（Ｘ_ｔ｜ｙ）の和集合であるｐ’_１（Ｘ_ｔ｜ｙ）∪ｐ’_２（Ｘ_ｔ｜ｙ）を、入力分布として入力分布記憶部３４に格納する。

さらに、例えば事後分布重み決定部３３が、第１のシステムモデル２１−１及び第２のシステムモデル２１−２が入力として読み出す分布を、以下に示すように入力分布からのランダムサンプリングによって生成する。事後分布重み決定部３３は、生成した分布を入力分布記憶部３４に格納する。第１のシステムモデル２１−１及び第２のシステムモデル２１−２は、決定した分布を入力分布記憶部３４から読み出す。第１のシステムモデル２１−１及び第２のシステムモデル２１−２は、読み出した分布を、次の時間ステップにおける計算に使用する。以下の説明では、第１のシステムモデル２１−１によって入力として読み出される分布を、「第１分布」表記する。同様に、第２のシステムモデル２１−２によって入力として読み出される分布を、「第２分布」表記する。

具体的には、事後分布重み決定部３３は、第１分布の積分値が、生成された新しい事後分布ｐ’_１（Ｘ_ｔ｜ｙ）の積分値と同じになるように、入力分布からのランダムサンプリングによって、第１分布を生成する。同様に、事後分布重み決定部３３は、第２分布の積分値が、事後分布重み決定部３３によって生成された新しい事後分布ｐ’_２（Ｘ_ｔ｜ｙ）の積分値と同じになるように、入力分布からのランダムサンプリングによって、第２分布を生成する。アンサンブル近似が使用されている場合、事後分布重み決定部３３は、アンサンブルメンバの総和がｐ’_１（Ｘ_ｔ｜ｙ）となるように、和集合ｐ’_１（Ｘ_ｔ｜ｙ）∪ｐ’_２（Ｘ_ｔ｜ｙ）からランダムサンプリングによってアンサンブルメンバを選択する。事後分布重み決定部３３は、選択したアンサンブルメンバを、第１分布として入力分布記憶部３４に格納する。第１のシステムモデル２１−１は、選択されたアンサンブルメンバ（すなわち第１分布）を、入力として、入力分布記憶部３４から読み出す。同様に、事後分布重み決定部３３は、アンサンブルメンバの総和がｐ’_２（Ｘ_ｔ｜ｙ）となるように、和集合ｐ’_１（Ｘ_ｔ｜ｙ）∪ｐ’_２（Ｘ_ｔ｜ｙ）からランダムサンプリングによってアンサンブルメンバを選択する。事後分布重み決定部３３は、選択したアンサンブルメンバを、第２分布として入力分布記憶部３４に格納する。第２のシステムモデル２１−２は、選択したアンサンブルメンバ（すなわち第２分布）を、入力として、入力分布記憶部３４から読み出す。

図２に示す例では、第２のシステムモデル２１−２のモデル尤度が第１のシステムモデル２１−１のモデル尤度より高い。その場合、第１のシステムモデル２１−１は、第１のシステムモデル２１−１へ入力されるアンサンブルメンバの数が、モデル尤度が高い第２のシステムモデル２１−２へ入力されるアンサンブルメンバの数より少なくなるように、ランダムサンプリングを行う。

事後分布重み決定部３３は、データ出力部４０に、式（９）に示す、時間ステップごとの新しい事後分布の重ね合せを出力する。

すなわち、事後分布重み決定部３３は、式（９）に示す、時間ステップごとの新しい事後分布の重ね合せを、時系列状態ベクトル記憶部４１に格納する。以上により、時系列状態ベクトル記憶部４１には、シミュレーション開始から終了までの時間ステップごとの式（９）によって表される新しい事後分布の重ね合せが、シミュレーション結果として保存される。出力部４２は、例えば端末装置（図示されない）からの要求に応じて、時系列状態ベクトル記憶部４１に格納されているシミュレーション結果を読み出し、読み出したシミュレーション結果を、その端末装置に出力する。

ここで、図１に示す第１の実施の形態のシミュレーション装置１００の動作について、図面を参照して詳細に説明する。

図３は、本実施形態のシミュレーション装置１００の動作の例を表すフローチャートである。図３に示す動作の前に、例えば、シミュレーション装置１００のユーザは、シミュレーション装置１００に実行させるシミュレーション（対象とするシミュレーション）の実施に必要な条件を、例えば端末装置等（図示されない）を使用して入力する。取得部１４は、入力された、対象とするシミュレーションの実施に必要な条件を受信し、その対象とするシミュレーションの実施に必要な条件をシミュレーション条件記憶部１１に格納する。シミュレーションの実施に必要な条件については後述する。そのユーザは、同様に、そのシミュレーションにおいて使用される観測データを、端末装置等を使用して、シミュレーション装置１００に入力する。取得部１４は、入力された観測データを、観測データ記憶部１２に格納する。シミュレーション装置１００は、例えば、そのユーザによる端末装置等を使用して行われた指示に従って、シミュレーションを開始する。

シミュレーション装置１００がシミュレーションを開始すると、まず、ｍ個のシステムモデル２１のそれぞれは、対象とするシミュレーションの実施に必要な条件を、例えば初期設定部１３を介して、シミュレーション条件記憶部１１から読み出す。シミュレーションの実施に必要な条件は、例えば、時間ステップ、終了までの規定時間、及び時間発展が求められる状態変数である。ｍ個のシステムモデル２１のそれぞれは、読み出した、対象とするシミュレーションの実施に必要な条件を、シミュレーションに使用する条件として設定する（ステップＳ１０１）。システムモデル２１は、設定した時間ステップ、終了までの規定時間、及び状態変数を、シミュレーション条件記憶部１１に格納する。次に、初期設定部１３が、シミュレーション条件記憶部１１に格納された情報を元に、初期状態のアンサンブルとパラメータとを生成する。なお、例えばあらかじめ選択された、ｍ個のシステムモデル２１のいずれかが、初期設定部１３として動作してもよい。初期設定部１３は、さらに、初期状態のアンサンブル及びパラメータを、ｍ個のシステムモデル２１（第１のシステムモデル２１−１〜第ｍのシステムモデル２１−ｍ）のそれぞれへ入力する（ステップＳ１０２）。そして、それぞれのシステムモデル２１が、次ステップのアンサンブル、すなわち事前分布を計算する（ステップＳ１０３）。このステップの時刻において、ｍ個の事後分布生成部３１（第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍ）は、観測データが得られているか否かを判定する（ステップＳ１０４）。観測データが無い場合（ステップＳ１０４においてＮＯ）、シミュレーション装置１００の動作は、ステップＳ１０３に戻る。ステップＳ１０３において、事前分布がそれぞれのシステムモデル２１に戻され、時間ステップは１増加する。言い換えると、ステップＳ１０３において、それぞれのシステムモデル２１は、例えば時間ステップを１増加させ、その時間ステップにおけるアンサンブルを計算する。観測データがある場合（ステップＳ１０４においてＹＥＳ）、シミュレーション装置１００の動作は、次ステップ（ステップＳ１０５）に進む。

そして、ｍ個の事後分布生成部３１は、それぞれ独立に、ｍ個のシステムモデル２１が計算した事前分布から事後分布を生成する（ステップＳ１０５）。上述のように、ｍ個の事後分布生成部３１は、第１の事後分布生成部３１−１〜第ｍの事後分布生成部３１−ｍである。ｍ個のシステムモデル２１は、第１のシステムモデル２１−１〜第ｍのシステムモデル２１−ｍである。その後、ｍ個の事後分布生成部３１は、さらに、モデル尤度を計算する（ステップＳ１０６）。事後分布重み決定部３３は、計算されたモデル尤度に基づき、事後分布重みを決定する（ステップＳ１０７）。そして、事後分布重み決定部３３は、次の時間ステップにおける入力である入力分布をサンプリングする（ステップＳ１０８）。事後分布重み決定部３３は、サンプリングによって生成した、次の時間ステップにおける入力である入力分布を、入力分布記憶部３４に格納する。また、事後分布重み決定部３３は、決定された重みを反映した事後分布の重ね合せにより生成したシミュレーション結果を、時系列状態ベクトル記憶部４１に格納する。データ出力部４０は、例えば、上述の端末装置（図示されない）を使用して入力されたユーザからの指示に従って、時系列状態ベクトル記憶部４１に格納されているシミュレーション結果をその端末装置に送信してもよい。

ここで、シミュレーションが開始してから規定の時間が経過している場合（ステップＳ１０９においてＹＥＳ）、例えば時間ステップが規定のステップ数に達している場合、図３に示す動作（すなわちシミュレーション）は終了する。シミュレーションが開始してから規定の時間が経過していない場合（ステップＳ１０９においてＮＯ）、すなわち、例えば時間ステップが規定のステップ数に達していない場合、シミュレーション装置１００の動作は、ステップＳ１０３に戻る。ステップＳ１０３において、入力分布記憶部３４から、ｍ個のシステムモデル２１（第１のシステムモデル２１−１〜第ｍのシステムモデル２１−ｍ）へ、それぞれのシステムモデル２１に応じた入力分布が入力される。例えば、ｍ個のシステムモデル２１は、それぞれ、計算した事後分布にもとに計算された入力分布を、入力分布記憶部３４から読み出す。そして、ｍ個のシステムモデル２１は、次ステップにおける計算を開始する。

以上で説明した本実施形態には、予測能力及び計算リソース効率の双方を向上させることができるという効果がある。

その理由は、事後分布重み決定部３３が、少なくとも２つの異なるモデルについて、事後分布とモデル尤度とに基づいて事後分布の重みを決定するからである。事後分布は、モデルアンサンブル部２０によって２つの異なるモデル毎に算出された、入力としての状態ベクトルの時間発展の、観測データに基づく事後分布である。モデル尤度は、観測データに基づくモデルの尤もらしさを表す値である。事後分布重み決定部３３は、さらに、決定した重みと事後分布とに基づいて、次の時間ステップにおける、モデルアンサンブル部２０への入力としての状態ベクトルを表す入力分布を算出する。本実施形態では、モデル毎に独立に、状態ベクトルの時間発展、事後分布、モデル尤度、及び、次の入力分布を算出することができる。従って、処理の並列化が可能である。さらに、モデル尤度が大きい場合、事後分布重み決定部３３で決定される重み係数の総和を減少させることによって、すなわちアンサンブルメンバの総数を減少させることによって、計算量を減らすことができる。このように、リソース効率を向上させることができる。また、本実施形態では、事後分布重み決定部３３が、モデル毎の、モデル尤度に基づく重みと事後分布とに基づいて、入力分布として表される状態ベクトルを算出する。例えば、モデル尤度が大きいほどモデルの重みを大きくなり、そして、モデル尤度が小さいほどモデルの重みが小さくなるように重みを決定することによって、予測能力を向上させることができる。

（第２の実施の形態）
次に、本発明の第２の実施の形態について図面を参照して説明する。本実施形態のシミュレーション装置２００は、第１の実施形態のシミュレーション装置１００を、作物の生育（成長）のシミュレーションに適用することによって実現される。作物生育（成長）についてシミュレーションを行うことによって、農業関係機関や農家をターゲットに、営農支援情報を提供することが出来る。なお、以下の本発明の第２の実施の形態の説明において参照される図面において、本発明の第１の実施の形態と同様の構成およびステップには同一の符号を付して、本実施の形態における詳細な説明を省略する。

まず、本発明の第２の実施の形態のシミュレーション装置２００の構成を、図４に示す。図４は、本実施形態のシミュレーション装置２００の構成を表すブロック図である。図４と第１の実施形態のシミュレーション装置１００の構成を示す図１とを比較すると、本実施形態のシミュレーション装置２００は、データ入力部１０及びモデルアンサンブル部２０の代わりに、データ入力部２１０及びモデルアンサンブル部２０Ａを含む。モデルアンサンブル部２０Ａとモデルアンサンブル部２０とを比較すると、モデルアンサンブル部２０Ａは、システムモデル２１の代わりに作物生育モデル２０１を含む。作物生育モデル２０１は、ｍ個の作物生育モデル（第１の作物生育モデル２０１−１〜第ｍの作物生育モデル２０１−ｍ）を含む。データ入力部２１０とデータ入力部１０とを比較すると、データ入力部２１０は、入力データ記憶部２１１を更に含む。

前述のように、本発明の第１の実施の形態のシミュレーション装置１００は、ｍ個のシステムモデル２１（第１のシステムモデル２１−１〜第ｍのシステムモデル２１−ｍ）を含む。図４に示すように、シミュレーション装置２００は、ｍ個のシステムモデル２１の代わりに、ｍ個の作物生育モデル２０１（第１の作物生育モデル２０１−１〜第ｍの作物生育モデル２０１−ｍ）を含む。本実施形態のｍ個の作物生育モデル２０１は、作物生育（成長）のシミュレーションに特化した、ｍ個のシステムモデル２１に相当する。言い換えると、本実施形態のｍ個の作物生育モデル２０１は、作物生育（成長）のシミュレーションに適用されたｍ個のシステムモデル２１である。また、データ入力部２１０のシミュレーション条件記憶部１１には、具体的なシミュレーション条件として、例えば、土壌の初期状態、地形、及び、作物を特徴付けるパラメータ等が格納されている。さらに、本実施の形態では、シミュレーション装置２００のデータ入力部２１０は、入力データ記憶部２１１を更に含む。入力データ記憶部２１１には、各種フィールドセンサやパブリックな情報として得られる、時系列の気象データや土壌水分量データなどが格納されている。なお、必要なシミュレーション条件、及び入力データは、個々の作物生育モデルの仕様に依存する。そのため、シミュレーション条件記憶部１１には、利用する作物生育モデルに応じて、その作物生育モデルに適したシミュレーション条件等が格納されていればよい。また観測データ記憶部１２には、具体的な観測データとして、作物の生育状態を表す、衛星や航空機によるリモートセンシングデータや、カメラ画像などが格納される。その他の部分の構成は第１の実施の形態のシミュレーション装置１００と同様であるので説明を省略する。

本実施形態では、作物の生育状態を示す観測データとして、植生指標として一般的に用いられているＮＤＶＩ（Ｎｏｒｍａｌｉｚｅｄ Ｄｉｆｆｅｒｅｎｃｅ Ｖｅｇｅｔａｔｉｏｎ Ｉｎｄｅｘ：正規化植生指標）が利用できる。この値は、可視赤バンド及び近赤外バンドの、２つのバンドの反射率から算出できる。さらに本実施の形態では、作物生育モデル２０１における状態ベクトルとして、葉面積指数（Ｌｅａｆ Ａｒｅａ Ｉｎｄｅｘ：ＬＡＩ）を用いることができる。ＬＡＩは、植生指数ＮＤＶＩと相関があることが知られている。このようなＬＡＩは、例えば、土壌の初期状態、地形及び作物のパラメータ、及び気象データ等のシミュレーション条件及びパラメータをもとに、作物モデル２０１によって算出される。ＬＡＩの算出に使用されるシミュレーション条件及びパラメータは、あらかじめシミュレーション条件記憶部１１又は入力データ記憶部２１１に格納され、作物生育モデル２０１によって読み出されればよい。ただし観測データ、及び状態ベクトルとして用いられる情報は、これらの量に限定されない。

観測データとしてのＮＤＶＩは、例えば、Ｔｅｒｒａ衛星またはＡｑｕａ衛星搭載のセンサＭＯＤＩＳ（Ｍｏｄｅｒａｔｅ Ｒｅｓｏｌｕｔｉｏｎ Ｉｍａｇｉｎｇ Ｓｐｅｃｔｒｏｒａｄｉｏｍｅｔｅｒ；Ｔｅｒｒａ／Ａｑｕａ ＭＯＤＩＳ）から得られるデータから算出できる。詳細には、Ｔｅｒｒａ／Ａｑｕａ ＭＯＤＩＳによる可視赤バンド（波長０．５８−０．８６μｍ）および近赤外バンド（波長０．７２５−１．１００μｍ）における太陽光に対する反射光の強度のデータが入手可能である。このデータは基本的に毎日取得可能であるが、地上での空間分解能は約２５０ｍと低い。またＬａｎｄｓａｔ衛星、Ｐｌｅａｄｅｓ衛星、及び、ＡＳＮＡＲＯ（Ａｄｖａｎｃｅｄ Ｓａｔｅｌｌｉｔｅ ｗｉｔｈ Ｎｅｗ ｓｙｓｔｅｍ Ａｒｃｈｉｔｅｃｔｕｒｅ ｆｏｒ Ｏｂｓｅｒｖａｔｉｏｎ）衛星などから得られるデータも利用できる。これらにより取得される波長域も、上述の波長域とほぼ同一である。ただし取得頻度と地上分解能は、Ｌａｎｄｓａｔ衛星の場合、８〜１６日間隔で約３０ｍであり、Ｐｌｅａｄｅｓ衛星、及び、ＡＳＮＡＲＯ衛星の場合、２〜３日間隔で約２ｍである。カメラ画像は、上述の可視赤バンド及び近赤外バンドにおける強度を表す画像であれば良い。ただし、本実施形態における観測データが観測される波長域は、必ずしもこれらのバンドに限定されるものではない。

複数の作物生育モデル２０１が従うモデルとして、例えば、ＤＳＳＡＴ（Ｄｅｃｉｓｉｏｎ Ｓｕｐｐｏｒｔ Ｓｙｓｔｅｍ ｆｏｒ Ａｇｒｏｔｅｃｈｎｏｌｏｇｙ Ｔｒａｎｓｆｅｒ）、ＡＰＳＩＭ（ｔｈｅ Ａｇｒｉｃｕｌｔｕｒａｌ Ｐｒｏｄｕｃｔｉｏｎ Ｓｙｓｔｅｍｓ ｓｉＭｕｌａｔｏｒ）、及び、ＷＯＦＯＳＴ（ＷＯｒｌｄ ＦＯｏｄ ＳＴｕｄｉｅｓ）などを用いることが出来る。そして、これらモデルの実環境への当てはまり度合が、モデル尤度という定量的かつ客観的な指標として、事後分布生成部３１によって算出され、モデル尤度記憶部３５に保存される。農業、特に露地栽培では、１つの作物が生育する回数は、１年に１回、多くても２〜３回程度である。また、土壌や作物品種などの、シミュレーションに必要な条件は、多様かつ多数であるため、取得が困難なことが多い。従って、このモデル尤度に関して蓄積された情報は、例えば、翌年度栽培開始の初期段階で、作物の生育状態を表す観測データが不十分な場合などに用いることが出来る。また、地域及び作物の組み合わせについてのモデル尤度が蓄積される。そのため、さらに他の地域及び作物へ展開するたびにデータを蓄積していくことによって、作物生育モデルごとの、地域及び作物に対する得手不得手を、定量的かつ客観的に判断することができるデータベースを作成できる。

さらに、複数の作物生育モデル２０１を用いたモデルアンサンブル部２０Ａについて、地域及び作物に対する得手不得手の克服以外の効果を述べる。近年、地球温暖化やゲリラ豪雨などの、気候及び気象の変動が顕在化してきている。特に露地栽培の農業は気象の影響を大きく受ける。従って、シミュレーションにより営農支援情報を提供する際、気象が作物生育に与える影響（特に悪影響）を精度よく示す必要がある。気象が作物生育に与える悪影響は、具体的には、日照りの継続による渇水ストレスや高温ストレス、ゲリラ豪雨による過剰水ストレスなどである。しかし、気象が作物生育に与える悪影響の算出精度は、作物生育モデルに大きく依存しているので、モデルの中で想定されていない気象状態は、算出されたシミュレーションの結果に反映されない。本実施形態の複数の作物生育モデル２０１は、異なるモデルに従って計算を行うことができる。例えば、現実の気象状態に応じた、渇水ストレスの計算精度が高いモデル、高温ストレスの計算精度が高いモデル、及び、過剰水ストレスの計算精度が高いモデルにそれぞれ従って影響の計算を行う、３つの異なる作物生育モデル２０１を用いることができる。従って、本実施形態では、複数の異なるシステムモデルを用いることによって、単独のモデルを用いる場合と比較して、予測能力を向上できる可能性がある。このように、現実世界の変化に適応できるように、複数のシステムモデルの組み合わせを設計できることも、本発明の特徴である。

また、モデルの尤度が低くても、すなわちモデルの重みが小さくても、そのモデルをモデルアンサンブルの構成要素として残しておくことによる効果について述べる。一般的には、尤度が高い、すなわち重みが大きいモデルさえ残しておけば、予測能力は維持できると考えられる。しかし、本実施形態において例示した作物生育モデルなどのように、環境や個体差の影響を受けるために不確実性が高い適用先では、異なる複数のシステムモデルを用いることによって、次の２つの効果が発揮される。第１の効果は、複数のシステムモデルを本実施形態のモデルアンサンブルとして扱うことによって、個々のシステムモデルに内在する不確実性が相殺されうる点である。第２の効果は、複数のシステムモデルに予測値が異なるシステムモデルが含まれていることによって、本実施形態のモデルアンサンブルにより分散とバイアスを低減しつつも、予測値にある程度の幅、すなわち分布を持たせることができる点である。従って、確定論的な予測のみならず、確率的なシナリオの予測も可能となる。

なお、本発明の第２の実施の形態のシミュレーション装置２００は、作物の育成のシミュレーションに適用された、第１の実施の形態のシミュレーション装置１００である。具体的には、作物生育モデル２０１は、作物の育成のシミュレーションに適用されたシステムモデル２１である。そして、第２の実施形態における観測データは、作物の生育状態を表すデータである。これらはあくまでも例である。本発明はこれらの例に制限されない。システムモデル２１が従うモデルは、適用する対象に応じて適切に選択されれば良い。本実施の形態は、シミュレーションモデル及び観測データの他の組み合わせにも適用することができる。シミュレーションモデル及び観測データの組み合わせは、例えば、水動態及び流体モデルと河川の水位センサデータ及び衛星データとの組み合わせ、生体及び人体モデルと医療データ又はヘルスケアデータとの組み合わせ、気象モデルと気象センサデータや衛星データとの組み合わせ等である。

以上の説明では、システムモデル２１が従うモデルは、物理法則、理論、機能等に基づいてすでにモデル式が与えられている演繹的なモデルである。そのようなモデルとして、例えば、第１の実施の形態で記載した、物理法則に基づく連続時間及び空間の偏微分方程式、第２の実施の形態で記載した作物生育モデル、水動態及び流体モデル、生体及び人体モデル、及び、気象モデルなどがある。しかし、システムモデル２１が従うモデルは、例えば、何らかのデータに基づき回帰的に生成されたモデル、または、機械学習などによって生成された帰納的なモデルであってもよい。例えば人間が関与する購買行動や、市場価格の推移などについて、物理法則や理論に基づくモデル式を立てることは容易ではない。しかし、帰納的なモデルを用いることによって、物理法則や理論に基づくモデル式を立てることが容易ではない対象を、モデルによって表すことができる。例えば、人間が関与する購買行動や、市場価格の推移などの対象を、購買モデルや、市場価格の推移モデルなどの帰納的なモデルによって表すことができる。すなわち、帰納的なモデルを用いることによって、本発明の各実施形態に係るシミュレーション装置を、物理法則や理論に基づくモデル式を立てることが容易ではない対象に適用することができる。さらに、演繹的なモデルと帰納的なモデルを組み合わせて使用しても良い。演繹的なモデル及び帰納的なモデルの組み合わせによって、いずれかのモデル単独では精度よく表すことができない対象を精度よく表せるようになることが期待できる。

以上で説明した本実施形態には、第１の実施形態の効果と同じ効果がある。その理由は、第１の実施形態の効果が生じる理由と同じである。

次に、本発明の第３の実施形態について、図面を参照して詳細に説明する。

図５は、本実施形態のシミュレーション装置３００の構成の例を表すブロック図である。図５を参照すると、本実施形態のシミュレーション装置３００は、モデルアンサンブル部２０と、事後分布生成部３１と、事後分布重み決定部３３と、を備える。モデルアンサンブル部２０は、入力としての状態ベクトルを基に、少なくとも２つの異なるモデルについて、当該モデルに従ったシミュレーションによって、前記状態ベクトルの時間発展を算出する。事後分布生成部３１は、前記状態ベクトルの時間発展と観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とを生成する。事後分布重み決定部３３は、前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する。

以上で説明した本実施形態には、第１の実施形態の効果と同じ効果がある。その理由は、第１の実施形態の効果が生じる理由と同じである。

シミュレーション装置１００、２００及び３００は、それぞれ、コンピュータ及びコンピュータを制御するプログラム、専用のハードウェア、又は、コンピュータ及びコンピュータを制御するプログラムと専用のハードウェアの組合せにより実現することができる。

図６は、シミュレーション装置１００、２００及び３００を実現することができる、コンピュータ１０００のハードウェア構成の一例を表す図である。図６を参照すると、コンピュータ１０００は、プロセッサ１００１と、メモリ１００２と、記憶装置１００３と、Ｉ／Ｏ（Ｉｎｐｕｔ／Ｏｕｔｐｕｔ）インタフェース１００４とを含む。また、コンピュータ１０００は、記録媒体１００５にアクセスすることができる。メモリ１００２と記憶装置１００３は、例えば、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ハードディスクなどの記憶装置である。記録媒体１００５は、例えば、ＲＡＭ、ハードディスクなどの記憶装置、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、可搬記録媒体である。記憶装置１００３が記録媒体１００５であってもよい。プロセッサ１００１は、メモリ１００２と、記憶装置１００３に対して、データやプログラムの読み出しと書き込みを行うことができる。プロセッサ１００１は、Ｉ／Ｏインタフェース１００４を介して、例えば、サーバ装置及び端末装置など（図示されない）にアクセスすることができる。プロセッサ１００１は、記録媒体１００５にアクセスすることができる。記録媒体１００５には、コンピュータ１０００を、シミュレーション装置１００、２００又は３００として動作させるプログラムが格納されている。

プロセッサ１００１は、記録媒体１００５に格納されている、コンピュータ１０００を、シミュレーション装置１００、２００又は３００として動作させるプログラムを、メモリ１００２にロードする。そして、プロセッサ１００１が、メモリ１００２にロードされたプログラムを実行することにより、コンピュータ１０００は、シミュレーション装置１００、２００又は３００として動作する。

第１のグループに含まれる部は、例えば、プログラムを記憶する記録媒体１００５からメモリ１００２に読み込まれた、各部の機能を実現することができる専用のプログラムと、そのプログラムを実行するプロセッサ１００１により実現することができる。第１のグループは、初期設定部１３、取得部１４、システムモデル２１、事後分布生成部３１、事後分布重み決定部３３、出力部４２、作物モデル２０１である。また、第２のグループに含まれる部は、コンピュータ１０００が含むメモリ１００２やハードディスク装置等の記憶装置１００３により実現することができる。第２のグループは、シミュレーション条件記憶部１１、観測データ記憶部１２、事前分布記憶部２２、事後分布記憶部３２、入力分布記憶部３４、モデル尤度記憶部３５、時系列状態ベクトル記憶部４１、入力データ記憶部２１１である。あるいは、第１のグループに含まれる部及び第２のグループに含まれる部の一部又は全部を、各部の機能を実現する専用の回路によって実現することもできる。

また、上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。

（付記１）
入力としての状態ベクトルを基に、少なくとも２つの異なるモデルについて、当該モデルに従ったシミュレーションによって、前記状態ベクトルの時間発展を算出するモデルアンサンブル手段と、
前記状態ベクトルの時間発展と観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とを生成する事後分布生成手段と、
前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する事後分布重み決定手段と、
を備えるシミュレーション装置。

（付記２）
前記事後分布生成手段は、前記事後分布が生成されるまでに得られている観測データに基づいて、または、新たに得られた観測データに基づいて逐次的に、モデル尤度を算出する
付記１に記載のシミュレーション装置。

（付記３）
前記事後分布重み決定手段は、前記モデル尤度に比例するように前記重みを決定する
付記１又は２に記載のシミュレーション装置。

（付記４）
前記モデルアンサンブル手段は、入力としての前記状態ベクトルを近似する確率分布である入力分布を基に、当該状態ベクトルの時間発展としての確率分布を算出し、
前記事後分布重み決定手段は、前記事後分布と前記重みとに基づいて、次の前記状態ベクトルとしての前記入力分布を算出する
付記１乃至３のいずれか１項に記載のシミュレーション装置。

（付記５）
前記事後分布重み決定手段は、前記重みを使用して前記事後分布を更新し、更新された前記事後分布の和集合を生成し、当該和集合からサンプリングされた分布として、前記入力分布を生成する
付記４に記載のシミュレーション装置。

（付記６）
前記モデルは、理論的に生成された演繹的モデル、又は、データに基づいて生成された帰納的モデルである
付記１乃至５のいずれか１項に記載のシミュレーション装置。

（付記７）
前記事後分布重み決定手段は、前記重みの総和が１以下になるよう、前記重みを決定する
付記１乃至６のいずれか１項に記載のシミュレーション装置。

（付記８）
前記シミュレーションにおける条件であるシミュレーション条件と、前記状態ベクトルの初期状態と、前記観測データとを取得する取得手段と、
前記状態ベクトルの時系列を出力する出力手段と、
を備える付記１乃至７のいずれか１項に記載のシミュレーション装置。

（付記９）
入力としての状態ベクトルを基に、少なくとも２つの異なるモデルについて、当該モデルに従ったシミュレーションによって、前記状態ベクトルの時間発展を算出し、
前記状態ベクトルの時間発展と観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とを生成し、
前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する、
シミュレーション方法。

（付記１０）
前記事後分布が生成されるまでに得られている観測データに基づいて、または、新たに得られた観測データに基づいて逐次的に、モデル尤度を算出する
付記９に記載のシミュレーション方法。

（付記１１）
前記モデル尤度に比例するように前記重みを決定する
付記０又は１０に記載のシミュレーション方法。

（付記１２）
入力としての前記状態ベクトルを近似する確率分布である入力分布を基に、当該状態ベクトルの時間発展としての確率分布を算出し、
前記事後分布と前記重みとに基づいて、次の前記状態ベクトルとしての前記入力分布を算出する
付記９乃至１１のいずれか１項に記載のシミュレーション方法。

（付記１３）
前記重みを使用して前記事後分布を更新し、更新された前記事後分布の和集合を生成し、当該和集合からサンプリングされた分布として、前記入力分布を生成する
付記１２に記載のシミュレーション方法。

（付記１４）
前記モデルは、理論的に生成された演繹的モデル、又は、データに基づいて生成された帰納的モデルである
付記９乃至１３のいずれか１項に記載のシミュレーション方法。

（付記１５）
前記重みの総和が１以下になるよう、前記重みを決定する
付記９乃至１４のいずれか１項に記載のシミュレーション方法。

（付記１６）
前記シミュレーションにおける条件であるシミュレーション条件と、前記状態ベクトルの初期状態と、前記観測データとを取得し、
前記状態ベクトルの時系列を出力する、
付記９乃至１５のいずれか１項に記載のシミュレーション方法。

（付記１７）
コンピュータを、
入力としての状態ベクトルを基に、少なくとも２つの異なるモデルについて、当該モデルに従ったシミュレーションによって、前記状態ベクトルの時間発展を算出するモデルアンサンブル手段と、
前記状態ベクトルの時間発展と観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とを生成する事後分布生成手段と、
前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する事後分布重み決定手段と、
して動作させるプログラム。

（付記１８）
コンピュータを、
前記事後分布が生成されるまでに得られている観測データに基づいて、または、新たに得られた観測データに基づいて逐次的に、モデル尤度を算出する前記事後分布生成手段と、
して動作させる付記１７に記載のプログラム。

（付記１９）
コンピュータを、
前記モデル尤度に比例するように前記重みを決定する前記事後分布重み決定手段と、
して動作させる付記１７又は１８に記載のプログラム。

（付記２０）
コンピュータを、
入力としての前記状態ベクトルを近似する確率分布である入力分布を基に、当該状態ベクトルの時間発展としての確率分布を算出する前記モデルアンサンブル手段と、
前記事後分布と前記重みとに基づいて、次の前記状態ベクトルとしての前記入力分布を算出する前記事後分布重み決定手段と、
して動作させる付記１７乃至１９のいずれか１項に記載のプログラム。

（付記２１）
コンピュータを、
前記重みを使用して前記事後分布を更新し、更新された前記事後分布の和集合を生成し、当該和集合からサンプリングされた分布として、前記入力分布を生成する前記事後分布重み決定手段と、
して動作させる付記２０に記載のプログラム。

（付記２２）
前記モデルは、理論的に生成された演繹的モデル、又は、データに基づいて生成された帰納的モデルである
付記１７乃至２１のいずれか１項に記載のプログラム。

（付記２３）
コンピュータを、
前記重みの総和が１以下になるよう、前記重みを決定する前記事後分布重み決定手段と、
して動作させる付記１７乃至２２のいずれか１項に記載のプログラム。

（付記２４）
コンピュータを、
前記シミュレーションにおける条件であるシミュレーション条件と、前記状態ベクトルの初期状態と、前記観測データとを取得する取得手段と、
前記状態ベクトルの時系列を出力する出力手段と、
して動作させる付記１７乃至２３のいずれか１項に記載のプログラム。

以上、実施形態を参照して本発明を説明したが、本発明は上記実施形態に限定されるものではない。本発明の構成や詳細には、本発明のスコープ内で当業者が理解し得る様々な変更をすることができる。

この出願は、２０１５年６月２日に出願された日本出願特願２０１５−１１１８５５を基礎とする優先権を主張し、その開示の全てをここに取り込む。

１０ データ入力部
１１ シミュレーション条件記憶部
１２ 観測データ記憶部
１３ 初期設定部
１４ 取得部
２０ モデルアンサンブル部
２１ システムモデル
２１−１ 第１のシステムモデル
２１−２ 第２のシステムモデル
２１−ｍ 第ｍのシステムモデル
２２ 事前分布記憶部
２２−１ 第１の事前分布記憶部
２２−２ 第２の事前分布記憶部
２２−ｍ 第ｍの事前分布記憶部
３０ アンサンブル更新部
３１ 事後分布生成部
３１−１ 第１の事後分布生成部
３１−２ 第２の事後分布生成部
３１−ｍ 第ｍの事後分布生成部
３２ 事後分布記憶部
３２−１ 第１の事後分布記憶部
３２−２ 第２の事後分布記憶部
３２−ｍ 第ｍの事後分布記憶部
３３ 事後分布重み決定部
３４ 入力分布記憶部
３５ モデル尤度記憶部
４０ データ出力部
４１ 時系列状態ベクトル記憶部
４２ 出力部
１００ シミュレーション装置
２００ シミュレーション装置
２０１ 作物生育モデル
２０１−１ 第１の作物生育モデル
２０１−１ 第１の作物生育モデル
２０１−ｍ 第ｍの作物生育モデル
２１０ データ入力部
２１１ 入力データ記憶部
３００ シミュレーション装置
１０００ コンピュータ
１００１ プロセッサ
１００２ メモリ
１００３ 記憶装置
１００４ Ｉ／Ｏインタフェース
１００５ 記録媒体

Claims (10)

1. 入力としての状態ベクトルを基に、少なくとも２つの異なる独立なモデルについて、当該モデルの各々に従った独立なシミュレーションによって、前記状態ベクトルの時間発展を、それぞれ異なる状態ベクトルとして独立に算出するモデルアンサンブル手段と、
   前記それぞれ異なる状態ベクトルの時間発展と、同一の観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とをそれぞれ独立に生成する事後分布生成手段と、
   前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する事後分布重み決定手段と、
   を備えるシミュレーション装置。
2. 前記事後分布生成手段は、前記事後分布が生成されるまでに得られている観測データに基づいて、または、新たに得られた観測データに基づいて、前記モデルごとに独立かつ逐次的に、モデル尤度を算出する
   請求項１に記載のシミュレーション装置。
3. 前記事後分布重み決定手段は、前記モデル尤度に比例するように、前記モデルごとに独立に前記重みを決定する
   請求項１又は２に記載のシミュレーション装置。
4. 前記モデルアンサンブル手段は、入力としての前記状態ベクトルを近似する確率分布である入力分布を基に、当該状態ベクトルの時間発展としての確率分布を前記モデルごとに独立に算出し、
   前記事後分布重み決定手段は、前記事後分布と前記重みとに基づいて、次の前記状態ベクトルとしての前記入力分布を前記モデルごとに独立に算出する
   請求項１乃至３のいずれか１項に記載のシミュレーション装置。
5. 前記事後分布重み決定手段は、前記重みを使用して前記モデルごとに独立に前記事後分布を更新し、更新されたそれぞれの異なる前記事後分布から和集合を生成し、当該和集合から、前記モデルの各々の重みであるアンサンブル数を設定してサンプリングされた分布として、前記入力分布を生成する
   請求項４に記載のシミュレーション装置。
6. 前記モデルは、理論的に生成された演繹的モデル、又は、データに基づいて生成された帰納的モデルである
   請求項１乃至５のいずれか１項に記載のシミュレーション装置。
7. 前記事後分布重み決定手段は、前記モデルの重みの総和が１にありアンサンブル数が保存されるように、又は、前記モデル尤度に基づいて前記モデルの重みのうちいずれかを小さくすることによって前記重みの総和が１未満になりアンサンブル数が減少するように、前記重みを決定する
   請求項１乃至６のいずれか１項に記載のシミュレーション装置。
8. 前記シミュレーションにおける条件であるシミュレーション条件と、前記状態ベクトルの初期状態と、前記観測データとを取得する取得手段と、
   前記状態ベクトルの時系列を出力する出力手段と、
   を備える請求項１乃至７のいずれか１項に記載のシミュレーション装置。
9. 入力としての状態ベクトルを基に、少なくとも２つの異なる独立なモデルについて、当該モデルの各々に従った独立なシミュレーションによって、前記状態ベクトルの時間発展を、それぞれ異なる状態ベクトルとして独立に算出し、
   前記それぞれ異なる状態ベクトルの時間発展と、同一の観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とをそれぞれ独立に生成し、
   前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する、
   シミュレーション方法。
10. コンピュータを、
    入力としての状態ベクトルを基に、少なくとも２つの異なる独立なモデルについて、当該モデルの各々に従った独立なシミュレーションによって、前記状態ベクトルの時間発展を、それぞれ異なる状態ベクトルとして独立に算出するモデルアンサンブル手段と、
    前記それぞれ異なる状態ベクトルの時間発展と、同一の観測データとに基づいて、前記状態ベクトルの前記時間発展の事後分布と、前記モデルの前記観測データに基づく尤もらしさであるモデル尤度とをそれぞれ独立に生成する事後分布生成手段と、
    前記事後分布と前記モデル尤度とに基づいて前記事後分布の重みを決定し、当該重みと前記事後分布とに基づいて次の前記状態ベクトルを算出する事後分布重み決定手段と、
    して動作させるプログラム。

Images