JP4542395B2 - 非定常時系列データ分類方法、装置、プログラム及びそのプログラムを記録した記録媒体 - Google Patents

Description

本発明は、非定常な時系列データの分類を行う技術に係り、特に、音声基本周波数Ｆ０のような欠測値を含む非定常な時系列データを分析する方法、装置、プログラム及びそのプログラムを記録した記録媒体に関する。

発達心理学の分野では、親と子の音声の基本周波数（Ｆ０）パタンについて類似性を分類し、子の月齢増加に伴いパタン類似性がどのように変化するかという先行研究がある。音声のＦ０は、発話のイントネーションやアクセント等の韻律情報を含み、発声器官の発達、性差、加齢に伴い変化するため、その特性は人間の言語獲得のメカニズムを明らかにする情報を含んでいると考えられる。Ｆ０パタン解析はそれらの情報を抽出するために必要で、その分析過程の一環として、客観的なパタン分類手法が要求される。

一方、データの時系列分類の研究は様々な手法が提案されている。例えば、データマイニングでは、シーケンシャルな記号列や、欠測値のない時系列データに対して、全区間あるいは小区間ごとで時系列間の距離を計算して行うクラスタリング等が提案されている（例えば、非特許文献１参照。）。また、時系列解析の分野では、判別分析を基礎とした時系列分類がよく議論されている（例えば、非特許文献２参照。）。
Keogh, E., Lin, J. and Truppel, W. (2003) "Clustering of time series subsequences is meaningless: implications for previous and future research", ICDM. Taniguchi, M. and Kakizawa, Y. (2000) Asymptotic Theory of Statistical Inference for Time Series, 477-534, Springer in statistics. springer-verlag New York.

しかし、従来、Ｆ０パタンのような欠測値を含む非定常時系列データを客観的な測度に基づき分類する方法は開示されていなかった。
すなわち、Ｆ０データは、有声音と無声音とから成り、この無声音を欠測値する時系列データとみなすことができる。ところが、従来の時系列分類は、データ全区間に対し定常性を仮定したスペクトル密度関数か、或いは局所定常性を仮定することによって得られた時変スペクトル密度関数間の測度により分類する手法が主流であり、欠測値を含む非定常時系列データに対するアプローチはなされていない。そのため、従来、Ｆ０パタンのような欠測値を含む非定常時系列データに客観的な測度に基づく統計解析を適用することはできず、Ｆ０パタンの分類は手作業で行うしかなかった。

この発明はこのような点に鑑みてなされたものであり、欠測値を含む非定常時系列データを客観的な測度に基づき分類することが可能な技術を提供することを目的とする。

この発明では上記課題を解決するために、欠測値を含む時系列データの変動傾向をトレンド推定値として求めるために、多項式回帰モデルを適用する。
すなわち、まず多項式回帰モデルを適用し、時系列データｙ^(s) _ｔの変動傾向を示すトレンド関数の推定関数Ｔ＾^(s)（ｔ）をサンプルｓごとに算出する。次に、この推定関数Ｔ＾^(s)（ｔ）から判別関数を生成し、その判別関数を用いて推定関数Ｔ＾^(s)（ｔ）のクラスター分析を行い、その分類結果を出力する。

ここで、好ましくは、推定関数Ｔ＾^(s)（ｔ）の算出は、最小二乗法によって、時系列データｙ^(s) _ｔからトレンド関数の回帰係数推定値β＾^(s) _ＬＳＥを算出し、
Ｔ＾^(s)（ｔ）＝（１，ｔ，ｔ^２，...，ｔ^ｐ−１）・β＾^(s) _ＬＳＥ
によって行う。なお、この次数ｐはモデルｙ^（ｓ） _ｔ＝Ｔ＾^(s)（ｔ）＋ｕ^(s) _ｔの情報量基準ＡＩＣ^(s)（ｐ）を最小値化する値である。また、判別関数は、好ましくは二群判別分析によって生成される。

この発明では、多項式回帰モデルを適用して、欠測値を含む時系列データの変動傾向を示すトレンド関数の推定関数Ｔ＾^(s)（ｔ）を算出し、この推定関数Ｔ＾^(s)（ｔ）を用いてパタン分類を行うこととしたため、Ｆ０パタンのような欠測値を含む非定常時系列データを客観的な測度に基づき分類することが可能となる。

以下、この発明の実施の形態を説明する。なお、以下ではまず本形態の原理について説明を行った後、その装置及び処理について説明を行う。
〔本形態の原理〕
＜基本的枠組み＞
本形態の基本的枠組みは、１．音声の基本周波数（Ｆ０）データのような欠測値を含む時系列データに対して多項式回帰モデルを適用し、その変動パタンを示すトレンド関数を推定し、２．推定されたトレンド関数に対して２群判別分析を適用して判別関数を生成し、３．判別関数をもとに最近隣接法でパタンを分類することである。

＜欠測値を含む時系列データに対する多項式回帰モデルの適用＞
Ｆ０データのような欠測値を含む時系列データに対し、
ｙ_ｔ＝Ｔ（ｔ）＋ｕ_ｔ …(1)
という時系列回帰モデルを導入する。
ここで、｛ｔ｝はサンプリング時刻を示す離散値であり、｛ｙ_ｔ｝は各サンプリング時刻ｔにおける観測値を示す時系列データである。なお、各｛ｙ_ｔ｝は、Ｆ０データのような欠測値を含む時系列データであり、あるサンプリング時刻ｔにおける時系列データは存在しない（欠測値）。例えば、欠測値の区間がサンプリング時刻ｔ_１＜ｔ≦ｔ_２の区間であった場合、式(1)は
ｙ_ｔ＝Ｔ（ｔ）＋ｕ_ｔ ，ｔ∈｛１，２，...，ｔ_１｝∪｛ｔ_２＋１，ｔ_２＋２，...，ｎ｝…(2)
と記述される。なお、ｎはサンプリング時刻ｔの最大値である。

また｛Ｔ（ｔ）｝は、時系列データの傾向を示すトレンド関数であり、回帰係数
β≡（β_０，β_１，...，β_ｐ−１）^Ｔ∈Ｒ^ｐ
に対して、
Ｔ（ｔ）＝β_０＋β_１ｔ＋...＋β_ｐ−１ｔ^ｐ−１
と表される。ここで、＊^Ｔは＊の転置行列を示し、Ｒ^ｐはｐ次の実数を示す。また｛≡｝は左辺を右辺のように定義することを意味する。なお、トレンド関数は全区間（ｔ＝｛１，２，...，ｎ｝）で定義される。
また、｛ｕ_ｔ｝はノイズ値を示すノイズ項であり、定常過程とする。

本形態では、このように定義された時系列回帰モデルを用い、欠測値を有する時系列データから、回帰係数βの最小二乗推定値である回帰係数推定値β＾^（ｓ） _ＬＳＥを求め、全区間（ｔ＝｛１，２，...，ｎ｝）におけるトレンド関数の推定関数
Ｔ＾^(s)（ｔ）＝（１，ｔ，ｔ^２，...，ｔ^ｐ−１）・β＾^(s) _ＬＳＥ …(3)
求める。なお、この例における次数ｐの決定は、赤池情報量規準
ＡＩＣ^(s)（ｐ）＝−２ｌｏｇ（最大対数尤度）＋２（ｐ＋１）
を適用して行う（例えば、中溝高好 (1988)「信号解析とシステム同定」コロナ社 201-205.参照）。またｓは各サンプルを示すパラメータであり、ここでいう「サンプル」とは、１つの観測結果に対応する時系列データ列を意味する。

＜推定関数に対する２群判別分析の適用＞
式(1)で示された時系列データが、以下のような２つのカテゴリのどちらかに分類されると考える：
Π_１：トレンド関数Ｔ_１（ｔ）をもつ式（１）のモデル
Π_２：トレンド関数Ｔ_２（ｔ）をもつ式（１）のモデル

そして、トレンド関数Ｔ（ｔ）をもつ新しい観測値が得られた場合に、

のような距離測度を基準にして分類することにする。実際の問題では、トレンド関数Ｔ（ｔ）は未知である。そこで、式（３）で定義された推定関数Ｔ＾^(s)（ｔ）を式（４）のＴ（ｔ）に代入したものを距離測度とする。故に、分類手順では、Ｌ（Ｔ＾^(s)（ｔ）：Ｔ_２）＞Ｌ（Ｔ＾^(s)（ｔ）：Ｔ_１）という関係を満たすときにカテゴリΠ_１に属すると判断し、それ以外の場合にはカテゴリΠ_２に属すると判断する。

また、上述した２カテゴリに分類する判別手法はシンプルであるが、これをｋカテゴリの判別に拡張することができる。この場合、式(1)で示された時系列データが、以下のようなｋ個のカテゴリのいずれかに分類されると考える：
Π_ｒ：トレンド関数Ｔ_ｒ（ｔ）をもつ式（１）のモデル ，ｒ＝{１，２，...，ｋ}
また、その分類手順としては、すべてのｒ_２（ｒ_１，ｒ_２∈ｒ，ｒ_２≠ｒ_１）に対して、Ｌ（Ｔ＾^(s)（ｔ）：Ｔ_ｒ２）＞Ｌ（Ｔ＾^(s)（ｔ）：Ｔ_ｒ１）という関係を満たすときに、その推定関数Ｔ＾^(s)（ｔ）はカテゴリΠ_ｒ１に属すると判断する。

〔装置及び処理〕
＜概要＞
図１は、一部が欠測した観測値の時系列データの傾向を分類する非定常時系列データ分類方法を例示した図である。
本形態の方法では、前処理において各サンプルｓのサンプリング時刻ｔに対応する時系列データｙ^（ｓ） _ｔがメモリに格納されている。
分類処理が開始されると、まずこのメモリから時系列データｙ^(s) _ｔが読み込まれる（ステップＳ１）。そして、上述の多項式回帰モデルを適用して、この時系列データｙ^(s) _ｔのサンプルｓごとの変動傾向を示すトレンド関数の推定関数Ｔ＾^(s)（ｔ）＝（１，ｔ，ｔ^２，...，ｔ^ｐ−１）・β＾^(s) _ＬＳＥが生成される（ステップＳ２）。ここで、次数ｐは情報量基準ＡＩＣ^(s)（ｐ）を最小値化するように決定され、回帰係数推定値β＾^(s) _ＬＳＥは最小二乗法によって算出される。

次に、生成された推定関数Ｔ＾^(s)（ｔ）に二群判別分析を適用して判別関数が生成される（ステップＳ３）。この例の場合、二群判別分析の判別基準Ｔ＾^(ｒ)（ｔ）は、ステップＳ２で生成された推定関数Ｔ＾^(s)（ｔ）から選択される。そして、生成された判別関数を用い、最近隣接法等により推定関数Ｔ＾^(s)（ｔ）のクラスター分析が行われ（ステップＳ４）、その分析結果が出力される（ステップＳ５）。

＜詳細＞
次に、本形態における装置及び処理の詳細について説明する。
図２は、本形態における非定常時系列データ分類装置１全体の機能構成を例示したブロック図である。なお、本形態の非定常時系列データ分類装置１は、例えば、ＣＰＵ（central processing unit）、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）等がバスで接続されたノイマン型コンピュータに所定のプログラムを実行させることにより構築されるものである。また、図３は、図２における推定関数生成手段１３の詳細構成を例示したブロック図であり、図４（ａ）は、次数選択手段１４の詳細構成を例示したブロック図であり、図４（ｂ）は、判別分析手段１５の詳細構成を例示したブロック図である。なお各ブロック図中の矢印はデータの流れを示し、そのデータの出力元と出力先とが情報のやり取り可能に接続されていることを示しているが、制御手段１１から出入りするデータの流れ及び接続に関する記載は省略してある。

また、図５及び図６は、制御手段１１の制御のもと、非定常時系列データ分類装置１が実行するデータ分類処理を説明するための図である。
以下、これらの図を用いて非定常時系列データ分類装置１の構成及び処理について説明する。

＜前処理＞
まず前処理として、サンプルｓ（ｓ＝｛１，２，...，ｓ_ｍａｘ｝）のサンプリング時刻ｔにおける時系列データｙ^(s) _ｔがメモリ１２の観測データ領域１２ａ（図２）に格納される。
この時系列データｙ^(s) _ｔは、一部が欠測した観測値の時系列データであり、例えば、マイクロホン等により観測した周波数領域の音声信号を短時間離散フーリエ変換等によって時間領域に変換し、その時間領域データの逆数をとること等によって生成されるＦ０データ等である。図７（ａ）はこのようなＦ０データを例示したグラフである。ここで、図７（ａ）のグラフの縦軸はｌｏｇＦ０／（Ｆ０の平均値）を示し、横軸はサンプリング時刻ｔを示している。図７（ａ）に例示するように、このＦ０データはｔ＝４１０ｍｓｅｃ〜４９０ｍｓｅｃ付近で欠測している。なおＦ０データの場合、この欠測部分は無声音部分に相当する。

図７（ｂ）は、メモリ１２の観測データ領域１２ａに格納される時系列データｙ^（ｓ） _ｔのデータ構造の例示である。この例の時系列データｙ^（ｓ） _ｔはサンプリング時刻ｔ_１+1≦ｔ≦ｔ_２の区間が欠測した
｛ｙ_ｔ；ｔ＝１，...，ｔ_１，ｔ_２＋１，...，ｎ｝ …(5)
というデータであり、各サンプリング時刻ｔとその時系列データｙ^（ｓ） _ｔとが対応付けられて格納されている。また、欠測区間のサンプリング時刻ｔ_１+1≦ｔ≦ｔ_２に対応する時系列データは存在せずｎｕｌｌ値となっている。なお、ｙの下付き添え字の（ｔ１）はｔ_１を示し、（ｔ２）はｔ_２を示す。また、図７（ｂ）はサンプリング時刻ｔ_１+1≦ｔ≦ｔ_２の一箇所のみが欠測している場合の例であるが、複数箇所が欠測しているデータの場合も同様である。すなわち、時系列データが存在するサンプリング時刻ｔにその時系列データｙ^(s) _ｔが対応付けられて格納され、時系列データが存在しないサンプリング時刻ｔにはｎｕｌｌ値が対応付けられる。なお、メモリ１２の観測データ領域１２ａに格納されるサンプリング時刻ｔと時系列データｙ^(s) _ｔとの組をＹ^(s) _０と表現する。

＜データ分類処理＞
次に、前処理によってメモリ１２観測データ領域１２ａに格納された時系列データｙ^(s) _ｔの分類処理について例示する。
まず、推定関数生成手段１３（図２）が、多項式回帰モデルを適用し、メモリ１２の観測データ領域１２ａに格納された時系列データｙ^(s) _ｔの変動傾向を示すトレンド関数の推定関数Ｔ＾^(s)（ｔ）をサンプルｓごとに算出する。本形態では、最小二乗法によって、時系列データｙ^(s) _ｔからトレンド関数の回帰係数推定値β＾^(s) _ＬＳＥを算出し、Ｔ＾^(s)（ｔ）＝（１，ｔ，ｔ^２，...，ｔ^ｐ−１）・β＾^(s) _ＬＳＥの演算によって推定関数Ｔ＾^(s)（ｔ）を算出する。

まず、制御手段１１が、サンプルを特定するパラメータｓに１を代入し（ステップＳ１１）、次数を特定するパラメータｐに１を代入する（ステップＳ１２）。なお、これらの代入結果はレジスタ１１ａ（図２）に格納される。
次に、分離ベクトル生成手段１３ａ（図３）が、メモリ１２の観測データ領域１２ａ（図２）にアクセスし、サンプルｓの各サンプリング時刻ｔにおける時系列データｙ^(s) _ｔの列を読み出す。そして、分離ベクトル生成手段１３ａは、これを欠測値前後の時系列データ列である分離ベクトルＳＹ^(s) _ｉ（欠測値を含まない時系列データ）に分離する（ステップＳ１３）。例えば、式（５）の例の場合、欠測区間はサンプリング時刻ｔ_１+1≦ｔ≦ｔ_２の一箇所のみであり、この各サンプリング時刻ｔにおける時系列データｙ^(s) _ｔの列は、以下のような２つの分離ベクトルＳＹ^(s) _１，ＳＹ^(s) _２に分離される。
ＳＹ^(s) _１＝（ｙ_１，ｙ_２，...，ｙ_（ｔ１）），ＳＹ^(s) _２＝（ｙ_{（ｔ２）＋１}，ｙ_{（ｔ２）＋２}，...，ｙ_ｎ） …(6)

なお、この例では、欠測区間が１箇所のみであるため２つの分離ベクトルＳＹ^(s) _１，ＳＹ^(s) _２に分離されるが、欠測区間がｘ箇所あった場合には、ｘ＋１個の分離ベクトルＳＹ^(s) _ｉ（ｉ＝｛１，２，...，ｘ＋１｝）に分離されることになる。
このように生成された各分離ベクトルＳＹ^(s) _ｉはベクトル結合手段１３ｂに送られ、ベクトル結合手段１３ｂは、これらの分離ベクトルＳＹ^(s) _ｉを結合した結合ベクトルＹ^(s)を生成する（ステップＳ１４）。例えば、式（６）の例の場合、ベクトル結合手段１３ｂは、結合ベクトル
Ｙ^(s)＝（ｙ_１，ｙ_２，...，ｙ_（ｔ１），ｙ_{（ｔ２）＋１}，ｙ_{（ｔ２）＋２}，...，ｙ_ｎ） …(7)
を生成する。

生成された結合ベクトルＹ^(s)は、転置ベクトル生成手段１３ｃに送られ、転置ベクトル生成手段１３ｃは、この結合ベクトルＹ^(s)の転置ベクトルＴＹ^(s)を生成して、レジスタ１３ｇに格納する（ステップＳ１５）。
また、行列生成手段１３ｄが、メモリ１２の観測データ領域１２ａ（図２）にアクセスし、サンプルｓに係る欠測値前後のサンプリング時刻ｔの列を読み出す。例えば、式（５）の例の場合、ｔ＝１，...，ｔ_１というサンプリング時刻列と、ｔ＝ｔ_２＋１，...，ｎというサンプリング時刻列とが読み出される。さらに、行列生成手段１３ｄはレジスタ１１ａ（図２）からパラメータｐを読み出す。これらを読み出した行列生成手段１３ｄは、読み出した欠測値前後のサンプリング時刻列ごとに、このサンプリング時刻列を構成する各サンプリング時刻ｔに対応する従属変数（１，ｔ，ｔ^２，…，ｔ^ｐ−１）を行とする従属変数行列Ｚ^(s) _ｉを生成する（ステップＳ１６）。例えば、ｔ＝１，...，ｔ_１というサンプリング時刻列と、ｔ＝ｔ_２＋１，...，ｎというサンプリング時刻列とが読み出されていた場合、行列生成手段１３ｄは、以下のような従属変数行列Ｚ^(s) _１，Ｚ^(s) _２を生成する。

生成された従属変数行列Ｚ^(s) _ｉは行列結合手段１３ｅに送られ、行列結合手段１３ｅは、送られた従属変数行列Ｚ^(s) _ｉを列方向（縦）に結合した結合行列Ｚ^(s)を生成し、レジスタ１３ｇに格納する（ステップＳ１７）。例えば、式（８）の例の場合には、

という結合行列が生成・格納される。
生成された結合行列Ｚ^(s)はさらに転置行列生成手段１３ｆにも送られ、転置行列生成手段１３ｆは、結合行列Ｚ^(s)の転置行列ＴＺ^(s)を生成し、レジスタ１３ｇに格納する（ステップＳ１８）。

次に、係数推定値算出手段１３ｈが、レジスタ１３ｇから転置ベクトルＴＹ^(s)、結合行列Ｚ^(s)、転置行列ＴＺ^(s)を読み出し、これらを用いてトレンド関数の回帰係数推定値
β＾^(s) _ＬＳＥ＝（Ｚ^(s)・ＴＺ^(s)）^−１・ＴＺ^(s)・ＴＹ^(s) …(10)
を算出する（ステップＳ１９）。例えば、式（７）〜（９）の例の場合、式（１０）の右辺は、
β＾^(s) _ＬＳＥ＝｛（Ｚ^(s) _１）^Ｔ・Ｚ^(s) _１＋（Ｚ^(s) _２）^Ｔ・Ｚ^(s) _２｝^−１・｛（Ｚ^(s) _１）^Ｔ・（ＳＹ^(s) _１）^Ｔ＋（Ｚ^(s) _２）^Ｔ・（ＳＹ^(s) _２）^Ｔ｝ …(11)
となる。なお、式（１０）の演算によって回帰係数の推定値が求まることは、例えば「中溝高好 (1988) “信号解析とシステム同定”コロナ社 10-11.」に記載されている。

算出された回帰係数推定値β＾^(s) _ＬＳＥは係数乗算手段１３ｉに送られ、係数乗算手段１３ｉは、さらにレジスタ１１ａ（図２）からパラメータｐを読み出し、トレンド関数の推定関数
Ｔ＾^(s)（ｔ）＝（１，ｔ，ｔ^２，...，ｔ^ｐ−１）・β＾^(s) _ＬＳＥ …(12)
を算出する（ステップＳ２０）。なお、推定関数Ｔ＾^(s)（ｔ）はすべてのサンプリング時刻ｔ＝｛１，２，.....，ｎ｝上の関数である。
算出された推定関数Ｔ＾^(s)（ｔ）はパラメータｐとともにメモリ１２のテーブル領域１２ｂに送られ、そこで相互に対応付けられて格納される。また、推定関数Ｔ＾^(s)（ｔ）は次数選択手段１４（図２）にも送られる。

次数選択手段１４は、モデルｙ^(s) _ｔ＝Ｔ＾^(s)（ｔ）＋ｕ^(s) _ｔの情報量基準ＡＩＣ^(s)（ｐ）を最小値化する値を次数ｐとして選択する。すなわちまず、制御手段１１においてサンプリング時刻ｔの１つを示すパラメータｔ’に１を代入し、このパラメータｔ’をレジスタ１１ａに格納する（ステップＳ２１）。
次に次数選択手段１４のノイズ項算出手段１４ａ（図４（ａ））が、メモリ１２の観測データ領域１２ａ（図２）からサンプリング時刻ｔ＝ｔ’に対応する時系列データｙ^(s) _ｔ’を読み出し、これとこれに対応する推定関数Ｔ＾^(s)（ｔ’）を用い、ノイズ項
ｕ^(s) _ｔ＝ｙ^(s) _ｔ−Ｔ＾^(s)（ｔ） …(13)
を算出してレジスタ１４ｂ（図４（ａ））に格納する（ステップＳ２２）。

次に制御手段１１（図２）が、レジスタ１１ａからパラメータｔ’を読み出し、これに１を加算した値を新たなパラメータｔ’としてレジスタ１１ａに格納する（ステップＳ２３）。次に、制御手段１１は、レジスタ１１ａに格納され最新のｔ’がｔ’＞ｎであるか否かを判断し、ｔ’＞ｎとなるまでステップＳ２２〜２４までの処理を繰り返す（ステップＳ２４）。これは、サンプルｓに関するすべてのサンプリング時刻ｔにおけるノイズ項ｕ^(s) _ｔを算出することに相当する。
ｔ’＞ｎとなると、分散算出手段１４ｃ（図４（ａ））が、レジスタ１４ｂからノイズ項ｕ^(s) _ｔ（ｔ＝｛１，２，...，ｎ｝）を読み出し、これを用いてノイズ項ｕ^(s) _ｔの分散σ_（ｓ） ^２を算出する（ステップＳ２５）。算出された分散σ_（ｓ） ^２は最大対数尤度算出手段１４ｄに送られ、最大対数尤度算出手段１４ｄは、メモリ１２の観測データ領域１２ａからサンプリング時刻ｔの最大値ｎを抽出する。そして、最大対数尤度算出手段１４ｄは、これらｔの最大値ｎ及び分散σ_（ｓ） ^２を用い、モデルｙ^（ｓ） _ｔ＝Ｔ＾^(s)（ｔ）＋ｕ^(s) _ｔの最大対数尤度
ＭＬ^(s)＝ｃｏｎｓｔ−ｎ・（ｌｎσ_（ｓ） ^２）／２ …(14)
を算出する（ステップＳ２６）。なお「ｃｏｎｓｔ」は定数を示す。算出された最大対数尤度ＭＬ^(s)は情報量基準算出手段１４ｅに送られ、情報量基準算出手段１４ｅは、情報量基準
ＡＩＣ^(s)（ｐ）＝−２ｌｏｇ（ＭＬ^(s)）＋２（ｐ＋１） …(15)
を算出し、それをメモリ１２のテーブル領域１２ｂ（図２）にパラメータｐと対応付けて格納する（ステップＳ２７）。

次に、制御手段１１は、レジスタ１１ａに格納されている最新のパラメータｐを抽出し、抽出したパラメータｐがｐ_ｍａｘ（定められたｐの最大値）であるか否か判断し（ステップＳ２８）、ｐ＝ｐ_ｍａｘとなるまで、パラメータｐに１を加えた値を新たなパラメータｐとしてレジスタ１１ａに格納し（ステップＳ２９）ステップＳ１６に戻る処理を繰り返す。すなわち、従属変数行列Ｚ^(s) _ｉを生成する手順（ステップＳ１６）から情報量基準ＡＩＣ^(s)（ｐ）を算出してメモリに格納する手順（ステップＳ２７）までは、パラメータｐの値を更新しながら所定回数繰り返される。これはｐ＝｛１，２，...，ｐ_ｍａｘ｝におけるすべての情報量基準ＡＩＣ^(s)（ｐ）を算出することに相当する。なお、ｐ_ｍａｘの値が小さすぎると正確な分類ができず、逆に大きすぎると演算速度が遅くなる。そのため、ｐ_ｍａｘの値は例えば５〜６程度が望ましい。

その後、推定関数選択手段１４ｆ（図４（ａ））が、メモリ１２のテーブル領域１２ｂ（図２）にアクセスし、公知の並び替えアルゴリズム等を用い、サンプルｓに対するＡＩＣ^(s)（ｐ）の最小値を求める。そして、推定関数選択手段１４ｆは、このＡＩＣ^(s)（ｐ）を最小値化するｐに対応する推定関数Ｔ＾^(s)（ｔ）（T^^(s)(t) for AIC^(s)(p)_min）をメモリ１２のテーブル領域１２ｂ（図２）から選択し、選択関数領域１２ｃに格納する（ステップＳ３０）。次に、制御手段１１が、レジスタ１１ａに格納されている最新のパラメータｓを抽出し、パラメータｓがｓ_ｍａｘ（ｓの最大値であり、メモリ１２の観測データ領域１２ａに格納されたサンプル数）であるか否か判断する（ステップＳ３１）。そして、制御手段１１は、ｓ＝ｓ_ｍａｘと判断するまで、レジスタ１１ａに格納されている最新のパラメータｓに１を加算した値を新たなパラメータｓとして（ステップＳ３２）レジスタ１１ａに格納し、ステップＳ１２に戻る処理を繰り返す。これは、すべてのサンプルについて、ＡＩＣ^(s)（ｐ）を最小値化するｐに対応する推定関数Ｔ＾^(s)（ｔ）（T^^(s)(t) for AIC^(s)(p)_min）を求めることに相当する。

ステップＳ３１でパラメータｓ＝ｓ_ｍａｘと判断されると、次に判別分析手段１５が、それぞれのパラメータｓ＝｛１，２，...，ｓ_ｍａｘ｝に対する判別関数を生成する。なお、判別関数の生成は、例えば情報量基準ＡＩＣ（ｐ）を最小値化するｐに対応する推定関数Ｔ＾^(s)（ｔ）を用い、二群判別分析によって行われる。

この例ではまず、基準関数設定手段１５ａ（図４（ｂ））が、２個以上の判別基準Ｔ＾^（ｒ）（ｔ）（ｒ∈ｓ）を設定してレジスタ１５ｂに格納する（ステップＳ３３）。この例の場合、基準関数設定手段１５ａは、メモリ１２の選択関数領域１２ｃに格納されている推定関数Ｔ＾^(s)（ｔ）（T^^(s)(t) for AIC^(s)(p)_min）のいずれかを判別基準Ｔ＾^（ｒ）（ｔ）として抽出し、レジスタ１５ｂに格納する。なお、ノイズ項ｕ^（ｓ） _ｔが発散しない場合、判別基準のＴ＾^（ｒ）（ｔ）間の距離が離れていれば離れているほど誤分類確率が低くなるため、できるだけ相互に距離が離れた推定関数Ｔ＾^(s)（ｔ）を判別基準Ｔ＾^（ｒ）（ｔ）として選択することが望ましい。ここでＴ＾^（ｓ１）（ｔ）とＴ＾^（ｓ２）（ｔ）との距離とは、

を意味する。

［Ｔ＾^（ｒ）（ｔ）間の距離が離れているほど誤分類確率が低くなることの証明］
まず、欠測値の区間を示すサンプリング時刻ｔ_１，ｔ_２及び｛ｕ_ｔ｝について以下のような仮定を設定する。
[仮定１]
（i）ｔ_１＝［αｎ］（０＜α＜１）、ｔ_２＝［（１−γ）ｎ］（０＜γ＜１）そしてα＋γ＜１、ここで［・］はガウス記号を示す。
（ii）｛ｕ_ｔ｝は以下に定義される定常過程である。

を満たす。

のようなスペクトル密度関数を持つ０平均定常過程となる。

また、以下のような仮定を設定する。
[仮定２]

[定理]
上述した仮定１、２が成立する場合、次の理由により次式が成立する。

を基礎とした分類が基本的に良好であることと、その判別関数が漸近正規性をもつことを意味する。

[定理の証明]
ノイズ項からなる転置ベクトルＵ_１＝（ｕ_１，ｕ_２，...，ｕ_（ｔ１））^Ｔと、Ｕ_２＝（ｕ_{（ｔ２）＋１}，ｕ_{（ｔ２）＋２}，...，ｕ_ｎ）^Ｔとを用いて、Ｕ≡（ｕ_１，ｕ_２，...，ｕ_{（ｔ１），}ｕ_{（ｔ２）＋１}，ｕ_{（ｔ２）＋２}，...，ｕ_ｎ）^Ｔ＝（Ｕ_１ ^Ｔ，Ｕ_２ ^Ｔ）^Ｔとする。
前述の式（１０）（１１）から回帰係数の推定値の誤差は、
β^^(s) _LSE-β=(Z^(s)・TZ^(s))^-1・TZ^(s)・U={(Z^(s) ₁)^T・Z^(s) ₁+(Z^(s) ₂)^T・Z^(s) ₂}^-1・{(Z^(s) ₁)^T・U₁+(Z^(s) ₂)^T・U₂}
となる。また、Ｓ_ｎ≡ｄｉａｇ｛ｓ_１，...，ｓ_ｐ｝が成り立つとする。ここで、S_j ²=z_1j ^T・z_1j+z_2j ^T・z_2jである。なお、以下の証明では、Z^(s) ₁ ，(Z^(s) ₁)^T，Z^(s) ₂ ，(Z^(s) ₂)^Tを、それぞれZ₁ ，Z₁'，Z₂ ，Z₂'と表す。

が成り立つとする。

式（１６）と仮定１（i）から

となる。また、式（１７）と同様に、以下の式を得ることができる。

の漸近分散行列を評価する。

まず、

ここで式（１８）から、ｐ×ｐ行列Ａ＝｛Ａ_ｊｋ｝のとき、ｎ→∞でＡ_ｎ→Ａとなる。σ（ｊ）＝Ｅ｛ｕ_ｔ＋１ｕ_ｔ｝とすると（Ｅ｛＊｝は＊の平均値を示す）、式（１８）と「E J Hannan, "Multiple time series" ,New York, Wiley [1970] , p.216」に記載された定理８から、以下の式を示す事ができる。

ということがいえる。これは時間差が開くとノイズ項の相関も落ちてくるということを意味する。

そして、これと式（１８）から以下の関係が成り立つ。

従って、式（２１）を式（１９）に代入し、ｐ×ｐ行列Ａ＝｛Ａ_ｊｋ｝のとき、ｎ→∞でＡ_ｎ→Ａとなることと、式（２０）とにより、

が成り立つことが分かる。

さて、判別問題に戻って、

が成り立つことがいえる。ここで、ｗ_ｔ＝（１，ｔ，ｔ^２，...，ｔ^ｐ−１）^Ｔであり、Ｔ_１（ｔ），Ｔ_２（ｔ）はそれぞれＴ＾^(１)（ｔ），Ｔ＾^（２）（ｔ）に相当し、＊’は、＊の転置行列を意味する。

ここで式（２２）から

についても同様に証明できる（「Ｔ＾^（ｒ）（ｔ）間の距離が離れているほど誤分類確率が低くなることの証明」終わり）。

次に、判別関数生成手段１５ｃが、レジスタ１５ｂから２つの判別基準Ｔ＾^（ｒ１）（ｔ），Ｔ＾^（ｒ２）（ｔ）（ｒ１，ｒ２∈ｒ，ｒ１≠ｒ２）を抽出し、メモリ１２の選択関数領域１２ｃからいずれかのパラメータｓに対する推定関数Ｔ＾^(s)（ｔ）（T^^(s)(t) for AIC^(s)(p)_min）を読み出し、判別関数

を生成してメモリ１２の判別関数領域１２ｄ（図２）に格納する（ステップＳ３４）。なお、このステップＳ３４の処理は、ｒ１，ｒ２∈ｒ、ｒ１≠ｒ２を満たすｒ１，ｒ２の組と、ｒ１，ｒ２以外のパラメータｓ＝｛１，２，...，ｓ_ｍａｘ｝との組み合わせすべてについて行われる。

その後、クラスター分析手段１６が、メモリ１２の判別関数領域１２ｄから判別関数Ｄ＾^{（ｓ，ｒ１，ｒ２）}を抽出し、それを用いて推定関数Ｔ＾^(s)（ｔ）のクラスター分析を行う（ステップＳ３５）。この例では、最近隣接クラスター分析（「統計学辞典」東洋経済新報社 ｐ．３８５ 参照）を適用してクラスター分析を行い、その分析結果を樹形型で出力する。

〔シミュレーションによる手法の妥当性の検討〕
次に本形態における非定常時系列データ分類手法が妥当であるかどうかを、モンテカルロシミュレーション実験（津田孝夫（１９９２）「モンテカルロ法とシミュレーション」培風館 参照）により検討する。
＜シミュレーションに使用するサンプル＞
このシミュレーションでは、実際の観測データを用いるのではなく、以下のように人工的に生成したサンプルを用いてデータ傾向の分類を行った。

・次数ｐ：１から４の間の整数からランダムに決定（擬似乱数を用いて生成）
・回帰係数β＝（β_０，β_１，...，β_ｐ−１）^Ｔ：−０．８から０.８の間の擬似乱数を採用
・トレンド関数：Ｔ（ｔ）＝β_０＋β_１ｔ＋...＋β_ｐ−１ｔ^ｐ−１
（ただし、所定の値を超える項については、その値以下となるような補正を行った）
・ノイズ項ｕ_ｔ：平均０、分散１の白色雑音系列を採用
・サンプリング点の数：ｎ＝１００
・サンプルデータ値：ｙ_ｔ＝Ｔ（ｔ）＋ｕ_ｔ
・欠測値の位置：任意に決定
・サンプル数：２０
図８（ａ）は、このように生成された２０個のサンプルを示したグラフである。なお、これらのグラフにおける横軸はサンプリング時刻ｔであり、縦軸はサンプルデータ値ｙ_ｔである。また、各グラフの上部に付された番号は各サンプルを示すパラメータｓの値である。

＜本形態の非定常時系列データ分類手法の適用＞
このシミュレーションでは、図８（ａ）のサンプルに対し、前述した非定常時系列データ分類手法を適用する。
図８（ｂ）は、図８（ａ）のサンプルに対し、本形態の非定常時系列データ分類手法を適用した場合に生成されるトレンド関数の推定関数Ｔ＾^(s)（ｔ）（T^^(s)(t) for AIC^(s)(p)_min）を示したグラフである。なお、これらのグラフにおける横軸はサンプリング時刻ｔであり、縦軸は推定関数Ｔ＾^(s)（ｔ）の出力値である。また、各グラフの上部に付された番号はｓの値である。

また、このシミュレーションでは、３番目と１８番目のサンプルに対応するトレンド関数の推定関数Ｔ＾^（３）（ｔ），Ｔ＾^（１８）（ｔ）を判別基準Ｔ＾^（ｒ）（ｔ）とし、判別関数Ｄ＾^{（ｓ，３，１８）}（ｓ＝｛１，...，２０｝）を算出した（式（１６））。そして、これらの判別関数Ｄ＾^{（ｓ，３，１８）}に対して、最近隣接クラスター分析を適用し、樹形図描かせた。
図９は、このように生成された樹形図である。なお、この図における横軸の数字は各サンプルを示すパラメータｓの値を示し、縦軸は高さ（判別基準からの距離）を示している。

この図に示すように、この推定関数Ｔ＾^(s)（ｔ）は、右下がり(falling)と右上がり（rising）との２つの大きなカテゴリに分かれた。また、右下がり(falling)は、さらに直線的な右下がり（Ｃ１）、若干凸に近いカーブ（Ｃ２）及び波打った形状（Ｃ３）に分かれ、右上がり（rising）は、さらに指数関数的な右上がり（Ｃ４）と、直線的な右上がり（Ｃ６）と、波打った形状（Ｃ５）と分かれた。この結果より、全体的にデータの変動傾向を示すパタンが類似したもの同士で分類されたことが分かる。これにより、本形態の手法が、統計的に妥当で、実データに適用し得る手法であることを確認ことができた。

〔実データへの本形態の適用〕
次に、本形態の手法の実データへの適用例を示す。
＜実データ＞
乳幼児音声データベースから、２４ヶ月齢の子とその父母のＦ０データに対し、本形態のデータ分類手法を適用した。このデータベースは、５組の日本人の乳幼児とその両親の音声を長期的に収録したものである。乳幼児は東京近郊に生まれ育ち、その両親は標準語を日常語として話す。デジタルオーディオテープレコーダで、彼らの自宅で静かな環境の中、日常会話を録音した。その際、特別なタスクは設けていない。子の誕生から５年間に渡り、月に１から４時間程度の頻度でおこなった。そのデータベースは、音声ファイルに変換され、Ｆ０、発話の特性、時間などと共にデータベースに収められている。データベースに関する詳しい情報は、「Amano, S. , Nakatani, T. and Kondo, T. (2003) “Fundamental frequency analysis of longitudinal recording in a Japanese infant speech database”, In: Proceedings of ICPhS-2003, 1983-1986, Aug.」参照。

＜データ解析＞
Ｆ０データを解析する際に問題となるのは、Ｆ０データの長さ、すなわち発話長に依存したデータ長は、当然発話ごとに異なる。そこで最も頻度が高い、すなわち最も頻度高く発話している発話長に相当するデータ長付近のデータを対象とすることにする。乳児、母親は３００ｍｓ、６００ｍｓのデータ長付近のデータを取り扱った。父親はデータ数が全体的に少ないため、２００ｍｓから４００ｍｓの間のデータ長を取り上げた。母親の発話のＦ０パタンを図１０に示す。なお、図１０の各グラフの上部にふされた数字は各サンプルを示すパラメータｓの値である。このデータに対し、前述した多項式回帰モデルを適用し、トレンド関数の推定を行う。その後におこなう判別分析の際、データ長が揃っていることが望ましいので、トレンド関数の推定関数Ｔ＾^(s)（ｔ）をリサンプリングし、データ長を揃えた。具体的には式（１６）を計算する際、推定関数Ｔ＾^(s)（ｔ）に対し、各発話を通して最大振幅と最小振幅の差で規格化した。また、判別関数の算出は、図１０における３番目と２２番目のサンプルに対応するトレンド関数の推定関数Ｔ＾^（３）（ｔ），Ｔ＾^（２２）（ｔ）を判別基準Ｔ＾^（ｒ）（ｔ）とし、判別関数Ｄ＾^{（ｓ，３，２２）}（ｓ＝｛１，...，４６｝）を算出した（式（１６））。

そして、これらの判別関数Ｄ＾^{（ｓ，３，２２）}に対して、最近隣接クラスター分析を適用し、樹形図描かせた。図１１（ａ）は、このように生成された樹形図である。なお、この図の横軸は各サンプルに対応するパラメータｓを示し、縦軸は高さを示す。図１１（ａ）に示すように、本形態の手法を適用することにより、平坦、右下がり、右上がり、ベル型といった適切なＦ０の形状に分類することができた。また、図１１（ｂ）（ｃ）は、それぞれ乳児、父親に関する分類結果を示した図である。
以上説明した通り、本形態の手法により、Ｆ０データの変動傾向を特定のパタンに分類することに成功した。これにより、従来マニュアルで行われていたＦ０パタン分析を客観的な統計的測度に基づいて行うことが可能となった。

なお、この発明は上述の実施の形態に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
また、上述の構成をコンピュータによって実現する場合、上述の装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

本発明は、例えば、発達心理学で乳幼児音声を用いた言語獲得の変遷を分析する際に適用が可能である。更に、ロボットとの対話システムにおける、発話のイントネーション、文末表現のマッチング等認識に関連した分野に代表される、欠測値を含み平均値が推移する変動を示すデータの分類が要求される分野にも応用できる。

本形態の非定常時系列データ分類方法を例示した図。 本形態における非定常時系列データ分類装置全体の機能構成を例示したブロック図。 図２における推定関数生成手段の詳細構成を例示したブロック図。 （ａ）は、次数選択手段の詳細構成を例示したブロック図であり、（ｂ）は、判別分析手段の詳細構成を例示したブロック図。 非定常時系列データ分類装置が実行するデータ分類処理を説明するための図。 非定常時系列データ分類装置が実行するデータ分類処理を説明するための図。 （ａ）はＦ０データを例示したグラフ。（ｂ）は、メモリの観測データ領域に格納される時系列データｙ^（ｓ） _ｔのデータ構造の例示。 （ａ）は生成された２０個のサンプルを示したグラフ。（ｂ）は、（ａ）のサンプルに対し、本形態の非定常時系列データ分類手法を適用した場合に生成されるトレンド関数の推定関数Ｔ＾^(s)（ｔ）を示したグラフ。 樹形図によるシミュレーションパタン分類結果を示した図。 母親の発話のＦ０パタンの例示。 （ａ）は樹形図による母親実データのパタン分類結果を示した図。（ｂ）（ｃ）は、それぞれ乳児、父親に関する分類結果を示した図。

符号の説明

１ 非定常時系列データ分類装置
１２ メモリ
１３ 推定関数生成手段
１４ 次数選択手段
１５ 判別分析手段
１６ クラスター分析手段

Claims (8)

1. コンピュータを、推定関数生成手段と、基準関数設定手段と、判別関数生成手段と、クラスター分析手段と、を有する非定常時系列データ分類装置として機能させて実行し、一部が欠測した観測値の時系列データの傾向を分類する非定常時系列データ分類方法であって、
   サンプルｓのサンプリング時刻ｔにおける上記時系列データｙ^（ｓ） _ｔがメモリに格納されており、
   上記推定関数生成手段が、ｙ_ｔをサンプリング時刻ｔにおける観測値を示す時系列データとし、*^Ｔを*の転置とし、Ｔ（ｔ）を時系列データｙ_ｔの変動傾向を示すトレンド関数Ｔ（ｔ）＝β_０＋β_１ｔ＋...＋β_ｐ−１ｔ^ｐ−１とし、（β_０，β_１，...，β_ｐ−１）^Ｔを回帰係数βとし、ｕ_ｔをノイズ値とした場合における、多項式回帰モデルｙ_ｔ＝Ｔ（ｔ）＋ｕ_ｔに対し、メモリに格納された上記時系列データｙ^（ｓ） _ｔ のうち、欠測値以外の時系列データを結合して得られた欠測値を含まない時系列データｙ ^（ｓ） _ｔ ’を適用し、ｙ_ｔ＝ｙ ^（ｓ） _ｔ ’とした場合における多項式回帰モデルｙ_ｔ＝Ｔ（ｔ）＋ｕ_ｔのトレンド関数Ｔ（ｔ）の推定関数Ｔ＾^（ｓ）（ｔ）をサンプルｓごとに算出し、算出した推定関数Ｔ＾^（ｓ）（ｔ）をメモリに格納する手順と、
   上記基準関数設定手段が、メモリに格納されている上記推定関数Ｔ＾^（ｓ）（ｔ）から、２個以上の推定関数Ｔ＾^（ｓ）（ｔ）を選択し、選択した推定関数Ｔ＾^（ｓ）（ｔ）を判別基準Ｔ＾^（ｒ）（ｔ）（r∈s）としてメモリに格納する手順と、
   上記判別関数生成手段が、メモリから２つの判別基準Ｔ＾^（ｒ１）（ｔ），Ｔ＾^（ｒ２）（ｔ）(ｒ１，ｒ２∈ｒ、ｒ１≠ｒ２)を抽出し、上記推定関数Ｔ＾^（ｓ）（ｔ）に対する判別関数
   

   

   
   (ｎはｔの最大値)
   を生成してメモリに格納する手順と、
   上記クラスター分析手段が、メモリから上記判別関数Ｄ＾^{（ｓ，ｒ１，ｒ２）}を抽出し、当該判別関数Ｄ＾^{（ｓ，ｒ１，ｒ２）}を用い、メモリから読み出した上記推定関数Ｔ＾^（ｓ）（ｔ）のクラスター分析を行い、その分析結果を出力する手順と、
   を実行する非定常時系列データ分類方法。
2. 請求項１記載の非定常時系列データ分類方法であって、
   上記推定関数Ｔ＾^（ｓ）（ｔ）を算出する手順は、
   上記推定関数生成手段が、
   最小二乗法によって、上記時系列データｙ^（ｓ） _ｔから上記回帰係数βの推定値である回帰係数推定値β＾^（ｓ） _ＬＳＥを算出し、
   Ｔ＾^（ｓ）（ｔ）＝（１，ｔ，ｔ^２，...，ｔ^ｐ−１）・β＾^（ｓ） _ＬＳＥ
   の演算によって上記推定関数Ｔ＾^（ｓ）（ｔ）を算出する手順であり、
   上記ｐは、モデルｙ^（ｓ） _ｔ＝Ｔ＾^（ｓ）（ｔ）＋ｕ^（ｓ） _ｔの情報量基準ＡＩＣ^（ｓ）（ｐ）を最小値化する値である、
   ことを特徴とする非定常時系列データ分類方法。
3. 請求項１又は２記載の非定常時系列データ分類方法であって、
   上記推定関数生成手段は、分離ベクトル生成手段と、ベクトル結合手段と、転置ベクトル生成手段と、行列生成手段と、行列結合手段と、転置行列生成手段と、係数推定値算出手段と、係数乗算手段とを有し、
   上記推定関数Ｔ＾^（ｓ）（ｔ）を算出する手順は、
   上記分離ベクトル生成手段が、上記サンプルｓごとの各時系列データｙ^（ｓ） _ｔを、欠測値前後の時系列データ列である分離ベクトルＳＹ^（ｓ） _ｉに分離する手順と、
   上記ベクトル結合手段が、上記分離ベクトルＳＹ^（ｓ） _ｉを結合した結合ベクトルＹ^（ｓ）を生成する手順と、
   上記転置ベクトル生成手段が、上記結合ベクトルＹ^（ｓ）の転置ベクトルＴＹ^（ｓ）を生成する手順と、
   上記行列生成手段が、上記欠測値前後のサンプリング時刻列ごとに、従属変数（１，ｔ，ｔ^２，…，ｔ^ｐ−１）を行とする従属変数行列Ｚ^（ｓ） _ｉを生成する手順と、
   上記行列結合手段が、上記従属変数行列Ｚ^（ｓ） _ｉを結合した結合行列Ｚ^（ｓ）を生成する手順と、
   上記転置行列生成手段が、上記結合行列Ｚ^（ｓ）の転置行列ＴＺ^（ｓ）を生成する手順と、
   上記係数推定値算出手段が、上記トレンド関数の回帰係数推定値
   β＾^（ｓ） _ＬＳＥ＝（Ｚ^（ｓ）・ＴＺ^（ｓ））^−１・ＴＺ^（ｓ）・ＴＹ^（ｓ）
   を算出する手順と、
   上記係数乗算手段が、上記推定関数
   Ｔ＾^（ｓ）（ｔ）＝（１，ｔ，ｔ^２，...，ｔ^ｐ−１）・β＾^（ｓ） _ＬＳＥ
   を算出する手順と、を有する、
   ことを特徴とする非定常時系列データ分類方法。
4. 請求項３記載の非定常時系列データ分類方法であって、
   上記推定関数Ｔ＾^（ｓ）（ｔ）を算出する手順の後に行われる、
   ノイズ項算出手段が、上記時系列データｙ^（ｓ） _ｔ及びこれに対応する上記推定関数Ｔ＾^（ｓ）（ｔ）を用い、ノイズ項ｕ^（ｓ） _ｔ＝ｙ^（ｓ） _ｔ−Ｔ＾^（ｓ）（ｔ）を算出する手順と、
   分散算出手段が、上記ノイズ項ｕ^（ｓ） _ｔの分散σ_（ｓ） ^２を算出する手順と、
   最大対数尤度算出手段が、ｔの最大値ｎ及び上記分散σ_（ｓ） ^２を用い、モデルｙ^（ｓ） _ｔ＝Ｔ＾^（ｓ）（ｔ）＋ｕ^（ｓ） _ｔの最大対数尤度
   ＭＬ^（ｓ）＝ｃｏｎｓｔ−ｎ・（ｌｎ σ_（ｓ） ^２）／２
   を算出する手順と、
   情報量基準算出手段が、情報量基準ＡＩＣ^（ｓ）（ｐ）＝−２ｌｏｇ（ＭＬ^（ｓ））＋２（ｐ＋１）を算出してメモリに格納する手順と、を有し、
   上記従属変数行列Ｚ^（ｓ） _ｉを生成する手順から上記情報量基準ＡＩＣ^（ｓ）（ｐ）を算出してメモリに格納する手順までは、ｐの値を更新しながら所定回数繰り返され、
   上記判別関数を生成する手順は、
   上記情報量基準ＡＩＣ^（ｓ）（ｐ）を最小値化するｐに対応する上記推定関数Ｔ＾^（ｓ）（ｔ）から判別関数を生成する手順である、
   ことを特徴とする非定常時系列データ分類方法。
5. 請求項１から４の何れかに記載の非定常時系列データ分類方法であって、
   上記時系列データｙ^（ｓ） _ｔは、
   音声の基本周波数（Ｆ０）データである、
   ことを特徴とする非定常時系列データ分類方法。
6. 一部が欠測した観測値の時系列データの傾向を分類する非定常時系列データ分類装置であって、
   サンプルｓのサンプリング時刻ｔにおける上記時系列データｙ^（ｓ） _ｔを格納するメモリと、
   ｙ_ｔをサンプリング時刻ｔにおける観測値を示す時系列データとし、*^Ｔを*の転置とし、Ｔ（ｔ）を時系列データｙ_ｔの変動傾向を示すトレンド関数Ｔ（ｔ）＝β_０＋β_１ｔ＋...＋β_ｐ−１ｔ^ｐ−１とし、（β_０，β_１，...，β_ｐ−１）^Ｔを回帰係数βとし、ｕ_ｔをノイズ値とした場合における、多項式回帰モデルｙ_ｔ＝Ｔ（ｔ）＋ｕ_ｔに対し、メモリに格納された上記時系列データｙ^（ｓ） _ｔ のうち、欠測値以外の時系列データを結合して得られた欠測値を含まない時系列データｙ ^（ｓ） _ｔ ’を適用し、ｙ_ｔ＝ｙ^（ｓ） _ｔ ’とした場合における多項式回帰モデルｙ_ｔ＝Ｔ（ｔ）＋ｕ_ｔのトレンド関数Ｔ（ｔ）の推定関数Ｔ＾^（ｓ）（ｔ）をサンプルｓごとに算出し、算出した推定関数Ｔ＾^（ｓ）（ｔ）をメモリに格納する推定関数生成手段と、
   メモリに格納されている上記推定関数Ｔ＾^（ｓ）（ｔ）から、２個以上の推定関数Ｔ＾^（ｓ）（ｔ）を選択し、選択した推定関数Ｔ＾^（ｓ）（ｔ）を判別基準Ｔ＾^（ｒ）（ｔ）（r∈s）としてメモリに格納する基準関数設定手段と、
   メモリから２つの判別基準Ｔ＾^（ｒ１）（ｔ），Ｔ＾^（ｒ２）（ｔ）(ｒ１，ｒ２∈ｒ、ｒ１≠ｒ２)を抽出し、上記推定関数Ｔ＾^（ｓ）（ｔ）に対する判別関数
   

   

   
   (ｎはｔの最大値)
   を生成してメモリに格納する判別関数生成手段と、
   メモリから上記判別関数Ｄ＾^{（ｓ，ｒ１，ｒ２）}を抽出し、当該判別関数Ｄ＾^{（ｓ，ｒ１，ｒ２）}を用い、メモリから読み出した上記推定関数Ｔ＾^（ｓ）（ｔ）のクラスター分析を行い、その分析結果を出力するクラスター分析手段と、
   を有することを特徴とする非定常時系列データ分類装置。
7. 請求項１から５の何れかに記載の非定常時系列データ分類方法の各手順をコンピュータに実行させるためのプログラム。
8. 請求項７記載のプログラムを格納したコンピュータ読取可能な記録媒体。

Images