"測定されたデータを受け入れる受入手段と、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、を含むデータ処理装置。請求項1記載のデータ処理装置であって、前記受入手段は、測定されたデータとして時間変化するデータを受け入れ、前記推定手段は、時刻tにおけるデータの値mtから、時刻tの信号の値を含む情報の確率密度を与える演算式として、当該情報がxtのときに測定されるデータの値mtの確率密度である事後確率と、前記信号の値を含む情報の初期値x0、及び当該情報の確率密度の時間発展を表す模式的な演算式の積により表現される事前確率とを用いたベイズ推定により前記信号パラメータを推定する推定手段であるデータ処理装置。請求項2記載のデータ処理装置であって、複数の演算装置要素を含み、前記推定手段は、前記ベイズ推定における、事前確率と事後確率との積を被積分関数とする積分を行うにあたり、積分範囲内からランダムに選択される積分変数の値を代入したときの前記被積分関数の値の演算を、複数の演算装置要素に分散的に行わせ、モンテカルロ法により前記積分を遂行するデータ処理装置。請求項1から3のいずれか一項に記載のデータ処理装置であって、前記ノイズパラメータを決定する手段は、前記ノイズ演算式に含まれる少なくとも一つのノイズパラメータについては、データの実際の測定に先立って、前記受入手段が受け入れるデータを用いて定められる基準に基づいて決定するデータ処理装置。コンピュータを、測定されたデータを受け入れる受入手段と、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、として機能させるプログラム。しかしながら、上記加算平均を用いるノイズの除去方法では、測定回数をなるべく多くしなければならない。ところが測定の対象によっては多数回の測定が困難である場合もある。また脳磁計等の場合、測定回数を増大させるほど被験者の負担が大きくなってしまう。本発明は上記実情に鑑みて為されたもので、単一の測定結果からでも信号の推定を可能とするデータ処理装置を提供することを、その目的の一つとする。上記従来例の問題点を解決するための本発明は、データ処理装置であって、測定されたデータを受け入れる受入手段と、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、前記ノイズパラメータを決定する手段と、前記測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、前記受入手段が受け入れる測定されたデータと、前記ノイズ演算式とに基づいて、前記信号パラメータを推定する推定手段と、前記推定された信号パラメータと、前記信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する手段と、を含むこととしたものである。本発明によれば、ノイズ及び信号を、パラメータを含んだ模式的な演算式として表現し、そのパラメータを推定するようにしたことで、単一の測定結果からでも信号の推定を可能とした。本発明の実施の形態に係るデータ処理装置を図面を参照しつつ説明する。本実施の形態のデータ処理装置1は、図1に例示するように、制御部11、記憶部12、操作部13、出力部14、及び入力インタフェース部15を含んで構成される。以下では、具体的に、脳磁計MEGの出力するデータ時間変化する一連のデータを入力インタフェース部15を介して受け入れる例について説明するが、本実施の形態のデータ処理装置1は、脳磁計のデータを処理するものに限られない。脳磁計のデータには、被験者のいない状態で脳磁計が出力するデータノイズのみのデータと、被験者が安静な状態にあるときに脳磁計が出力するデータ外的刺激を受けていない状態:安静時データと呼ぶと、被験者に対して外的刺激を与えたときに脳磁計が出力するデータノイズ、自発脳磁場及び、刺激による誘発脳磁場が累算されたデータ:刺激時データと呼ぶとがある。制御部11は、CPUCentralProcessingUnit等のプログラム制御デバイスであり、記憶部12に格納されたプログラムに従って動作する。この制御部11は、入力インタフェース部15から入力される、外部の測定装置ここでは例えば脳磁計とするにより測定されたデータを受け入れる。また、この制御部11は、少なくとも一つのノイズパラメータを含み、ノイズを模式的に表すノイズ演算式について、過去に得られたデータを用いて定められる基準に基づいて、ノイズパラメータを決定する。さらに制御部11は、測定されたデータに含まれる信号を模式的に表す信号演算式であって、少なくとも一つの信号パラメータを含む信号演算式と、測定されたデータと、ノイズ演算式とに基づいて、信号パラメータを推定し、当該推定した信号パラメータと、信号演算式とを用いて、測定されたデータに含まれる信号の推定結果を生成し、出力する。この制御部11の詳しい処理の内容は後に説明する。記憶部12は、メモリデバイス等であり、制御部11によって実行されるプログラムを保持している。このプログラムは、例えばDVD-ROMDigitalVersatileDiscReadOnlyMemory等のコンピュータ可読な記録媒体に格納されて提供され、この記憶部12に複写されたものであってもよい。また、この記憶部12は制御部11のワークメモリとしても動作する。操作部13は、キーボードやマウスなどを含み、利用者の指示操作を受け入れて、制御部11に当該指示操作の内容を出力する。出力部14は、ディスプレイや、プリンタ、ネットワークインタフェース等であり、制御部11から入力される指示に従って、情報を出力する。入力インタフェース部15は、USBUniversalSerialBusやパラレルポート等のデータの入力を受け入れるデバイスであり、ここでの例では脳磁計に接続されて、脳磁計が出力するデータを受け入れ、当該受け入れたデータを制御部11に出力する。本実施の形態では、脳磁計から入力されるデータは時々刻々と変化する時間変化する。そして時刻tにおけるデータmtが、測定の対象たる信号stとノイズとの和により表現されているものとする。すなわち、ここにノイズは、要因ごとに複数の成分の和として表現されており、1式ではτtがベースラインドリフト成分、qtが交流電源成分、etがベースラインドリフトや交流電源成分以外の定常ノイズ成分、ntは、このようにモデル化したことによるデータからの差を表すモデル化誤差成分であるものとしている。また本実施形態では、不定期に生じる測定値のとび例えばグリッチノイズ、眼球運動や心拍等の生体から生じるノイズはないものとする。ノイズの各成分は、具体的に次のような形の演算式で表されるものとする。ベースラインドリフト成分τtの時間発展を、平滑化事前分布を用いた確率的トレンドモデルを採用し、と表す。ここでvτtはトレンドの変化を起こす外乱であり、平均ゼロ分散στのガウス分布と仮定しておく。なお、この2式、及び以下の説明において時刻tjは、tから予め定めた単位時間例えばTとするを用いてj×Tだけ後の時刻を意味するものである。さらに1式におけるqtの時間発展を、サンプリング周波数と交流周波数とから決まる自己回帰モデルARモデルを用いて、と表す。ここでも同様に、vqtは平均ゼロ分散σqのガウス分布と仮定する。またΔTはサンプリング周波数と交流電源の周波数から決定される、離散化された周期を表す。定常ノイズ成分etの時間発展は次数lのARモデルを仮定する。すなわち、とする。これにおいても、vetは平均ゼロ分散σeのガウス分布と仮定する。モデル化誤差成分ntは、平均ゼロ分散σnのガウス分布と仮定しておく。本実施の形態の制御部11は、ノイズ演算式として、1から4式の和、を用いる。このノイズ演算式に含まれるノイズパラメータを成分とするベクトルθモデルパラメータθベクトルは、従って、となる。本実施の形態の制御部11は、機能的には、図2に例示するように、データ記録部21、モデル学習部22、ノイズパラメータ管理部23、信号推定部24、及び信号出力部25を含んで構成される。ここにデータ記録部21は、入力インタフェース部15を介して受け入れられるデータを記憶部12に蓄積して保持する。本実施の形態では、このデータ記録部21は、操作部13から入力される指示に従い、受け入れたデータを、ノイズのみのデータ、ノイズと自発脳磁場とを含んだデータ、ノイズ、自発脳磁場、及び誘発脳磁場を含んだデータのいずれかに分類して記憶部12に格納する。モデル学習部22は、データ記録部21によって記憶部12に蓄積されたデータを参照し、過去に得られたデータ被験者のいない状態で得られたノイズのみのデータを用いて定められる基準に基づいて、6式のモデルパラメータθベクトルの成分各ノイズパラメータを決定する。具体的に、このモデル学習部22は、過去に得られたノイズのみの時系列のデータm1,m2,...,mtを記憶部12から読み出す。そしてベクトルθが与えられたときにデータがこの順に観測される確率pm1,m2,...,mt|θを最大とするベクトルθを求める。すなわち、現実の観測変数mtを生成する真の確率モデルは未知であるものの、当該真の確率モデルを、近似的に5式で与えられるものとし、真の確率モデルと近似的なモデルとの擬距離であるKLカルバックライブラーダイバージェンスKL情報量を最小とする5式のパラメータベクトルθの成分を求めることとするのである。例えば上記KL情報量を最小とする赤池情報量基準を用いてベクトルθを決定すればよい。このためモデル学習部22は、上記確率pm1,m2,...,mt|θの対数logpm1,m2,...,mt|θを対数尤度とし、赤池情報量基準AICを、AIC=-2×pm1,m2,...,mt|θ2×Nとして演算する。ここにNは自由に設定できるモデルパラメータの数を表し、ここではθの次元である、l4となる。モデル学習部22は、このAICの値を、複数の予め定めたθである値の組ごとに演算する。具体的にここで予め定めたθの値の組は、a1,a2,...,al,στ,σq,σe,σnで張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したAICの値が最小となるθを、上記複数の予め定めたθである値の組のうちから選択する。つまりモデル学習部22は、AICの値を最小とするベクトルθを、まずはグリッドサーチにより大域的に粗く探索する。モデル学習部22は、さらにグリッドサーチにより求めたベクトルθθ0とするを初期値として、非線形最適化により最適化されたベクトルを求める。ここで非線形最適化の一例としては、上記対数尤度logpm1,m2,...,mt|θをカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθ0を初期値として、対数尤度を最大とするベクトルθを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθを、とする。ここで終了条件は、更新前のθと更新後のθとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。ノイズパラメータ管理部23は、モデル学習部22が演算により求めたを、ノイズパラメータとして記憶部12に格納する。また、このノイズパラメータ管理部23は、信号推定部24からノイズパラメータの要求を受けて、記憶部12に格納したノイズパラメータを読み出し、信号推定部24に対して出力する。本実施の形態の信号推定部24は、自発脳磁場推定部24aと、誘発脳磁場推定部24bとを含んで構成される。自発脳磁場推定部24aは、記憶部12に格納されている、ノイズと自発脳磁場とを含んだ安静時データ安静状態にある被験者を測定した脳磁場のデータを読み出す。ここで自発脳磁場推定部24aは、自発脳磁場がによって表される確率モデルで近似的に表されるものとする。ここにrtは時刻tの自発脳磁場成分を意味し、δ波,θ波,α波,β波,lowγ波,highγ波の成分から構成されるものとし、各成分を、と表したものである。なお、ここでのθは、θ波に対応する自発脳磁場成分であることを意味するものとして便宜的に利用しているもので、パラメータのベクトルθとは異なる。また、これらの成分の各々はいずれも、次の準周期振動モデルで表現されるものとする。ここでωは、δ,θ,α,β,γl,γhのそれぞれをとる。また、fωは、成分ωδ,θ,α,β,γl,γhのそれぞれをとるの中心周波数であり、vtは、平均が0、分散がσω2であるガウス分布に属するとの条件、を満足する確率的入力である。以上より、自発脳磁場の成分を表すモデルパラメータであるθベクトルは、となる。ここにおいて、左辺のθpのθはθベクトルを表すθであり、右辺中、σθ2とあるθは、θ波に対応する自発脳磁場成分を表す便宜的なものである。自発脳磁場推定部24aは、記憶部12に格納されているノイズと自発脳磁場とを含んだ時系列の安静時データm1,m2,...,mtを記憶部12から読み出す。また自発脳磁場推定部24aは、ノイズパラメータ管理部23に対してノイズパラメータを要求し、この要求に応答してノイズパラメータ管理部23が出力するノイズパラメータと5式とを用いて、ノイズの時系列推定データν1,ν2,...を得る。自発脳磁場推定部24aは、ノイズと自発脳磁場とを含んだ時系列の安静時データm1,m2,...,mtのそれぞれ対応する時刻のデータからノイズの時系列推定データの対応する時刻のデータを差引きして、m′i=mi-νiを求め、ノイズがないと仮定される状態のデータ、m′1,m′2,...,m′tを得る。そして自発脳磁場推定部24aは、10式で表されるベクトルθpが与えられたときに安静時データが実際に測定された時系列順に観測される確率pm′1,m′2,...,m′t|θpを最大とするベクトルθpを求める。この場合も、現実の観測変数m′tを生成する真の確率モデルは未知であるものの、当該真の確率モデルを、近似的に7式で与えられるものとし、真の確率モデルと近似的なモデルとの擬距離であるKLカルバックライブラーダイバージェンスKL情報量を最小とする7式のパラメータ10式のベクトルθpの成分を求めることとするのである。例えば上記KL情報量を最小とする赤池情報量基準を用いて10式のベクトルθを決定すればよい。このため自発脳磁場推定部24aは、上記確率pm′1,m′2,...,m′t|θpの対数logpm′1,m′2,...,m′t|θpを対数尤度とし、赤池情報量基準AICを、AIC=-2×pm′1,m′2,...,m′t|θp2×Nとして演算する。ここにNは自由に設定できるモデルパラメータの数を表し、ここでは10式のベクトルθpの次元である、6となる。自発脳磁場推定部24aは、このAICの値を、複数の予め定めたθpである値の組ごとに演算する。具体的にここで予め定めたθpの値の組は、σω2ただしωはδ,θ,α,β,γl,γhのそれぞれをとるで張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したAICの値が最小となるθの値の組を、上記複数の予め定めたθpとしての値の組のうちから選択する。つまり自発脳磁場推定部24aは、AICの値を最小とするベクトルθpを、まずはグリッドサーチにより大域的に粗く探索する。自発脳磁場推定部24aは、さらにグリッドサーチにより求めたベクトルθpθp0とするを初期値として、非線形最適化により最適化された自発脳磁場のモデルのためのパラメータのベクトルを求める。ここで非線形最適化の一例としては、上記対数尤度logpm′1,m′2,...,m′t|θpをカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθp0を初期値として、対数尤度を最大とするベクトルθpを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθpを、とする。ここで終了条件は、更新前のθpと更新後のθpとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。自発脳磁場推定部24aは、こうして求めたを、誘発脳磁場推定部24bに出力する。誘発脳磁場推定部24bは、記憶部12に格納されている、ノイズと自発脳磁場とを含んだ刺激時データ刺激を受けた状態にある被験者を測定した脳磁場のデータを読み出す。ここで誘発脳磁場推定部24bは、誘発脳磁場を次の2次の線型モデルのインパルス応答で近似できるものとする。本実施の形態の誘発脳磁場推定部24bは、誘発脳磁場を、確率的入力を持つ2次の定常線形モデルとして11式で表現しておく。この11式においてモデルパラメータはシステム入力強度と、システム入力に対するモデルの応答を決めるシステムパラメータb1,b2である。ここでは単一の線型モデルのインパルス応答で誘発脳磁場を近似できるものとして、この刺激システム入力の確率分布を次のようにして与える。すなわち、誘発脳磁場推定部24bは、システムパラメータとしてb1,b2をもつシステムのインパルス応答hの時系列h1,h2,...と、得られている刺激時データm1,m2,...との相互相関関数を求める。そして誘発脳磁場推定部24bは、刺激の確率分布が、上記相互相関の絶対値が最大となる時間的ラグt誘発脳磁場が現れる時刻tを中心とする二項分布B2t,0.5であるとする。ここでは時間が離散化されているので二項分布とした。また、入力強度はラグtでの相互相関関数の値を平均値Iとするガウス分布NI,σI2としておく。誘発脳磁場推定部24bは、ノイズパラメータ管理部23に対してノイズパラメータを要求し、この要求に応答してノイズパラメータ管理部23が出力するノイズパラメータと5式とを用いて、ノイズの時系列推定データν1,ν2,...を得る。誘発脳磁場推定部24bは、また自発脳磁場推定部24aから入力される自発脳磁場に関するパラメータと7式とを用いて、自発脳磁場の時系列推定データg1,g2,...を得る。そして誘発脳磁場推定部24bは、刺激時データm1,m2,...,mtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、m′′i=mi-νi-giを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、m′′1,m′′2,...,m′′tを得る。11式では、モデルから直接にデータが得られるものではないので、上述のように刺激の確率分布を、ガウス分布を用いた演算式これには上記パラメータb1,b2が含まれるで表現しておき、この演算式における時刻tでの期待値μtと分散σ2とを用いた確率差分方程式により、刺激に対する推定される誘発脳磁場の時刻tでの信号st及び時刻t-1での信号st1を、その成分に含む状態ベクトルxの時間発展の演算式を得ておく。この演算式は、時刻t,t-1において推定される誘発脳磁場stから、時刻t1において推定される誘発脳磁場st1が得られる時間発展の確率運動のモデルpst1|st,st1を表す。誘発脳磁場推定部24bは、システムパラメータb1,b2を用いて表される信号st,st1を成分として含む状態ベクトルxtについて、初期状態px0これは誘発脳磁場が発生していないものとして既知であると、上記時間発展の確率とを用いて、時刻tで観測されるデータがm′′tであったとき、状態ベクトルがxtである確率pxt|m′′tを演算する。本実施の形態の誘発脳磁場推定部24bは、この確率pxt|m′′tを、次のように求める。すなわち、ベイズの定理より、が成立するので、これを変形して、を得る。ここでpm′′tは、と表すことができ、またpxtは、時間発展の確率を用いて、と表すことができるから、結局13式は、とすることができる。この16式は、時刻tで観測されるデータがm′′tであったとき、状態ベクトルがxtである確率密度pxt|m′′tが、観測のモデルpm′′t|xtと、運動のモデルpxt1|xtと、初期状態px0とから推定できることを表す。誘発脳磁場推定部24bは、システムパラメータb1,b2を含むベクトルθsが与えられたときに確率pxt|m′′tを最大にするベクトルθsを求める。具体的に誘発脳磁場推定部24bは、KL情報量を最小とする赤池情報量基準を用いてベクトルθsを決定する。このため誘発脳磁場推定部24bは、上記確率pxt|m′′tの対数logpxt|m′′tを対数尤度とし、赤池情報量基準AICを、AIC=-2×pxt|m′′t2×Nとして演算する。ここにNは自由に設定できるモデルパラメータの数を表し、ここではθsの次元である、2となる。誘発脳磁場推定部24bは、このAICの値を、複数の予め定めたθsである値の組ごとに演算する。具体的にここで予め定めたθsの値の組は、b1,b2で張られるベクトル空間上に設定された複数の格子点とすることができる。そして演算したAICの値が最小となるθsを、上記複数の予め定めたθsとしての値の組のうちから選択する。つまり誘発脳磁場推定部24bは、AICの値を最小とするベクトルθsを、まずはグリッドサーチにより大域的に粗く探索する。誘発脳磁場推定部24bは、さらにグリッドサーチにより求めたベクトルθsθs0とするを初期値として、非線形最適化により最適化されたベクトルを求める。ここで非線形最適化の一例としては、上記対数尤度logpxt|m′′tをカルマンフィルタにより求める方法がある。上記グリッドサーチにより求めたベクトルθs0を初期値として、対数尤度を最大とするベクトルθsを、例えば最急降下法により予め定めた終了条件を満足するまで繰り返し逐次的に更新して求め、当該求めたθsを、とする。ここで終了条件は、更新前のθsと更新後のθsとの距離が予め定めたしきい値を下回っているなどの広く使われる条件を適用すればよい。信号出力部25は、誘発脳磁場推定部24bが演算により得たと、11式で表されるモデルの式とを用いて、誘発脳磁場の時間変化波形を演算して出力する。さらに本実施の形態の制御部11は、以上の演算を、複数のセンサ脳磁計の場合はSQUIDsから得られた各データについて行い、センサごとにのパラメータを得てもよい。このときには、複数のセンサの一つを注目センサとするとき、注目センサに係るパラメータを得るためには、ノイズのみのデータや安静時データ、刺激時データは、当該注目センサによって得られたものを用いてもよい。そして制御部11は、複数のセンサの各々に対応して得られたパラメータと11式とを用い、各センサの位置における誘発脳磁場の時間変化波形を演算し、等磁界線図を作成して出力してもよい。この等磁界線図の作成処理は広く知られているので、ここでの詳細な説明を省略する。[モンテカルロ法による演算]ところで制御部11における16式の演算では、右辺分母の積分ベイズ推定における、事前確率と事後確率との積を被積分関数とする積分が解析的な解が得られないので困難になる。そこで本実施の形態の制御部11は、複数の演算素子を含み、次のように演算を行うこととしてもよい。ここで複数の演算素子は、CPUでなく、GPUGraphicsProcessingUnit等の演算素子を用いてもよい。このように複数の演算素子を利用できる場合、制御部11は、各演算素子に対し、積分範囲内からランダムにxtとxt1との値を選択させ、当該ランダムに選択した積分変数の値複数ある場合はそれぞれ積分範囲からランダムに選択するを被積分関数に代入して得た被積分関数値を、並列的に演算させる。そして制御部11は、各演算素子で得られた複数の被積分関数値を累算してモンテカルロ法により積分を遂行する。[ノイズのみのデータを得るタイミング]また、以上の説明において、ノイズパラメータの決定に用いた、ノイズのみのデータは、演算を行う時点より過去であれば、いつ得られたものであってもよかった。例えば刺激データを得る直前刺激データを実際に測定する被験者が入る直前に、測定に先立って得られたデータであってもよいし、刺激データを測定した後、被験者が退席してから得られたデータであってもよい。[自発脳磁場をベイズ推定で求める例]ここまでの説明において、自発脳磁場推定部24aは、10式で表されるベクトルθが与えられたときに、安静時データが、測定された時系列順に観測される確率pm′1,m′2,...,m′t|θを最大とするベクトルθを求めていた。しかしながら、本実施の形態はこれに限られない。すなわち、自発脳磁場推定部24aは、時刻tでの状態ベクトルxtとして、時刻t以前の自発脳磁場を表すデータgt,gt1,...と、時刻t以前の各ノイズを表すデータτt,qtqt1...qtT2et...etl1,nt...とを含めたxtτt,qtqt1...qtT2et...etl1,nt,gt,gt1,...Tを生成する。なお、肩のTは、転置を表す。ここで、各ノイズのデータや自発脳磁場のデータを、時刻tよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻t1における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なり、例えば自発脳磁場でいえば、gt一つだけの場合もあり得る。自発脳磁場推定部24aは、誘発脳磁場推定部24bが行う16式の演算と同様の演算により、10式のベクトルθが与えられたときに確率pxt|m′tを最大にするベクトルθを求める。具体的に自発脳磁場推定部24aは、KL情報量を最小とする赤池情報量基準を用いてベクトルθを決定する。このため自発脳磁場推定部24aは、上記確率pxt|m′tの対数logpxt|m′tを対数尤度とし、赤池情報量基準AICを、AIC=-2×pxt|m′t2×Nとして演算する。ここにNは自由に設定できるモデルパラメータの数を表し、ここではθの次元である、6となる。そして自発脳磁場推定部24aは、このAICを最小とするθを、既に述べたグリッドサーチや非線形最適化を組み合わせた方法で見出すこととしてもよい。この際の16式の演算においても、モンテカルロ法による演算を行うことができる。[自発脳磁場もベイズ推定に含める例]さらに、ここまでの説明では、誘発脳磁場推定部24bは、自発脳磁場推定部24aから入力される自発脳磁場に関するパラメータと7式とを用いて、自発脳磁場の時系列推定データg1,g2,...を得て、刺激時データm1,m2,...,mtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、m′′i=mi-νi-giを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、m′′1,m′′2,...,m′′tを得ていた。しかしながら本実施の形態はこれに限られない。すなわち誘発脳磁場推定部24bは、刺激時データm1,m2,...,mtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータを差引きして、m′′i=mi-νiを求め、ノイズがないと仮定される状態のデータ、m′′1,m′′2,...,m′′tを得ておき、次のように処理を遂行してもよい。この場合、誘発脳磁場推定部24bは、時刻tでの状態ベクトルxtに、時刻t以前の自発脳磁場を表すデータgt,gt1,...を含めて、xtstst1,gt,gt1,...Tとする。なお、肩のTは、転置を表す。ここで、自発脳磁場を時刻tよりどれだけ前のデータから含めるかは、gt1を予測するに十分な数とする。これは採用するモデルによっても異なり、例えばgt一つだけの場合もあり得る。誘発脳磁場推定部24bが行う16式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部24bは、既に述べたと同様にしてを求めることができる。[ノイズもベイズ推定に含める例]さらに、ここまでの説明では、誘発脳磁場推定部24bは、自発脳磁場の時系列推定データg1,g2,...を得て、刺激時データm1,m2,...,mtのそれぞれ対応する時刻のデータから、ノイズの時系列推定データの対応する時刻のデータ及び、自発脳磁場の時系列推定データの対応する時刻のデータを差引きして、m′′i=mi-νi-giを求め、ノイズ及び自発脳磁場がないと仮定される状態のデータ、m′′1,m′′2,...,m′′tを得ていた。しかしながら本実施の形態はこれに限られない。すなわち誘発脳磁場推定部24bは、刺激時データm1,m2,...,mtのそれぞれ対応する時刻のデータから、自発脳磁場の時系列推定データg1,g2,...の対応する時刻のデータを差引きして、m′′i=mi-giを求め、自発脳磁場がないと仮定される状態のデータ、m′′1,m′′2,...,m′′tを得ておき、次のように処理を遂行してもよい。この場合、誘発脳磁場推定部24bは、時刻tでの状態ベクトルxtに、時刻t以前の各ノイズを表すデータτt,qtqt1...qtT2et...etl1,nt...を含め、xtstst1,τt,qtqt1...qtT2et...etl1,ntTとする。なお、肩のTは、転置を表す。ここで、各ノイズのデータを、時刻tよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻t1における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なる。誘発脳磁場推定部24bが行う16式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部24bは、既に述べたと同様にしてを求めることができる。[ノイズも自発脳磁場もベイズ推定に含める例]さらに誘発脳磁場推定部24bは、刺激時データm1,m2,...,mtのそれぞれを、そのまま上記m′′1,m′′2,...,m′′tとして、次のように処理を遂行してもよい。すなわち誘発脳磁場推定部24bは、時刻tでの状態ベクトルxtに、時刻t以前の各ノイズを表すデータτt,qtqt1...qtT2et...etl1,nt...を含め、xtstst1,τt,qtqt1...qtT2et...etl1,nt,gt,gt1,...Tとする。なお、肩のTは、転置を表す。ここで、各ノイズのデータを、時刻tよりどれだけ前のデータから含めるかは、ノイズによって異なり、それぞれ時刻t1における対応するノイズを予測するに十分な数とする。これは採用するモデルによっても異なる。また誘発脳磁場推定部24bは、時刻t以前の自発脳磁場を表すデータgt,gt1,...も時刻tでの状態ベクトルxtに含める。ここで、自発脳磁場を時刻tよりどれだけ前のデータから含めるかは、gt1を予測するに十分な数とする。これは採用するモデルによっても異なり、例えばgt一つだけの場合もあり得る。誘発脳磁場推定部24bが行う16式の演算は、このように状態ベクトルの成分が増えても、既に説明した方法と同様に行うことができ、従って誘発脳磁場推定部24bは、既に述べたと同様にしてを求めることができる。[グリッドサーチにおけるARモデルの変数変換]さらにここまでの説明では定常ノイズ成分etの時間発展を表すARモデルにおいては、安定性の条件が見にくいので、次の方法により繰り返し変数変換を行って、1次のARモデルに変換してもよい。すなわち、L次のARモデルをL-1次のARモデルに変換する方法として、を用いる。ただし、i=1,2,...,l-1である。ここにaの肩にある添字のlは、l次の変数であることを意味する。この17式による変換を、1次となるまで繰り返して行うと、各変換の段階で、ajjが得られるここにj=1,2,...l。すると、安定なARモデルとなるための必要十分条件が、と得られる。そこで、グリッドサーチにおいては、この18式で表される範囲を定義域として、この定義域を均等に分割した格子点を用いればよい。[他のモデル]なお、一般に活動部位や刺激の種類に応じて、誘発脳磁場モデルの入力の確率分布は異なると考えられるうえ、線形モデルに限らず非線形なモデルが適当な場合が想定される。また本実施の形態は脳磁計のデータに限らず、一般的な測定データの処理に用いることができるものである。そこで利用者がそれぞれの仮説に従い、測定対象となるデータのモデルここでの例では誘発脳磁場のモデルの確率分布を作成してもよい。本実施の形態の処理は、種々の確率分布を仮定した場合でも実行可能なものである。また、自発脳磁場のモデルについても、ここでは8式によって表されるものとしたが、振幅の時間変化が大きい場合には、確率的入力項のパワーが時間的に変化するものと仮定し、当該パワーをパラメータを含む何らかの関数形で記述し、当該パラメータをモデルパラメータに含めて処理を行ってもよい。さらに、被験者の各自発脳磁場帯域のピーク周波数が典型的なものと大きく異なる場合は、中心周波数もモデルパラメータに含めてもよい。これらの場合も、上述の処理と同様にモデルパラメータの一つとして推定を行うことができる。[動作]本実施の形態は、以上の構成を有し、次のように動作する。すなわち、本実施の形態のデータ処理装置1は、ノイズのみのデータ、安静時データ、刺激時データを受け入れて記憶部12に保持しており、図3に例示するように、まずノイズのみのデータを用いて予め定めたノイズのモデルに含まれるノイズパラメータを推定するS1。またデータ処理装置1は、安静時データに含まれるノイズを、上記推定したノイズパラメータと予め定めたノイズのモデルとを用いて除去し、ノイズを除去した安静時データを参照して、予め定めた安静時データのモデルに含まれるパラメータを推定するS2。そしてデータ処理装置1は、刺激時データから、処理S1で推定したノイズパラメータと予め定めたノイズのモデルとを用いて除去し、刺激時データに含まれる自発脳磁場を、処理S2で推定したパラメータと、予め定めた安静時データのモデルとを用いて除去する。そしてデータ処理装置1は、こうしてノイズ及び自発脳磁場を除去したデータを参照して、予め定めた刺激時データのモデル誘発脳磁場のモデルに含まれるパラメータを推定するS3。データ処理装置1は、処理S3で推定したパラメータと誘発脳磁場のモデルとを用いて、誘発脳磁場のデータを推定して出力するS4。本発明の一実施例について説明する。ここでは被験者に、ボタンを押してもらうという刺激を与えたときの誘発脳磁場を得たときの例について述べる。図4は、LPFLowPassFilterを経た169回の測定データを加算平均した、従来の方法で得られた誘発脳磁場の波形複数のセンサの信号を重ね合わせて表している、以下も同様と、これらの波形から得られた等磁界線図である。図5は、ノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図である。この図5のようなデータを169回得て加算平均したものが図4に示した従来の結果である。図6は、この図5に示したデータ単一測定のデータに基づき、上記本実施の形態で説明した処理を経て本発明の実施例に係るデータ処理装置1が出力する誘発脳磁場の波形複数のセンサの信号を重ね合わせて表していると、これらの波形から得られた等磁界線図である。図4に示した従来の解析図と、図6に示した本実施例の解析図とを比較すると、略一致していることが理解される。図7は、自発脳磁場の強度が大きい場合のノイズを含めて実際に得られる時系列データの例と、当該時系列データを基に得られた等磁界線図である。この図7に例示したようなデータを複数回得て加算平均した、従来の方法によって推定された等磁界線図と、この図7に例示したような複数のデータのそれぞれから、本実施例のデータ処理装置1で推定して得た複数の等磁界線図図8に例示するとのコヒーレンスを演算して表したものを図9に示す。また、従来の方法で推定された等磁界線図を図10に示す。図9,図10から見られるように、図10で得られた等磁界線図のうち、信号の強度が大きい領域でのコヒーレンスは、図9に示されるように0.8を超えている。つまり、本実施例のデータ処理装置1での推定は、ロバストネスの面からも十分であることが理解される。"