(57)【要約】【課題】時系列データが有する周期の成分組成を精度よく推定することができるようにする。【解決手段】位相計算部３２が、計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、単位時間毎に位相を計算する。ピーク判定部３４が、単位時間毎に計算された位相と、予め定められた閾値とを比較して、単位時間毎にピークのオンオフを判定する。周期推定部３６が、単位時間毎のピークのオンオフの判定結果に基づいて、単位時間の数で表される周期の各々について、時系列データが有する周期成分組成を推定する。【選択図】図４
2020135639.tif 000002
【特許請求の範囲】
【請求項１】
計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算する位相計算部と、
前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定するピーク判定部と、
前記単位時間毎のピークのオンオフの判定結果に基づいて、単位時間の数で表される周期の各々について、前記時系列データが前記周期の成分を有する割合を推定する周期推定部と、
を含む周期推定装置。
【請求項２】
前記位相計算部は、計測対象の個体の集合に含まれる個体の各々について、前記個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算し、
前記ピーク判定部は、前記計測対象の個体の集合に含まれる個体の各々について、前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定し、前記単位時間毎にピークのオンオフの代表値を求め、
前記周期推定部は、前記単位時間毎のピークのオンオフの代表値に基づいて、前記周期の各々について、前記計測対象の個体の集合の前記時系列データが前記周期の成分を有する割合を推定する請求項１記載の周期推定装置。
【請求項３】
前記位相計算部は、計測対象の個体の集合に含まれる個体の各々について、前記個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算し、
前記ピーク判定部は、前記計測対象の個体の集合に含まれる個体の各々について、前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定し、
前記周期推定部は、前記計測対象の個体の集合に含まれる個体の各々について、前記単位時間毎のピークのオンオフの判定結果に基づいて、前記周期の各々について、前記時系列データが前記周期の成分を有する割合を推定し、前記周期の各々について、前記時系列データが前記周期の成分を有する割合の代表値を求める請求項１記載の周期推定装置。
【請求項４】
前記単位時間は、１年であり、
前記計測対象の個体は、樹木作物であり、
前記計測値は、前記樹木作物の結実量である請求項１〜請求項３の何れか１項記載の周期推定装置。
【請求項５】
位相計算部が、計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算し、
ピーク判定部が、前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定し、
周期推定部が、前記単位時間毎のピークのオンオフの判定結果に基づいて、単位時間の数で表される周期の各々について、前記時系列データが前記周期の成分を有する割合を推定する、
周期推定方法。
【請求項６】
コンピュータを、
計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算する位相計算部、
前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定するピーク判定部、及び
前記単位時間毎のピークのオンオフの判定結果に基づいて、単位時間の数で表される周期の各々について、前記時系列データが前記周期の成分を有する割合を推定する周期推定部
として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、周期推定装置、方法、及びプログラムに係り、特に、時系列データが有する周期を推定するための周期推定装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
従来より、ドングリのMastingの数値実験結果に対して、離散的な時系列の自己相関関数が得られることが知られている（非特許文献１、２）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】井鷺 裕司、「植物個体の物質収支モデルでmastingを考える」、個体群生態学会会報 52、 49-54、 1995-06-30
【非特許文献２】Y. ISAGI, al.,"How Does Masting Happen and Synchronize?",J. theor. Biol. (1997)187,231-239
【発明の概要】
【発明が解決しようとする課題】
【０００４】
上記の非特許文献１、２に記載の方法では、離散的な時系列の自己相関関数を求めているため、原理的に、正しい周期成分以外の周期成分がアーチファクトとして現れてしまう。
また、離散的フーリエスペクトル（離散的フーリエ変換）により、時系列の平均的な周期成分を検出する手法も考えられるが、この手法でも、原理的に、正しい周期成分以外の周期成分がアーチファクトとして現れてしまう。
【０００５】
本発明は、上記の事情を鑑みてなされたもので、時系列データが有する周期の成分組成を精度よく推定することができる周期推定装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００６】
上記の目的を達成するために本発明に係る周期推定装置は、計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算する位相計算部と、前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定するピーク判定部と、前記単位時間毎のピークのオンオフの判定結果に基づいて、単位時間の数で表される周期の各々について、前記時系列データが前記周期の成分を有する割合を推定する周期推定部と、を含んで構成されている。
【０００７】
本発明に係る周期推定方法は、位相計算部が、計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算し、ピーク判定部が、前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定し、周期推定部が、前記単位時間毎のピークのオンオフの判定結果に基づいて、単位時間の数で表される周期の各々について、前記時系列データが前記周期の成分を有する割合を推定する。
【０００８】
本発明に係るプログラムは、コンピュータを、計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算する位相計算部、前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定するピーク判定部、及び前記単位時間毎のピークのオンオフの判定結果に基づいて、単位時間の数で表される周期の各々について、前記時系列データが前記周期の成分を有する割合を推定する周期推定部として機能させるためのプログラムである。
【０００９】
本発明によれば、位相計算部が、計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、前記単位時間毎に位相を計算する。ピーク判定部が、前記単位時間毎に計算された位相と、予め定められた閾値とを比較して、前記単位時間毎にピークのオンオフを判定する。そして、周期推定部が、前記単位時間毎のピークのオンオフの判定結果に基づいて、単位時間の数で表される周期の各々について、前記時系列データが前記周期の成分を有する割合を推定する。
【００１０】
このように、計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、単位時間毎に位相を計算し、予め定められた閾値と比較して、単位時間毎にピークのオンオフを判定し、単位時間の数で表される周期の各々について、前記時系列データが前記周期の成分を有する割合を推定することにより、時系列データが各周期の成分を有する割合を精度よく推定することができる。
【発明の効果】
【００１１】
以上説明したように、本発明の周期推定装置、方法、及びプログラムによれば、時系列データが有する周期の成分組成を精度よく推定することができる、という効果が得られる。
【図面の簡単な説明】
【００１２】
【図１】（Ａ）時系列データの一例を示す図、（Ｂ）時間毎の位相の一例を示す図、及び（Ｃ）位相の分布の一例を示す図である。
【図２】各周期を説明するための図である。
【図３】本発明の第１の実施の形態に係る周期推定装置を示すブロック図である。
【図４】本発明の第１の実施の形態に係る周期推定装置を示す機能ブロック図である。
【図５】本発明の第１の実施の形態に係る周期推定装置の周期推定処理ルーチンの内容を示すフローチャートである。
【図６】計測対象の個体の集合についての時系列データの一例を示す図である。
【図７】本発明の第２の実施の形態に係る周期推定装置の周期推定処理ルーチンの内容を示すフローチャートである。
【図８】本発明の第３の実施の形態に係る周期推定装置の周期推定処理ルーチンの内容を示すフローチャートである。
【図９】実施例における各周期の組成割合の推定結果を示す図である。
【図１０】（Ａ）アンサンブル時系列データの一例を示す図、（Ｂ）個別の各周期の組成割合の推定結果を示す図、（Ｃ）各周期の組成割合の推定結果の平均値を示す図、（Ｄ）各周期の組成割合の推定結果の中央値を示す図、（Ｅ）各周期の組成割合の推定結果の最頻値を示す図、及び（Ｆ）アンサンブル時系列データから得られる平均時系列データに対する、各周期の組成割合の推定結果を示す図である。
【図１１】（Ａ）時系列データの一例を示す図、（Ｂ）実施例における各周期の組成割合の推定結果を示す図、（Ｃ）比較例における各周期の組成割合の推定結果を示す図、及び（Ｄ）比較例における各周期の組成割合の推定結果を示す図である。
【図１２】（Ａ）時系列データの一例を示す図、（Ｂ）実施例における各周期の組成割合の推定結果を示す図、（Ｃ）比較例における各周期の組成割合の推定結果を示す図、及び（Ｄ）比較例における各周期の組成割合の推定結果を示す図である。
【発明を実施するための形態】
【００１３】
以下、図面を参照して本発明の実施の形態を詳細に説明する。以下では、樹木作物の毎年の結実量を表す時系列データが有する周期の成分組成を推定する周期推定装置に本発明を適用した場合を例に説明する。
【００１４】
＜本発明の実施の形態の概要＞
振動子集団の同期現象は様々な分野において遍在する。多年生植物の結実同期現象も同様である。これを隔年結果、マスティングとよぶ。樹木作物としては、柑橘（ウンシュウミカン、レモン、ポンカン）、ナッツ類（ピスタチオ、ピーカン）などの隔年結果現象が見られ、広葉果樹（リンゴ，ナシなど）においても隔年結果現象が見られる。また、樹木堅果（ブナ，ミズナラ，コナラ，クリなど）のマスティング現象は同期現象に相当し、スギ花粉なども同期現象が見られる。
【００１５】
これらの同期現象を解明し、その豊凶予測を行うためには、集団レベルと個体レベルの双方において変動の位相を求めることが必要である。
【００１６】
【数１】
2020135639.tif 000003
(1)
【００１７】
という時系列集団データに対して、下記の演算を行う。
【００１８】
i番目の時系列{xi(t), t=1,2,3,…,T}に対して、式（２）の演算を行う。
【００１９】
【数２】
2020135639.tif 000004
(2)
【００２０】
ここで、HTはヒルベルト変換、angle(W)は複素数Wの偏角、
2020135639.tif 000005
【００２１】
【数３】
2020135639.tif 000006
(3)
【００２２】
以下の式を満たす場合に、ONとし、それ以外をOFFと判定する。
【００２３】
|θ(t)|<π/2
【００２４】
判定するための関数をON(t)とし、個体i毎に式(4)で定める。
【００２５】
【数４】
2020135639.tif 000007
(4)
【００２６】
１個体の時系列データが図１（Ａ）のようであった場合（即ち式(1)においてn=1）、式(4)によってON年を黒四角，OFF年を白丸のように判定することができる。この時、図１（Ｂ）ではON年（黒四角）は全て、式(4)の条件を満たしていることが分かる。当然のことながら、複素平面において第１象限・第４象限にON年の状態がプロットされ、OFF年は第2象限・第3象限にある。
【００２７】
そして、各周期の成分を有する割合を以下のように計算する。
【００２８】
まず、周期Ｑを図２のように定義する。
【００２９】
図２（Ａ）のようにON-OFF-ONのパターンでは周期２であり枠で囲った２年間が周期２であると判定する。同様に、図２（Ｂ）では、ON-OFF-OFF-ONで、枠で囲った３年間が周期３、図２（Ｃ）では、ON-OFF-OFF-OFF-ONで、枠で囲った４年間が周期４となる。図２（Ｄ）も枠で囲った３年間が周期３となる。よって、ONの次にＱ−１回のOFFが連続して発生し、その後ONとなったパターンにおいて最初のＱ年間が周期Ｑであると判定する。
【００３０】
上記の手続きを以下の式（５）又は式（５）’で定式化する。
【００３１】
【数５】
2020135639.tif 000008
(5)
2020135639.tif 000009
(5)'
【００３２】
ここで、uは単位ステップ関数である。
【００３３】
【数６】
2020135639.tif 000010
(6)
【００３４】
式(5)又は式(5)’によって、時系列{x(t), t=1,2,…,T}に含まれるQ周期の度合いFP(Q)を概算できる。Tが大きくなるほど精度は高くなる。Q=1,2,…P<T-1であり、Pは適用する用途に応じて定める。これによって時系列データが有する各周期の組成割合が求まる。
【００３５】
［第１の実施の形態］
＜本発明の第１の実施の形態の周期推定装置の構成＞
図３に示すように、本発明の第１の実施の形態に係る周期推定装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６、ＨＤＤ１８、通信インタフェース２１、及びこれらを相互に接続するためのバス２２を備えている。
【００３６】
ＣＰＵ１２は、各種プログラムを実行する。ＲＯＭ１４には、各種プログラムやパラメータ等が記憶されている。ＲＡＭ１６は、ＣＰＵ１２による各種プログラムの実行時におけるワークエリア等として用いられる。記録媒体としてのＨＤＤ１８には、後述する周期推定処理ルーチンを実行するためのプログラムを含む各種プログラムや各種データが記憶されている。
【００３７】
本実施の形態における周期推定装置１０を、周期推定処理ルーチンを実行するためのプログラムに沿って、機能ブロックで表すと、図４に示すようになる。周期推定装置１０は、入力部２０、演算部３０、及び出力部４０を備えている。
【００３８】
入力部２０は、計測対象の個体の単位時間毎の計測値である時系列データを受け付ける。例えば、樹木作物１本の１年毎の結実量を表す時系列データを受け付ける。
【００３９】
演算部３０は、位相計算部３２、ピーク判定部３４、及び周期推定部３６を備えている。
【００４０】
位相計算部３２は、計測対象の個体の単位時間毎の計測値である時系列データに対して、ヒルベルト変換を用いた上記式（２）に従って、単位時間ｔ毎に位相θ（ｔ）を計算する。
【００４１】
ピーク判定部３４は、上記式（４）に従って、単位時間ｔ毎に計算された位相θ（ｔ）と、予め定められた閾値であるπ／２とを比較して、単位時間ｔ毎にピークのオンオフを判定する。
【００４２】
周期推定部３６は、単位時間ｔ毎のピークのオンオフの判定結果に基づいて、上記式（５）又は式（５）’に従って、単位時間の数で表される周期Ｑの各々について、時系列データが周期Ｑの成分を有する割合ＦＰ（Ｑ）を推定する。
【００４３】
出力部４０は、周期Ｑの各々について推定された、時系列データが周期Ｑの成分を有する割合ＦＰ（Ｑ）を出力する。例えば、これにより、樹木作物の１年毎の結実量を表す時系列データにおける豊作凶作の周期成分組成を知ることができる。
【００４４】
＜周期推定装置の動作＞
次に、本発明の実施の形態に係る周期推定装置１０の動作について説明する。
【００４５】
入力部２０によって、計測対象の個体の単位時間毎の計測値である時系列データを受け付けると、周期推定装置１０によって、図５に示す周期推定処理ルーチンが実行される。
【００４６】
まず、ステップＳ１００において、位相計算部３２は、計測対象の個体の単位時間毎の計測値である時系列データに対して、ヒルベルト変換を用いた上記式（２）に従って、単位時間ｔ毎に位相θ（ｔ）を計算する。
【００４７】
ステップＳ１０２において、ピーク判定部３４は、上記式（４）に従って、単位時間ｔ毎に計算された位相θ（ｔ）と、予め定められた閾値であるπ／２とを比較して、単位時間ｔ毎にピークのオンオフを判定する。
【００４８】
ステップＳ１０４において、周期推定部３６は、単位時間ｔ毎のピークのオンオフの判定結果に基づいて、上記式（５）又は式（５）’に従って、単位時間の数で表される周期Ｑの各々について、時系列データが周期Ｑの成分を有する割合ＦＰ（Ｑ）を推定し、出力部４０により出力して、周期推定処理ルーチンを終了する。
【００４９】
以上説明したように、本発明の第１の実施の形態に係る周期推定装置によれば、計測対象の個体の単位時間毎の計測値である時系列データに対してヒルベルト変換を用いて、単位時間毎に位相を計算し、予め定められた閾値と比較して、単位時間毎にピークのオンオフを判定し、周期の各々について、時系列データが当該周期の成分を有する割合を推定することにより、時系列データが有する周期の成分組成を精度よく推定することができる。
【００５０】
また、従来の同期現象を記述する方法は、定性的、外観的なものにとどまるため、世界中の多くの異なる研究分野で膨大なデータが取得され蓄積されているにもかかわらず、それらの相互関係、相違、同一性などを定量的に統計的に比較することができない。また、樹木作物（柑橘、広葉果樹、ナッツ類）は日本でも樹木堅果の結実量だけではなく、世界中でデータが蓄積されており、また、気象データや各種経済指標データの関係性を同期と周期組成の観点から定量的に解析することができない。これに対し、本実施の形態に係る周期推定装置によれば、時系列データが各周期の成分を有する割合を推定することができるため、種々の時系列データについて、定量的に解析することができる。
【００５１】
［第２の実施の形態］
次に、第２の実施の形態に係る周期推定装置について説明する。なお、第１の実施の形態に係る周期推定装置と同様の構成であるため、同一符号を付して説明を省略する。
【００５２】
＜概要＞
第２の実施の形態では、計測対象の個体の集合について、時系列データが有する周期成分組成を推定する点が、第１の実施の形態と異なっている。
【００５３】
上記の第１の実施の形態で説明したようにして、個体毎に各年ｔのON,OFFを判定できるが、第２の実施の形態では、n個の集団のON,OFFの判定を、式(７)のようにして行う。
【００５４】
【数７】
2020135639.tif 000011
(7)
【００５５】
式(７)では、t年のON,OFFの判定を{ＯＮi(t),i=1,2,…,n}の中央値(メディアン)として行っている。
【００５６】
それ以外に、modなど、集合データの代表値によって行うこともできる。さらに、任意の方法で代表値を定めることができるが、中央値(メディアン)を用いるのが好ましい。
【００５７】
図６は、式（１）においてn=10,T=30としたときの時系列データの集合（アンサンブル時系列）に適用した例である。個体毎のON、OFFは、式(4)より判定している。
【００５８】
n=１０個体からなる集団の時間tにおける平均値
2020135639.tif 000012
を太線で示している。式(7)を用いることに因って、各年の集団のON,OFFを判定し、ONを黒四角で示し、OFFを白丸で示した。図中AやBの場合はONと判定するのは自明である。また、ＣやＤのように、n=１０個体からなる集団の時間tにおける平均値（集合平均）
2020135639.tif 000013
がゼロの場合もOFFと判定することは自明と言える。しかし、点線枠線で示したような中間的な値をとる場合の判定がこれまでは恣意的であった。これに対して、本発明においては、式(7)によってON,OFFが判定できている。注目すべき点は、（２）はＯＦＦ年と判定されており、（１）、（３）、（４）はＯＮと判定されている。後者におけるn=１０個体からなる集団の時間tにおける平均値
2020135639.tif 000014
は前者におけるn=１０個体からなる集団の時間tにおける平均値
2020135639.tif 000015
より小さいにも関わらず、このような判定結果となっている。これは、１０個体それぞれでON、OFF判定を行い、各年においてON、OFFの占める割合が大きい方を当該年における集団のON、OFFとして判定しているためである（式(7)）。
【００５９】
これは、個体には個体間変異があり、すなわち、同じ樹種であっても、個体毎に大きさが大きく異なり、果実やドングリの生産量は個体毎に大きく異なる。よって、個体毎にON、OFFを判定し（式(4)）、それを式(7)で用いることによって、同期する集団のON,OFF判定を行うという方法は客観的で合理的となる。
【００６０】
＜本発明の第２の実施の形態の周期推定装置の構成＞
第２の実施の形態に係る周期推定装置１０の入力部２０は、計測対象の個体の集合に含まれる個体の各々についての単位時間毎の計測値である時系列データを受け付ける。例えば、複数本の樹木作物からなる集合に含まれる１本１本の１年毎の結実量を表す時系列データを受け付ける。
【００６１】
位相計算部３２は、計測対象の個体の集合に含まれる個体の各々について、当該個体の単位時間毎の計測値である時系列データに対して、ヒルベルト変換を用いた上記式（2）に従って、単位時間ｔ毎に位相θ（ｔ）を計算する。
【００６２】
ピーク判定部３４は、計測対象の個体の集合に含まれる個体の各々について、上記式（4）に従って、単位時間ｔ毎に計算された位相θ（ｔ）と、予め定められた閾値であるπ／２とを比較して、単位時間ｔ毎にピークのオンオフを判定し、単位時間ｔ毎にピークのオンオフの代表値を求める。
【００６３】
周期推定部３６は、単位時間ｔ毎のピークのオンオフの代表値に基づいて、上記式（5）に従って、単位時間の数で表される周期Ｑの各々について、時系列データが周期Ｑの成分を有する割合ＦＰ（Ｑ）を推定する。
【００６４】
出力部４０は、周期Ｑの各々について推定された、時系列データが周期Ｑの成分を有する割合ＦＰ（Ｑ）を出力する。例えば、単位時間毎の各個体の計測値ｘi（ｔ）と、周期Ｑの各々について推定された、時系列データが周期Ｑの成分を有する割合ＦＰ（Ｑ）とを表示する。これにより、例えば、複数本の樹木作物からなる集合に含まれる１本１本の１年毎の結実量を表す時系列データにおける豊作凶作の周期を知ることができる。
【００６５】
＜周期推定装置の動作＞
次に、本発明の第２の実施の形態に係る周期推定装置１０の動作について説明する。
【００６６】
入力部２０によって、計測対象の個体の集合に含まれる個体の各々についての単位時間毎の計測値である時系列データを受け付けると、周期推定装置１０によって、図７に示す周期推定処理ルーチンが実行される。
【００６７】
まず、ステップＳ２００において、位相計算部３２は、計測対象の個体の各々について、単位時間毎の計測値である時系列データに対して、ヒルベルト変換を用いた上記式（2）に従って、単位時間ｔ毎に位相θ（ｔ）を計算する。
【００６８】
ステップＳ２０２において、ピーク判定部３４は、計測対象の個体の各々について、上記式（4）に従って、単位時間ｔ毎に計算された位相θ（ｔ）と、予め定められた閾値であるπ／２とを比較して、単位時間ｔ毎にピークのオンオフを判定する。
【００６９】
ステップＳ２０４において、上記ステップＳ２０２の計測対象の個体の各々についての判定結果に基づいて、単位時間ｔ毎に、ピークのオンオフの代表値を求める。
【００７０】
ステップＳ２０６において、周期推定部３６は、単位時間ｔ毎のピークのオンオフの代表値に基づいて、上記式（5）又は式（5）'に従って、単位時間の数で表される周期Ｑの各々について、時系列データが周期Ｑの成分を有する割合ＦＰ（Ｑ）を推定し、出力部４０により出力して、周期推定処理ルーチンを終了する。
【００７１】
以上説明したように、第２の実施の形態に係る周期推定装置によれば、計測対象の個体の集合について、時系列データが各周期の成分を有する割合を推定することにより、計測対象の個体の集合について、時系列データが有する周期の成分組成を精度よく推定することができる。
【００７２】
また、樹木作物（柑橘、落葉果樹、ナッツ類）について、時系列データが各周期の成分を有する割合の推定結果を利用して、世界の果樹生産量の品種別、地域別予測を行うことができる。また、世界の果樹生産量変動の品種間、地域間の相関関係、因果関係の解明を行うことができる。また、気象、気候と果樹生産量変動の解明を行うことも可能である。
【００７３】
［第３の実施の形態］
次に、第３の実施の形態に係る周期推定装置について説明する。なお、第１の実施の形態に係る周期推定装置と同様の構成であるため、同一符号を付して説明を省略する。
【００７４】
＜概要＞
第３の実施の形態では、計測対象の個体の集合について、周期を推定する際に、個体の各々について推定された周期の代表値を求める点が、第１の実施の形態と異なっている。
【００７５】
上記式(1)で示した時系列データの集合{xi(t),i=1,2,…,n, t=1,2,…,T}について周期を推定する場合には、以下の式(8-1)〜式(8-2)、式(8-1)'〜式(8-2)'の何れかに従って、個体毎に各周期Ｑの成分の割合を推定し、式（9）に従って、各周期Ｑの成分の割合の代表値として中央値を求める。
【００７６】
【数８】
2020135639.tif 000016
(8-1)
2020135639.tif 000017
(8-2)
2020135639.tif 000018
(8-1)'
2020135639.tif 000019
(8-2)'
【００７７】
【数９】
2020135639.tif 000020
(9)
【００７８】
なお、代表値としてはmedian（中央値）の他、modやmean（平均）を採用することも可能である．
【００７９】
＜本発明の第３の実施の形態の周期推定装置の構成＞
第３の実施の形態に係る周期推定装置１０の位相計算部３２は、計測対象の個体の集合に含まれる個体の各々について、当該個体の単位時間毎の計測値である時系列データに対して、ヒルベルト変換を用いた上記式（2）に従って、単位時間ｔ毎に位相θ（ｔ）を計算する。
【００８０】
ピーク判定部３４は、計測対象の個体の集合に含まれる個体の各々について、上記式（4）に従って、単位時間ｔ毎に計算された位相θ（ｔ）と、予め定められた閾値であるπ／２とを比較して、単位時間ｔ毎にピークのオンオフを判定する。
【００８１】
周期推定部３６は、計測対象の個体の集合に含まれる個体の各々について、単位時間ｔ毎のピークのオンオフの判定結果に基づいて、上記式(8-1)〜式(8-2)、式(8-1)'〜式(8-2)'の何れかに従って、単位時間の数で表される周期Ｑの各々について、時系列データが周期Ｑの成分を有する割合ＦＰi（Ｑ）を推定する。そして、周期推定部３６は、周期Ｑの各々について、上記式（9）に従って、時系列データが周期Ｑの成分を有する割合ＦＰi（Ｑ）の代表値ＦＰ（Ｑ）を求める。
【００８２】
出力部４０は、周期Ｑの各々について推定された、時系列データが周期Ｑの成分を有する割合の代表値ＦＰ（Ｑ）を出力する。例えば、上記図７に示すように、単位時間毎の各個体の計測値ｘi（ｔ）と、周期Ｑの各々について推定された、時系列データが周期Ｑの成分を有する割合の代表値ＦＰ（Ｑ）とを表示する。
【００８３】
＜周期推定装置の動作＞
次に、本発明の第３の実施の形態に係る周期推定装置１０の動作について説明する。
【００８４】
入力部２０によって、計測対象の個体の集合に含まれる個体の各々についての単位時間毎の計測値である時系列データを受け付けると、周期推定装置１０によって、図８に示す周期推定処理ルーチンが実行される。
【００８５】
まず、ステップＳ３００において、位相計算部３２は、計測対象の個体の各々について、単位時間毎の計測値である時系列データに対して、ヒルベルト変換を用いた上記式（2）に従って、単位時間ｔ毎に位相θ（ｔ）を計算する。
【００８６】
ステップＳ３０２において、ピーク判定部３４は、計測対象の個体の各々について、上記式（4）に従って、単位時間ｔ毎に計算された位相θ（ｔ）と、予め定められた閾値であるπ／２とを比較して、単位時間ｔ毎にピークのオンオフを判定する。
【００８７】
ステップＳ３０４において、計測対象の個体の各々について、上記ステップＳ３０２の判定結果に基づいて、上記式(8-1)〜式(8-2)、式(8-1)'〜式(8-2)'に従って、単位時間の数で表される周期Ｑの各々について、時系列データが周期Ｑの成分を有する割合ＦＰi（Ｑ）を推定する。
【００８８】
ステップＳ３０６では、上記ステップＳ３０４で計測対象の個体の各々について推定された、周期Ｑの成分を有する割合ＦＰi（Ｑ）に基づいて、周期Ｑの各々について、上記式（9）に従って、時系列データが周期Ｑの成分を有する割合ＦＰi（Ｑ）の代表値ＦＰ（Ｑ）を求め、出力部４０により出力して、周期推定処理ルーチンを終了する。
【００８９】
以上説明したように、第３の実施の形態に係る周期推定装置によれば、計測対象の個体の集合について、時系列データが各周期の成分を有する割合を推定することにより、計測対象の個体の集合について、時系列データが有する周期の成分組成を精度よく推定することができる。
【００９０】
＜実施例＞
コンピュータシミュレーションによって、１００個体、２４年間の結実パターンの同期現象のデータセット（アンサンブル時系列データ（時系列データ集合））を生成し、上記第２の実施の形態で説明した方法を用いて、各周期Ｑの成分を有する割合を算出した。周期２から周期７までの算出結果を図９右に示した。周期２は36.36％、周期３は42.86%、周期４は20%であり、周期２から周期４までの合計で99.22％を占めることが分かる。
【００９１】
また、図１０（Ａ）に示した、コンピュータシミュレーションにより生成したアンサンブル時系列データに対して、上記の実施の形態で説明した方法を適用した。個体数は１０６、時間は１５ステップである。図１０（Ｂ）に１０６個体毎の周期成分を示した。周期３が卓越しており、周期２も有意な成分があり、また周期５の成分を持つ個体がわずかであるが存在していることが分かる。集合の代表値を平均値、中央値、最頻値にとったものを図１０（Ｃ）、（Ｄ）、（Ｅ）に示す。また、上記図１０（Ａ）のアンサンブル時系列データから得られる平均時系列に対する適用結果を図１０（Ｆ）に示す。時系列データ集合の周期成分（周期組成）に対しては、目的に応じて、図１０（Ｂ）〜（Ｆ）の何れかの表示方法を採用することができる。
【００９２】
＜実験例＞
上記の第１の実施の形態で説明した方法を、比較手法（コレログラム（自己相関関数）および離散的フーリエ変換の各々を用いた手法）と比較した結果について説明する。図１１（Ａ）に示す周期３の時系列データに対して、３つの方法を適用した結果をそれぞれ図１１（Ｂ）〜（Ｄ）に示す。図１１（Ｂ）に示すように、上記の第１の実施の形態で説明した方法では、周期３の成分のみで生成されていることが明確である。図１１（Ｃ）に示すように、コレログラムを用いた手法では、周期３の整数倍の周期（周期６、周期９）が現れる。一方、その他の周期成分も現れる。図１１（Ｄ）に示すように、離散的フーリエ変換を用いた手法では、周期３の成分が卓越的であるが、それ以外に存在しない周期成分がアーチファクトとして現れる。
【００９３】
また、図１２（Ａ）に示すような、周期３、周期６、周期２で構成された時系列データについて、３つの方法を適用した結果をそれぞれ図１２（Ｂ）〜（Ｄ）に示す。図１２（Ｂ）に示すように、上記の第１の実施の形態で説明した方法では、明確に周期成分（周期組成）を表示できているが、図１２（Ｃ）、（Ｄ）に示すように、コレログラムを用いた手法や離散的フーリエ変換を用いた手法では、存在しない周期成分がアーチファクトとして現れる。このように、アーチファクト成分の発生は、コレログラムおよび離散的フーリエ変換において原理的に発生するものである。
【００９４】
＜変形例＞
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００９５】
例えば、上記の実施の形態において、個体毎の計測値のｎ個体の平均値（集合平均）の時系列
2020135639.tif 000021
に対して上記式(2)と同様に、ｎ個体の集合を代表する位相の時系列を（１２）式で求める。
【数１０】
2020135639.tif 000022
（１２）
ただし、
2020135639.tif 000023
これに式(4)を適用して、ON(t)を求めて、それを式(5)又は式(5)'に適用することにより、計測対象の個体の集合について、時系列データが各周期の成分を有する割合を推定するようにしてもよい。
【００９６】
また、ON(t)を求める際の閾値がπ／２である場合を例に説明したが、これに限定されるものではなく、π／２以外を閾値として用いてよい。
【００９７】
また、樹木作物の結実量以外の時系列データの周期推定に、本発明を適用してもよい。例えば、森林に関する計測値（スギ花粉の発生量、ブナ、ミズナラ、コナラなどのドングリの発生量）の時系列データの周期成分組成の推定に、本発明を適用してもよい。また、野生動物の発生数の時系列データの周期成分組成の推定に、本発明を適用してもよい。これにより、場所ごとの野生動物の発生数の予測や、場所ごとの野生動物の密度の予測が可能となる。
【００９８】
また、野生動物による被害に関する計測値の時系列データの周期成分組成の推定に、本発明を適用してもよい。これにより、野生動物被害予測が可能となり、また、野生動物被害対策の最適化が可能となる。
【００９９】
また、サーカディアンリズムに関する計測値の時系列データや、脳波に関するに関する計測値の時系列データの周期成分組成の推定に、本発明を適用してもよい。また、心拍、脈波、生体信号に関する時系列データの周期成分組成の推定に、本発明を適用してもよい。また、身体運動の各部に関する計測値の時系列データの周期成分組成の推定、電子回路に関する計測値の時系列データの周期成分組成の推定、心筋細胞に関する計測値の時系列データの周期成分組成の推定に、本発明を適用してもよい。これにより、種々の同期現象の解析や分析が可能となる。
【０１００】
また、本発明のプログラムは、記憶媒体に格納して提供するようにしてもよい。
【符号の説明】
【０１０１】
１０   周期推定装置
２０   入力部
３０   演算部
３２   位相計算部
３４   ピーク判定部
３６   周期推定部
４０   出力部
【図１】
2020135639.tif 000024
【図２】
2020135639.tif 000025
【図３】
2020135639.tif 000026
【図４】
2020135639.tif 000027
【図５】
2020135639.tif 000028
【図６】
2020135639.tif 000029
【図７】
2020135639.tif 000030
【図８】
2020135639.tif 000031
【図９】
2020135639.tif 000032
【図１０】
2020135639.tif 000033
【図１１】
2020135639.tif 000034
【図１２】
2020135639.tif 000035