JPWO2018207524A1 - 識別方法、分類分析方法、識別装置、分類分析装置および記憶媒体 - Google Patents

Abstract

［課題］計測データ集合から適切に不適合データを識別して、例えば、先端センシングデバイスによる計測結果の信頼性向上に寄与する識別方法、計測データに対する高精度の分類分析を行うことのできる分類分析方法、識別装置、分類分析装置、識別用記憶媒体および分類分析用記憶媒体を提供することである。［解決手段］パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、あらかじめ求めた特徴量を機械学習のための学習データとし、ＰＵ分類手法に基づく分類器により高精度に識別した不適合データが取り除かれた被分析データから得られる特徴量を変数にして、分類分析プログラムを実行することによって分析物に関する分類分析が行われる。高精度に識別した不適合データを取り除いた被分析データにより分析物に関する分類分析が高精度に実行可能である。［選択図］図４

Description

本発明は、計測系から得られた計測データに含まれる不適合データを識別する識別方法、該不適合データを取り除いたデータによる分類分析を行う分類分析方法、識別装置、分類分析装置および記憶媒体に関する。

例えば、非特許文献１に記載されているように、ナノセンシング、微量計測、量子計測など先端センシングデバイス開発分野において、微細・微量な対象を計測するためのデバイスが次々と開発されている。

ＷＯ２０１３−１３７２０９号公報

「Ｒｏｓｅｎｓｔｅｉｎ，Ｊ．Ｋ．，Ｗａｎｕｎｕａ， Ｍ．，Ｍｅｒｃｈａｎｔ，Ｃ．Ａ．，Ｄｒｎｄｉｃ，Ｍ．，ａｎｄ Ｓｈｅｐａｒｄ，Ｋ．Ｌ．： Ｉｎｔｅｇｒａｔｅｄ ｎａｎｏｐｏｒｅ ｓｅｎｓｉｎｇ ｐｌａｔｆｏｒｍ ｗｉｔｈ ｓｕｂ−ｍｉｃｒｏｓｅｃｏｎｄ ｔｅｍｐｏｒａｌ ｒｅｓｏｌｕｔｉｏｎ，Ｎａｔｕｒｅ Ｍｅｔｈｏｄｓ，ｐｐ．４８７−４９２（２０１２）」 「Ｗｅｋａ３：Ｄａｔａ Ｍｉｎｉｎｇ Ｓｏｆｔｗａｒｅ ｉｎ Ｊａｖa」、Ｍａｃｈｉｎｅ Ｌｅａｒｎｉｎｇ Ｇｒｏｕｐ ａｔ ｔｈｅ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｗａｉｋａｔｏ、インターネット<ＵＲＬ:ｈｔｔｐ：／／ｗｗｗ．ｃｓ．ｗａｉｋａｔｏ．ａｃ．ｎｚ／ｍｌ／ｗｅｋａ／> 「Ｅｌｋａｎ，Ｃ．ａｎｄ Ｎｏｔｏ，Ｋ．： Ｌｅａｒｎｉｎｇ Ｃｌａｓｓｉｆｉｅｒｓ ｆｒｏｍ Ｏｎｌｙ Ｐｏｓｉｔｉｖｅ ａｎｄ Ｕｎｌａｂｅｌｅｄ Ｄａｔａ，ｉｎ ＫＤＤ '０８ Ｐｒｏｃｅｅｄｉｎｇｓ ｏｆ ｔｈｅ １４ｔｈ ＡＣＭ ＳＩＧＫＤＤ ｉｎｔｅｒｎａｔｉｏｎａｌ ｃｏｎｆｅｒｅｎｃｅ ｏｎ Ｋｎｏｗｌｅｄｇｅ ｄｉｓｃｏｖｅｒｙ ａｎｄ ｄａｔａ ｍｉｎｉｎｇ，ｐｐ．２１３−２２０， Ｌａｓ Ｖｅｇａｓ，Ｎｅｖａｄａ，ＵＳＡ（２００８），ＡＣＭ Ｎｅｗ Ｙｏｒｋ，ＮＹ，ＵＳＡ」 「Ｔｓｕｔｓｕｉ，Ｍ．，Ｔａｎｉｇｕｃｈｉ，Ｍ．，Ｙｏｋｏｔａ，Ｋ．，ａｎｄ ＫａｗａｉＴ．： Ｉｄｅｎｔｉｆｙｉｎｇ Ｓｉｎｇｌｅ Ｎｕｃｌｅｏｔｉｄｅｓ ｂｙ Ｔｕｎｎｅｌｉｎｇ Ｃｕｒｒｅｎｔ，Ｎａｔｕｒｅ Ｎａｎｏｔｅｃｈｎｏｌｏｇｙ，Ｖｏｌ．５，ｐｐ．２８６−２９０（２０１０）」

しかしながら、上記の先端センシングデバイスの多くは、計測系や計測対象が微小であるが故に該対象の部分的な情報のみを出力するものであり、出力が熱雑音や量子ノイズなどの影響を受けることが多い。このため、対象信号よりもノイズ信号レベルの大きい場合、すなわち、ＳＮ比が非常に悪いという場合が多く、一次的な計測段階では計測精度が低すぎて実用化に適さないという問題を生じていた。また、このような計測状況下では、小さい信号がノイズで大きい信号が対象を表すと仮定して、信号強度でノイズ成分を除去するという一般的なノイズ除去手法を採用する余地はない。さらに、対象に関する知識や問題固有の性質を用いた各種ノイズフィルタリングを適用する場面も、係る知識や性質が明らかでない場合が多く適用しがたい。特に、１分子計測技術を用いた次々世代ＤＮＡシークエンサーでは、対象分子や計測系の性質に未知な部分が多くかつノイズ信号が大きいため、ノイズの影響が深刻な課題となっている。

本発明の目的は、計測データ集合から適切に不適合データを識別して、例えば、先端センシングデバイスによる計測結果の信頼性向上に寄与する識別方法、計測データに対する高精度の分類分析を行うことのできる分類分析方法、識別装置、分類分析装置、識別用記憶媒体および分類分析用記憶媒体を提供することである。

本発明は、上記課題に鑑み、正例集合と未知集合から分類器を学習する機械学習手法に着目し、例えば、正例／負例の２値分類に適したＰＵ分類手法（Ｃｌａｓｓｉｆｉｃａｉｏｎ ｏｆ Ｐｏｓｉｔｉｖｅ ａｎｄ Ｕｌａｂｅｌｅｄ Ｅｘａｍｐｌｅｓ）により構成した分類器を用いることによって、計測パターンから不適合データを高精度に識別し得るという知見にも基づいてなされた発明である。ＰＵ分類手法の詳細は、非特許文献３に記載されている。

本発明に係る第１の形態は、
計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された不適合データの識別をコンピュータ制御プログラムの実行によって行う識別方法であって、
前記コンピュータ制御プログラムは、正例集合の正例データと、正例負例のいずれかが不明である未知集合の未知データとから正負例を分類する分類器を学習する機械学習を用いた識別分析プログラムを有し、
前記計測空間に分析物を含まない試料を導入して計測する第１計測条件の下で得られるパルス状信号の第１種データと、前記計測空間に分析物を含む試料を導入して計測する第２計測条件の下で得られるパルス状信号の第２種データとを記憶する記憶手段を有し、
前記第１種データを前記正例データとし、前記第２種データを前記未知データとして、前記識別分析プログラムを実行することによって、前記第２種データに含まれる前記不適合データを識別することを特徴とする識別方法である。

本発明に係る第２の形態は、
第１の形態に係る識別方法により識別した不適合データを記憶する不適合データ記憶手段を有し、
計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された前記不適合データを取り除いた被分析データの分類分析をコンピュータ制御プログラムの実行によって行う分類分析方法であって、
前記コンピュータ制御プログラムは、機械学習を用いた分類分析を行う分類分析プログラムを有し、
前記パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、
あらかじめ求めた特徴量を前記機械学習のための学習データとし、前記不適合データを取り除いた被分析データのパルス状信号から得られる特徴量を変数にして、前記分類分析プログラムを実行することによって前記分析物に関する分類分析を行うことを特徴とする分類分析方法である。

本発明に係る第３の形態は、
前記特徴量は、
所定の時間幅内における波形の波高値、
パルス波長ｔ_a、
パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
該波形の鋭さを表す尖度、
パルス開始からパルスピークに至る傾きを表す俯角、
波形を所定の時間毎に区分した時間区分面積の総和を表す面積、
パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比、
パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値をベクトルの成分とする平均値ベクトル、
前記平均値ベクトルに対し波長が基準値になるように規格化した場合の規格化された平均値ベクトル、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメント、
のいずれか１または２以上である分類分析方法である。

本発明に係る第４の形態は、
計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された不適合データの識別をコンピュータ制御プログラムの実行によって行う識別装置であって、
前記コンピュータ制御プログラムは、正例集合の正例データと、正例負例のいずれかが不明である未知集合の未知データとから正負例を分類する分類器を学習する機械学習を用いた識別分析プログラムを有し、
前記計測空間に分析物を含まない試料を導入して計測する第１計測条件の下で得られるパルス状信号の第１種データと、前記計測空間に分析物を含む試料を導入して計測する第２計測条件の下で得られるパルス状信号の第２種データとを記憶する記憶手段を有し、
前記第１種データを前記正例データとし、前記第２種データを前記未知データとして、前記識別分析プログラムを実行することによって、前記第２種データに含まれる前記不適合データを識別することを特徴とする識別装置である。

本発明に係る第５の形態は、
第４の形態に係る識別装置により識別した不適合データを記憶する不適合データ記憶手段を有し、
計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された前記不適合データを取り除いた被分析データの分類分析をコンピュータ制御プログラムの実行によって行う分類分析装置であって、
前記コンピュータ制御プログラムは、機械学習を用いた分類分析を行う分類分析プログラムを有し、
前記パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、
あらかじめ求めた特徴量を前記機械学習のための学習データとし、前記不適合データを取り除いた被分析データのパルス状信号から得られる特徴量を変数にして、前記分類分析プログラムを実行することによって前記分析物に関する分類分析を行うことを特徴とする分類分析装置。

本発明に係る第６の形態は、
前記特徴量は、
所定の時間幅内における波形の波高値、
パルス波長ｔ_a、
パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
該波形の鋭さを表す尖度、
パルス開始からパルスピークに至る傾きを表す俯角、
波形を所定の時間毎に区分した時間区分面積の総和を表す面積、
パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比、
パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値をベクトルの成分とする平均値ベクトル、
前記平均値ベクトルに対し波長が基準値になるように規格化した場合の規格化された平均値ベクトル、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメント、
のいずれか１または２以上である分類分析装置である。

本発明に係る第７の形態は、第１の形態に係るコンピュータ制御プログラムを記憶したことを特徴とする識別用記憶媒体である。

本発明に係る第８の形態は、第２の形態コンピュータ制御プログラムを記憶したことを特徴とする分類分析用記憶媒体である。

第１の形態によれば、計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された不適合データの識別をコンピュータ制御プログラムの実行によって行う識別方法であって、前記コンピュータ制御プログラムは、正例集合の正例データと、正例負例のいずれかが不明である未知集合の未知データとから正負例を分類する分類器を学習する機械学習を用いた識別分析プログラムを有し、前記計測空間に分析物を含まない試料を導入して計測する第１計測条件の下で得られるパルス状信号の第１種データと、前記計測空間に分析物を含む試料を導入して計測する第２計測条件の下で得られるパルス状信号の第２種データとを記憶する記憶手段を有し、前記第１種データを前記正例データとし、前記第２種データを前記未知データとして、前記識別分析プログラムを実行することによって、前記第２種データに含まれる前記不適合データを識別することができる。したがって、本形態においては、ＰＵ分類手法に基づく分類器を構成して、計測の結果得られたパルス状信号に含まれる不適合データを高精度に識別でき、例えば、先端センシングデバイスによる計測結果の信頼性向上に寄与することができる。

特には、本形態における分類器は、対象に関する知識や問題固有の性質を用いることなく、過去に収集された不適合データ集合と、正負不明の実測データ集合の各データで構成することができるので、単純な信号強度で識別する従来手法では達成できない優れた不適合データの除去性能を具備し、種々の計測データの解析への広汎な適用可能性を有している。

第２の形態によれば、第１の形態に係る識別方法により識別した不適合データを記憶する不適合データ記憶手段を有し、計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された前記不適合データを取り除いた被分析データの分類分析をコンピュータ制御プログラムの実行によって行う分類分析方法であって、前記コンピュータ制御プログラムは、機械学習を用いた分類分析を行う分類分析プログラムを有し、前記パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、あらかじめ求めた特徴量を前記機械学習のための学習データとし、前記不適合データを取り除いた被分析データのパルス状信号から得られる特徴量を変数にして、前記分類分析プログラムを実行することによって前記分析物に関する分類分析を行うことができる。したがって、本形態においては、第１の形態に係るＰＵ分類手法に基づく分類器により高精度に識別した不適合データが取り除かれた被分析データにより前記分析物に関する分類分析を高精度に行うことができる。

第３の形態によれば、前掲の各特徴量は、パルス状信号の波形形態由来の特徴量であり、これらの特徴量群のいずれかのうち１または２以上の特徴量を使用することにより、機械学習による分類分析をより一層高精度に行うことができる。

本形態においては、上記特徴量群のうち少なくとも１つ以上の特徴量を使用して分類分析する場合に限らず、上記特徴量群のうち２つ以上の組み合わせを使用して分類分析を行うことができる。

第４の形態によれば、計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された不適合データの識別をコンピュータ制御プログラムの実行によって行う識別装置であって、前記コンピュータ制御プログラムは、正例集合の正例データと、正例負例のいずれかが不明である未知集合の未知データとから正負例を分類する分類器を学習する機械学習を用いた識別分析プログラムを有し、前記計測空間に分析物を含まない試料を導入して計測する第１計測条件の下で得られるパルス状信号の第１種データと、前記計測空間に分析物を含む試料を導入して計測する第２計測条件の下で得られるパルス状信号の第２種データとを記憶する記憶手段を有し、前記第１種データを前記正例データとし、前記第２種データを前記未知データとして、前記識別分析プログラムを実行することによって、前記第２種データに含まれる前記不適合データを識別することができる。したがって、本形態においては、ＰＵ分類手法に基づく分類器を構成して、計測の結果得られたパルス状信号に含まれる不適合データを高精度に識別でき、例えば、先端センシングデバイスによる計測結果の信頼性向上に寄与し得る識別装置を提供することができる。

特には、本形態に係る分類器は、対象に関する知識や問題固有の性質を用いることなく、過去に収集された不適合データ集合と、正負不明の実測データ集合の各データで構成することができるので、本形態は、単純な信号強度で識別する従来手法では達成できない優れた不適合データの除去性能を具備し、種々の計測データの解析への広汎な適用可能性を有する識別装置を実現することができる。

第５の形態によれば、第４の形態に係る識別装置により識別した不適合データを記憶する不適合データ記憶手段を有し、計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された前記不適合データを取り除いた被分析データの分類分析をコンピュータ制御プログラムの実行によって行う分類分析装置であって、前記コンピュータ制御プログラムは、機械学習を用いた分類分析を行う分類分析プログラムを有し、前記パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、あらかじめ求めた特徴量を前記機械学習のための学習データとし、前記不適合データを取り除いた被分析データのパルス状信号から得られる特徴量を変数にして、前記分類分析プログラムを実行することによって前記分析物に関する分類分析を行うことができる。したがって、本形態においては、第４の形態に係るＰＵ分類手法に基づく分類器により高精度に識別した不適合データが取り除かれた被分析データにより前記分析物に関する分類分析を高精度に行える分類分析装置を提供することができる。

第６の形態によれば、前掲の各特徴量は、パルス状信号の波形形態由来の特徴量であり、これらの特徴量群のいずれかのうち１または２以上の特徴量を使用することにより、機械学習による分類分析をより一層高精度に行える分類分析装置を提供することができる。

本形態においては、上記特徴量群のうち少なくとも１つ以上の特徴量を使用して分類分析する場合に限らず、上記特徴量群のうち２つ以上の組み合わせを使用して分類分析を行える分類分析装置を実現することができる。

第７の形態によれば、第１の形態に係るコンピュータ制御プログラムを記憶した識別用記憶媒体を提供することができる。したがって、本形態に係る記憶媒体は、第１の形態で説明した前記コンピュータ制御プログラムによる効果を有するので、識別用記憶媒体に記憶したコンピュータ制御プログラムをコンピュータにインストールして該コンピュータに識別分析動作させることによって高精度な識別分析を行うことができる。

第８の形態によれば、第２の形態に係るコンピュータ制御プログラムを記憶した分類分析用記憶媒体を提供することができる。したがって、本形態に係る記憶媒体は、第２の形態で説明した前記コンピュータ制御プログラムによる効果を有するので、分類分析用記憶媒体に記憶したコンピュータ制御プログラムをコンピュータにインストールして該コンピュータに分類分析動作させることによって高精度な分類分析を行うことができる。

第７および第８の形態に係る記憶媒体としては、フレキシブルディスク、磁気ディスク、光ディスク、ＣＤ、ＭＯ、ＤＶＤ、ハードディスク、モバイル端末等、コンピュータにより読み取り可能な記憶媒体のいずれかを選択することができる。

本発明によれば、コンピュータ端末を利用して、ＤＮＡ記憶媒体の情報圧縮技術や人工塩基対を用いた医薬品創薬、あるいは、計測試料に混入する微細な塵埃、あるいは体液などに含まれる分析物質を計測対象とする場合における、赤血球、白血球、血小板等の微小物質などに起因する不適合データの識別・除去技術等の分野におけるデータ分析を高精度に行うことができる。

本発明に係る実施形態における分析対象となる計測データを測定するための計測系を模式的に示す概要図、該計測系により計測したパルス状信号の波形例を示す図である。 前記計測系によりＤＮＡ構成分子について計測したパルス状信号の波形例を示す図である。 本発明の一実施形態である分類分析装置の概略構成を示す概略ブロック図である。 前記実施形態に用いるＰＣ１の識別・分類分析プログラムにより実行可能な処理内容の概要を示す図である。 ＰＣ１による識別処理を示すフローチャートである。 本発明による識別精度の検証に使用した２２種の分類器用ソフトウエアのリストを示す表である。 前記実施形態に用いる波高ベクトルを示す図である。 前記実施形態に用いる波長方向時間ベクトルを示す図である。 ＰＵ法の学習アルゴリズムの処理手順の概要を説明するための図である。 ＰＵ法における主要な解析内容を示す図である。 ＰＵ法による分類器の処理内容をまとめた概要説明図である。 前記識別処理におけるＰＵ法による２値分類器の識別処理を示すフローチャートである。 本発明に係る識別方法の識別精度を検証するための識別実験から得られたパルスピーク波高のヒストグラムを示す図である。 前記識別実験から得られたＦ−尺度（Ｆ−Ｍｅａｓｕｒｅのヒストグラムを示す図である。 マイクロ・ナノポアデバイスの概略構成を示す概略側断面図である。 ＰＣ１による分析処理の説明に必要な処理プログラム構成を示す図である。 実施例の大腸菌と枯草菌につき実測した粒子通過によるパルス波形例を示す図である。 本発明に係る各種特徴量を説明するためのパルス波形図である。 カルマンフィルターを説明するための図である。 カルマンフィルターの各因子を実際の計測電流データで説明するための図である。 カルマンフィルターの予測（８Ａ）と更新（８Ｂ）の繰返しの詳細を示す図である。 ＢＬ推定処理プログラムに基づくＢＬ推定処理を示すフローチャートである。 カルマンフィルターの因子調整に使用したビーズモデルの波形図である。 大腸菌２２と枯草菌２３が電解質溶液２４中に混在する様子を模式的に示した貫通孔１２周辺の拡大図である。 調整因子のｍ、ｋ、αの組合せに応じてビーズモデルの波形から拾われたパルスの数を示す表である。 特徴量抽出プログラムの実行処理内容の概要を示すフローチャートである。 粒子種推定処理を示すフローチャートである。 １つの波形データに関する各特徴量（１５Ａ）と、大腸菌と枯草菌の粒子種における確率密度関数のイメージ図（１５Ｂ）とを示す図である。 大腸菌と枯草菌の粒子種の個々より得られた確率密度分布の重ね合わせのイメージ図である。 ｋ個の粒子種別の粒子総数と、粒子種別の出現確率と、データ全体の出現頻度の期待値との関係を示すイメージ図である。 ラグランジュ未定乗数法により最適化を行う制約付き対数尤度最大化式の導出過程を説明するための図である。 データファイル作成処理を示すフローチャートである。 確率密度関数の推定処理を示すフローチャートである。 粒子数の推定処理を示すフローチャートである。 Ｈａｓｓｅｌｂｌａｄ反復法による粒子数推定処理を示すフローチャートである。 ＥＭアルゴリズムによる処理手順を示すフローチャートである。 本実施形態に係る個数分析機能により分析した結果の一例を示す図である。 特徴量としてパルス波長、波高を使用した検証例と、特徴量としてパルス波長、ピーク位置比を使用した検証例の各推定結果データを示す表である。 特徴量としてピーク付近波形の広がり、パルス波長を使用した検証例と、特徴量としてピーク付近波形の広がり、波高を使用した検証例の各推定結果データを示す表である。 特徴量として尖度と、パルス波高を使用した場合における個数推結果を示す図である。 ＢＬ推定処理プログラムに基づくＢＬ推定処理を示すフローチャートである。 大腸菌と枯草菌の混合比をそれぞれ、１：１０、２：１０、３：１０、３５：１００とした場合における各個数推定結果を示すヒストグラムである。 大腸菌と枯草菌の混合比をそれぞれ、４：１０、４５：１００、１：２とした場合における各個数推定結果を示すヒストグラムである。 特徴量としてパルス波長、パルス波高を使用した場合における各粒子の散布状態を合成した図である。 特徴量としてピーク付近波形の広がり、パルス波長を使用した場合、特徴量としてピーク付近波形の広がり、ピーク位置比を使用した場合、ピーク付近波形の広がり、パルス波高を使用した場合における各粒子の散布状態を合成した図である。 マイクロ・ナノポアデバイス８を用いて、３種の粒子３３ａ、３３ｂ、３３ｃが貫通孔１２を通過して得られる検出信号の波形例と、特徴量に基づいて得られる確率密度関数の導出例を示す図である。 俯角および面積の特徴量を説明するためのパルス波形図である。 波高ベクトルの取得の仕方を説明するための図である。 ｄ次元の波高ベクトルとデータサンプリングとの関係を説明するための図である。 時間（波長）および波幅に関する第２類型の特徴量を説明するためのパルス波形図である。 ｄ_w次元の波幅ベクトルとデータサンプリングとの関係を説明するための図である。 波幅に関する慣性モーメントを波幅ベクトルにより取得する取得過程を説明するための図である。 複数の方向に分割した場合の特徴量作成用波形ベクトルの一例を説明するための図である。 特徴量抽出の処理内容を示すフローチャートである。 １ＭＨｚ、５００ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である ２５０ｋＨｚ、１２５ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である。 ６３ｋＨｚ、３２ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である。 １６ｋＨｚ、８ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である。 ４ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である。 全サンプリングデータに対する各特徴量組合せに関する推定評価表である。 １ＭＨｚ〜１２５ｋＨｚでの高密度サンプリングしたときの各特徴量組合せに関する推定評価表である。 ６３ｋＨｚ〜４ｋＨｚでの低密度サンプリングしたときの各特徴量組合せに関する推定評価表である。 全サンプリングデータを使用したとき（５０Ａ）および高密度にサンプリングしたとき（５０Ｂ）に高い個数推定精度が得られる上位５種の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフである。 低密度にサンプリングしたときに高い個数推定精度が得られる上位５種の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフ（５１Ａ）と、全サンプリングデータを使用したときの４種類の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフ（５１Ｂ）である。 ４種類の各特徴量組合せに対する、特徴量作成に要する計算時間とＨａｓｓｅｌｂｌａｄ法による反復計算に要する計算時間との合計計算時間を示すサンプリング周波数（ｋＨｚ）−所要計算時間（秒）のグラフ（５２Ａ）と、各特徴量組合せに対する特徴量作成に要する計算時間を示すサンプリング周波数（ｋＨｚ）−所要計算時間（秒）のグラフ（５２Ｂ）である。 ４種類の各特徴量組合せに対する、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する計算時間を示すサンプリング周波数−所要計算時間（秒）のグラフである。 本発明に係る分類分析方法の概要を説明するための概要図である。 本実施形態における主な制御処理を示す図である。 本実施形態における分類分析処理を示すフローチャートである。 分類分析処理の検証により評価結果と、該検証における分析試料の詳細を示す表である。 Ｆ−尺度（Ｆ−Ｍｅａｓｕｒｅ）の説明図である。

本発明の一実施形態に係る分類分析装置を図面を参照して以下に説明する。本実施形態においては、分析物の一例としてＤＮＡ構成分子を分類分析する塩基種分析形態で説明する。

図１の（１Ａ）は、本実施形態における分析対象となる計測データを測定するための計測系を模式的に示す概要図である。

計測系は、塩基分子を含む溶液を収容する収容容器で構成された計測空間ＭＳと、計測空間ＭＳ内に対向して配置された1対の微細形状の電極Ｄ１、Ｄ２とを有する。電極Ｄ１、Ｄ２は、金（Ａｕ）元素で形成されたナノギャップ電極であり、互いに微細距離を隔てて配設されている。微細距離は、約１ｎｍに形成されている。計測空間ＭＳには、測定試料は、溶媒（純水）と、溶媒に混入されたＤＮＡ構成分子とを含む溶液サンプルである。

非特許文献４に記載されているように、ナノギャップ電極は、次々世代ＤＮＡシークエンサーとして期待されるデバイスである。この電極は、機械的破断接合と呼ばれる手法を用いて作成されたごく微細な隙間をもつ電極ギャップである。この電極ギャップに一定の電圧をかけると、ギャップ付近を物質が通過する際に量子力学的トンネル効果による電流（トンネル電流：図１の破線参照）が流れる。このトンネル電流が、物質が通過した瞬間のパルス電流として電流計測器ＭＥにより計測される。このナノ ギャップ電極によるトンネル電流パルスを計測することにより、ＤＮＡ塩基分子の種類を１分子単位で識別することが可能になり、既存技術では困難であったペプチドのアミノ酸配列や疾病マーカーとなる修飾アミノ分子の識別などが可能になってきている 。図１の測定系において、約１ｎｍの電極間隙を有するナノギャップ電極（Ｄ１、Ｄ２）を用いて、電極付近を通り過ぎる１分子に流れるトンネル電流パルスを計測して検出したパルス状信号のデータを分析対象とする。

被計測分子には、人工核酸塩基であるジチオフェンウラシル誘導体（以下、ＢｉｔｈｉｏＵと略す。）とＴＴＦウラシル誘導体（以下、ＴＴＦと略す。）の２種類を用いた。これらの分子は、識別を容易にするためにエピジェネティック部位（ＤＮＡメチル化などが起こる後天修飾部位）を化学的に修飾したものである。矢印Ｆで示すように、ギャップ付近をＤＮＡ分子を通過させる駆動力源は、分子自体のブラウン運動の他に、電気泳動、電気浸透流、誘電泳動によるものを使用することができる。

図１の（１Ｂ）および図２は、図１の計測系により計測したパルス状信号の波形例を示す。これらの図において、横軸は計測時間（×１０^-4ｓｅｃ）、縦軸は計測電流値（ｎＡ）を示す。

（１Ｂ）に示すように、パルス状信号のパルス判定部分は、計測波形中央の１／３の部分であり、このパルス波形データを被分析データに使用する。

（２Ａ）の２Ａ１、２Ａ２は、塩基分子ＢｉｔｈｉｏＵを検出したときの波形例を示す。（２Ｂ）の２Ｂ１、２Ｂ２は、塩基分子ＴＴＦを検出したときの波形例を示す。（２Ａ）の２Ａ３、２Ａ４および（２Ｂ）の２Ｂ３、２Ｂ４は、塩基分子を検出したときに混在するノイズ波形例を示す。

図１の計測系において、ＤＮＡの１塩基分子を電流パルスとして計測して検出する。計測されたパルスには、塩基分子由来のものだけではなく電極表面の金属原子のゆらぎや不純物による電流パルスも含まれている（図２の（２Ａ）の２Ａ３、２Ａ４および（２Ｂ）の２Ｂ３、２Ｂ４参照）。これらのノイズパルスのために、本来は塩基由来であるパルスを見逃したり、逆にノイズパルスであったのに塩基分子パルスが計測されたと誤判定する可能性が起こり得るので、計測結果としてＤＮＡ塩基分子の識別が困難になる。本発明は、計測されたパルスの波形データ集合から適切にノイズパルスの不適語データを識別、除去して、高精度に塩基種類の分類分析を可能にすることができる。

図３は本実施形態の分類分析装置の概略構成を示す。この分類分析装置は、パーソナルコンピュータ（以下、ＰＣという。）１により構成され、ＰＣ１にはＣＰＵ２、ＲＯＭ３、ＲＡＭ４およびデータファイル記憶部５を有する。ＲＯＭ３には、コンピュータ制御プログラムが格納されている。コンピュータ制御プログラムには、機械学習を用いた不適合データの識別処理および分類分析を行うための識別・分類分析プログラムおよび識別・分類分析に必要な特徴量の作成用プログラム等の各種処理プログラムが含まれている。分類分析プログラム等の各種処理プログラムは、各プログラムを記憶した記憶媒体（ＣＤ、ＤＶＤ等）からインストールされて格納可能になっている。ＰＣ１にはキーボード等の入力手段６および液晶ディスプレイ等の表示手段７が入出力可能に接続されている。データファイル記憶部５には分析用データが格納可能になっている。

ＰＣ１は、不適合データの識別処理機能および分類分析処理機能を具備するが、本発明においては、識別処理機能および分類分析処理機能をそれぞれ別個に備えた専用端末で構成することができる。

図４は、ＰＣ１の識別・分類分析プログラム（コンピュータ制御プログラム）により実行可能な処理内容の概要を示す。図５は、ＰＣ１の識別処理のフローチャートである。

ＰＣ１の識別処理は、以下の識別方法に基づく処理手順により行われる。識別・分類分析プログラムには、正例集合の正例データと、、正例負例のいずれかが不明である未知集合の未知データとから正負例を分類する分類器を学習する機械学習を用いたＰＵ手法に基づいた識別処理プログラムが格納されている。
（処理１−１）計測空間ＭＳに分析物（ＤＮＡ構成分子）を含まない試料（溶媒のみ）を導入して計測する第１計測条件の下で得られるパルス状信号の第１種データを記憶手段のＲＡＭ４に取り込んで記憶させる。
（処理１−２）計測空間ＭＳに分析物（ＤＮＡ構成分子）を含む試料（溶媒＋ＤＮＡ構成分子）を導入して計測する第２計測条件の下で得られるパルス状信号の第２種データを記憶手段のＲＡＭ４に取り込んで記憶させる。
（処理１−３）識別処理プログラムの入力形式に合わせるために、第１種データおよび第２種データの属性ベクトルが作成される。
（処理１−４）第１種データを正例データとし、第２種データを未知データとして、識別処理プログラムを実行する。
（処理１−５）識別処理プログラムの実行により、確率ｐ（ｓ＝１｜ｘ）を抽出して求める。該確率データは、ＲＡＭ４の所定エリアに記憶、保存される。なお、以下のＰＵ法による解析において使用される属性ベクトルは、多次元データによるもので、ベクトル表記されるものであるが、以下の説明では、特にベクトル表記を省略している。
（処理１−６）該確率により、第２種データに含まれている、分析物以外の要素（前述の塩基分子由来のものでなく電極表面の金属原子のゆらぎや不純物）に起因して検出された不適合データを検出、識別する。検出した不適合データは、ＲＡＭ４の所定エリアに記憶、保存される。

識別処理プログラムには、非特許文献２に開示されている、機械学習プラットフォームフリーウェア Ｗｅｋａの分類器用ソフトウエアを用いることができる。

図６は、本発明による識別精度の検証に使用した２２種の分類器用ソフトウエアのリストである。識別処理プログラムとして、２２種のいずれも使用可能であり、ＲＯＭ３に格納可能になっている。ＰＵ法におけるｐ（ｓ＝１｜ｘ）の計算と、ＰＵ法によるノイズデータ除去後の識別処理のいずれに対してもＷｅｋａのプログラムを使用して検証を行った。

計測したパルス波形データは、波長も波高もまちまちであるため、機械学習分類器によって塩基種類を識別するためには次元のそろった属性ベクトルを入力として用いる必要があり、機械学習処理プログラムの実行の際には、入力形式に合わせる前処理として、一種の粗視化を施し、パルス波形を反映した属性ベクトルを作成する前処理が（処理１−１）および（処理１−２）において行われる。

計測したパルス波形データは、波長も波高もまちまちであるため、機械学習分類器によって塩基種類を識別するためには次元のそろった属性ベクトルを入力として用いる必要があり、識別処理プログラムの実行の際には、入力形式に合わせる前処理として、一種の粗視化を施し、パルス波形を反映した属性ベクトルを作成する前処理が（処理１−３）において行われる。

図７は、波高ベクトルを示す。図８は、波長方向時間ベクトルを示す。

図７に示すように、計測パルス波形について、波長方向にｄ_h分割し各分割区分ごとに計測電流値の平均値を計算して、これを成分にしたｄ_h次元の属性ベクトルを波高ベクトルとする。この属性ベクトルには、波高方向に規格化したものとしないものの２種類が作成される。

図８に示すように、パルスのピーク前後で計測電流値を２つのグループに分けた上で、波高方向にｄ_w 分割すると、パルスの計測電流値は、２ｄ_wのグループに分割される。この分割区分ごとにパルス開始時点からのステップ数の平均値を 算出して、これらの値を成分としてもつ２ｄ_w次元の波長方向時間ベクトルが作成される。また、パルス開始時点から終了時点までの時間を「 1 」とする規格化を施した規格化波長方向時間ベクトルも作成される。以上の波高ベクトルと波長方向時間ベクトルに加え、これらを単に連結した属性ベクトルも作成される。これらのベクトルデータは、ＲＡＭ４の所定エリアに記憶される。

本実施形態において２値分類器を構成するために、波高と波長の２つの特徴量を使用している。本実施形態に係る不適合データの識別精度を検証する検証実験においては、１つのパルス波形データから作成した、下記のＶ１〜Ｖ８の８通りの属性ベクトルを使用して検証した。

（Ｖ１）パルスピーク値を「１」に規格化した波高ベクトル（ｈｖＮｒｍｄ）
（Ｖ２）規格化しない波高ベクトル（ｈｖＲａｗ）
（Ｖ３）パルス波長時間を「１」に規格化した波長方向時間ベクトル（ｗｖＮｒｍｄ）
（Ｖ４）規格化しない波長方向時間ベクトル（ｗｖＲａｗ）
（Ｖ５）Ｖ１とＶ２を連結した（ｄ_h＋２ｄ_w）次元ベクトル
（Ｖ６）Ｖ１とＶ４を連結した（ｄ_h＋２ｄ_w）次元ベクトル
（Ｖ７）Ｖ２とＶ３を連結した（ｄ_h＋２ｄ_w）次元ベクトル
（Ｖ８）Ｖ２とＶ４を連結した（ｄ_h＋２ｄ_w）次元ベクトル

検証実験では、上記８通りの属性ベクトルを作成し、これらの識別精度の比較を行った。属性ベクトル作成時の分割数は予備解析を行った上で、一律にｄ_h ＝１０、ｄ_w＝５とした。

通常の２値分類器の場合、正例と負例が与えられたデータから学習して分類器が生成される。これに対し、本実施形態においては、計測データに不適合データが混在する場合であるから、ＰＵ法による分類器を使用している。本実施形態に使用されるＰＵ法は、非特許文献３に詳述されているように、正例とラベルなしデータから学習し、正例／負例の２値分類をするための半教師あり学習アルゴリズムの一種である。ＲＯＭ３に格納したＰＵ法の学習アルゴリズムの処理手順の概要は以下の通りである。

図９は、ＰＵ法の学習アルゴリズムの処理手順の概要を説明するための図である。同図（９Ａ）は、学習に使用する変数およびラベルフラグを示し、（９Ｂ）は、（９Ａ）の前提条件の詳細を示す。図１０は、ＰＵ法における主要な解析内容を示す図である。図１１は、以下に説明する、ＰＵ法による分類器の処理内容をまとめた概要説明図である。図１１において、正例集合、負例集合をそれぞれ、Ｐ、Ｎとし、Ｐには、ラベル付き部分集合Ｌとラベル無し部分集合Ｕが含まれ、ＮにはＵのみを含むとしている。

事例ｘ（入力データ）をパルス波形に関する属性ベクトルとし、ｙをそのクラスラベル、事例にクラスラベルが付けられているか否かを示すフラグをｓとする。入力事例の集合のなかで、正例（ｙ＝１）の一部のみがラベルされており（ｓ＝ １）、他の正例と全ての負例（ｙ＝０）はラベルされていない（ｓ＝０）。すなわち、サンプルが負例であるならばラベルされている確率はゼロであり、p(ｓ＝ １｜ｘ，ｙ＝０)＝０である。このような事例集合を２値分類器の学習アルゴリズムの入力とし、サンプルがラベルされている確率 ｇ（ｘ）＝ｐ（ｓ＝１｜ｘ）を求めることができる。さらに、本来求めたいものはｇ（ｘ）ではなくｐ（ｙ＝１｜ｘ）であるから、次の補正を加えてｐ（ｙ＝１｜ｘ）が抽出される。

全事例集合において、サンプルがラベルされている確率である、ｇ（ｘ）＝ｐ（ｓ＝１｜ｘ）は、図１０の（１０ａ）に示す導出過程により、ｇ（ｘ）＝ｐ（ｙ＝１｜ｘ）ｐ（ｓ＝１｜ｙ＝１）の関係式に導出される。ｃ＝ｐ（ｓ＝１｜ｙ＝１）とすると、サンプルが正例である確率は、ｐ（ｙ＝１｜ｘ）＝ｇ（ｘ）／ｃと与えられる。

ここで、正事例集合中でラベル付けされる確率が一様ランダム、すなわちｘによらずｐ（ｓ＝１｜ｙ＝１，ｘ） ＝ｐ（ｓ＝１｜ｙ＝１）＝ｃは一定値であると仮定している。これは、取り扱う計測データが意図的に偏った恣意的なデータではないことによる。

ここで、ｃは、次のようにして推定することが可能である。正事例集合中で一 様ランダムにラベル付けされているならば、ｇ（ｘ）は、ｘが正例である場合には正例中に含まれるラベル付き事例集合の割合 に一致し、ｇ（ｘ）＝ｐ（ｓ＝１｜ｙ＝１）＝ｃとなる。そこで、ＰＵ法によらない通常の２値分類器で求めたｇ（ｘ）を用いて、正事例であるラベル付き事例集合Ｌ中の平均（下式の数１としてｃを推定することができる。

図１２は、（処理１−５）におけるＰＵ法による２値分類器の識別処理を示す。

処理Ｐ１−１において、ｇ（ｘ）を学習データ集合により学習する処理が行われる。次に、処理Ｐ１−２において、数１に基づいてｃを検証用データ集合により推定する処理が行われる。処理Ｐ１−３において、ｇ（ｙ＝１｜ｘ）＝ｇ（ｘ）／ｃの関係から確定したｇ（ｙ＝１｜ｘ）によってテストデータに対する正例／負例の識別を行う処理が行われる。この場合の判断基準は、ｇ（ｙ＝１｜ｘ）＞０．５とすることができる。

本発明は、図１１に示すＰＵ法による分類器の構成において、ラベルなし事例が正例である確率を抽出することができる。以下に、ラベルなし事例が正例である確率を抽出する抽出手順を説明する。

ラベル付き事例は全て正例であるが、ラベルなし事例は正例、負例のいずれの可能性もある。ラベルなし事例が正例である確率をｗ（ｘ）とすると、その負例である確率は、 １−ｗ（ｘ）である。そこで、ラベルなし事例をすべて２倍に複製し、一方を正例として扱い、もう一方を負例として扱う。正例として扱うラベルなし事例ｘには重みｗ（ｘ）を与え、負例として扱うラベルなし事例ｘ には重み「１−ｗ（ｘ）」を与える。ラベル付き事例は、すべて正例であるから、重み「１」で正例として扱う。これら重み付き事例集合を学習データとして分類器を作成する。

ここで、ｃとｇ（ｘ）＝ｐ（ｓ＝１｜ｘ）は、図９〜図１１に示した手法により得られているとした場合、ラベルなし事例が正例である確率ｗ（ｘ）は、図１０の（１０ｂ）に示す導出過程により、ｗ（ｘ）＝（１−ｃ）ｇ（ｘ）／（ｃ（１−ｇ（ｘ）））となるので、ｃおよびｇ（ｘ）の抽出によりラベルなし事例が正例である確率ｗ（ｘ）を求めることができる。

本実施形態による識別精度の検証実験を以下に説明する。

ＤＮＡ構成分子を分析物とし、ナノギャップ電極を用いて計測した計測パルス集合から、１）まず前処理としてＰＵ法による分類器を構成し、ノイズ由来のパルス（不適合データ）を識別し（図１２参照）、抽出した不適合データを第２種データから除去して塩基由来のパルスのみのデータ集合を取得した。２）そのようにして得た塩基由来のパルス集合に対して塩基種別の識別精度を評価した。

ノイズ除去は、あらかじめ塩基（ＢｉｔｈｉｏＵ、ＴＴＦ）を含んでいない溶 媒のみに対して、ナノギャップ電極により計測したトンネル電流パルスを取得しておく。このパルス集合は塩基とは関係のないノイズ由来のパルスであり、「ノイズパルス集合」と呼ぶことにする。次に、溶媒に塩基ＢｉｔｈｉｏＵを混入したものについて計測した電流パルスを取得する。ＴＴＦについても同様に取得する。このパルス集合には、塩基由来の「塩基パルス」とノイズパルスの双方が含まれている。そこでこれを「塩基＋ノイズパルス集合」と呼ぶことにする。

ノイズパルス集合中のパルスは必ずノイズパルスであるので、それを正事例集合（第１データのデータ集合）とみなし、塩基＋ノイズパルス集合中のパルスは、いずれのパルスであるか不明なので、それをラベルなし事例集合とみなして、図１２に示したＰＵ分類器処理によって、ノイズパルス（正例）と塩基パルス（負例）の識別を行うことができ、正例であるノイズデータを除去することにより、ほぼ塩基パルスのみからなる集合（塩基パルス集合）を得ることができる。

ＰＵ分類器処理は、このノイズパルスと塩基パルスの正負例分類のために 1 度使用するだけであり過学習による問題は起こらないので、全パルス集合を学習用データとして用いてＰＵ分類器を作成し、それにより全パルス集合をノイズパルスと塩基パルスに分離する分類分析を行った。このようにして、ＢｉｔｈｉｏＵの塩基＋ノイズパ ルス集合からＰＵ分類器によりＢｉｔｈｉｏＵの塩基パルス集合を取得、ＴＴＦの塩基＋ノイズパルス集合からＰＵ分類器によりＴＴＦの塩基パルス集合を取得した。

塩基パルスと塩基＋ノイズパルスの識別精度評価のために、ノイズデータを除去、分離したＢｉｔｈｉｏＵとＴＴＦの塩基パルス集合に対し、通常の２値分類器による塩基種類の識別実験を行った。識別実験では、２種塩基の塩基パルス数のいずれかが１０に満たない場合は、学習用事例が少なすぎるために実験対象から除外した。識別精度の指標にはＦ−ｍｅａｓｕｒｅ（後述の図７１に示す、Ｆ−尺度）を用い、１０倍交差検定（１０−Ｆｏｌｄ ｃｒｏｓｓｖａｌｉｄａｔｉｏｎ：以下、１０ＣＶと略す。）により精度評価を行った。１０ＣＶの際には、ＢｉｔｈｉｏＵとＴＴＦの塩基パルス数は同数とした．すなわち、ＰＵ分類器により得たＢｉｔｈｉｏＵとＴＴＦの塩基パルス数がそれぞれＮＢ、ＮＴであるとき、１０ＣＶに用いる塩基パルス数をＢｉｔｈｉｏＵ、ＴＴＦともにＮ＝ｍｉｎ（ＮＢ，ＮＴ） に揃えた。パルス数がＮより大きい塩基パルス集合については、Ｎ個の塩基パルスをランダム抽出した。また、ＰＵ分類器によるノイズ除去の効果を見るために、ノイズデータ除去を施す以前の塩基＋ノイズパルス集合に対しても、ＢｉｔｈｉｏＵとＴＴＦの識別実験を行った。ＢｉｔｈｉｏＵ、ＴＴＦそれぞれの塩基＋ノイズパルス集合からＮ個ずつランダム抽出して得たパルス集合に対し、塩基パルス同様に １０ＣＶで精度評価を行った。

以下に検証実験の実験条件を説明する。

図６に示した機械学習ソフトウエアによる、２２種の分類器で種々の分析条件の下で識別精度を調べた。

パルス抽出パラメータとして、パルス抽出の際には、計測電流値のベースラインからどれだけ外れたらパルス開始と判定するかという波高閾値αと、何ステップ以上波高閾値を超えたらパルスであると判定するかという波長閾ｋ値の２つのパラメータ（後述の図２２に示す調整因子）を用いている。これらのパラメータを様々試し、「波高閾値αについて４通り×波長閾ｋ値について４通りの計１６通り」に対して実験を行った。属性ベクトルとしては、Ｖ１〜Ｖ８の８種類の属性ベクトルについて試した。

分類分析用の分類器としてアンサンブル学習「Ｒｏｔａｔｉｏｎ Ｆｏｒｅｓｔ」を採用し、その内部で用いるベース分類器として、Ｗｅｋａに実装されているもののうち、入力事例連続値ベクトルの２値分類を行える、図６の２２種類の分類器を使用した。ＰＵ手法としては、図１０に示したｇ（ｘ）およびｗ（ｘ）の２種の手法を使用した。

上記の実験条件下で行う不適合データの識別実験は、上記の実験条件下で行う不適合データの識別実験は、パルス抽出パラメータ１６通り×属性ベクトル８通り×Ｗｅｋａに実装された２２分類器×ＰＵ手法２通りの全ての組合せのうちで、塩基パルス数が２塩基とも１０以上であった３２７２ケースについて行った。この識別実験では、単純化のため、抽出したパ ルスに対して、１）ＢｉｔｈｉｏＵのノイズ除去、２）ＴＴＦのノイズ除去、３）ノイズ除去後の２塩基識別のこれら３者に用いた条件（パルス抽出パラメータ、属性ベクトル、分類器、ＰＵ分類器手法）は全て共通とした。同様な条件でＰＵ分類器によるノイズ除去を用いないで塩基＋ノイズパ ルス集合に対しても識別実験を行った。

図１３は、３２７２ケースから得られた、Ｆ−ｍｅａｓｕｒｅ＞０．９であった１事例（Ｆ−ｍｅａｓｕｒｅが０．９３であった解析条件で使用した計測パルス集合）で、塩基と判定したパルスと、ノイズと判定したパルスについてのパルスピーク波高のヒストグラムを示す。横軸は、パルスピーク波高（ｎＡ）、縦軸は、パルス数を示す。（１３Ａ）は、ＢｉｔｈｉｏＵに関するノイズパルスと塩基パルスのヒストグラムを示し、（１３Ｂ）は、ＴＴＦに関するノイズパルスと塩基パルスのヒストグラムを示す。（１３Ａ）において、ノイズパルスと塩基パルスとでは、それぞれ、０〜０．３、０．０２〜０．４の波高範囲で分布している。（１３Ｂ）において、ノイズパルスと塩基パルスとでは、それぞれ、０〜０．２、０〜１．２の波高範囲で分布している。

図１３からわかるように、パルスピーク波高のヒストグラムは、ノイズパルスと塩基パルス間で重なり合う部分が多く、パルスピーク波高だけではノイズパルスと塩基パルスとの識別が困難であることがわかる。

図１４は、ノイズ除去あり／なしのそれぞれ３２７２ケースについて得られた、Ｆ−尺度（Ｆ−Ｍｅａｓｕｒｅ）のヒストグラムを示す。横軸は塩基識別精度、縦軸は、ノイズ除去あり／なし別の各種条件下における解析事例数を示す。ノイズ除去なしの場合と、ノ
イズ除去ありの場合とでは、それぞれ、０．３〜０．６、０．５〜１．０の精度範囲で分布している。解析事例総数は、パルス抽出パラメータと属性ベクトルと分類器の各種組合せの３２７２である。

図１４からわかるように、ノイズ除去なしの場合と、ノイズ除去ありの場合とでは、重なり部分が大きいものの、１００％ないし１００％に近い精度まで識別精度が向上している。したがって、ＰＣ１による不適合データの識別処理性能は、パルスピーク波高だけでは塩基／ノイズの判定が困難な場合であっても、パルス波形特徴を的確に把握した特徴量の属性ベクトルを用いることにより、適切にノイズパルスを除去して高い塩基分類精度を得ることができる。

ＰＣ１による分類分析装置は、上記の識別処理により識別された不適合データを除去した被分析データに対する高精度の分類分析機能を有する。この分類分析機能は、以下の分析手順によって構成されている。

（Ｃ１）上記の識別処理により識別した、分析物以外の要素に起因して検出された不適合データを不適合データ記憶手段のＲＡＭ４の所定エリアに記憶する。不適合データは、ＰＣ１で検出して記憶するだけでなく、あらかじめ外部端末に記憶した不適合データファイルをＰＣ１に導入して記憶させてもよい。

（Ｃ２）計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、記憶した不適合データを取り除いたデータ群を被分析データとしてＲＡＭ４の所定エリアに記憶する。あらかじめ不適合データを取り除いたデータ群を被分析データとしてＰＣ１に導入して記憶させてもよい。

（Ｃ３）ＰＣ１に搭載されたコンピュータ制御プログラムには、機械学習を用いた分類分析を行う分類分析プログラムが含まれ、ＲＯＭ３に格納されている。

（Ｃ４）パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、あらかじめ求めた特徴量を機械学習のための学習データとし、不適合データを取り除いた被分析データのパルス状信号から得られる特徴量を変数にして、分類分析プログラムを実行することによって分析物に関する分類分析が実行可能になっている。

ＰＣ１による分類分析装置によれば、パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、あらかじめ求めた特徴量を機械学習のための学習データとし、ＰＵ分類手法に基づく分類器により高精度に識別した不適合データが取り除かれた被分析データから得られる特徴量を変数にして、分類分析プログラムを実行することによって分析物に関する分類分析を行うので、該被分析データにより分析物に関する分類分析を高精度に行うことができる。

本実施形態においては、高精度の不適合データの識別と分類分析を行えるので、例えば、人工塩基の識別可能となる途を開き、ＤＮＡ記憶媒体の情報圧縮技術や、人工塩基対を用いた医薬品創薬などへの応用展開を図ることができる。

本発明は、上記の実施形態における電流出力波形に限らず、広範囲の出力波形、例えば、電圧波形、インピーダンス波形などに対する不適合データの識別および分類分析に適用することができる。

本発明は、ナノギャップ電極による測定系による検出波形に限らず、試料対象が通過する試料対象相当の微小構造、例えば、貫通孔、ウェル（凹型），ピラー（凸型）、流路等の測定系による検出波形に適用することができる。本発明における計測データの適応範囲は、時系列の計測データ中のすべてであって，電気計測に限らず光計測、音等の物理現象の検出データを含むことが可能である。

本発明における分析物以外の要素に起因する除去対象には、上記の実施形態における量子レベルの要素に限らず、例えば、計測器、計測デバイス、溶液中に存在する分析物以外で混入しているシグナルに適用可能である。すなわち、本発明は、例えば、計測試料に混入する微細な塵埃、あるいは血液などに含まれる分析物質を計測対象とする場合における、赤血球、白血球、血小板等の微小物質などに起因する不適合データの識別・除去技術に適用することができる。

本発明に係る識別処理および分類分析に使用される特徴量は、上記の波高・波長に限らず、波形形態由来の各種特徴量を使用することができる。本発明者ら、マイクロ・ナノポアデバイスを用いた粒子検出技術で登場する波形データの解析から有効な特徴量を把握するに至った。マイクロ・ナノポアデバイスを用いた粒子検出技術は、特許文献１などに開示されている。以下、本発明に有効な特徴量およびそれを使用したときの分類分析処理を詳述する。

図１５は、マイクロ・ナノポアデバイス８を用いた粒子検出装置の概略構成を示す。

粒子検出装置は、マイクロ・ナノポアデバイス８およびイオン電流検出部により構成されている。マイクロ・ナノポアデバイス８は、チャンバー９と、チャンバー９を上下の収容空間に区画する隔壁１１と、隔壁１１の表裏側に配置された一対の電極１３、１４とを有する。隔壁１１は基板１０上に形成されている。隔壁１１の中央付近には、微小径の貫通孔１２が穿設されている。貫通孔１２の下方には、基板１０の一部を下向きに凹状に取り除いた凹部１８を設けている。マイクロ・ナノポアデバイス８は半導体デバイス等の製造技術（例えば、電子線描画法やフォトリソグラフィ）を用いて作成される。すなわち、基板１０はＳｉ材で構成され、表面上にＳｉ₃Ｎ₄膜による隔壁１１が薄膜形成されている。凹部１８は基板１０の一部をエッチングにより除去して形成されている。

隔壁１１は、大きさ１０ｍｍ角で、厚さ０．６ｍｍのＳｉ基板に５０ｎｍのＳｉＮ膜を積層して形成されている。Ｓｉ₃Ｎ₄膜にレジストを塗布し、電子線描画法により、直径３μｍの円形の開口パターンを形成して、貫通孔１２が穿設されている。貫通孔１２の裏側においては、ＫＯＨによるウェットエッチングを行い５０μｍ角の開口を形成して凹部１８を設けている。凹部１８の形成は、ウェットエッチングに限らず、例えば、ＣＦ₄系ガスによるドライエッチングなどによる等方的なエッチング等により行うことができる。

隔壁１１用の膜には、ＳｉＮ膜の他に、ＳｉＯ₂膜、Ａｌ₂Ｏ₃膜、ガラス、サファイア、セラミック、樹脂、ゴム、エラストマー等の絶縁性膜を使用することができる。基板１０の基板材料もＳｉに限定されるものではなく、ガラス、サファイア、セラミック、樹脂、ゴム、エラストマー、ＳｉＯ₂、ＳｉＮ、Ａｌ₂Ｏ₃等を使用することができる。

貫通孔１２は、上記の基板上の薄膜に形成する場合に限らず、例えば、貫通孔１２を形成した薄膜状シートを基板上に接着することによって、貫通孔を有する隔壁を形成するようにしてもよい。

イオン電流検出部は、電極１３、１４の電極対と、電源１５、増幅器１６および電圧計２０により構成されている。電極１３、１４は貫通孔１２を介して対向配置されている。増幅器１６はオペアンプ１７と帰還抵抗１９とにより構成されている。オペアンプ１７の（−）入力端子と電極１３は接続されている。オペアンプ１７の（＋）入力端子は接地されている。オペアンプ１７の出力側と電源１５の間に電圧計２０が接続配置されている。電源１５により電極１３、１４間には、０．０５〜１Ｖの印加電圧を使用できるが、本実施例では０．０５Ｖが印加されるようになっている。増幅器１６は電極間に流れる電流を増幅して電圧計２０側に出力する。電極１３、１４の電極材料としては、例えば、Ａｇ／ＡｇＣｌ電極、Ｐｔ電極、Ａｕ電極等を使用でき、好ましくはＡｇ／ＡｇＣｌ電極である。

チャンバー９は、マイクロ・ナノポアデバイス８周囲を密閉状に囲む流動性物質収容容器であり、電気的および化学的に不活性な材質、例えば、ガラス、サファイア、セラミック、樹脂、ゴム、エラストマー、ＳｉＯ₂、ＳｉＮ、Ａｌ₂Ｏ₃等により形成されてよい。

チャンバー９内には注入口（図示せず）から検体２１を含む電解質溶液２４が充填される。検体２１は、例えば、細菌、微小粒子状物質、分子状物質等の分析物である。検体２１を、流動性物質である電解質溶液２４に混入して、マイクロ・ナノポアデバイス８による計測が行われる。イオン電流検出部による計測終了時には排出口（図示せず）より充填溶液は排出可能になっている。電解質溶液には、例えば、リン酸緩衝生理食塩水（ＰＢＳ）、Ｔｒｉｓ−ＥＤＴＡ（ＴＥ）バッファーやそれらの希釈液の他、これらと同様なすべての電解質溶液剤を使用することができる。計測は、検体含有電解質溶液をチャンバー９内に導入、充填するごとに行う場合に限らず、溶液溜から簡易ポンプ装置により検体含有電解質溶液（流動性物質）を汲み出して注入口よりチャンバー９内に充填し、計測後に排出口から排出させ、さらに、別の溶液溜あるいは新たな溶液を溶液溜に貯留し、新たに汲み出して次の計測を行う連続計測システムを構成するようにしてもよい。

電解質溶液２４をチャンバー９内に充填した状態で、貫通孔１２の上下の電極１３、１４間に電源１５の電圧印加を行うと、貫通孔１２に比例した一定のイオン電流が電極間に流れる。電解質溶液２４中の細菌等の検体が貫通孔１２を通過する際には、一部のイオン電流が検体により阻害されるため、電圧計２０によりパルス状のイオン電流減少を計測することができる。したがって、マイクロ・ナノポアデバイス８を用いた粒子検出装置によれば、計測電流の波形変化を検出することにより、検体（例えば、粒子）毎の貫通孔１２通過による流動性物質中に含まれる粒子個々の存在を高精度に検出することができる。計測態様には、流動性物質を強制的に流動させながら計測する場合に限らず、流動性物質を非強制的に流動させながら計測する場合を含むことができる。

電圧計２０によるイオン電流の計測出力は外部出力可能になっている。この外部出力は、変換回路装置（図示せず）によりデジタル信号データ（計測電流データ）に変換されて記憶装置（図示せず）に一旦保存された後、データファイル記憶部５に格納される。データファイル記憶部５には、マイクロ・ナノポアデバイス８を用いた粒子検出装置によりあらかじめ取得した計測電流データを外部入力することができる。

図６８は、ＰＣ１による分析物（例えば、大腸菌Ｅｃや枯草菌Ｂｓ）に対する分類分析処理の概要を説明するための概要図を示す。

図６８の分類分析処理は、以下の分析手順（ａ）〜（ｄ）によって構成されている。
（ａ）所定の分析物（例えば、大腸菌Ｅｃや枯草菌Ｂｓ）を含む流動性物質に対しナノポアデバイス８ａによる計測の結果、各種別毎の検出信号として得られた貫通孔８ｂの分析物通過に対応するパルス状信号Ｄｅ、Ｄｂの波形形態の特徴を示す特徴量をあらかじめ求める。パルス状信号Ｄｅ、Ｄｂは、それぞれ大腸菌Ｅｃ、枯草菌Ｂｓの貫通孔８ｂ通過によって得られた信号である。
（ｂ）コンピュータ解析部１ａには、機械学習による分類分析を行う分類分析プログラムが内蔵されている。（ａ）においてあらかじめ求めた特徴量は、大腸菌Ｅｃ、枯草菌Ｂｓの既知データから得られた特徴量であり、機械学習のための学習データとしてコンピュータ解析部１ａにおいて使用される。
（ｃ）例えば、大腸菌Ｅｃおよび枯草菌Ｂｓの含有比ないし含有数が不明の状態で流動性物質中に混入された混合物を被分類分析物Ｍｂとした場合、（ａ）の既知データ取得の場合と同様に、ナノポアデバイス８ｃによる計測を行う。この計測により、被分類分析物Ｍｂの貫通孔８ｄ通過によって被分析データとしてパルス状信号Ｄｍが得られる。
（ｄ）既知データによる特徴量を学習データとし、被分析データのパルス状信号Ｄｍから得られる特徴量を変数にして、分類分析プログラムを実行することによって、該被分析データにおける所定の分析物に関する分類分析を行うことができる。

上記の分類分析により、特徴量に基づいて機械学習による分類分析を行って、種別の分からない被分析データを大腸菌Ｅｃまたは枯草菌Ｂｓの通過に由来するもの１ｂと由来しないものに分類することができる。なお、本発明に係る特徴量は、コンピュータ解析部１ａにおいて作成してもよいし、別の特徴量作成プログラムを使用して作成してコンピュータ解析部１ａに与えるようにしてもよい。

図６９は、ＰＣ１による主な制御処理を示す。

主な制御処理には、入力処理（ステップＳ１００）、入力データから特徴量を取得する特徴量取得処理（ステップＳ１０１）、分類分析処理（ステップＳ１０４）、個数分析処理（ステップＳ１０５）および出力処理（ステップＳ１０６）が含まれている。入力処理（ステップＳ１００）において、ＰＣ操作に必要な各種入力、内蔵プログラムの起動入力、各種分析の実行指示入力、計測電流データおよび／または特徴量データの入力、出力態様の設定入力、分析時に特徴量を指定する場合の指定特徴量の入力等が行われる。この入力処理には、不適合データの除去処理も含まれている。入力手段６による分析種別の指定操作を行うことによって、分類分析処理（ステップＳ１０４）または個数分析処理（ステップＳ１０５）が実行可能になっている（ステップＳ１０２、Ｓ１０３）。分類分析処理は、特徴量取得処理（ステップＳ１０１）において入力データから取得した特徴量のベクトル量データを使用して分類分析可能になっている。個数分析処理は、特徴量取得処理において入力データから取得した特徴量のスカラーデータを使用して個数分析可能になっている。本実施形態は、分類分析処理機能に加え、個数分析処理機能を具備する実施形態であるが、本発明は、分類分析処理機能のみを具備した実施形態により実施することができる。

本実施形態に係るコンピュータ制御プログラムには、粒子種別の個数ないし個数分布を分析するための個数分析プログラムが含まれている。個数分析処理（ステップＳ１０５）において、個数分析プログラムの実行が可能になっている。出力処理（ステップＳ１０６）において、分類分析処理（ステップＳ１０４）および個数分析処理（ステップＳ１０５）における分析結果データの出力が可能であり、例えば、表示手段７に各種分析結果データが表示出力される。ＰＣ１に出力手段としてプリンタ（図示せず）を接続した場合には、各種分析結果データのプリント出力が可能になっている。

＜個数分析処理について＞
本実施形態に係る分類分析装置は、個数分析プログラムの実行によって、分析対象として例えば、１種または２種以上の粒子（分析物の一例）を含む流動性物質（電解質溶液２４）を隔壁１１上側の一面側に供給し、粒子が貫通孔１２を通過することにより生ずる電極１３、１４間の通電変化を検出した検出信号のデータ（計測電流データ）に基づいて粒子種別の個数ないし個数分布を分析する個数分析機能を有する。すなわち、ＰＣ１は、ＣＰＵ２の制御によりＲＯＭ３に格納した個数分析プログラムを実行することにより、データファイル記憶部５に格納、記憶した計測電流データに対する個数分析処理を行うことができる。個数分析処理は、検出信号に含まれ粒子通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、粒子種別の個数を導出する個数分析方法に基づいて、粒子種別個数の自動分析を行うことができる。

図１６は、ＰＣ１の分析処理の説明に必要な処理プログラム構成を示す。各処理プログラムはＲＯＭ３に格納されている。分析対象の一データ実施例として、分析物としての２種の粒子（大腸菌と枯草菌）を含む電解質溶液２４を用いて抽出した計測電流データ（各粒子のパルス抽出データ）を元データに使用している。

個数分析用処理プログラム（個数分析プログラム）には、検出信号として得られた、貫通孔１２の粒子通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度関数を求める確率密度関数モジュールプログラムと、確率密度推定の結果から粒子種別の個数を導出する粒子種分布推定プログラムとが含まれている。分類分析および個数分析に使用される処理プログラムには、データ群から抽出したベースラインを基準にして、パルス状信号の波形形態の特徴を示す特徴量を抽出する特徴量抽出プログラムと、抽出した特徴量に基づいて得られる粒子毎のパルス特徴量データによるデータファイルを作成するデータファイル作成プログラムとが含まれている。分類分析処理および個数分析処理は、データファイル作成プログラムにより作成されたデータに対して実行される。特徴量抽出プログラムには、元の計測電流データから該ベースラインを抽出するベースライン推定処理プログラムを含む。特徴量取得処理（ステップＳ１０１）においては、特徴量抽出プログラムおよびデータファイル作成プログラムを実行して、入力処理（ステップＳ１００）で入力されたデータから特徴量を作成して、ＲＡＭ４の特徴量記憶用データファイルに記憶させる処理が行われる。分類分析用の入力データは、学習データに供される特徴量の作成に必要な既知データと、被分析用のデータ（分析データ）である。既知データから作成された特徴量データは、既知データによる特徴量記憶用データファイルＤＡに記憶され、分析データから作成された特徴量データは、分析データによる特徴量記憶用データファイルＤＢに記憶される。分類分析を行う場合、これらのデータファイルＤＡ、ＤＢから特徴量のベクトル量データを取り込んで分析処理が実行可能になっている。個数分析用の入力データは、被分析用のデータ（分析データ）のみである。個数分析用の入力データから作成された特徴量データは、個数分析用データファイルＤＣに記憶され、個数分析を行う場合、該データファイルＤＣから特徴量のスカラーデータを取り込んで分析処理が実行可能になっている。

粒子種分布推定の前提として真の確率密度関数の形式が未知であるから、確率密度関数モジュールプログラムの実行により、カーネル法と呼ばれるノンパラメトリック（関数形式を指定しない）確率密度推定が行われる。推定対象の元データは、パルス状信号から得られた、例えば、波高ｈ・時間幅Δｔ・出現数等を含むパルス出現分布データである。元の計測データ分布の各データを計測誤差不確定性を導入したガウス分布で表し、各ガウス分布の重ね合わせにより確率密度関数が得られる。確率密度関数モジュールプログラムの実行により確率密度推定処理を行い、元データを該元データに基づいた未知の複雑な確率密度関数（例えば、特徴量のパルス波高・パルス幅・出現確率）で表すことができる。

図４６は、マイクロ・ナノポアデバイス８を用いて、３種の粒子３３ａ、３３ｂ、３３ｃが貫通孔１２を通過して得られる検出信号の波形例と、特徴量に基づいて得られる確率密度関数の導出例を示す。同図（４６Ａ）は、マイクロ・ナノポアデバイス８を用いた粒子検出装置を模式的に示す。同図（４６Ｂ）〜（４６Ｄ）は、各検出信号の波形データを示す。同図（４６Ｅ）〜（４６Ｇ）は、各波形データから得られた確率密度関数の３次元分布図を示す。（４６Ｅ）〜（４６Ｇ）における、ｘ軸、ｙ軸、ｚ軸は、それぞれ、特徴量のパルス波高、パルス幅および確率密度推定により得られた確率密度を示す。

上記のように、ノンパラメトリックな密度関数の推定法の一つであるカーネル法に基づいて確率密度推定処理が行われる。カーネル法は、１つのデータ点にある関数（カーネル関数）を当てはめ、これを全てのデータ点について行い、配置された関数を重ね合わせる推定法であり、滑らかな推定値を得るに適する。

確率密度関数モジュールプログラムの実行により、計測電流波形のパルス波高、パルス幅等のデータから多変数多次元確率密度とみなして２次元以上に拡張して加重の最適推定を行い粒子種別個数分布の推定処理が行われる。加重の最適推定にはＨａｓｓｅｌｂａｌｄ反復法に基づいて実行されるＥＭアルゴリズムソフトウエアが使用される。ＥＭアルゴリズムは、ＰＣ１にあらかじめインストールされている。粒子種別個数分布の推定処理により得られた粒子種別個数分布結果は、表示手段７に粒子種別に対する出現頻度（粒子個数）のヒストグラムで表示出力可能になっている。

本発明に係る特徴量は、パルス状信号由来のパラメータとして、該パルス状信号の波形の局所的特徴を示す第１類型に属するものと、該パルス状信号の波形の全体的特徴を示す第２類型に属するもののいずれかである。これらのうちの１または２以上の特徴量を使用して個数分析を行うことによって、粒子種等の分析物種別に応じた個数ないし個数分布を高精度に分析することができる。

図２４は、大腸菌２２と枯草菌２３の２種の粒子が電解質溶液２４中に混在する様子を模式的に示した貫通孔１２周辺の拡大図である。
＜特徴量について＞
図１７は、実施例の大腸菌と枯草菌につき実測した粒子通過によるパルス波形例を示す。図１７の（４−１）〜（４−９）は、大腸菌の実測パルス波形例（９種類）を示し、（４−１０）〜（４−１８）は、枯草菌の実測パルス波形例（９種類）を示す。両者を外観で比較すると、両者間に波高や波長には差異はあまりないが、ピーク位置や波形尖度等の粒子通過パルス波形形態の属性に顕著な相違がみられる。例えば、大腸菌の場合、ピークが時間経過に伴い前倒し傾向にあり、全体的に波形が尖っている（波形尖度が大きい）。枯草菌の場合、ピークが時間経過に伴い後倒し傾向にあり、波形尖度が小さい。

上記の粒子通過パルス波形形態の属性の違いに基づいて、確率分布作成のベースに用いる特徴量をパルス波形データから粒子種（大腸菌と枯草菌）別に抽出することができる。

図１８は、本発明に係る各種特徴量を説明するためのパルス波形図である。図１８において、横軸は時間、縦軸はパルス波高を示す。

第１類型の特徴量は、
所定の時間幅内における波形の波高値、
パルス波長ｔ_a、
パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
該波形の鋭さ（ピーク波形の広がり）を表す尖度、
パルス開始からパルスピークに至る傾きを表す俯角、
波形を所定の時間毎に区分した時間区分面積の総和を表す面積、および
パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比のいずれかである。

図１８の５ａ〜５ｄは、それぞれ、パルス波長、波高値、ピーク位置比、尖度を示す。図５のＢＬは、パルス波形データから抽出（後述のＢＬ抽出処理参照）した基準ライン（以下、ベースラインという。）を示す。これら４種類のパルス特徴量は、図１９に基づいて示すと以下の（１）〜（４）で定義される。

（１）波長（パルス幅）Δｔ： Δｔ＝ｔ_e−ｔ_s（ｔ_sはパルス波形の開始時間、ｔ_eはパルス波形の終了時間、Δｔ＝ｔ_a）
（２）波高｜ｈ｜： ｈ＝ｘ_p−ｘ_o（ＢＬのｘ_oを基準にしてパルスピークＰＰのｘ_pまでのパルス波形の高さ）
（３）ピーク位置比ｒ： ｒ＝（ｔ_p−ｔ_s）／（ｔ_e−ｔ_s）（パルス波長（＝Δｔ）と、パルス開始からパルスピークｐｐに至るまでの時間ｔ_b（＝ｔ_p−ｔ_s）との比）
（４）ピーク尖度κ： 波高｜ｈ｜＝１、ｔ_s＝０、ｔ_e＝１となるように正規化し、パルスピークＰＰから波高３０％の水平線と交差する時刻の時刻集合[Ｔ]＝[[ｔ_i]｜ｉ＝１，・・・，ｍ]を収集して、下記数２に示すように、時刻集合[Ｔ]のデータの分散をパルス波形広がりとしてκが求められる。

図４７は、俯角、面積および面積比の特徴量を説明するためのパルス波形図である。同図において横軸は時間、縦軸はパルス波高を示す。これら３種類のパルス特徴量は、図に基づいて示すと以下の（５）、（６）、（７）で定義される。

（５）俯角θは、（３４Ａ）に示すように、パルス開始からパルスピークに至る傾きであり、下記数３により定義される。

（６）面積ｍは、下記数４に示すように、単位べクトル[ｕ]と波高ベクトル[ｐ]との内積による面積[ｍ]で定義される。なお、以下の説明において、変数Ａのベクトル表記は[Ａ]で示される。例えば、（３４Ｂ）の１０分割例に示すように、面積ｍは、一つの波形を所定の時間毎に１０分割したときの時間区分面積ｈ_i（幅ｈ_x、高さｈ_yとしたとき、ｈ_i＝ｈ_x×ｈ_y、i＝１〜１０）の総和を表す面積である。

ここで、特徴量計算の準備として、以下に定義するｄ次元波高ベクトル[ｐ]（＝（ｈ₁，ｈ₂，・・・，ｈ_d））をあらかじめ計算して求めておく必要がある。

図４８は、波高ベクトルの取得の仕方を説明するための図である。

（３５Ａ）に示すように、一つの波形データにつき、波長をｄ等分してｄ個のデータグループの分化が行われる。ついで、（３５Ｂ）に示すように、各グループ（各分割区間）ごとに波高の値を平均化して、例えば、１０等分するときには平均値Ａ１〜Ａ１０が求められる。この平均化には、波高値を規格化しない場合と、波高値を規格化する場合とを含むことができる。数４で表記した面積[ｍ]は、規格化しない場合を示す。このようにして求めた平均値を成分とするｄ次元ベクトルが「波高ベクトル」と定義される。

図４９は、ｄ次元の波高ベクトルとデータサンプリングとの関係を説明するための図である。

（３６Ａ）に示すように、パルスデータの取得に関わるサンプリングレートが大きい場合、パルス部分におけるステップ数（データ数）Ｔがベクトルの次元数ｄを上回るので、上述の取得手順により各区分の平均値を成分した波高ベクトルを得ることができる。一方、サンプリングレートを下げていくと、パルス部分におけるステップ数Ｔがベクトルの次元数ｄ（＞Ｔ）を下回る事態が生ずる。Ｔ＜ｄの場合、上述の取得手順により各区分の平均値を取得できないので、３次スプライン補間によってｄ次元の波高ベクトルを取得することができる。

特徴量抽出プログラムには、波高ベクトルデータを取得するための波高ベクトル取得プログラムが含まれている。波高ベクトル取得プログラムの実行により、パルスステップ数Ｔがベクトルの次元数ｄを上回るか（Ｔ＞ｄ）同等の場合（Ｔ＝ｄ）、時間方向にｄ等分した各区分の平均値を求め、該平均値を成分とするｄ次元波高ベクトルを取得し、パルスステップ数Ｔがベクトルの次元数ｄを下回る場合（Ｔ＜ｄ）、３次スプライン補間を実行してｄ次元波高ベクトルを取得するようになっている。すなわち、３次スプライン補間法を用いた補間処理を行うことにより、パルスステップ数が少ない場合にもベクトルの次元数を一定にすることができる。

（７）面積比ｒ_mは、（３４Ｂ）に図示した時間区分面積ｈ_iをパルス開始からパルスピークに至るまでの区間での和の、全波形面積に対する面積比で定義される。下記数５は、面積比ｒ_mを示す。

第１類型の特徴量は、パルス波高、パルス波長、パルス面積等のパルス状信号の波形に一義的に由来し、局所的特徴を示す特徴量である。第２類型の特徴量は、第１類型の局所的特徴に対し全体的特徴を示す特徴量である。

第２類型の特徴量は、
パルス開始時点を中心にして時間区分面積を質量に、かつ該中心から時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値をベクトルの成分とする平均値ベクトル、
前記平均値ベクトルに対し波長が基準値になるように規格化した場合の規格化された平均値ベクトル、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメントのいずれかである。

図５０は、時間（波長）および波幅に関する第２類型の特徴量を説明するためのパルス波形図である。同図において横軸は時間、縦軸はパルス波高を示す。これらのパルス特徴量は、同図に基づいて示すと、以下の（８）〜（１５）で定義される。

（８）時間慣性モーメントは、（３４Ｂ）と同様に一つの波形を所定の時間毎にｉ次元で等分割したときの時間区分面積ｈ_iを質量に、かつ該中心から時間区分面積ｈ_iに至る時間を回転半径に擬制したときに定まる特徴量である。すなわち、時間慣性モーメントの特徴量は、下記数６に示すように、ベクトル[ｖ]と波高ベクトル[ｐ]との内積による[Ｉ]で定義される。ここで、ベクトルの次元をｎとしたとき、[ｖ]＝（１²，２²，３²，・・・ｎ²）および[ｐ]＝（ｈ₁，ｈ₂，・・・，ｈ_d）である。例えば、時間慣性モーメントは、（３７Ａ）の１０分割例に示すように、（３４Ｂ）と同様に一つの波形を所定の時間毎に１０分割したとき、時間区分面積ｈ_i（幅ｈ_x、高さｈ_yとしたとき、ｈ_i＝ｈ_x×ｈ_y、i＝１〜１０）を質量に、かつ該中心から時間区分面積ｈ_iに至る時間を回転半径に擬制したときに定まる特徴量であり、（６）の面積ｍと同様に、波高ベクトルにより求めることができる。

（９）規格化された時間慣性モーメントは、（８）で示した時間区面積を作成した波形に対して波高が基準値の「１」になるように波高方向に規格化した波形を用いて、（８）と同様にして作成した波高ベクトルｈ_iにより数６で定義される特徴量である。

（１０）平均値ベクトルは、（３７Ｂ）の１０分割例に示すように、一つの波形を波高方向にｉ次元で等分割し、パルスピーク前後それぞれにおいて各分割単位(分割領域ｗ_i）毎に時刻値の平均値を算出し、分割領域ｗ_iの同一波高位置の平均値をベクトルの成分とする特徴量である。

（１１）規格化された平均値ベクトルは、（１０）の平均値ベクトルに対し波長が基準値になるように規格化した場合の特徴量である。

（１２）波幅平均値慣性モーメントは、（３７Ｂ）の１０分割例に示すように、一つの波形を波高方向にｉ次元で等分割し、パルスピーク前後それぞれにおいて各分割単位(分割領域ｗ_i）毎に時刻値の平均値を算出し、分割領域ｗ_iの同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布ｈ_i（ベクトルの次元数をｎとしたとき、ｉ＝１〜ｎ）と擬制して波形裾野の時間軸Ａｔを回転中心にした場合の慣性モーメントとして定義される特徴量である。定義式は、数６と同じであり、（１２）の特徴量は、ベクトル[ｖ]と質量分布ｈ_iとの内積により求めることができる。

（１３）規格化された波幅平均値慣性モーメントは、上記の分割領域ｗ_iを作成した波形に対して波長が基準値の「１」になるように波長方向に規格化した波形を用いて、（１２）と同様にして作成した質量分布ｈ_iにより数６で定義される特徴量である。

（１４）波幅分散慣性モーメントは、波幅平均値慣性モーメントと同様に、一つの波形を波高方向にｉ次元で等分割し、パルスピーク前後それぞれにおいて各分割単位(分割領域ｗ_i）毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布ｈ_i（ベクトルの次元数をｎとしたとき、ｉ＝１〜ｎ）と擬制して波形裾野の時間軸Ａｔを回転中心にした場合の慣性モーメントとして定義される特徴量であり、波幅平均値慣性モーメントと同様に、数６で定義される。

（１５）規格化された波幅分散慣性モーメントは、上記の分割領域ｗ_iを作成した波形に対して波長が基準値の「１」になるように波長方向に規格化した波形を用いて、（１４）と同様にして作成した質量分布ｈ_iにより数６で定義される特徴量である。

波幅平均値慣性モーメントおよび波幅分散慣性モーメントは、上記のように、数６で定義される特徴量であり、該定義中のベクトル[ｐ]は、波幅平均値慣性モーメントの場合、時刻値の平均値の差のベクトルであり、波幅分散慣性モーメントの場合、時刻値の分散ベクトルである。以下の説明で、（１２）〜（１５）の波幅に関する慣性モーメントにおけるベクトル[ｐ]を[ｐ_w]と表す。

（１２）〜（１５）の波幅に関する慣性モーメントのデータ作成演算には、図３６に示した波高ベクトルの縦横軸を交換した波幅ベクトル[ｐ_w]（＝[ｐ₁，ｐ₂，・・・，ｐ_dw]）を用いて行われる。波幅ベクトルは、（１２）〜（１５）の特徴量の定義で示した平均値の差ベクトルまたは分散ベクトルである。波幅ベクトルを密度分布として捉えることによって、（１２）と（１３）の波幅平均値慣性モーメントおよび（１４）と（１５）の波幅分散慣性モーメントを求めることができる。波幅ベクトルは、パルス波形データを波高方向にｄ_w等分し、各区分毎に求めた波高値の平均値の差または分散を成分にもつｄ_w次元ベクトルである。（３７Ｂ）の場合、波幅ベクトルの次元は１０次元である。（３７Ｂ）に示した時間軸Ａｔは、ベースラインＢＬと異なり、波幅ベクトルから得られるパルス裾野周りの回転軸ラインである。

図５１は、ｄ_w次元の波幅ベクトルとデータサンプリングとの関係を説明するための図である。

特徴量抽出プログラムには、以下のｄ_w次元の波幅ベクトルの作成演算処理により波幅ベクトルを取得するための波幅ベクトル取得プログラムが含まれている。

パルス波形データは、波高方向に様々な間隔で分布しているため、波高方向に分割した区画においてデータ点が存在しない不存在領域Ｂｄを１または２以上含む場合が生ずる。（３８Ａ）においては、不存在領域Ｂｄの１例を矢印で示している。不存在領域Ｂｄは、データ間隔が粗くなってデータ点が存在せず、上記の数６で定義される波幅に関する慣性モーメントの成分を得ることができない。そこで、前述のパルス波形広がりの場合と同様に、パルスピークまでの波高をｄ_w等分したときの各波高の時刻集合[Ｔｋ]＝[[ｔ_i]｜ｉ＝１，・・・，ｍ]を収集することにより波幅ベクトルの成分が作成される。このとき、データ点が存在しない不存在領域Ｂｄにおいては、線形補間によって成分データの取得が行われる。該線形補間は、パルスピークの（１０ｋ＋５）％（ｋ＝０，１，２，３、・・・）の値をまたぐ連続した２つのデータに対して行われる。（３８Ｂ）は、データ点ｔ_iとｔ_i+1の間に生じた不存在領域Ｂｄの高さｋについての線形補間点ｔ_kの一例を示している。なお、波幅ベクトルの作成に際して、（３８Ｃ）に示すように、パルス波形データの裾野領域ＵＲに波高データの食い違いが生じている場合、パルスピークに近い方に揃えるように、パルスピークに遠い側の波高データＤｕは切り捨てられる。波幅ベクトル取得プログラムの実行処理には、不存在領域Ｂｄに対する線形補間処理と、波高データの食い違いに対する波高データＤｕの切り捨て処理が含まれている。

図５２は、波幅に関する慣性モーメントを波幅ベクトルにより取得する取得過程を説明するための図である。

（３９Ａ）は、波形を波高方向に１０等分した例であり、一つの波形３９ａに対して、上記の線形補間処理および切り捨て処理を行うことにより得られた、波幅ベクトルの分割領域３９ｂおよび回転軸ライン３９ｃを示す。

（３９Ｂ）に示すように、各分割単位ごとに、パルスピーク前後それぞれにおいて時刻値の平均値を算出し、分割領域の同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルの波幅ベクトルを取得することができる。この平均値の差ベクトルを質量分布と擬制して回転軸ライン３９ｃ（時間軸）を回転中心にした（１２）の波幅平均値慣性モーメントを作成することができる。また、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを取得することができる。この分散ベクトルを質量分布と擬制して回転軸ライン３９ｃ（時間軸）を回転中心にした（１４）の波幅分散慣性モーメントを作成することができる。なお、（１０）、（１１）の平均値ベクトルは、
時刻値の平均値を算出して分割領域の同一波高位置の各平均値を成分としたベクトルであり、ｄ_w等分した場合、２ｄ_w次元の時刻ベクトルで表わされる。

特徴量の作成に用いた波高ベクトルおよび波幅ベクトルのベクトル次元数は、分割数に拘泥する必要はなく任意に設定することができる。波高ベクトルおよび波幅ベクトルは、波長ないし波高の一方向に細分化した場合であるが、特徴量の作成には、複数方向に細分化したベクトルを使用することができる。

図５３は、複数の方向に分割した場合の特徴量作成用波形ベクトルの一例を説明するための図である。

（４０Ａ）は、一つの波形データをメッシュ状に分割したデータマップ４０ａを示す。データマップ４０ａは、波形データを横軸の時間軸方向にｄ_n分割し、縦軸の波高方向にｄ_w分割してマトリクス状にデータ点の数の分布状態を示す。（４０Ｂ）は、マトリクス状の区画（格子）の一部を拡大した分布状態を示す。（４０Ｂ）の分布状態において、１１×１３個の格子に０〜６個のデータ点数が分布している。このマトリクス分割によって、各格子内のデータ点数／総データ点数をｄ_n×ｄ_w次元ベクトルの成分とした波形ベクトルを、マトリクス配列のデータ群を走査状に並べ替えたベクトルに変換することにより波高ベクトルや波幅ベクトルに代えて特徴量の作成に使用することができる。
＜ベースラインの推定について＞
一般に細菌等は微細に異なる形態を有する微小物体である。例えば、平均的な大腸菌の場合、２〜４μｍの体長で外径が０．４〜０．７μｍである。平均的な枯草菌の場合、２〜３μｍの体長で外径が０．７〜０．８μｍである。さらに、大腸菌などには２０〜３０ｎｍの鞭毛が付随している。

検体粒子として細菌等を使用する場合には、パルス波形データから僅かの違いを見逃すと、個数判定精度の低下をもたらしてしまう。このため、特徴量を正確に算出して確率分布の推定基礎とするために、粒子通過パルス波高を正確に把握する必要があり、これには計測信号のベースラインの推定を行う必要がある。しかし、計測信号の元データのベースラインには、ノイズデータや微弱な計測電流による揺らぎが含まれているために、この揺らぎ成分等を除いたベースラインを確定してからパルス波高等を検出する必要がある。ベースラインの推定（以下、ＢＬの推定という。）は実用上、コンピュータによりオンラインで（即時的に）行うのが好ましい。

ＢＬの推定をコンピュータ上で行う手法として、離散的な誤差のある観測から時々刻々変化する量を推定するに好適なカルマンフィルターを使用すれば、これにより外乱（システム雑音や観測雑音）を取り除いてベースラインＢＬを推定することができる。

カルマンフィルターとは、離散的制御過程が図１９の（６Ａ）に示す線形差分方程式により定義され、更新可能な状態ベクトル[ｘ]の時刻[ｔ]における値を推定する手法である。カルマンフィルターにおいては、状態ベクトル[ｘ]およびシステム制御入力[ｕ_t]の値は直接観測できないものとされる。

状態ベクトル[ｘ]は、図１９の（６Ｂ）に示す観測モデルにより、間接的に推定されるものとする。システム制御入力[ｕ_t]については、その統計的変動幅[σ_u,t]のみをパラメータとして仮定する。

計測電流データ[Ｘ]はベクトルではなくスカラーであり、さらに各種行列もスカラーであり、[Ｆ]＝[Ｇ]＝[Ｈ]＝[１]とみなすことができる。したがって、時刻ｔの実際の電流値のベースラインレベル、時刻ｔで計測された電流、時刻ｔの観測ノイズをそれぞれ、[ｘ_t]、[ｙ_t]、[ν_t]とすると、[ｘ_t]および[ｙ_t]は図１９の（６Ｃ）に示すように表される。[ｘ_t]、[ｕ_t]、[ν_t]は観測不可能な因子であり、[ｙ_t]は観測可能な因子である。イオン電流検出部による計測周波数をｆ（Ｈｚ）とすると時刻データは１／ｆ（秒）刻みとなる。システム制御入力[ｕ_t]の影響は実際上非常に小さいものと仮定してベースラインの推定を行うことができる。

図２０は、上記の各因子を実際の計測電流データで示した図である。イオン電流検出部による実際の計測の際には、粒子が貫通孔１２に詰まったりして、ベースラインの歪みが生ずるが、計測時には歪みの発生時点で中断して、歪み原因を除去してから計測が行われるので、元のデータ集合には歪みのないベースラインを含むデータのみが収集されている。

カルマンフィルターによる推定は予測と更新の繰返しにより行われる。ベースラインの推定にもカルマンフィルターによる予測と更新を繰返して実行される。

図２１は、カルマンフィルターにおける予測（８Ａ）と更新（８Ｂ）の繰返しの詳細を示す図である。図８において、ベクトル表記に付加した「ハット」記号は推定値を示している。添え字の「ｔ｜ｔ−１」は（ｔ−１）時点の値に基づく、ｔ時点の値の推定値であることを示している。

図２２は、ＢＬ推定処理プログラムに基づくＢＬ推定処理を示す。ＢＬ推定処理にいては、ＢＬの推定と、それに基づくパルス波高値の抽出が行われる。

ＢＬ推定処理の実行に際しては、カルマンフィルターにおける予測と更新の処理に必要な調整因子の開始時刻ｍ、定数ｋ、αの値は推定対象のデータ属性に応じて適切な値にあらかじめ調整（チューニング）して決めておく必要がある。αの値はベースラインの推定値の分散を調整するための値である。ｋの値は図２１に示したカルマンフィルターにおける更新Ａの実行回数に関係する値である（ステップＳ５７、Ｓ６２参照）。開始時刻ｍは計測サンプリングの１個分を１ステップとして計算されたステップ数分の時間データである。

図２３は、該調整に使用したビーズモデルの波形図を示す。図１５においては、粒子として細菌等と同程度の大きさの微小ビーズ玉を混入させた場合（ビーズモデル）の溶液状態を示している。図２３の（１０Ａ）はイオン電流検出部によってサンプリング周波数９０００００Ｈｚで取得した波形データである。（１０Ａ）に示すビーズモデルの波形はなだらかに減衰していく波形を示している。（１０Ａ）の右端部分に激しい落ち込みが生じており、それを拡大して（１０Ｂ）に示している。

ビーズモデルの波形から（１０Ｂ）に示すベースラインの段差部分（１０Ｃ）が検出された場合に、その直前期間が初期値計算期間となる。例えば、ｍ＝１０００００とした場合、当該初期値計算期間を除いた期間において有意性のあるパルスが１１〜１２個目視で確認することができる。

図２５は、調整因子のｍ、ｋ、αの組合せに応じてビーズモデルの波形から拾われたパルスの数を示す表である。

図２５の（１２Ａ）は、ｍ＝１００００の場合のｋ値（１０、３０、５０、７０、９０）、α値（２、３、４、６）の組合せによるパルス数を示す。同図（１２Ｂ）は、ｍ＝５００００の場合のｋ値（１０、３０、５０、７０、９０）、α値（２、３、４、６）の組合せによるパルス数を示す。（１２Ｃ）は、ｍ＝１０００００の場合のｋ値（１０、３０、５０、７０、９０）、α値（２、３、４、６）の組合せによるパルス数を示す。

図２５の３種のシミュレーション結果を比較すると、（１２Ａ）と（１２Ｂ）の場合、計測されるべきパルス数は１２になり、（１２Ｃ）では１１となっている。したがって、実施例では、パルス数の最大値の最も小さい（１２Ｃ）を採用して、ｍ＝１０００００、ｋ＝５０、α＝６のチューニング設定を行っている。これらのチューニング設定データはあらかじめＲＡＭ２３の設定エリアに記憶、設定されている。

図２２のＢＬ推定処理は上記チューニング設定下、図２１に示したカルマンフィルターによるＢＬ推定で行われる。まず、ステップＳ５１にて、時刻ｍにおけるカルマンフィルターの初期値がＲＡＭ２３のワークエリアに設定される。このとき、データファイル記憶部５に格納したパルス波形データはＲＡＭ２３のワークエリアに読み込まれる。ついで、時刻（ｍ＋１）におけるカルマンフィルターの予測と更新（図２１のＡおよびＢ）が実行される（ステップＳ５２）。予測と更新においては、図２１に示したカルマンフィルターの各演算が実行され、ＲＡＭ２３に記憶される。以降、所定の単位時間毎に予測と更新（ＡおよびＢ）が繰り返し実行され、時刻ｔにおけるカルマンフィルターの予測と更新Ａが
行われたとき、下記数６の条件が満たされたか否かが判断される（ステップＳ５３、Ｓ５４）。単位時間は元データのサンプリング周波数により定まる値であり、あらかじめＲＡＭ２３にセットされている。

数７の条件が満たされない場合、時刻ｔにおけるカルマンフィルターの更新Ｂが実行され、単位時間経過したデータ毎にステップＳ５３〜Ｓ５５の処理が繰り返される。上記数７の条件が満たされた場合、その回数値が１回ごとにＲＡＭ２３のカウントエリアに累積記憶される（ステップＳ５４、Ｓ５６）。ついで、該カウント値に基づき、数７の条件が時刻ｓを起点としてｋ回連続して満たされたか否かが判断される（ステップＳ５７）。ｋ回連続していない場合はステップＳ５５に進み、更新Ｂが行われる。

ｋ回連続した場合はステップＳ５８に進み、ＢＬ確定のためのホールド必要期間が開始したと判定される。このとき、ホールド必要期間のホールド開始時刻をｓとしてＲＡＭ２３に記憶されるとともに、時刻（ｓ＋１）〜時刻（ｓ＋ｋ−１）の間のカルマンフィルターの演算結果は記憶されずに捨てられる。

ホールド必要期間の開始により、時刻ｔにおけるパルスの落ち込み最大値がＲＡＭ２３に更新可能に記憶される（ステップＳ５９）。ついで、ステップＳ５４と同様に、ホールド必要期間における、下記数８の条件が満たされているか否かの判断が行われる（ステップＳ６０）。

上記数８の条件が満たされない場合、パルスの落ち込み最大値の更新が行われる（ステップＳ５９、Ｓ６０）。数８の条件が満たされた場合、その回数値が１回ごとにＲＡＭ２３のカウントエリアに累積記憶される（ステップＳ６０、Ｓ６１）。ついで、該カウント値に基づき、数８の条件が時刻ｓ２を起点としてｋ回連続して満たされたか否かが判断される（ステップＳ６２）。ｋ回連続していない場合はステップＳ５９に戻る。

ｋ回連続した場合はステップＳ６３に進み、このとき更新記憶されたパルスの落ち込み最大値がパルス波高値の推定値としてＲＡＭ２３に記憶される。パルス波高値の推定値はパルス開始時刻およびパルス終了時刻のデータとともに記憶される。パルス波高値の推定を終えると、ホールド必要期間は終了と判定される。この終了によりホールド必要期間のホールド終了時刻がｓ２としてＲＡＭ２３に記憶される（ステップＳ６４）。次に、ステップＳ６５に進み、時刻ｓの値はカルマンフィルターの演算処理の再開時の初期値として、時刻ｓ２〜時刻（ｓ＋ｋ−１）の期間について遡及してカルマンフィルターの演算が実行される。ステップＳ６５の後は、全パルス波形データのＢＬ推定処理を行ったか否かが判断されて（ステップＳ６６）、全パルス波形データの推定完了で終了し、残データがあるときはステップＳ５３に移る。
＜特徴量抽出について＞
図２６は、特徴量抽出プログラムの実行処理内容の概要を示す。

特徴量抽出処理は、図２２の上記ＢＬ推定処理の実行によってパルス波高値（波高｜ｈ｜）の抽出データがあることを条件に実行可能になる（ステップＳ４１）。パルス波高値の抽出データがある場合、前述の波高ベクトル取得プログラムおよび波幅ベクトル取得プログラムが実行されて、各種ベクトルのデータ作成演算が実行される（ステップＳ４２）。波高ベクトルおよび波幅ベクトルの全てのデータ取得を終えると、該ベクトルデータが保存される（ステップＳ４３、Ｓ４４）。ついで、各種の特徴量の抽出処理が実行される（ステップＳ４５）。波高ベクトルおよび波幅ベクトルのデータ取得に際しては、３次スプライン補間法を用いた補間処理、線形補間処理および切り捨て処理が随時行われる。

図５４は、特徴量の抽出処理（ステップＳ４５）の実行処理内容を示す。ステップＳ７１〜Ｓ８３は、それぞれ上記（１）〜（１３）で定義された第１類型および第２類型の特徴量の算出と、算出された特徴量の記憶、保存の処理を示す。

第１類型の特徴量は、ステップＳ７１〜Ｓ７６において算出される。波長（パルス幅）Δｔは、パルス波高値の抽出データ群に対し時系列的に順次算出されて記憶される（ステップＳ７１）。算出された特徴量は、ＲＡＭ４の特徴量記憶用メモリエリアに記憶される。パルス幅は、Δｔ（＝ｔ_e−ｔ_s;ｔ_sはパルス波形の開始時刻、ｔ_eはパルス波形の終了時刻）を演算して求められる。ピーク位置比ｒは、パルス波高値の抽出データ群に対し時系列的に順次算出されて記憶される（ステップＳ７２）。ピーク位置比ｒは、ｒ＝（ｔ_p−ｔ_s）／（ｔ_e−ｔ_s）（パルス幅Δｔと、パルス開始からパルスピークｐｐに至るまでの時間（＝ｔ_p−ｔ_s）との比）を演算して求められる。

ピーク尖度κは、パルス波高値の抽出データ群に対し時系列的に順次算出されて記憶される（ステップＳ７３）。パルス波高値｜ｈ｜＝１、ｔ_s＝０、ｔ_e＝１となるように正規化し、パルスピークＰＰから波高３０％の水平線と交差する時刻の時刻集合Ｔ＝[[ｔ_i]｜ｉ＝１，・・・，ｍ]を収集して、時刻集合Ｔのデータの分散を演算してパルス波形広がりとしてκが求められる。

俯角θは、パルス開始からパルスピークまでの時刻と波高のデータと、前掲の数２の演算とに基づいて求められる（ステップＳ７４）。面積ｍは、波高ベクトルのデータにより求められ、時間区分面積ｈ_iを区分数に応じて求め、それらの総和を求めることにより算出、記憶される（ステップＳ７５）。該区分数は、任意に設定可能であり、例えば、１０である。面積比ｒ_mは、全波形面積と、時間区分面積ｈ_iをパルス開始からパルスピークに至るまでの区間での部分和とをそれぞれ求め、部分和の全波形面積に対する面積比を算出して記憶される（ステップＳ７６）。

第２類型の特徴量は、ステップＳ７７〜Ｓ８２において算出される。時間慣性モーメントは、波高ベクトルのデータにより求められ、区分数に応じて求めた時間区分面積ｈ_iと、前掲の数６の演算とに基づいて算出され、記憶される（ステップＳ７７）。（９）の規格化された時間慣性モーメントは、ステップＳ７７において得られた時間慣性モーメントに対し波高が基準値の「１」になるように波高方向に規格化処理（波高ベクトルと規格化ベクトルの内積）した規格化データとして記憶される（ステップＳ７８）。波幅平均値慣性モーメントは、ステップＳ４２〜Ｓ４４で求めた波幅ベクトル（平均値の差ベクトル）のデータから、パルスピーク前後それぞれにおいて分割単位(あらかじめ設定された分割数：１０）毎に算出した時刻値の平均値の差と、前掲の数７の演算とに基づいて算出され、記憶される（ステップＳ７９）。（１１）の規格化された波幅平均値慣性モーメントは、ステップＳ７９において得られた波幅平均値慣性モーメントに対し波長が基準値「１」になるように波長方向に規格化処理（平均値の差ベクトルと規格化ベクトルの内積）した規格化データとして記憶される（ステップＳ８０）。波幅分散慣性モーメントは、波幅ベクトル（分散ベクトル）のデータから、分割単位毎に算出した時刻値の分散と、前掲の数７の演算とに基づいて算出され、記憶される（ステップＳ８１）。（１３）の規格化された波幅分散慣性モーメントは、ステップＳ８１において得られた波幅分散慣性モーメントに対し波長が基準値「１」になるように波長方向に規格化処理（分散ベクトルと規格化ベクトルの内積）した規格化データとして記憶される（ステップＳ８２）。

全データからの特徴量の抽出を終了すると、各データのファイル保存が行われ、別のデータ群があるか否かが判断される（ステップＳ８３、Ｓ８４）。引き続き別ファイルのデータ群があれば、上記処理（ステップＳ７１〜Ｓ８２）が繰返し実行可能になっている。処理すべきデータがなくなれば特徴量の抽出処理は終了する（ステップＳ８５）。上記の抽出処理においては、第１類型および第２類型の特徴量を全て求めるようにしているが、入力手段６の指定入力により所望の特徴量を指定可能にし、該指定による特徴量だけを抽出可能にすることができる。

図２７は、粒子種分布推定プログラムに基づいて実行される粒子種推定処理を示す。＜確率密度関数の推定について＞ 同種の粒子であっても計測されるパルス波形が一定とは限らないので、粒子種分布推定のための準備として、テストデータから予め粒子種別のパルス波形の確率密度関数の推定が行われる。確率密度関数の推定により導出される確率密度関数によって各パルスの出現確率を表すことができる。

図２８の（１５Ｂ）は大腸菌と枯草菌の粒子種において、パルス波形の特徴量としてパルス幅とパルス波高が用いて得られたパルス波形に対する確率密度関数のイメージ図であり、図中の濃淡によりパルスの出現確率を表している。図２８の（１５Ａ）は１つの波形データに関する第１類型の特徴量の一部を示す。

パルス幅Δｔとパルス波高ｈの真の密度関数は未知であるので、ノンパラメトリックな確率密度関数の推定を行う必要がある。本実施形態では、カーネル関数としてガウス関数を採用したカーネル密度推定を用いている。

カーネル密度推定とは、計測データにカーネル関数で与えられた確率密度分布を想定し、それらの分布を重ね合わせた分布を確率密度関数とみなす手法である。カーネル関数としてガウス関数を用いた場合は、各データに対して正規分布を想定し、それらを重ね合わせた分布を確率密度関数とみなすことができる。

図２９は、大腸菌と枯草菌の粒子種の個々より得られた確率密度分布の重ね合わせのイメージ図である。同図（１６Ｃ）は、パルス幅Δｔとパルス波高ｈの特徴量データ（１６Ａ）から、各粒子につき求めた確率密度分布（１６Ｂ）を重ね合わせた状態を示している。

入力データ[ｘ]に対する確率密度関数[ｐ（ｘ）]は教師データ数[Ｎ]、教師データ[μ_i]、分散共分散行列[Σ]を用いて、下記数９で表される。

さらに、確率密度関数[ｐ（ｘ）]は下記数１０に示すように、各次元のガウス関数の積で表すことができる。

数１０からわかるように、各パルス属性が正規分布に従う独立な確率変数であると仮定していることに相当し、これは３次元以上にも同様に拡張可能である。したがって、本実施形態においては２種以上の粒子種個数の分析が可能である。

確率密度関数モジュールプログラムは、２種の特徴量に対する確率密度関数を演算して求める機能を有する。すなわち、２つの特徴量[（β，γ）]による推定対象データを使う場合、カーネル関数としてガウス関数を採用したカーネル密度推定における確率密度関数[ｐ（β，γ）]は下記数１１で表される。

数１１に基づいて確率密度関数モジュールプログラムにより実行される確率密度関数推定処理は後述の図３３により詳述するように、２つの特徴量における確率密度関数の推定処理を行う。

図３０は、ｋ個の粒子種別の粒子総数と、粒子種別の出現確率と、データ全体の出現頻度の期待値との関係を示すイメージ図である。同図（１７Ａ）はデータ全体の出現頻度を示す。同図（１７−１）〜（１７−ｋ）は粒子種別の出現頻度を示す。パルス[ｘ]が計測される出現頻度の期待値は、粒子種別の確率密度関数に従ってパルス[ｘ]が計測される出現頻度の期待値の和となる。図３０に示すように、粒子種別の粒子総数[ｎ_i]と、粒子種別出現確率[ｐ_i（ｘ）]とから粒子種別の期待値の和として下記数１２で表すことができる。

本実施形態においては、あらかじめ求められた粒子種別の確率密度関数の推定を行った確率密度関数データ（数１０参照）が分析参照データとしてＲＡＭ２３に記憶されている。粒子種別個数分析は、数１２に基づいて、分析対象の全体データの出現頻度を各分析データから適合する粒子種別の個数を割り出すことにより行われる。個数分析は異なる粒子種のヒストグラム（粒子種に対する出現頻度（粒子数））を推定することにより行われる。

図２７の粒子種推定処理においては、データの編集によって特徴量によるデータファイルを作成するデータファイル作成処理（ステップＳ１）と、粒子数の推定処理（ステップＳ２）と、推定粒子種分布の算出処理（ヒストグラム作成処理）（ステップＳ３）が行われる。粒子数の推定処理においては、最尤法、ラグランジュ未定乗数法およびＨａｓｓｅｌｂｌａｄ反復法による推定手法が用いることができる。
＜最尤法（ｍａｘｉｕｍ ｌｉｋｅｌｉｈｏｏｄ ｅｓｔｉｍａｔｉｏｎ：統計学において与えられたデータからそれが従う確率分布の母数を点推定する方法である。）について＞
今、実際のパルス推定結果として、データセット[Ｄ]＝[ｘ₁，ｘ₂，ｘ₃，・・・ｘ_N]が得られているとする。推定されたｊ番目のパルス波高データが出現する尤度（尤もらしさ）は下記数１３で表わされる。

するとデータセットＤが出現する尤度は下記数１４で表される。

数１４の尤度を最大化する様な粒子種分布の値セット[ｎ]＝[ｎ₁，・・・，ｎ_k]^Tが最も尤もらしい粒子種分布である。
＜ラグランジュ未定乗数法(束縛条件のもとで最適化を行う解析学的方法であり、各束縛条件に対して未定乗数を用意し、これらを係数にする線形結合を新しい関数（未定乗数も新たな変数とする）として捉えることによって束縛問題を通常の極値問題として解く方法である。）について＞
データセットＤが出現する尤度を最大化することは、データセット[Ｄ]が出現する対数尤度を最大化することに等しい。下記数１５はラグランジュ未定乗数法の適否を調べるための対数尤度を導出する過程を示す。

数１５において、途中の係数１／Ｎ^Nは最終式では省略している。

ここで、粒径個数分布の値セットｎ＝[ｎ₁，・・・，ｎ_k]^Tには、「合計がＮである」という制約（下記数１６参照）がある。

したがって、最も尤もらしい粒子種分布を得るという命題は、制約付き対数尤度最大化の問題になるので、ラグランジュ未定乗数法により最適化を行うことが可能である。ラグランジュ未定乗数法により最適化を行う制約付き対数尤度最大化式を下記数１７で表すことができる。

数１７に示す制約付き対数尤度最大化式からは、図３１に示す数学的導出過程を経て下記数１８に示す[ｋ]個の連立方程式を導き出すことができる。

数１８に示す連立方程式を数値的に解くには、Ｈａｓｓｅｌｂｌａｄが提唱した反復法を用いて行うことができる。Ｈａｓｓｅｌｂｌａｄ反復法によれば、下記数１９の反復計算を行えばよい。この反復法の詳細は提唱論文（Ｈａｓｓｅｌｂｌａｄ Ｖ ．，１９６６，Ｅｓｔｉｍａｔｉｏｎ ｏｆ ｐａｒａｍｅｔｅｒｓ ｆｏｒ ａ ｍｉｘｔｕｒｅ ｏｆ ｎｏｒｍａｌ ｄｉｓｔｒｉｂｕｔｉｏｎｓ． Ｔｅｃｈｎｏｍｅｒｉｃｓ，８，ｐｐ．４３１−４４４）に記述されている。

数１９の反復計算には、市販のＥＭアルゴリズムのソフトウエアを利用して行われる。ＥＭアルゴリズムは、命名の由来から明らかなように、確率分布のパラメータを、尤度関数を最大化することで計算する方法、つまり尤度関数である確率分布の期待値（Ｅｘｐｅｃｔａｔｉｏｎ）を最大化（Ｍａｘｉｍｉｚａｔｉｏｎ）することができるアルゴリズムである。ＥＭアルゴリズムによれば、求めたいパラメータの初期値を設定して、その値から尤度(期待値)を計算して、多くの場合、尤度関数の偏微分が０になる条件を使って、繰り返し計算で最大尤度のパラメータを計算することができる。ＥＭアルゴリズムを使用して行うＨａｓｓｅｌｂｌａｄ反復法の演算処理は、求めるパラメータの初期値を設定して、その値から尤度(期待値)を計算し、さらに尤度関数の偏微分が０になる条件を使って繰り返し計算を行って最大尤度のパラメータを計算する工程を有する。

＜粒子種推定処理について＞
図２７に示した粒子種推定処理において実行可能な、データファイル作成処理（ステップＳ１）、確率密度関数の推定処理（ステップＳ２）、粒子数の推定処理（ステップＳ３）および推定粒子種分布の算出処理（ステップＳ４）を以下に詳述する。

図３２はデータファイル作成プログラムにより実行されるデータファイル作成処理（ステップＳ１）を示す。

ＰＣ１の入力手段６を使用して、データファイルを作成するｋ個（実施例では２個）ずつの特徴量の指定操作を行うことができる。指定された特徴量の組合せ入力がＲＡＭ２３に設定される（ステップＳ３０）。特徴量の設定毎の特徴量データファイルのデータがＲＡＭ２３のワークエリアに読み込まれる（ステップＳ３１）。特徴量データファイルは、図２２のＢＬ推定処理および図２６の特徴量抽出処理で抽出され、ファイル保存されている特徴量（パルス波高値等）データである。

個数推定に使う特徴量をｋ個指定することにより、Ｎ行ｋ列の行列データが作成される（ステップＳ３２）。作成された行列データは、粒子種分布推定用データファイルに出力され、指定特徴量別に保存される（ステップＳ３３）。指定特徴量に対するすべてのデータファイルの生成を終えると終了する（ステップＳ３４）。

図３３は、確率密度関数モジュールプログラムにより実行される確率密度関数の推定処理（ステップＳ２）を示す。確率密度関数推定処理は、数６に基づいて２つの特徴量における確率密度関数の推定処理が行われる。

データファイル作成処理（ステップＳ１）において作成された、確率密度関数推定対象のデータファイルのデータを読み込んで、Ｎ行２列の行列[Ｄ]が作成される（ステップＳ２０、Ｓ２１）。行列[Ｄ]の列ごとの下記数２０に示す分散

が算出される（ステップＳ２２）。ついで、下記数２１に示す分散パラメータが標準偏差係数ｃを用いて下記数２２に示すように設定される（ステップＳ２３）。

分散パラメータおよび行列[Ｄ]の各行を下記数２３に示す教師データとして代入されて確率密度関数が求められ、ＲＡＭ２３の所定エリアに記憶される（ステップＳ２４、Ｓ２５）。上記のステップＳ２０〜Ｓ２５の処理は全ての処理対象データからの確率密度関数の導出を行うまで行われる（ステップＳ２６）。

図３４は粒子数の推定処理（ステップＳ３）を示す。

まず、上記のステップＳ２０、Ｓ２１と同様に、データファイル作成処理において作成された、粒子数推定対象のデータファイルのデータを読み込んで、Ｎ行２列の行列[Ｄ]が作成される（ステップＳ１０、Ｓ１１）。行列[Ｄ]データに対して、Ｈａｓｓｅｌｂｌａｄ反復法による推定処理が実行される（ステップＳ１２）。

図３５は、ＥＭアルゴリズムによって実行される、Ｈａｓｓｅｌｂｌａｄ反復法による粒子数推定処理を示す。図２３はＥＭアルゴリズムによる処理手順を示す。

まず、初期値の設定（処理１９Ａ）を行った後、確率密度関数に基づく個数計算が順次実行される（処理１９Ｂ）（ステップＳ１２ａ、Ｓ１２ｂ）。個数計算の反復は（１９Ｃ）に示す収束条件を満たすまで実行される（ステップＳ１２ｃ）。ＥＭアルゴリズムの実行結果（粒子種ごとの推定個数データ）はＲＡＭ２３の所定エリアに格納される（ステップＳ１２ｄ）。

粒子数推定処理により得られた粒子種ごとの推定個数データは、ステップ４において、粒子種別の個数分布データに編集され、表示指定に応じて表示手段７にヒストグラム表示出力可能になる。図２７では省略しているが、本実施形態においては、分散図出力の指定を受けた場合には、特徴量データによる粒子種別の分散図を表示出力可能になっている。

図３７は、本実施形態の粒子種個数分析により分析した結果の一例を示す。同図（２４Ａ）および（２５Ｂ）は分析対象の粒子種である大腸菌、枯草菌の顕微鏡拡大写真である。（２４Ｃ）および（２５Ｄ）は特徴量としてパルス波高およびパルス尖度を注して粒子数推定処理の実行により得られた粒子種ごとの推定個数データのヒストグラム、分散図を示す。
＜特徴量による粒子種個数の分析精度の検証１について＞
本発明者らは、上記実施例の大腸菌と枯草菌の計測電流データを用いて下記の評価条件下で粒子種個数の分析性能の検証１を行った。

検証１の評価条件は以下の通りである。

（１）大腸菌と枯草菌の１０００ｋＨｚ実験測定データで評価を行う。

（２）特徴量には、波長Δｔ、波高ｈ、ピーク位置比ｒ、ピーク尖度ｋの４つの第１類型の特徴量を算出して用いる。

（３）各特徴量の組合せについて個数推定処理を実施する。

（４）大腸菌と枯草菌の実測データをランダムに学習用とテスト用に分けて推定評価する。この推定評価を１０回反復して実施し、それらの平均精度と標準偏差を算出する。この場合、実際に近い精度を評価する交差検定法（ｃｒｏｓｓ ｖａｌｉｄａｔｉｏｎ）により行う。

（５）検証粒子（大腸菌と枯草菌）の実測データの一部を個別に個数分析し、残りを所定の混合比δによりランダムに混合して検証用として、個数分析結果を比較する。ランダムデータ混合用データ混合プログラムをＲＯＭ３に格納しＰＣ１を利用してデータのランダム混合を実行して、そのランダム混合したデータに対する個数推定を行う。すなわち、図３２のステップＳ３２における行列データには、データ混合プログラムにより作成されたＮ行ｋ列のランダム置換行列データを使用する。混合比δには、大腸菌の混合率として１０、２０、３０、３５、４０、４５、５０％の７種類を使用する。ＢＬ推定用のパラメータ(調整因子）ｍ、ｋ、αの値はそれぞれ、１０００００、４００、６を使用し、確率密度関数の推定用の標準偏差係数ｃには０．１を設定する。粒子種個数推定時の収束条件αは０．１に設定する。なお、評価に使用した上記調整因子の値には、図２５で示したシミュレーション例と同様にしてより厳密な調整を行って得られた値が使用されている。

図３８の（２５Ａ）および（２５Ｂ）は、特徴量としてパルス波長、波高を使用した検証例と、特徴量としてパルス波長、ピーク位置比を使用した検証例の各推定結果データを示す。

本検証により得られた全パルスの数は大腸菌で１４６個、枯草菌で４０５個であった。

図３９の（２６Ａ）および（２６Ｂ）は特徴量としてピーク付近波形の広がり、パルス波長を使用した検証例と、特徴量としてピーク付近波形の広がり、波高を使用した検証例の各推定結果データを示す。

粒子種別個数の評価は、図４０の（２７Ｂ）に示す数式で表わされる「重み付き平均相対誤差」により行うことができる。「重み付き平均相対誤差」は各粒径の相対誤差にその粒径の真の個数割合を掛けたものを、全粒径について足した数値である。

図４０の（２７Ａ）は特徴量として尖度と、パルス波高を使用した場合における個数推定結果を示す。

図４１の（２８Ａ）および（２８Ｂ）は特徴量としてパルス波長、パルス波高を使用した場合における各混合比δ別の個数推定結果と、特徴量としてパルス波長、ピーク位置比を使用した場合における各混合比δ別の個数推定結果を示す。

図４２の（２９Ａ）〜（２９Ｄ）は大腸菌と枯草菌の混合比をそれぞれ、１：１０、２：１０、３：１０、３５：１００とした場合における各個数推定結果を示すヒストグラムである。

図４３の（３０Ａ）〜（３０Ｃ）は、大腸菌と枯草菌の混合比をそれぞれ、４：１０、４５：１００、１：２とした場合における各個数推定結果を示すヒストグラムである。

図４４の（３１Ａ）および（３１Ｂ）は、特徴量としてパルス波長、パルス波高を使用した場合における各粒子の散布状態を合成した図である。

図４５の（３２Ａ）、（３２Ｂ）および（３２Ｃ）は、特徴量としてピーク付近波形の広がり、パルス波長を使用した場合、特徴量としてピーク付近波形の広がり、ピーク位置比を使用した場合、ピーク付近波形の広がり、パルス波高を使用した場合における各粒子の散布状態を合成した図である。

上記の性能評価実験から、以下の評価結果が得られた。

（１）図４４および図４５のデータ散布図では、４つの特徴量に関して、大腸菌と枯草菌の特徴は大きく重なるが、明らかな相違があることが認められる。

（２）図４０（２７Ａ）等に示す種別個数分布の推定結果からは、この評価検証の特徴量の中でパルス波高とピーク尖度の特徴量を組み合わせた場合が最も精度がよく、重み付き平均相対誤差の評価で４〜１２％の分析精度を得ることができる。上記実施形態においては４種類全部の特徴量を抽出しているが、上記検証結果を踏まえて一部の特徴量（例えば、パルス波高とピーク尖度）だけを抽出して個数分析するようにしてもよい。
＜特徴量による粒子種個数の分析精度の検証２について＞
上記実施例の大腸菌と枯草菌の計測電流データを用いて、検証１とは別の粒子種個数の分析性能の検証２を行った。検証２においては、検証１とは異なり、第１類型および第２類型の特徴量（（１）〜（１３）の１３種類）を算出して用い、これらの組合せに係る特徴量とサンプリングデータ数との関連性および各組合せの分析性能を検証した。

図５５の（４２Ａ）および（４２Ｂ）は、それぞれ、全データのうち、１ＭＨｚ、５００ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図４３の（４３Ａ）および（４３Ｂ）は、それぞれ、全データのうち、２５０ｋＨｚ、１２５ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図４４の（４４Ａ）および（４４Ｂ）は、それぞれ、全データのうち、６３ｋＨｚ、３２ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図５８の（４５Ａ）および（４５Ｂ）は、それぞれ、全データのうち、１６ｋＨｚ、８ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図５９は、４ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。これらの表中の各組合せごとの推定評価結果は、検証１の（４）と同様に交差検定法により得られた、上側に記載の平均精度と、下側に括弧書きで示した標準偏差を表す。表中の慣性Ｉ、慣性Ｉ（規格化）、慣性Ｉ ｗ、慣性Ｉ ｗｖ、慣性Ｉ ｗ（規格化）、慣性Ｉ ｗｖ（規格化）は、それぞれ、（８）の時間慣性モーメント、（９）の規格化された時間慣性モーメント、（１０）の波幅平均値慣性モーメント、（１２）の波幅分散慣性モーメント、（１１）の規格化された波幅平均値慣性モーメント、（１３）の規格化された波幅分散慣性モーメントの特徴量を示す。

図６０は、全サンプリングデータにおける各特徴量組合せに関する推定評価結果を示す。図６１は、全データのうち１ＭＨｚ〜１２５ｋＨｚでの高密度サンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図６２は、全データのうち６３ｋＨｚ〜４ｋＨｚでの低密度サンプリングしたときの各特徴量組合せに関する推定評価結果を示す。

図６３は、全サンプリングデータを使用したとき（５０Ａ）および高密度にサンプリングしたとき（５０Ｂ）に高い個数推定精度が得られる上位５種の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフである。図６３における上位５種の特徴量の組合せは、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、ピーク位置比ｒ−慣性Ｉ、俯角θ−慣性Ｉ、慣性Ｉ−慣性Ｉ ｗ（規格化）である。

図６４は、低密度にサンプリングしたときに高い個数推定精度が得られる上位５種の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフ（５１Ａ）と、全サンプリングデータを使用したときの４種類の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフ（５１Ｂ）である。図６３および図６４の縦軸の値は、５０回の交差検定を行って得られた重み付き平均相対誤差の平均値である。５１Ａにおける上位５種の特徴量の組合せは、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、ピーク位置比ｒ−面積ｍ、俯角θ−面積ｍ、面積ｍ−慣性Ｉ ｗｖ（規格化）である。５１Ｂにおける４種類の特徴量の組合せは、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、尖度ｋ−波高｜ｈ｜、尖度ｋ−ピーク位置比ｒである。

検証２から得られた結果は、以下の通りである。
（Ｒ１）図６０および図６３に示すように、全サンプリングデータを使用したときには、上位５種の組合せ、すなわち、波長Δｔ−慣性Ｉ、波長Δｔ−面積ｍ、ピーク位置比ｒ−慣性Ｉ、俯角θ−慣性Ｉ、慣性Ｉ−慣性Ｉ ｗ（規格化）の特徴量の場合、高い個数推定精度を得ることができる。これらの特徴量の組合せによる個数推定精度（重み付き平均相対誤差）は、例えば、波長Δｔ−慣性Ｉで２５０〜１０００ｋＨｚのサンプリング領域において約９〜１０％、波長Δｔ−面積ｍで１２５〜２５０ｋＨｚのサンプリング領域において約９〜１０％、波長Δｔ−慣性Ｉで１６〜６３ｋＨｚのサンプリング領域において約１３〜１５％である。
（Ｒ２）図６１に示すように、全サンプリングデータより少ないが高密度なサンプリングデータを使用したときに高い個数推定精度が得られる特徴量は、上位５種の組合せで示せば、波長Δｔ−慣性Ｉ、波長Δｔ−面積ｍ、ピーク位置比ｒ−慣性Ｉ、慣性Ｉ−慣性Ｉ ｗ、俯角θ−慣性Ｉの５種である。これらの特徴量の組合せによる個数推定精度（重み付き平均相対誤差）は、例えば、波長Δｔ−慣性Ｉで２５０〜１０００ｋＨｚのサンプリング領域において約９〜１０％、波長Δｔ−面積ｍで１２５〜２５０ｋＨｚのサンプリング領域において約９〜１０％、波長Δｔ−慣性Ｉで１６〜６３ｋＨｚのサンプリング領域において約１３〜１５％である。
（Ｒ３）図６２に示すように、高密度サンプリングデータと比べてさらに少ない低密度サンプリングデータを使用した場合に高い個数推定精度が得られる特徴量は、上位５種の組合せで示せば、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、俯角θ−面積ｍ、面積ｍ−慣性Ｉ ｗｖ（規格化）、ピーク位置比ｒ−面積ｍの５種である。これらの特徴量の組合せによる個数推定精度（重み付き平均相対誤差）は、波長Δｔ−慣性Ｉで２５０〜１０００ｋＨｚのサンプリング領域において約９〜１０％、波長Δｔ−面積ｍで１２５〜２５０ｋＨｚのサンプリング領域において約９〜１０％、、波長Δｔ−慣性Ｉで１６〜６３ｋＨｚのサンプリング領域において約１３〜１６％である。
（Ｒ４）（Ｒ１）〜（Ｒ３）からわかるように、第１類型と第２類型の特徴量の組合せを使用しても高精度の個数推定を行うことができる。さらに、本発明に係る個数分析方法によれば、サンプリング数が十分に多くなくても、所定のサンプリング数が得られれば十分あるときと同程度の精度で個数分析を行うことができる。例えば、検証１で調べた尖度ｋとピーク位置比ｒとの組合せでは、１２％の最大誤差を生じていたが、例えば、波長Δｔ−慣性Ｉの特徴量による場合には、全データを使用しなくとも１ＭＨｚ〜１２５ｋＨｚでの高密度サンプリングデータを使用して、つまり部分的なデータであっても個数推定処理を約９％の高精度で行うことができる。したがって、本実施形態に係る個数分析機能は、定常的な個数分析にとどまらず、例えば、緊急性を要する検疫検査や医療現場において、菌類等の粒子有無や個数の判別に即応的実施に好適な検査ツールとして使用することができる。
＜個数分析処理時間の検証３について＞
個数推定には、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する所要計算時間がかかるので、この所要計算時間とサンプリング周波数との関係について特徴量の比較検討を検証３で検証した。検証３の比較検討例には、図６４の（５１Ｂ）に示した、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、尖度ｋ−波高｜ｈ｜、尖度ｋ−ピーク位置比ｒの４種類の特徴量の組合せを使用した。これらの組合せは、他の組合せと比較して交差検定精度の良い組合せである。個数分析の計算に要する時間には、特徴量作成に要する時間と、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する計算時間とが含まれるので、特徴量作成に要する計算時間ＣＴ１、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する計算時間ＣＴ２およびそれらの合計計算時間ＣＴ３（＝ＣＴ１＋ＣＴ２）について比較検討を行った。この場合も、それぞれの所要計算時間は、５０回の交差検定を行って得られた各計算時間の平均値である。

図６５は、４種類の各特徴量組合せに対する合計計算時間ＣＴ３を示すサンプリング周波数（ｋＨｚ）−所要計算時間（秒）のグラフ（５２Ａ）と、各特徴量組合せに対する特徴量作成に要する計算時間ＣＴ１を示すサンプリング周波数（ｋＨｚ）−所要計算時間（秒）のグラフ（５２Ｂ）である。図６６は、各特徴量組合せに対する計算時間ＣＴ２を示すサンプリング周波数−所要計算時間（秒）のグラフである。

（５２Ａ）に示すように、波長Δｔ−面積ｍと波長Δｔ−慣性Ｉの特徴量組合せＧ１は、ほぼ同じ合計計算時間になっており、尖度ｋ−波高｜ｈ｜と尖度ｋ−ピーク位置比ｒの特徴量組合せＧ２は、ほぼ同じ合計計算時間になっている。（５２Ｂ）に示すように、特徴量組合せＧ１のそれぞれの特徴量作成に要する計算時間は同じであり、特徴量組合せＧ２のそれぞれの特徴量作成に要する計算時間は同じになっている。図５３に示すように、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する時間は、特徴量組合せＧ１、Ｇ２のいずれにおいても、１ＭＨｚ〜１６ｋＨｚでのサンプリング領域において約３，５秒以下の短時間で処理可能になっている。

検証３の特徴量組合せＧ１、Ｇ２の比較結果から明らかに、第１類型と第２類型の同一類型との組合せであっても異なる混合組合せであっても特徴量を使用して所要計算時間の短縮化を図ることができる。したがって、本実施形態に係る個数分析機能によれば、定常的な個数分析にとどまらず、例えば、緊急性を要する検疫検査や医療現場において、菌類等の粒子有無や個数の判別処理を迅速に行うことができる。

以上の性能評価からわかるように、ナノポアデバイス８により検出した検出信号のデータ群をベースにして、個数導出手段である粒子種分布推定プログラムの実行によって、該検出信号として得られた粒子通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、粒子種別の個数を導出することができる。したがって、ＰＣ１個数分析機能を用いることによって、例えば、細菌や微小粒子状物質等の分析物種別に応じた個数ないし個数分布を高精度に分析することができ、個数分析検査における簡易化および低コスト化を実現することができる。ナノポアデバイス８による検出信号を個数分析装置に直接的に取り込んでデータ保存可能にすることにより、検査・分析を統合した粒子種統合分析システムを構築するようにしてもよい。

特徴量に基づくデータ群から確率密度推定を行い、粒子種別の個数を導出した結果を、出力手段である表示手段７に表示出力あるいはプリンタにプリント出力することができる。したがって、本実施形態によれば、高精度の導出結果（粒子個数、粒子個数分布、推定精度等）を、例えば、ヒストグラムや散布図の出力形態で認知可能に即応的に報知することができるので、例えば、迅速な対応を要する医療現場や検疫場における有用な検査ツールとして本実施形態に係る個数分析機能を使用することができる。

本発明は、識別処理プログラムを搭載した特定のＰＣ等のコンピュータ端末に限らず、該識別処理プログラムの一部ないし全部を記憶した識別分析用記憶媒体に適用することができる。すなわち、所定のコンピュータ端末に該識別分析用記憶媒体に記憶した識別分析プログラムをインストールして所望のコンピュータに個数分析動作させることができるので、簡便かつ安価に個数分析を行うことができる。本発明の適用可能な記憶媒体には、フレキシブルディスク、磁気ディスク、光ディスク、ＣＤ、ＭＯ、ＤＶＤ、ハードディスク、モバイル端末等、コンピュータにより読み取り可能な記憶媒体のいずれかを選択して使用することができる。

図６９は、本実施形態に係る分類分析処理を示す。

図６７のコンピュータ解析部１ａは、本実施形態のＰＣ１に対応する。分析処理の準備作業として、入力処理（ステップＳ１００）において、不適合データの除去処理、特徴量の指定、既知データおよび被分析データのＰＣ１への入力が行われる。特徴量は、前述の（１）〜（１５）に示した第１類型および第２類型の一部または全部あるいは１以上の組合せの特徴量を該入力処理であらかじめ指定しておくことができる。例えば、大腸菌Ｅｃおよび枯草菌Ｂｓを、粒子種別が特定された分析物（特定分析物）とする場合、これらの特定分析物の個々につき、ナノポアデバイス８ａによる計測を行って、それぞれのパルス状信号のデータが既知データとしてＰＣ１に入力され、入力データは、ＲＡＭ４の既知データ記憶用メモリエリアに格納される。特定分析物の含有状態が不明の被分析対象に対するナノポアデバイス８ａによる計測を行って得られたパルス状信号のデータが分析データとしてＰＣ１に入力され、入力データは、ＲＡＭ４の分析データ記憶用メモリエリアに格納される。

分類分析処理が起動操作により起動されると、既知データの入力有無が判別される（ステップＳ１１０）。既知データの未入力の場合、表示手段７により既知データの入力を促すガイダンス表示が行われる。図６９において各種のガイダンス表示による報知処理ステップは省略している。既知データが入力されると、入力された既知データは、ＲＡＭ４の既知データ記憶用メモリエリアに格納され、特徴量の作成に供される（ステップＳ１００、Ｓ１０１）。

既知データ入力がある場合、特徴量の指定があるか否かが判断される（ステップＳ１１０、Ｓ１１１）。特徴量の指定がある場合、ＲＡＭ４の既知データによる特徴量記憶用データファイルＤＡから指定された特徴量のベクトル量データがＲＡＭ４の学習データ記憶エリアに取り込まれる（ステップＳ１１３）。特徴量の指定がない場合は、ＲＡＭ４の既知データによる特徴量記憶用データファイルＤＡからすべての特徴量のベクトル量データがＲＡＭ４の学習データ記憶エリアに取り込まれる（ステップＳ１１２）。

ついで、分析データの入力有無が判別される（ステップＳ１１４）。分析データ入力がない場合、表示手段７により分析データ入力を促すガイダンス表示が行われる。分析データが入力されると、取得した分析データはＲＡＭ４の分析データ記憶用メモリエリアに格納される（ステップＳ１００）。分析データが入力されると、記述のように、分析データに関する特徴量が作成され、ＲＡＭ４に記憶される（ステップＳ１０１）。分析データの入力がある場合、ＲＡＭ４の分析データによる特徴量記憶用データファイルＤＢから特徴量のベクトル量データがＲＡＭ４の変数データ記憶エリアに取り込まれる（ステップＳ１１５）。

既知データおよび分析データの入力済みで特徴量の取得状態において、分類分析の実行を促すガイダンス表示が行われる。該ガイダンス表示にしたがって所定の指示操作を行うことにより、分類分析プログラムが起動されて機械学習による分類分析の実行処理が行われる（ステップＳ１１６）。本実施形態において、例えば、ランダムフォレスト法に基づくアルゴリズムで構成された機械学習による分類分析プログラムがＲＯＭ３にあらかじめ格納されている。既知データによる特徴量を学習データとし、分析データから得られる特徴量を変数にして、該分類分析プログラムを実行することによって、該被分析データにおける特定分析物に関する分類分析を行うことができる。該分類分析プログラムの実行に際しては、パルス波形を同一次元の数値ベクトルに変換し、各ベクトルがどのように異なっているかを判別することにより個別のパルスを識別して分類分析が行われる。

本発明に係る機械学習による分類分析手法には、ランダムフォレスト法に限らず、例えば、ｋ近傍法、ナイーブベイズ分類器、決定木、ニューラルネットワーク、サポートベクターマシン、バギング法、ブースティング法等の集団学習によるものを使用することができる。

機械学習による分類分析の実行処理が、分析データによる特徴量のすべてについて実行されると、分類分析処理が終了し、分類分析結果の出力処理が行われる（ステップＳ１１７）。出力処理において、種別の分からない分析データの個々について、特定分析物の例示である大腸菌Ｅｃまたは枯草菌Ｂｓの通過に由来するものの割合である分類結果が表示手段７に表示可能になっている。出力可能な表示態様には、各分析データ毎の分類結果に限らず、分析物（例えば、大腸菌Ｅｃまたは枯草菌Ｂｓ）の該当総数、両者の該当比率等の表示態様を使用することができる。
＜分類分析処理の処理精度の検証＞
上記の分類分析処理の処理精度について、種々の機械学習による分析手法を適用して分類分析を試行して本実施形態による分類分析処理の精度を検証した。

図７０の（５７Ａ）は、同図（５７Ｂ）に示す分析試料を用いて、特徴量（Ｆｅａｔｕｒｅ）と機械学習による分析手法のアルゴリズム（以下、分類器という。）とを種々組み合わせた場合に、本発明に係る分類分析処理（図６９参照）を実行した評価結果を示す。

分析試料は、（５７Ｂ）に示すように、２種類の細菌種（大腸菌、枯草菌）である。各細菌種に対し、貫通孔１２の内径が４．５Φで、貫通孔１２の貫通距離（ポア深さ）が１５００ｍｍのマイクロ・ナノポアデバイス８を用いて通過波形を計測して得られたパルス状信号データを４２個（大腸菌の場合、計測パルスのすべてで、枯草菌の場合、計測パルス数２６５個のうちの４２個である。）を使用した。分類器の実行の際には、パルス状信号データのうち約９割を学習データとし、残りのデータを変数に振り分けた。

評価項目は、（５７Ａ）に示すように、Ｆ−尺度（Ｆ−Ｍｅａｓｕｒｅ）で表され、真陽性率（ＴＰＲａｔｅ）、偽陽性率（ＦＰＲａｔｅ）、適合率（Ｐｒｅｃｉｓｉｏｎ）、再現率（Ｒｅｃａｌｌ）、Ｆ値（ＦＭｅａｓｕｒｅ）、受信者操作特性曲線面積（ＲＯＣ（Ｒｅｃｅｉｖｅｒ Ｏｐｅｒａｔｉｎｇ Ｃｈａｒａｃｔｅｒｉｓｔｉｃ ）ｃｕｒｖｅ Ａｒｅａ）の項目からなる。

図７１は、Ｆ−尺度の説明図である。

Ｆ−尺度は、（５８Ａ）に示すように、２種類の細菌種の実数（大腸菌の実数：Ｐ、枯草菌の実数：Ｎ）に対し、各細菌種の予想値を割り付けた場合に、各組合せにおける真陽性（ＴＰ）、偽陽性（ＦＰ）、真陰性（ＦＮ）および偽陰性（ＴＮ）の総和が１であるとして、（５８Ｂ）に示すように、２ＴＰ／（２ＴＰ＋ＦＰ＋ＦＮ）で表される。

この検証において、アルゴリズムが異なる６７種類の分類器を用いて、各種特徴量ないし特徴量の組合せを用いて約４０００種のパターンに対して分類分析を試行した。この結果、６０種の特徴量の組合せに対して有意な分析結果が得られた。図７０の（５７Ａ）は、この検証で得られたＦ−尺度の優れた上位１０位内の分類結果を示す表である。

上位１０位内における特徴量には、（５７Ａ）に示すように、（１）〜（１１）、（１４）および（１５）の１３種類の特徴量を並べた１３次元の特徴量ベクトル（表中に「ｈｖ＆Ｆ」で略記）、波高ベクトル（表中に「ｈ」で略記）と（１０）の平均値ベクトル（表中に「ｗＶ」で略記）との組合せ（表中に「ｈ＆ｗＶ」で略記）、波高ベクトルと（１１）の規格化された平均値ベクトル（表中に「ｗＮｒｍｄＶ」で略記）との組合せ（表中に「ｈ＆ｗＮｒｍｄＶ」で略記）が含まれている。（５７Ａ）で最も優れた分類精度によるものは、特徴量としてｈ＆ｗＶの組合せを使用した、ランダムホレスト法による分類器（「４ｍｅｔａ．Ｒａｎｄｏｍ Ｃｏｍｍｉｔｅｅ」）による場合であり、その分類精度約９８．９％の高精度であった。

本発明は、分類分析プログラムを搭載した特定のＰＣ等のコンピュータ端末に限らず、該分類分析プログラムの一部ないし全部を記憶した分類分析用記憶媒体に適用することができる。すなわち、所定のコンピュータ端末に該分類分析用記憶媒体に記憶した分類分析プログラムをインストールして所望のコンピュータに分類分析動作させることができるので、簡便かつ安価に分析析を行うことができる。本発明の適用可能な記憶媒体には、フレキシブルディスク、磁気ディスク、光ディスク、ＣＤ、ＭＯ、ＤＶＤ、ハードディスク、モバイル端末等、コンピュータにより読み取り可能な記憶媒体のいずれかを選択して使用することができる。

尚、本発明は上記実施形態に限定されるものではなく、本発明の技術的思想を逸脱しない範囲における種々変形例、設計変更などをその技術的範囲内に包含するものであることは云うまでもない。

本発明によれば、高精度の不適合データの識別と分類分析を行えるので、例えば、ＤＮＡ記憶媒体の情報圧縮技術や人工塩基対を用いた医薬品創薬、あるいは、計測試料に混入する微細な塵埃、あるいは体液などに含まれる分析物質を計測対象とする場合における、赤血球、白血球、血小板等の微小物質などに起因する不適合データの識別・除去技術等の分野に広範囲に応用発展することができる。特には、本発明は、ＤＮＡやＲＮＡと夾雑物を含む試料を分析対象にして、例えば、下水中のＤＮＡ含有分析を行ってウイルスの発生を検知する検知技術におけるデータ分析に適用することができる。

１ パーソナルコンピュータ
２ ＣＰＵ
３ ＲＯＭ
４ ＲＡＭ
５ データファイル記憶部
６ 入力手段
７ 表示手段
８ マイクロ・ナノポアデバイス
９ チャンバー
１０ 基板
１１ 隔壁
１２ 貫通孔
１３ 電極
１４ 電極
１５ 電源
１６ 増幅器
１７ オペアンプ
１８ 凹部
１９ 帰還抵抗
２０ 電圧計
２１ 検体
２２ 大腸菌
２３ 枯草菌
２４ 電解質溶液
ＭＳ 計測空間
Ｄ１ 電極
Ｄ２ 電極
ＭＥ 電流計測器

Claims (8)

1. 計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された不適合データの識別をコンピュータ制御プログラムの実行によって行う識別方法であって、
   前記コンピュータ制御プログラムは、正例集合の正例データと、正例負例のいずれかが不明である未知集合の未知データとから正負例を分類する分類器を学習する機械学習を用いた識別分析プログラムを有し、
   前記計測空間に分析物を含まない試料を導入して計測する第１計測条件の下で得られるパルス状信号の第１種データと、前記計測空間に分析物を含む試料を導入して計測する第２計測条件の下で得られるパルス状信号の第２種データとを記憶する記憶手段を有し、
   前記第１種データを前記正例データとし、前記第２種データを前記未知データとして、前記識別分析プログラムを実行することによって、前記第２種データに含まれる前記不適合データを識別することを特徴とする識別方法。
2. 請求項１に記載の識別方法により識別した不適合データを記憶する不適合データ記憶手段を有し、
   計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された前記不適合データを取り除いた被分析データの分類分析をコンピュータ制御プログラムの実行によって行う分類分析方法であって、
   前記コンピュータ制御プログラムは、機械学習を用いた分類分析を行う分類分析プログラムを有し、
   前記パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、
   あらかじめ求めた特徴量を前記機械学習のための学習データとし、前記不適合データを取り除いた被分析データのパルス状信号から得られる特徴量を変数にして、前記分類分析プログラムを実行することによって前記分析物に関する分類分析を行うことを特徴とする分類分析方法。
3. 前記特徴量は、
   所定の時間幅内における波形の波高値、
   パルス波長ｔ_a、
   パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
   該波形の鋭さを表す尖度、
   パルス開始からパルスピークに至る傾きを表す俯角、
   波形を所定の時間毎に区分した時間区分面積の総和を表す面積、
   パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比、
   パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
   前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
   波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値をベクトルの成分とする平均値ベクトル、
   前記平均値ベクトルに対し波長が基準値になるように規格化した場合の規格化された平均値ベクトル、
   波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
   前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
   波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
   前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメント、
   のいずれか１または２以上である請求項２に記載の分類分析方法。
4. 計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された不適合データの識別をコンピュータ制御プログラムの実行によって行う識別装置であって、
   前記コンピュータ制御プログラムは、正例集合の正例データと、正例負例のいずれかが不明である未知集合の未知データとから正負例を分類する分類器を学習する機械学習を用いた識別分析プログラムを有し、
   前記計測空間に分析物を含まない試料を導入して計測する第１計測条件の下で得られるパルス状信号の第１種データと、前記計測空間に分析物を含む試料を導入して計測する第２計測条件の下で得られるパルス状信号の第２種データとを記憶する記憶手段を有し、
   前記第１種データを前記正例データとし、前記第２種データを前記未知データとして、前記識別分析プログラムを実行することによって、前記第２種データに含まれる前記不適合データを識別することを特徴とする識別装置。
5. 請求項４に記載の識別装置により識別した不適合データを記憶する不適合データ記憶手段を有し、
   計測空間に分析物を含む試料を導入して検出したパルス状信号のデータから、分析物以外の要素に起因して検出された前記不適合データを取り除いた被分析データの分類分析をコンピュータ制御プログラムの実行によって行う分類分析装置であって、
   前記コンピュータ制御プログラムは、機械学習を用いた分類分析を行う分類分析プログラムを有し、
   前記パルス状信号の波形形態の特徴を示す特徴量をあらかじめ求め、
   あらかじめ求めた特徴量を前記機械学習のための学習データとし、前記不適合データを取り除いた被分析データのパルス状信号から得られる特徴量を変数にして、前記分類分析プログラムを実行することによって前記分析物に関する分類分析を行うことを特徴とする分類分析装置。
6. 前記特徴量は、
   所定の時間幅内における波形の波高値、
   パルス波長ｔ_a、
   パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
   該波形の鋭さを表す尖度、
   パルス開始からパルスピークに至る傾きを表す俯角、
   波形を所定の時間毎に区分した時間区分面積の総和を表す面積、
   パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比、
   パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
   前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
   波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値をベクトルの成分とする平均値ベクトル、
   前記平均値ベクトルに対し波長が基準値になるように規格化した場合の規格化された平均値ベクトル、
   波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
   前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
   波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
   前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメント、
   のいずれか１または２以上である請求項５に記載の分類分析装置。
7. 請求項１に記載のコンピュータ制御プログラムを記憶したことを特徴とする識別用記憶媒体。
8. 請求項２に記載のコンピュータ制御プログラムを記憶したことを特徴とする分類分析用記憶媒体。

Images