(57)【要約】  モデル推定部（１１）は対象観測データと参照観測データとに基づいて分析システム（２０）のモデル関数を推定するが、その際に、まず対象観測データと参照観測データに基づいて対象システムと参照システムとのモデルの変動量を推定する。そして推定された変動量に基づいて観測データを補正したあとに対象システムのモデル関数を推定する。パラメータ決定部（１２）はモデル関数の平均及び分散に基づいて獲得関数を計算し、該獲得関数により次の観測の際のパラメータ値を決める。データ取得部（１３）はパラメータ値を分析システム（２０）に設定し、それに対応する観測値を取得する。観測データをフィードバックしたループを繰り返すことで最適なパラメータ値を求める。これにより、マルチタスクベイズ最適化法によりパラメータ探索を行う際に、対象観測データと参照観測データとの相関が低い場合であってもパラメータ探索の効率を改善することができる。
2019244474.tif 000002
【特許請求の範囲】
【請求項１】
与えられたパラメータ値を条件として観測を行って観測値を出力する対象システムにおいて該観測値又は該観測値から求められる指標値が所定の状態となるようなパラメータ値を探索するパラメータ探索方法であって、
前記対象システムで観測された対象観測データと、該対象システムと関連する参照システムで観測された参照観測データと、に基づいて、前記対象システムのモデル関数の事後分布を推定するステップであり、パラメータ空間における前記対象観測データと前記参照観測データとの間の変動量を該両観測データに基づいて推定する変動量推定ステップを含み、その推定された変動量に相当する補正がなされたモデル関数の事後分布を推定するモデル推定ステップと、
前記モデル推定ステップで得られたモデル関数の事後分布に基づいて、次の観測の際の条件であるパラメータ値を決定するパラメータ決定ステップと、
前記決定されたパラメータ値を前記対象システムに設定し、該対象システムから対象観測データを取得するデータ取得ステップと、
を有し、前記モデル推定ステップ、前記パラメータ決定ステップ、及び前記データ取得ステップ、の処理を繰り返すことによって、より適切なパラメータ値を探索する、パラメータ探索方法。
【請求項２】
前記変動量は、パラメータ空間における対象観測データと参照観測データとの間のシフト量である、請求項１に記載のパラメータ探索方法。
【請求項３】
前記変動量推定ステップでは、対象観測データと、参照観測データと、で独立にそれぞれのシステムのモデル関数の事後分布を推定し、いずれか一方のモデル関数にパラメータ空間における異なる変動量を複数与えてそれぞれ両モデルの相関性を求め、最も高い相関性を示す変動量を推定値とする、請求項１に記載のパラメータ探索方法。
【請求項４】
前記変動量推定ステップでは、対象観測データと参照観測データとのいずれか一方のデータにパラメータ空間における異なる変動量を複数与えてそれぞれモデル関数の事後分布を推定し、それぞれの事後分布の平均値モデルと観測データとの誤差が最も小さくなる変動量を推定値とする、請求項１に記載のパラメータ探索方法。
【請求項５】
前記変動量推定ステップでは、モデル関数がガウス過程に従うとの前提で該モデル関数の事後分布を算出し、ガウス過程のハイパーパラメータを最尤推定した際に最大の尤度を示す変動量を推定値とする、請求項４に記載のパラメータ探索方法。
【請求項６】
与えられたパラメータ値を条件として観測を行って観測値を出力する対象システムにおいて該観測値又は該観測値から求められる指標値が所定の状態となるようなパラメータ値を探索するパラメータ探索装置であって、
前記対象システムで観測された対象観測データと、該対象システムと関連する参照システムで観測された参照観測データと、に基づいて、前記対象システムのモデル関数の事後分布を推定するものであり、パラメータ空間における前記対象観測データと前記参照観測データとの間の変動量を該両観測データに基づいて推定する変動量推定部を含み、その推定された変動量に相当する補正がなされたモデル関数の事後分布を推定するモデル推定部と、
前記モデル推定部により得られたモデル関数の事後分布に基づいて、次の観測の際の条件であるパラメータ値を決定するパラメータ決定部と、
前記決定されたパラメータ値を前記対象システムに設定し、該対象システムから対象観測データを取得するデータ取得部と、
前記モデル推定部、前記パラメータ決定部、及び前記データ取得部による処理を繰り返し実行するようにそれら各部を制御する制御部と、
を備える、パラメータ探索装置。
【請求項７】
前記変動量は、パラメータ空間における対象観測データと参照観測データとの間のシフト量である、請求項６に記載のパラメータ探索装置。
【請求項８】
前記変動量推定部は、対象観測データと、参照観測データと、で独立にそれぞれのシステムのモデル関数の事後分布を推定し、いずれか一方のモデル関数にパラメータ空間における異なる変動量を複数与えてそれぞれ両モデルの相関性を求め、最も高い相関性を示す変動量を推定値とする、請求項６に記載のパラメータ探索装置。
【請求項９】
前記変動量推定部は、対象観測データと参照観測データとのいずれか一方のデータにパラメータ空間における異なる変動量を複数与えてそれぞれモデル関数の事後分布を推定し、それぞれの事後分布の平均値モデルと観測データとの誤差が最も小さくなる変動量を推定値とする、請求項６に記載のパラメータ探索装置。
【請求項１０】
前記変動量推定部は、モデル関数がガウス過程に従うとの前提で該モデル関数の事後分布を算出し、ガウス過程のハイパーパラメータを最尤推定した際に最大の尤度を示す変動量を推定値とする、請求項９に記載のパラメータ探索装置。
【請求項１１】
与えられたパラメータ値を条件として観測を行って観測値を出力する対象システムにおいて該観測値又は該観測値から求められる指標値が所定の状態となるようなパラメータ値を探索するためのパラメータ探索用プログラムであって、コンピュータを、
前記対象システムで観測された対象観測データと、該対象システムと関連する参照システムで観測された参照観測データと、に基づいて、前記対象システムのモデル関数の事後分布を推定するものであり、パラメータ空間における前記対象観測データと前記参照観測データとの間の変動量を該両観測データに基づいて推定する変動量推定機能部を含み、その推定された変動量に相当する補正がなされたモデル関数の事後分布を推定するモデル推定機能部と、
前記モデル推定機能部により得られたモデル関数の事後分布に基づいて、次の観測の際の条件であるパラメータ値を決定するパラメータ決定機能部と、
前記決定されたパラメータ値を前記対象システムに設定し、該対象システムから対象観測データを取得するデータ取得機能部と、
前記モデル推定機能部、前記パラメータ決定機能部、及び前記データ取得機能部による処理を繰り返し実行するようにそれら各機能部を制御する制御機能部と、
して動作させる、パラメータ探索用プログラム。
【請求項１２】
前記変動量は、パラメータ空間における対象観測データと参照観測データとの間のシフト量である、請求項１１に記載のパラメータ探索用プログラム。
【請求項１３】
前記変動量推定機能部は、対象観測データと、参照観測データと、で独立にそれぞれのシステムのモデル関数の事後分布を推定し、いずれか一方のモデル関数にパラメータ空間における異なる変動量を複数与えてそれぞれ両モデルの相関性を求め、最も高い相関性を示す変動量を推定値とする、請求項１１に記載のパラメータ探索用プログラム。
【請求項１４】
前記変動量推定機能部は、対象観測データと参照観測データとのいずれか一方のデータにパラメータ空間における異なる変動量を複数与えてそれぞれモデル関数の事後分布を推定し、それぞれの事後分布の平均値モデルと観測データとの誤差が最も小さくなる変動量を推定値とする、請求項１１に記載のパラメータ探索用プログラム。
【請求項１５】
前記変動量推定機能部は、モデル関数がガウス過程に従うとの前提で該モデル関数の事後分布を算出し、ガウス過程のハイパーパラメータを最尤推定した際に最大の尤度を示す変動量を推定値とする、請求項１４に記載のパラメータ探索用プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分析条件であるパラメータ値を変更したときに観測値が変化する分析装置などのシステムにおいて、最適な又はそれに近いパラメータ値を探索するパラメータ探索方法、パラメータ探索装置、及びパラメータ探索用のコンピュータプログラムに関する。
【背景技術】
【０００２】
例えば液体クロマトグラフ質量分析装置（ＬＣ−ＭＳ）では、質量分析装置に含まれるイオン源、各種のイオン輸送光学系、或いは、イオン検出器などにそれぞれ印加する電圧を変化させると、同じ試料に対する信号強度値が変化する。そのため、それら印加電圧をそれぞれパラメータとして測定感度ができるだけ高くなるように、各パラメータの値を調整する作業が測定に先立って行われる。そうした調整作業の際には、通常、一つのパラメータについてその値を所定の範囲で且つ所定のステップ幅で変化させながらそれぞれ測定を実行し、その各測定において得られる信号強度値を比較することで、最も高い信号強度を与える値を最適なパラメータ値として選択する方法が採られることが多い。【０００３】
上述したように或るパラメータの最適値を探索する際には、できるだけ少ない測定回数で以て最適な値を探索できることが望ましい。何故なら、測定回数が多いと時間も手間も掛かるし、パラメータ探索のための測定に標準試料を用いる場合には、測定回数が多いと使用する標準試料の量も多くなるからである。【０００４】
観測値が最も大きくなる（場合によっては最も小さくなる）パラメータの値を探索するような作業は、ＬＣ−ＭＳに限らず、様々な分析装置、測定装置、計測装置、さらにはこうした装置を含む各種システムにおいて一般的に行われる。また、実際の測定や分析を行うものではなく、分子動力学などに基づくコンピュータシミュレーションを実行するシステムや、深層学習などの機械学習アルゴリズムをコンピュータ上で実行するシステムなど、様々な装置やシステムにおいても、同様のパラメータ探索が広く行われている。【０００５】
上記問題は、一般化すれば、「繰返し実験において、できるだけ良好な観測値が得られるパラメータの値を探索する」という問題である。ここでいう良好な観測値とは、一般的には最大の観測値であるが、実験や観測の目的等によっては、ＳＮ比が最大である観測値、最小の観測値、などである場合もある。【０００６】
一例として、図１０（ａ）に示すように、或るパラメータ（例えば電圧）に関する異なる四つの値についての観測が終了しており、図中に黒丸点で示す四つの観測値が得られているものとする。いま、この四つの観測値よりもさらに高い観測値を得ることが期待できる次のパラメータ値、を探索したい場合を考える。【０００７】
作業者が自らの判断により、次に設定すべきパラメータ値を選択する場合には、図１０（ｂ）に点線及び一点鎖線でそれぞれ示すように、どのような実験対象モデルを想定するのかによって次のパラメータ値の選択がばらつく（図中のＡ、Ｂ）ことが避けられない。【０００８】
これに対し、パラメータ探索の一つの手法として広く知られているものとしてベイズ最適化（Bayesian Optimization）がある（特許文献１など参照）。ベイズ最適化法によれば、取得済みの観測データに基づいて、確率的に良い結果が期待できる次のパラメータ値（実験条件）を決定することができる。図１０（ｃ）は、上記四つの観測値からベイズ最適化法により実験対象のモデルを推定し、次に観測すべきパラメータ値を探索した結果を示す図である。図１０（ｃ）において、実線で示すカーブはベイズ最適化法により推定されたモデル関数の事後分布の平均値、斜線で塗りつぶした範囲はそのモデル関数事後分布の不確かさ（或いは分散）、「Next」と記載された縦の太実線は自動的に選択された次のパラメータ値を示している。【０００９】
より詳しく述べると、ベイズ最適化法では実験対象のモデルがガウス過程に従うという仮定の下で、取得済みの観測データに基づいてモデル関数の事後分布の平均値と分散値とを計算する。そして、それら計算値を基に、真のものにより近づくようにモデルを改善することが期待される次の実験条件（パラメータ値）を決定する。そして、決定された実験条件の下で次の実験を行って新たな観測データを取得し、このデータを既に取得済みの観測データに加えてモデル関数の事後分布を再度推定する。こうした処理を繰り返すことで、より高い観測値が得られるパラメータ値を得ることができる。【００１０】
図１１はベイズ最適化を用いたパラメータ探索の一例を示す図であり、図１１（ａ）は初期の四つの観測値（図１０（ａ）と同じ）にさらに一つの観測値が加わった状態の図、図１１（ｂ）は初期の四つの観測値（図１０（ａ）と同じ）にさらに七つの観測値が加わった状態の図である。図１１において、黒丸点は観測値、実線のカーブは推定されたモデル関数事後分布の平均値、点線のカーブは真のモデル、斜線で塗りつぶした範囲はモデル関数事後分布の不確かさ、「Next」と記載された縦の太実線は自動的に選択された次のパラメータ値を示している。図１１（ｂ）に示すように、パラメータ値を変更しながら実験を繰り返していくと、観測値が最も高くなるパラメータ値の付近の範囲で優先的にパラメータ探索が行われるようになる。これにより、観測値が最も高くなる又はそれに近くなるようなパラメータ値を最終的に見つけることができる。【００１１】
しかしながら、上述した一般的なベイズ最適化法では、観測点の数が少ない探索初期におけるモデル関数の事後分布の推定精度が低く、そのために効率的な探索が難しいという課題があった。この課題に対し、近年、ベイズ最適化法を改良したマルチタスクベイズ最適化（Multi-Task Bayesian Optimization）法が提案されている（非特許文献１参照）。マルチタスクベイズ最適化法では、観測対象である実験（以下、「対象実験」という）とそれに関連する実験（以下、「参照実験」という）においてそれぞれデータを取得済みであることを前提とし、対象実験による対象観測データと参照実験による参照観測データとに相関があるという仮定の下で対象実験のモデルを推定する。非特許文献１に示されているようにマルチタスクベイズ最適化では、観測点数が少ない探索初期の段階でも参照観測ータに基づいて対象実験のモデルをかなり高い精度で推定することができるため、特に探索初期における探索の効率を改善することができる、という利点がある。
【先行技術文献】
【特許文献】
【００１２】
【特許文献１】特開２０１８−７３３６０号公報
【非特許文献】
【００１３】
【非特許文献１】スウェルスキー（K. Swersky）ほか２名、「 マルチタスク・ベイジアン・オプティマイゼイション（Multi-Task Bayesian Optimization）」、［online］、［２０１８年６月１日検索］、NIPS、2013、インターネット＜https://papers.nips.cc/paper/5086-multi-task-bayesian-optimization.pdf＞
【非特許文献２】スノーク（J. Snoek）ほか２名、「プラクティカル・ベイジアン・オプティマイゼイション・オブ・マシン・ラーニング・アルゴリズムズ（Practical Bayesian Optimization of Machine Learning Algorithms）」、NIPS, 2012
【発明の概要】
【発明が解決しようとする課題】
【００１４】
従来提案されているマルチタスクベイズ最適化法は、対照実験のモデルと参照実験のモデルとの相関が高いことを前提としている。これは例えば分析装置でいえば、対照実験を行う分析装置と参照実験を行う分析装置とが個体は異なるものの同じ機種の分析装置である、或いは機種は異なるものの同じ方式の分析装置である、といった場合を想定していることによる。しかしながら、実際には、相関が十分に高い対照観測データと参照観測データとを常に用意できるとは限らない。【００１５】
いま一例として、図１２中に太実線で示すガウス関数に従うモデルを原モデルとして、そのゲイン、標準偏差（σ）、及び平均値（μ）をそれぞれ変更した３種類のモデルを考える。これら３種類のモデルと原モデルとの相関係数はそれぞれ、１、０．９６、及び０．３６である。変更後のモデルの関数の形状はいずれも原モデルに類似しているものの、特に横軸方向、つまりはパラメータ値の変化の方向にシフトさせるようにモデルを変化させた場合（図１２（ｃ））、原モデルと変更後のモデルとの相関計数は大きく低下している。【００１６】
本発明者らの検討によれば、上述したように相関が低い参照観測データと対象観測データとを用いてマルチタスクベイズ最適化を行うと、対象実験のモデルを正しく推定することができず、パラメータの探索効率が低下してしまう場合があることが分かった。この例を図１３に示す。図１３は、対象実験の真のモデルが参照実験の真のモデルに対してシフトしている場合にマルチタスクベイズ最適化を試行した結果を示す図である。図中、太点線は参照実験の真のモデル、太一点鎖線は対象実験の真のモデル、黒四角点は参照実験の観測点、黒丸点は対象実験の観測点、太実線は対象実験モデルの事後分布の平均値、斜線で示す範囲は対象実験モデルの事後分布の不確かさを示している。【００１７】
図１３（ａ）は対象観測データが１点である場合、図１３（ｂ）は対象観測データが５点に増加した場合の結果である。適切なパラメータ探索ができていれば、対象観測データの点数が増加すれば対象実験モデルの事後分布の不確かさは減少する筈である。ところが、この例では、正しいモデル推定ができていないために、誤って参照実験における最適条件の付近を優先的に探索してしまっており、対象実験による観測点が増加しても最適なパラメータ値が決まらない。このように従来のマルチタスクベイズ最適化法では、対照観測データと参照観測データとの相関が低いと、その参照観測データをモデル推定に有効に活用できなくなり、パラメータ探索効率の改善度合いが低下したり、最悪の場合、対象観測データのみを用いてパラメータ探索をする場合よりも探索効率が悪くなる場合さえある。【００１８】
本発明は上記課題を解決するためになされたものであり、その目的とするところは、参照実験で得られたデータが対象実験で得られたデータに対してパラメータ空間において変化しており両者の相関が低い場合であっても、マルチタスクベイズ最適化法を用いて効率よく最適なパラメータを探索することができるパラメータ探索方法、装置、及びプログラムを提供することにある。
【課題を解決するための手段】
【００１９】
上記課題を解決するために成された本発明に係るパラメータ探索方法は、与えられたパラメータ値を条件として観測を行って観測値を出力する対象システムにおいて該観測値又は該観測値から求められる指標値が所定の状態となるようなパラメータ値を探索するパラメータ探索方法であって、
前記対象システムで観測された対象観測データと、該対象システムと関連する参照システムで観測された参照観測データと、に基づいて、前記対象システムのモデル関数の事後分布を推定するステップであり、パラメータ空間における前記対象観測データと前記参照観測データとの間の変動量を該両観測データに基づいて推定する変動量推定ステップを含み、その推定された変動量に相当する補正がなされたモデル関数の事後分布を推定するモデル推定ステップと、
前記モデル推定ステップで得られたモデル関数の事後分布に基づいて、次の観測の際の条件であるパラメータ値を決定するパラメータ決定ステップと、
前記決定されたパラメータ値を前記対象システムに設定し、該対象システムから対象観測データを取得するデータ取得ステップと、
を有し、前記モデル推定ステップ、前記パラメータ決定ステップ、及び前記データ取得ステップ、の処理を繰り返すことによって、より適切なパラメータ値を探索するものである。【００２０】
上記課題を解決するために成された本発明に係るパラメータ探索装置は、本発明に係るパラメータ探索方法を実施するための装置であり、与えられたパラメータ値を条件として観測を行って観測値を出力する対象システムにおいて該観測値又は該観測値から求められる指標値が所定の状態となるようなパラメータ値を探索するパラメータ探索装置であって、
前記対象システムで観測された対象観測データと、該対象システムと関連する参照システムで観測された参照観測データと、に基づいて、前記対象システムのモデル関数の事後分布を推定するものであり、パラメータ空間における前記対象観測データと前記参照観測データとの間の変動量を該両観測データに基づいて推定する変動量推定部を含み、その推定された変動量に相当する補正がなされたモデル関数の事後分布を推定するモデル推定部と、
前記モデル推定部により得られたモデル関数の事後分布に基づいて、次の観測の際の条件であるパラメータ値を決定するパラメータ決定部と、
前記決定されたパラメータ値を前記対象システムに設定し、該対象システムから対象観測データを取得するデータ取得部と、
前記モデル推定部、前記パラメータ決定部、及び前記データ取得部による処理を繰り返し実行するようにそれら各部を制御する制御部と、
を備えるものである。【００２１】
上記課題を解決するために成された本発明に係るパラメータ探索用プログラムは、本発明に係るパラメータ探索方法をコンピュータ上で実行するためのコンピュータプログラムであり、与えられたパラメータ値を条件として観測を行って観測値を出力する対象システムにおいて該観測値又は該観測値から求められる指標値が所定の状態となるようなパラメータ値を探索するためのパラメータ探索用プログラムであって、コンピュータを、
前記対象システムで観測された対象観測データと、該対象システムと関連する参照システムで観測された参照観測データと、に基づいて、前記対象システムのモデル関数の事後分布を推定するものであり、パラメータ空間における前記対象観測データと前記参照観測データとの間の変動量を該両観測データに基づいて推定する変動量推定機能部を含み、その推定された変動量に相当する補正がなされたモデル関数の事後分布を推定するモデル推定機能部と、
前記モデル推定機能部により得られたモデル関数の事後分布に基づいて、次の観測の際の条件であるパラメータ値を決定するパラメータ決定機能部と、
前記決定されたパラメータ値を前記対象システムに設定し、該対象システムから対象観測データを取得するデータ取得機能部と、
前記モデル推定機能部、前記パラメータ決定機能部、及び前記データ取得機能部による処理を繰り返し実行するようにそれら各機能部を制御する制御機能部と、
して動作させるものである。【００２２】
本発明においてパラメータ値を探索する対象のシステムは、各種の分析装置や測定装置、計測装置、観察装置などの、物理的又は化学的な手法による測定や観察を行う装置やシステムを含むほか、分子動力学などに基づくシミュレーションや、深層学習などの機械学習アルゴリズムに基づく演算処理など、コンピュータ上での処理のみによって観測値を取得可能なシステムも含む。【００２３】
また本発明に係るパラメータ探索用プログラムは、汎用パーソナルコンピュータ、より高性能なコンピュータ、又は、各種のシステムに組み込まれた専用のコンピュータ、などを動作させるプログラムである。こうしたプログラムは、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、メモリカード、ＵＳＢメモリ（ドングル）などの、非一時的な記録媒体に収録されてユーザに提供されるようにすることができる。或いは、インターネットなどの通信回線を介したデータ転送の形式で、ユーザに提供されるようにすることもできる。もちろん、ユーザがシステムを新規に購入する場合、該システムに含まれるコンピュータに本発明に係るログラムを予め組み込んでおくこともできる。【００２４】
本発明に係るパラメータ探索方法において、パラメータ探索の基本的なアルゴリズムはマルチタスクベイズ最適化である。即ち、対象システムのモデル関数の事後分布を推定するために、その対象システムで観測された対象観測データのほかに、参照システムで観測された参照観測データを用いる。但し、与えられた参照観測データをそのままモデル関数の事後分布の推定に用いるのではなく、対象観測データと参照観測データとの間のパラメータ空間での変動量を、その両観測データに基づいて推定する。このパラメータ空間はパラメータの種類の数と同じ次元を持つ。したがって、パラメータが２種類であれば、パラメータ空間は２次元空間である。【００２５】
対象観測データと参照観測データとの間の変動、つまりは対象システムのモデルと参照システムのモデルとの間の変動の態様としては、パラメータ空間における回転、シフト（移動）、伸縮、などが考えられる。これら全ての変動の態様についてモデル関数形状の変動を考慮しても構わないが、多くの場合、両データの間の相関の低下に大きく寄与するのはシフトであるため、変動の態様としてシフトのみを考慮すれば十分であることが多い。即ち、本発明の一態様として、前記変動量は、パラメータ空間における対象観測データと参照観測データとの間のシフト量であるものとするとよい。【００２６】
変動量が上記シフト量であり、例えばパラメータ空間が１次元である場合、変動量推定ステップでは、その１次元方向のシフトのみを想定してシフト量を推定すればよい。変動量の推定方法としては、次の二つの方法のいずれかを採用すればよい。【００２７】
即ち、一つの変動量推定方法は、対象観測データと、参照観測データと、で独立にそれぞれのシステムのモデル関数の事後分布を推定し、いずれか一方のモデル関数にパラメータ空間における異なる変動量を複数与えてそれぞれ両モデルの相関性を求め、最も高い相関性を示す変動量を推定値とする、という方法である。【００２８】
この方法では、計算時間が掛かる、つまりは計算コストの高いモデル関数の事後分布の計算を２回のみに抑えることができるという利点がある。但し、モデル関数の事後分布を正確に求めるには或る程度のデータ点数が必要であるため、データ点数が少ない場合にはこの方法を採用しにくい。【００２９】
他の一つの変動量推定方法は、対象観測データと参照観測データとのいずれか一方のデータにパラメータ空間における異なる変動量を複数与えてそれぞれモデル関数の事後分布を推定し、それぞれの事後分布の平均値モデルと観測データとの誤差が最も小さくなる変動量を推定値とする、という方法である。【００３０】
この二番目の方法では一つ目の方法に比べると計算時間が掛かるものの、対象観測データ及び参照観測データの両方とモデル関数の誤差を基に変動値を推定するため、データ点数が少ない場合でも変動量を精度良く推定ができるという利点がある。【００３１】
特に二番目の変動量推定方法では、モデル関数がガウス過程に従うとの前提で該モデル関数の事後分布を算出し、ガウス過程のハイパーパラメータを最尤推定した際に最大の尤度を示す変動量を推定値とする、という方法を採ることができる。【００３２】
これによれば、ガウス過程のハイパーパラメータ最尤推定によるモデル関数の推定と、パラメータ空間での変動量の指標、例えば獲得関数の算出とを同時に行うため、計算時間を短縮することができる。【００３３】
上記いずれの方法でも、妥当な変動量を推定することができる。変動量が推定されたならば、例えばその変動量を補正値として、対象観測データと参照観測データとに基づいて推定される対象システムのモデル関数の事後分布を補正する。したがって、このときに得られるモデル関数の事後分布は、対象観測データと参照観測データとの間の変動量を反映したものとなる。パラメータ決定ステップでは、上述したように得られた補正済みのモデル関数の事後分布に基づいて、次の観測の際の条件であるパラメータ値を決定する。【００３４】
データ取得ステップでは、決定されたパラメータ値を対象システムにおける次の条件として設定し、該対象システムで観測される観測値を情報として含む対象観測データを取得する。これにより、対象観測データは増加する。そして、上述したモデル推定ステップ、パラメータ決定ステップ、及びデータ取得ステップにおける処理を、所定の終了条件を満たすまで繰り返し実行させることで、最適であると推定されるパラメータ値を取得する。この繰り返しの際に、モデル推定ステップでは対象観測データと参照観測データとの間の変動量が補正されたモデルが利用されるので、対象観測データと参照観測データとの相関が低い場合であっても、少ない繰返し回数で以て適切なパラメータ値を探索することができる。【００３５】
なお、考慮すべき変動の要素が多いと、その要素毎に複数の変動量を設定していずれが適当かを推定する必要がある。それに対し上述したように変動の要素を相関の低下への寄与が最も大きいシフトに限定することで、変動量の補正についての十分な効果を確保しながら演算時間を短縮することができる。
【発明の効果】
【００３６】
本発明によれば、対象観測データと参照観測データとの変動量を補正することで両データの相関が実質的に大きくなり、モデル関数の事後分布の推定精度が向上する。これにより、特に対象観測データの数が少ない場合であっても、参照観測データを有効に利用して効率的にパラメータ探索を行うことができる。即ち、より少ない繰返し回数で以て最適又はそれに近い状態のパラメータ値を見つけることができる。
【図面の簡単な説明】
【００３７】
【図１】本発明の一実施例であるパラメータ探索装置の概略ブロック構成図。
【図２】図１に示したパラメータ探索装置においてパラメータ探索を行う際の処理の流れを示すフローチャート。
【図３】本実施例のパラメータ探索装置におけるモデル推定部のブロック構成図。
【図４】本実施例のパラメータ探索装置における第１の変動量推定方法の処理の流れを示すフローチャート。
【図５】本実施例のパラメータ探索装置における第２の変動量推定方法の処理の流れを示すフローチャート。
【図６】参照システムのモデルと対象システムのモデルとが互いにシフトした状態であるときの各モデルの一例を示す図。
【図７】図６に示したシステムで得られるデータに対する、シフト補正なしの場合（従来法）とシフト補正ありの場合（本発明法）とのパラメータ探索結果の比較を示す図。
【図８】参照システムのモデルと対象システムのモデルとのモデル関数が異なる例を示す図。
【図９】図８に示したシステムで得られるデータに対する、シフト補正なしの場合（従来法）とシフト補正ありの場合（本発明法）とのパラメータ探索結果の比較を示す図。
【図１０】パラメータ探索における問題設定とその解法の一例を示す図。
【図１１】ベイズ最適化法によるパラメータ探索の一例を示す図。
【図１２】ガウス関数を原モデルとしてゲイン、標準偏差、及び平均値を変更した３種類のモデルを示す図。
【図１３】マルチタスクベイズ最適化法によりパラメータ探索が適切でない場合の一例を示す図。
【発明を実施するための形態】
【００３８】
以下、本発明に係るパラメータ検索方法、並びに、該方法を実施するための装置及びコンピュータプログラムの一実施例を、添付の図面を参照しつつ説明する。
図１は、本実施例によるパラメータ探索装置の概略ブロック構成図である。
【００３９】
このパラメータ探索装置１０は、機能ブロックとして、モデル推定部１１、パラメータ決定部１２、データ取得部１３、及び繰返し制御部１４、を備える。本発明の一実施例における対象システムに相当する分析システム２０はパラメータ探索装置１０には含まれず、与えられたパラメータ値を条件として所定の分析、例えば質量分析やクロマトグラフ分析を実行し、その分析結果を観測値として得るものである。例えば分析システム２０が質量分析装置である場合、パラメータ値は例えば一つのイオン輸送光学素子に印加される電圧の値であり、観測値は標準試料に対する質量分析により得られた特定の質量電荷比における信号強度値である。もちろん、これは一例であり、分析システム２０の分析手法はこれに限らず、また、後述するように別のシステムに置き換えることもできる。
【００４０】
モデル推定部１１は、参照観測データと対象観測データとに基づいて、対象システムである分析システム２０のモデル関数の事後分布を推定する。対象観測データは分析システム２０で得られた観測値を含むデータであり、参照観測データは対象システムとは異なるものの、それに類似した参照システムで得られた観測値を含むデータである。
【００４１】
パラメータ決定部１２は、推定された対象システムのモデル関数事後分布に基づいて、該対象システムにおいて次に観測する（つまりは分析する）際の条件であるパラメータ値を決定する。データ取得部１３は、この決定されたパラメータ値を分析条件として分析システム２０へと入力し、該分析システム２０においてそのパラメータ値を条件とする所定の分析により得られた観測値を受け取る。
【００４２】
繰返し制御部１４は、データ取得部１３により取得された観測値及び設定されたパラメータ値を含む観測データを対象観測データに追加して、モデル推定部１１、パラメータ決定部１２、及びデータ取得部１３が同じ処理を繰返し実行するようにそれら各部の動作を制御する。パラメータ探索装置１０は、このフィードバックのルーチンを多数繰り返し実行することによって、最適な又はそれに近いパラメータ値を探索する。
【００４３】
なお、パラメータ探索装置１０を構成する各部の機能は、例えばデジタルシグナルプロセッサを含むハードウェア回路で構成することも不可能でないが、実際には、その機能の一部又は分析システム２０とのインターフェイス部分などを除く大部分を、コンピュータにインストールしたプログラムを該コンピュータにより実行することで実現するものとするとよい。こうしたプログラムは、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、メモリカード、ＵＳＢメモリ（ドングル）などの、非一時的な記録媒体に収録されてユーザに提供されるようにすることができる。或いは、インターネットなどの通信回線を介したデータ転送の形式で、ユーザに提供されるようにすることもできる。
【００４４】
以下、パラメータ探索を行う際の処理の流れと各部の動作を、図２に示したフローチャートを参照して詳しく説明する。
なお、以下の説明から明らかであるが、基本となる処理の手順は、非特許文献１等に開示されているマルチタスクベイズ最適化法のアルゴリズムである。
【００４５】
まず、モデル推定部１１は、システムのモデル関数ｆ（ｘ）が例えばガウス過程に従うという仮定の下で、該モデル関数の事後分布を推定する。非特許文献２によれば、ガウス過程に従う関数ｆ（ｘ）は次の(1)式で示されるように、関数平均値ｍ（ｘ）と共分散行列ｋ（ｘ，ｘ’）とによって表される関数であることが知られている。
ｆ（ｘ）〜ＧＰ（ｍ（ｘ），ｋ（ｘ，ｘ’））      …(1)
【００４６】
上記共分散行列は、「二つの変数位置ｘ、ｘ’において出力ｆ（ｘ）、ｆ（ｘ’）がどの程度同じような挙動をとるか」を表す関数である。共分散行列は、通常、二つの変数位置ｘ、ｘ’が近いほど大きな値をとる連続的な単調関数で定義され、例えばRBFカーネルやMatern5/2カーネルなどが用いられる。RFBカーネルの例を(2)式に示す。
ｋRBF（ｘ，ｘ’）＝θ0 exp（−||ｘ−ｘ’||2 ／２θ12）      …（2）
計算を簡略化する場合には、(1)式中の関数平均値モデルｍ（ｘ）を０とすることもできるが、通常は、関数平均値モデルｍ（ｘ）として観測データの平均値を用いるとよい。
【００４７】
いまｎ個の観測データ（ｘn，ｙn）が取得済みであるとすると、モデル推定部１１は、この観測データに基づいて、モデル関数の事後分布の平均値μ及び分散値σ2を求める（ステップＳ１）。平均値μ及び分散値σ2は、変量正規分布（Multivariate normal distribution）の公式より以下の(3)、(4)式で計算される。
μ（ｘ；｛ｘn，ｙn｝）＝Ｋ（Ｘ，ｘ）TＫ（Ｘ，Ｘ）-1（ｙ−ｍ（Ｘ））      …(3)
σ2（ｘ；｛ｘn，ｙn｝）＝ｋ（ｘ，ｘ）−Ｋ（Ｘ，ｘ）TＫ（Ｘ，Ｘ）-1Ｋ（Ｘ，ｘ）      …(4)
ここでは、参照観測データと対象観測データとを共に利用するので、マルチタスクにおけるタスク間の共分散行列Ｋt（ｔ、ｔ’）を導入し、共分散関数を次のように拡張する。
Ｋmulti＝（（ｘ，ｔ），（ｘ’，ｔ’））＝Ｋt（ｔ、ｔ’）［×］Ｋx（ｘ，ｘ’）      …(5)
ここで、［×］はクロネッカー積を表す。
【００４８】
上記(2)式などで表されるカーネル関数中のハイパーパラメータθは、モデル関数に即した妥当な値に設定する必要がある。そこで、ハイパーパラメータθは例えば、そのθが与えられたときのｙの推定値の尤もらしさを表す尤度関数Ｌ（θ）が最大となるように、最尤推定により算出される。関数f（ｘ）の対数尤度関数logＬ（θ）は、次の(6)式で表される。
logＬ（θ）＝（−１／２）ｙTＫ（Ｘ，Ｘ）-1ｙ−（１／２）log｜Ｋ（Ｘ，Ｘ）｜−（ｎ／２）log２π      …(6)
ここで、argmaxθＬ（θ）は一般的な最適化問題であり、例えば準ニュートン法のような勾配法で求めることができる。
【００４９】
パラメータ決定部１２は、モデル推定部１１からモデル関数ｆ（ｘ）の事後分布の平均関数値μ（ｘ）と分散値σ2（ｘ）とを受け取り、それらに基づいて評価のための獲得関数ａ（ｘ）を計算する。そして、この獲得関数ａ（ｘ）が最大となるパラメータ値ｘを見つけ、それを観測値の改善が期待できる次の分析条件のためのパラメータ値として出力する（ステップＳ２）。獲得関数ａ（ｘ）としては、ベイズ最適化法において一般に用いられるProbability of Improvement（ＰＩ）、Expected Improvement（ＥＩ）、Upper Confidence Bound（ＵＣＢ）などのいずれかを用いることができる。
【００５０】
一例として、Expected Improvementの式を(7)式に示す（非特許文献２参照）。
ＥＩ（ｘ）＝（μ（ｘ）−ｆ（ｘ+））Φ（Ｚ）＋σ（ｘ）φ（Ｚ）    if σ（ｘ）＞０，ＥＩ（ｘ）＝０    if σ（ｘ）＝０      …（7）
但し、Ｚ＝（μ（ｘ）−ｆ（ｘ+））／σ（ｘ）
ｆ（ｘ+）は観測データの最大値を表し、Φ（＊）は標準正規分布の累積分布関数、φ（＊）は標準正規分布を表す演算子である。Expected Improvementでは、関数ｆ（ｘ）の事後確率モデルの中の観測データの最大値を超える関数の期待値を算出しており、大域探索と局所探索とをバランスよく実施する関数として知られている。
【００５１】
データ取得部１３は、上述したようにパラメータ決定部１２から与えられたパラメータ値を分析システム２０における分析条件として設定し、それに対する観測値を取得する（ステップＳ３）。分析システム２０は、分析・計測システムや制御システムといった実験系のシステムであるほか、分子動力学のようなシミュレーションのシステム、或いは、深層学習のような機械学習アルゴリズムを実行するシステムなど、入力されたパラメータ値の下で何らかの値を観測するシステムであればよい。
【００５２】
次に、繰返し制御部１４は、予め決められた繰返しの終了条件が満たされるか否かを判定する（ステップＳ４）。繰返しの終了条件は、例えばルーチンの繰返し回数が予め設定された値に到達したとき、或いは、観測値が予め定められた目標に到達したとき（目標の範囲に入ったとき）、などとすればよい。そして、繰返し終了条件を未だ満たしていなければ、新たに得られた観測値とそのときのパラメータ値とを含む対象観測データを観測データに加えたうえでステップＳ４からステップＳ１へと戻り、上述したステップＳ１〜Ｓ３の処理を繰り返す。一方、繰返し終了条件を満たしていれば、最終的なパラメータ値を探索結果である最適値として出力し（ステップＳ５）、処理を終了する。
【００５３】
以上説明した内容は、一般的なマルチタスクベイズ最適化法で行われている処理と同じである。これに対し本実施例のパラメータ探索装置１０において特徴的であるのは、参照観測データと対象観測データとの相関が低い場合に、モデル推定部１１においてデータ補正が行われる点である。これについて、さらに図３〜図５を参照して詳しく説明する。図３は、本実施例のパラメータ探索装置におけるモデル推定部１１の詳細な機能ブロック構成図である。
【００５４】
参照観測データと対象観測データとの相関が低い場合、タスク間のデータ（ｘ，ｔ）、（ｘ’，ｔ’）にそれぞれ対する出力ｆ（ｘ，ｔ）と出力ｆ（ｘ’，ｔ’）とが類似した挙動をとることが少ない。そのため、対象システムのモデル関数を推定するうえで参照観測データの有用性が低くなる。具体的には、共分散行列Ｋt（ｔ，ｔ’）を相関が強い共分散行列としてしまうと誤ったモデルを推定することになるし、一方、共分散行列Ｋt（ｔ，ｔ’）を相関が低い共分散行列としてしまうとマルチタスクが独立なタスクとして扱われ、実質的なマルチタスクでなくなる。
【００５５】
そこで、本実施例のパラメータ探索装置において、モデル推定部１１は機能ブロックとして、変動量推定部１１１、補正処理部１１２、及び、補正済みモデル推定部１１３を備える。変動量推定部１１１は、パラメータ空間での回転、シフト、伸縮などの変動を計算上で与え、参照観測データと対象観測データの変動量を予測する。補正処理部１１２は予測される変動量だけ観測データを補正し、補正済みモデル推定部１１３は補正された観測データに基づいて上述したようにモデル推定を行う。
【００５６】
２種類のパラメータに対して観測値が得られるシステムに適用する場合を例に挙げて、具体的に説明する。図６は、横軸、縦軸にそれぞれ異なるパラメータ（Parameter 1、Parameter 2）をとり、それら２種類のパラメータの値の組合せに対する観測値をグレイスケール（実際にはカラースケール）で示した図である。これは、どのようなパラメータ値の下でどのような観測値が得られるかというシステムのモデルを表している。図６（ｂ）は参照システムのモデルを示しており、ここでは規則的に分布している白丸点で示した１６点の位置に対応するパラメータ値の組合せの下での観測値が参照観測データとして取得済みであるものとする。一方、図６（ａ）はパラメータ探索対象である対象システムのモデルを示しており、ランダムな３点の位置に対応するパラメータ値の組合せの下での観測値のみが取得済みであるものとする。この状態から、最適なパラメータ値の組合せを効率的に探索するのが、ここで与えられている命題である。
【００５７】
なお、図６（ａ）に示した対象システムのモデルは、図６（ｂ）に示した参照システムのモデルにおける各参照観測データをParameter 1の方向（図６で横軸方向）に−１〜−４の範囲でずらすことで、対象システムと参照システムとの間の相関が低いモデルを人為的に作成したものである。結果として、対象システムのモデルは、参照システムのモデルをParameter 1の方向に−２だけずらしたものとなっている。
【００５８】
上記のような参照観測データ及び対象観測データが与えられた場合、変動量推定部１１１は、次の二つの方法のいずれかでパラメータ空間での変動量を推定することができる。
【００５９】
［第１の変動量推定方法］
この第１の変動量推定方法は、対象観測データと参照観測データとで独立にそれぞれのシステムのモデル関数の事後分布を推定し、予測される変動量をカバーするパラメータ範囲で一方のモデル関数に変動を与えた上で両モデル間の相関値を計算し、最も高い相関を示す変動量を推定値とする方法である。この方法における処理の流れを図４に示す。
【００６０】
まず、与えられた対象観測データと参照観測データとについてそれぞれ独立に、ガウス過程のハイパーパラメータ最尤推定を行ってモデル関数の事後分布を推定する（ステップＳ１１）。次に、推定された参照システムのモデルに基づくデータを、Parameter 1及びParameter 2の方向にそれぞれ−５〜＋５の範囲で１ずつシフトさせたシフトデータから成る変形モデルを生成する（ステップＳ１２）。これにより、二つのパラメータのシフト量の組合せの総数だけ、参照システムの変形モデルが生成される。
【００６１】
そのあと、参照システムの各変形モデルと対象システムのモデルとの間の相関値を計算する（ステップＳ１３）。この相関値としては例えば、同じパラメータ値の組合せにおける観測値同士の相関値の総和などを用いればよい。そして、得られた複数の相関値の中で最も大きな相関値を示すシフト量の組合せを、変動量の推定値として採用する（ステップＳ１４）。こうして変動量が推定されたならば、補正処理部１１２はその変動量に基づいて対象観測データを補正し、補正済みモデル推定部１１３はその補正された対象観測データと参照観測データとに基づいて、上述したように対象システムのモデル関数を推定すればよい。
【００６２】
この方法では、変動量を１回求める際にモデル推定を２回のみ実行すればよい。通常、モデル推定のための計算量は多いため、モデル推定の回数が多いほど計算時間が長くなるが、モデル推定が２回で済むために計算量が少なくて済むという利点がある。
【００６３】
［第２の変動量推定方法］
第２の変動量推定方法は、対象観測データと参照観測データとのいずれか一方に異なる変動をいくつか与えた上でモデル関数の事後分布を推定し、その事後分布の平均値モデルと観測データとの誤差が一番小さくなる変動量を推定値とする方法である。この方法における処理の流れを図５に示す。
【００６４】
まず、与えられた参照観測データを、Parameter 1及びParameter 2の二つの方向にそれぞれ−５〜＋５の範囲で１ずつシフトさせた変形データを生成する（ステップＳ２１）。これにより、二つのパラメータのシフト量の組合せの総数だけ、参照観測データに基づく変形データ群が生成される。そして、生成された各変形データと与えられた対象観測データとに基づいて、一般的なマルチタスクベイズ最適化法におけるガウス過程のハイパーパラメータ最尤推定を行って対象システムのモデル関数の事後分布を推定する（ステップＳ２２）。即ち、この場合には、二つのパラメータのシフト量の組合せの総数の回数だけ、モデルを推定する必要がある。このモデル推定の際に最尤推定の尤度関数が得られるから、その尤度関数が最も大きくなるシフト量の組合せを求め、これを変動量の推定値として採用する（ステップＳ２３）。
【００６５】
こうして変動量が推定されたならば、補正処理部１１２はその変動量に基づいて対象観測データを補正し、補正済みモデル推定部１１３はその補正された対象観測データと参照観測データとに基づいて上述したように対象システムのモデル関数を推定すればよい。
【００６６】
上記第２の変動量推定方法を用いて観測データの変動量を推定し、その変動を補正したうえでマルチタスクベイズ最適化法を実施して最適パラメータ値を求めたときの効果を図７に示す。図７において、横軸は探索回数（図２におけるステップＳ１〜Ｓ４のループの繰返し回数）、縦軸は求まるパラメータ値の下で得られる観測値の最大値を示しており、（ａ）はシフト量を補正しない場合（つまり従来のマルチタスクベイズ最適化法の場合）、（ｂ）はシフト量を補正した場合（つまり本発明の一実施例による方法の場合）の結果である。この結果から、本発明の一実施例による方法では、参照システムと対象システムとの相関が低い場合（図７ではシフト幅が３Ｖや５Ｖの場合）であっても、少ない探索回数で以て最適なパラメータ値を探索できていることが分かる。
【００６７】
上記具体例では、モデル又はデータの変動として単純なシフトのみを想定していた。これは、図１２でいえば、図１２（ｃ）に示したように平均値が原信号から変更された場合に相当する。但し、多くの場合、パラメータ空間は２次元以上の多次元であり、パラメータ空間での伸縮や回転等の変動も生じる。しかしながら、参照システムは対象システムと関連するものであり、それら両システムのモデルの相関性が低下する主たる要因はシフトであることから、変動としてシフトのみを考慮して補正を行ってもパラメータ探索の効率を上げるうえで効果は十分にある。
【００６８】
図８は図６と同様に参照システムと対象システムのモデルを表したものであるが、この例では、単純なシフトではなく、参照システムと対象システムとでモデル自体が異なるものとしている。図９は、図８に示したモデルについて一般的なベイズ最適化法、従来のマルチタスクベイズ最適化法（補正処理なし）、及び本発明の一実施例による方法（補正処理あり）を実施して最適パラメータ値を求めたときの効果を比較した結果である。但し、変動量推定方法は上記第２の変動量推定方法と同じであり、変動はシフトしか考慮していない。
【００６９】
図９から明らかであるように、本発明の一実施例による方法では探索回数が少ない状態でも最適パラメータが探索できている。このように、参照システムと対象システムとの間のモデルの変動がシフトだけでない場合でも、変動がシフトであるとみなして変動量の推定及びその補正を行うことで、パラメータ探索の効率改善が十分に可能であることが確認できる。
【００７０】
もちろん、参照システムと対象システムとの間のモデルの変動として回転や伸縮などを考慮しても構わないが、その変動要因を増やすほど上記ステップＳ１２で生成すべき変形モデル又はステップＳ２１で生成すべき変形データが多くなる。そのため、計算量が大幅に増加し計算時間が延びることになる。こうした負荷の増大に対しパラメータ探索の効率の改善の程度は小さいため、シフト以外の変動を考慮してもそれほど有効ではない。
【００７１】
図１に示した上記実施例において、分析システム２０としては質量分析装置やクロマトグラフ装置など、様々な装置（システム）が想定されるが、特に効率的なパラメータ探索が必要とされるものの具体例として、タンデム型質量分析装置を用いたＬＣ−ＭＳやＧＣ−ＭＳにおけるＭＲＭ（多重反応モニタリング）測定の際のパラメータ最適化がある。ＭＲＭ測定では、実測結果に基づいて高感度な測定が可能であるパラメータを探索するが、測定対象である試料をイオン化する部分のインターフェイスパラメータには、温度やガス流量など安定状態になるまでの待ち時間が必要なパラメータが多いうえに、化合物、移動相等の条件毎に最適なパラメータが異なるため、パラメータの探索に長い時間を要したり試料を多量に消費しまったりするという問題がある。このように、探索すべきパラメータの数が多く、しかも一つのパラメータを見つけるための測定に時間を要する場合に、本発明を用いると、少ない測定回数で以て効率的にパラメータを探索することができるので、時間削減と試料節約の効果が特に大きい。
【００７２】
また、上記実施例は本発明の一例にすぎず、本発明の趣旨の範囲で適宜変形、修正、追加等を行っても本願特許請求の範囲に包含されることは当然である。
【００７３】
例えば上述したように、本発明に係るパラメータ探索方法の対象となるシステムは分析システムに限らず、何らかのパラメータ条件の下で観測値が得られるシステム、或いは、得られた観測値からＳＮ比、ダイナミックレンジ等の何らかの指標値が計算されるシステムであれば、どのようなシステムにも適用可能である。したがって、そのパラメータの種類も数値化できるものであれば、特に限定されない。
【符号の説明】
【００７４】
１０…パラメータ探索装置
１１…モデル推定部
１１１…変動量推定部
１１２…補正処理部
１１３…モデル推定部
１２…パラメータ決定部
１３…データ取得部
１４…繰返し制御部
【図１】
2019244474.tif 000003
【図２】
2019244474.tif 000004
【図３】
2019244474.tif 000005
【図４】
2019244474.tif 000006
【図５】
2019244474.tif 000007
【図６】
2019244474.tif 000008
【図７】
2019244474.tif 000009
【図８】
2019244474.tif 000010
【図９】
2019244474.tif 000011
【図１０】
2019244474.tif 000012
【図１１】
2019244474.tif 000013
【図１２】
2019244474.tif 000014
【図１３】
2019244474.tif 000015
【国際調査報告】
2019244474.tif 000016
2019244474.tif 000017
2019244474.tif 000018
2019244474.tif 000019