JP2763021B2 - 電流ダイポール推定装置 - Google Patents
電流ダイポール推定装置Info
- Publication number
- JP2763021B2 JP2763021B2 JP7224973A JP22497395A JP2763021B2 JP 2763021 B2 JP2763021 B2 JP 2763021B2 JP 7224973 A JP7224973 A JP 7224973A JP 22497395 A JP22497395 A JP 22497395A JP 2763021 B2 JP2763021 B2 JP 2763021B2
- Authority
- JP
- Japan
- Prior art keywords
- dipoles
- current
- likelihood
- parameter
- distribution
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Description
【０００１】
【発明の属する技術分野】本発明は、生体の体表面頭皮
上で測定された電磁場分布に基づいて、生体内の電流ダ
イポールを推定する装置に関する。
上で測定された電磁場分布に基づいて、生体内の電流ダ
イポールを推定する装置に関する。
【０００２】
【従来の技術】生体内電流ダイポール推定は、従来、特
に人間の頭皮上で測定した電磁場分布から脳内の活動部
位を推定する手法として用いられており、この活動部位
を知ることによって、脳の高次機能や脳内疾患部位に関
する知見を得ることができる。
に人間の頭皮上で測定した電磁場分布から脳内の活動部
位を推定する手法として用いられており、この活動部位
を知ることによって、脳の高次機能や脳内疾患部位に関
する知見を得ることができる。
【０００３】被験者の頭皮上で測定した電磁場分布から
脳内電流ダイポールを推定するためには、電流ダイポー
ルによって頭部モデル表面に生ずる電磁場分布の理論値
と、測定した脳波または脳磁場データとの差異を最小に
するパラメータ（位置の３成分および方向の２成分、強
さ）を求める必要がある。これらのパラメータを推定す
る方法としては、従来、図７に示すような方法が知られ
ている。（例えば特開平７−１００１２１号公報）ここ
で説明する従来方法は、測定データを脳磁場に限定した
ものである。
脳内電流ダイポールを推定するためには、電流ダイポー
ルによって頭部モデル表面に生ずる電磁場分布の理論値
と、測定した脳波または脳磁場データとの差異を最小に
するパラメータ（位置の３成分および方向の２成分、強
さ）を求める必要がある。これらのパラメータを推定す
る方法としては、従来、図７に示すような方法が知られ
ている。（例えば特開平７−１００１２１号公報）ここ
で説明する従来方法は、測定データを脳磁場に限定した
ものである。
【０００４】ステップ１６では、電流ダイポール各パラ
メータの単位変化量ｄＰi （ｉ＝１，…，ｎ）の初期値
と各パラメータの収束条件εi を設定する。ステップ１
７では、電流ダイポールの数Ｎd と初期パラメータＰi0
を設定し、ステップ１８で評価関数の初期値Ｅ0 ＝Ｅ
（Ｐ10，…，Ｐn0）を計算する。ここで評価関数は、下
記の数１式で定義される。
メータの単位変化量ｄＰi （ｉ＝１，…，ｎ）の初期値
と各パラメータの収束条件εi を設定する。ステップ１
７では、電流ダイポールの数Ｎd と初期パラメータＰi0
を設定し、ステップ１８で評価関数の初期値Ｅ0 ＝Ｅ
（Ｐ10，…，Ｐn0）を計算する。ここで評価関数は、下
記の数１式で定義される。
【０００５】
【数１】 ただし、Ｂmkはｋ番目の測定点における磁場の測定値、
Ｂckはｋ番の測定点における磁場の計算値、Ｎm は測定
点の数である。
Ｂckはｋ番の測定点における磁場の計算値、Ｎm は測定
点の数である。
【０００６】ステップ１９では、変化させるパラメータ
Ｐi の指標ｉが集合Ｕの要素であるかどうかを調べる。
集合Ｕとは、電流ダイポールパラメータの方向と強さを
指定する成分（モーメント成分）のうち、値を連動して
変化させることが可能な指標の集合のことを言う。指標
ｉが集合Ｕの要素である場合はステップ２０〜２３に示
した処理を行い、指標ｉが集合Ｕの要素でない場合はス
テップ２４の処理を行う。これによって新たな評価関数
値Ｅijを得る。Ｅijはｉ番目のパラメータＰiを単位変
化量の整数倍だけ変化させて計算した評価関数値であ
り、下記の数２式で表わされる。
Ｐi の指標ｉが集合Ｕの要素であるかどうかを調べる。
集合Ｕとは、電流ダイポールパラメータの方向と強さを
指定する成分（モーメント成分）のうち、値を連動して
変化させることが可能な指標の集合のことを言う。指標
ｉが集合Ｕの要素である場合はステップ２０〜２３に示
した処理を行い、指標ｉが集合Ｕの要素でない場合はス
テップ２４の処理を行う。これによって新たな評価関数
値Ｅijを得る。Ｅijはｉ番目のパラメータＰiを単位変
化量の整数倍だけ変化させて計算した評価関数値であ
り、下記の数２式で表わされる。
【０００７】
【数２】 ステップ２５では初期値Ｅ0 と評価関数値Ｅijを比較
し、Ｅ0 ＜Ｅijならばステップ２９へ進み、ｄＰi をＮ
S で割ってこれを新たなｄＰi とし、ステップ３０へ進
む。Ｅ0 ＞Ｅijならばステップ２６へ進み、全てのパラ
メータ変化の中から評価関数を最小にするパラメータ変
化を選択する。例えば、ｘ成分をｙｄＰxだけ変化させ
た時に評価関数が最も小さくなる場合、Ｐx →Ｐx ＋ｙ
ｄＰx とする。さらにステップ２７においてステップ２
１と同様な方法でＥxy，Ｅrxy を計算し、ＥxyとＥrxy
を比較する。このときＥxy≠Ｅrxy ならばステップ１８
へ戻り、Ｅxy＝Ｅrxy ならばステップ２８でモーメント
qiを変更した後ステップ１８へ戻る。
し、Ｅ0 ＜Ｅijならばステップ２９へ進み、ｄＰi をＮ
S で割ってこれを新たなｄＰi とし、ステップ３０へ進
む。Ｅ0 ＞Ｅijならばステップ２６へ進み、全てのパラ
メータ変化の中から評価関数を最小にするパラメータ変
化を選択する。例えば、ｘ成分をｙｄＰxだけ変化させ
た時に評価関数が最も小さくなる場合、Ｐx →Ｐx ＋ｙ
ｄＰx とする。さらにステップ２７においてステップ２
１と同様な方法でＥxy，Ｅrxy を計算し、ＥxyとＥrxy
を比較する。このときＥxy≠Ｅrxy ならばステップ１８
へ戻り、Ｅxy＝Ｅrxy ならばステップ２８でモーメント
qiを変更した後ステップ１８へ戻る。
【０００８】ステップ３０では、ｄＰi が収束条件を満
たしているかを調べ、ｄＰi ＞εiならばステップ１８
へ戻り、ｄＰi ＜εi ならば処理を終了する。
たしているかを調べ、ｄＰi ＞εiならばステップ１８
へ戻り、ｄＰi ＜εi ならば処理を終了する。
【０００９】上に述べた従来手法を簡単にまとめると、
パラメータを単位変化量だけ変化させて残差二乗和の値
が減少したかどうかを調べ、なるべく残差二乗和の値を
小さくするようなパラメータを探索する方法であると言
える。
パラメータを単位変化量だけ変化させて残差二乗和の値
が減少したかどうかを調べ、なるべく残差二乗和の値を
小さくするようなパラメータを探索する方法であると言
える。
【００１０】
【発明が解決しようとする課題】しかしながら、上述の
方法には以下に述べるような問題点がある。従来手法で
は、推定したパラメータの評価関数として計算によって
求めた磁場と測定した磁場との残差二乗和を用いている
が、一般にパラメータの数が多くなるほど残差二乗和の
値は小さく出ることが知られている（例えば文献１、中
川ほか、「最小二乗法による実験データ解析」、ｐ．１
５０、東京大学出版会、１９８４）。パラメータの数は
電流ダイポールの数に比例するので、実際の電流ダイポ
ール数よりも多い数の電流ダイポールを用いて推定した
場合、正しい電流ダイポール数を用いて推定した場合よ
りも小さな残差二乗和が得られてしまう。
方法には以下に述べるような問題点がある。従来手法で
は、推定したパラメータの評価関数として計算によって
求めた磁場と測定した磁場との残差二乗和を用いている
が、一般にパラメータの数が多くなるほど残差二乗和の
値は小さく出ることが知られている（例えば文献１、中
川ほか、「最小二乗法による実験データ解析」、ｐ．１
５０、東京大学出版会、１９８４）。パラメータの数は
電流ダイポールの数に比例するので、実際の電流ダイポ
ール数よりも多い数の電流ダイポールを用いて推定した
場合、正しい電流ダイポール数を用いて推定した場合よ
りも小さな残差二乗和が得られてしまう。
【００１１】例えば実際の電流ダイポール数は２個であ
るのに、３個の電流ダイポールを仮定して推定を行って
も、２個と仮定した場合よりも小さな残差二乗和が得ら
れるので、ダイポールの数は３個であるという誤った判
断を下してしまう可能性がある。このような方法を、例
えば脳内の疾患部位推定に応用すると、正常な部位まで
をも異常部位と診断してしまう可能性があり、深刻な問
題を引き起こす。
るのに、３個の電流ダイポールを仮定して推定を行って
も、２個と仮定した場合よりも小さな残差二乗和が得ら
れるので、ダイポールの数は３個であるという誤った判
断を下してしまう可能性がある。このような方法を、例
えば脳内の疾患部位推定に応用すると、正常な部位まで
をも異常部位と診断してしまう可能性があり、深刻な問
題を引き起こす。
【００１２】本発明の課題は、上述の問題を克服するた
め、推定すべき電流ダイポールが未知の場合でも、電流
ダイポールの数とパラメータを正確に推定することが可
能な電流ダイポール推定装置を提供することにある。
め、推定すべき電流ダイポールが未知の場合でも、電流
ダイポールの数とパラメータを正確に推定することが可
能な電流ダイポール推定装置を提供することにある。
【００１３】
【課題を解決するための手段】本発明による電流ダイポ
ール推定装置は、生体体表面上の複数個所で測定した測
定電磁場分布を入力とし、前記測定電磁場分布の源泉と
して、１個または２個以上の電流ダイポールを生体内部
に仮定し、該電流ダイポールの個数・位置・向き・大き
さを指定するパラメータを推定する装置において、前記
電流ダイポールの数と前記電流ダイポールを表すパラメ
ータに比して十分多数の前記測定電磁場分布の測定点数
から前記電流ダイポール数の増加関数であるコンプレキ
シティを計算し、前記コンプレキシティを用いて前記電
流ダイポールの個数に関する事前分布を計算し、前記電
流ダイポールから計算される理論電磁場分布と、前記測
定電磁場分布から尤度を計算し、前記尤度を最大にする
ようなパラメータ変化を求めて元のパラメータに加える
演算を、前記パラメータ変化が所与の基準値よりも小さ
くなるまで繰返し、前記パラメータ変化が前記基準値よ
りも小さくなった時点で、評価関数を最大にする前記電
流ダイポールの個数とパラメータを出力し、前記評価関
数としては、前記事前分布と前記尤度を用いて計算し、
前記コンプレキシティの減少関数である事後分布を用い
ることを特徴とする。この装置は、第１の装置と呼ばれ
ても良い。
ール推定装置は、生体体表面上の複数個所で測定した測
定電磁場分布を入力とし、前記測定電磁場分布の源泉と
して、１個または２個以上の電流ダイポールを生体内部
に仮定し、該電流ダイポールの個数・位置・向き・大き
さを指定するパラメータを推定する装置において、前記
電流ダイポールの数と前記電流ダイポールを表すパラメ
ータに比して十分多数の前記測定電磁場分布の測定点数
から前記電流ダイポール数の増加関数であるコンプレキ
シティを計算し、前記コンプレキシティを用いて前記電
流ダイポールの個数に関する事前分布を計算し、前記電
流ダイポールから計算される理論電磁場分布と、前記測
定電磁場分布から尤度を計算し、前記尤度を最大にする
ようなパラメータ変化を求めて元のパラメータに加える
演算を、前記パラメータ変化が所与の基準値よりも小さ
くなるまで繰返し、前記パラメータ変化が前記基準値よ
りも小さくなった時点で、評価関数を最大にする前記電
流ダイポールの個数とパラメータを出力し、前記評価関
数としては、前記事前分布と前記尤度を用いて計算し、
前記コンプレキシティの減少関数である事後分布を用い
ることを特徴とする。この装置は、第１の装置と呼ばれ
ても良い。
【００１４】なお、前記評価関数として、前記尤度の対
数に前記コンプレキシティを掛けた関数を用いても良
い。この装置は、第２の装置と呼ばれても良い。
数に前記コンプレキシティを掛けた関数を用いても良
い。この装置は、第２の装置と呼ばれても良い。
【００１５】更に、前記尤度を最大化する演算におい
て、前記パラメータの初期値を複数個用意し、各初期値
に対して推定された複数のパラメータの組から推定パラ
メータの標準偏差または分散を計算し、前記尤度と前記
標準偏差または前記分散との比を取り、これを前記評価
関数として用いても良い。この装置は、第３の装置と呼
ばれても良い。
て、前記パラメータの初期値を複数個用意し、各初期値
に対して推定された複数のパラメータの組から推定パラ
メータの標準偏差または分散を計算し、前記尤度と前記
標準偏差または前記分散との比を取り、これを前記評価
関数として用いても良い。この装置は、第３の装置と呼
ばれても良い。
【００１６】
【発明の実施の形態】以下に、本発明で用いられる電流
ダイポール推定装置の原理について説明する。図２は電
流ダイポールモデルを説明するための図であり、被験者
の脳内活動部位１４とモデルの電流ダイポール１５とを
示している。生体体表面上のｎ点で測定された電磁場分
布をｙ＝（ｙ1 ，…，ｙn ）と書くと、この電磁場分布
ｙは、脳内神経細胞の活動（これを真のダイポールと呼
ぶ）によって頭皮上に生じる電磁場分布と、測定中ラン
ダムに混入する雑音から成る確率変数と見なせる。測定
された電磁場分布に対し、その源泉としてｓ個の電流ダ
イポールの存在を生体内に仮定する。これらの電流ダイ
ポールをモデル電流ダイポールと呼ぶが、以下では特に
必要でない限りモデル電流ダイポールのことを単にダイ
ポールと呼ぶことにする。
ダイポール推定装置の原理について説明する。図２は電
流ダイポールモデルを説明するための図であり、被験者
の脳内活動部位１４とモデルの電流ダイポール１５とを
示している。生体体表面上のｎ点で測定された電磁場分
布をｙ＝（ｙ1 ，…，ｙn ）と書くと、この電磁場分布
ｙは、脳内神経細胞の活動（これを真のダイポールと呼
ぶ）によって頭皮上に生じる電磁場分布と、測定中ラン
ダムに混入する雑音から成る確率変数と見なせる。測定
された電磁場分布に対し、その源泉としてｓ個の電流ダ
イポールの存在を生体内に仮定する。これらの電流ダイ
ポールをモデル電流ダイポールと呼ぶが、以下では特に
必要でない限りモデル電流ダイポールのことを単にダイ
ポールと呼ぶことにする。
【００１７】なお、ｓ個のダイポールを仮定する時、ｊ
番目のダイポールの位置ベクトルをｒj ＝（ｘj 1 ，ｘ
j 2 ，ｘj 3 ）と書き、方向・強さを表すモーメントベ
クトルをＭj ＝（Ｍj 1 ，Ｍj 2 ，Ｍj 3 ）と書くもの
とする。また、パラメータθ＝（θ1 ，…，θ6s）でｓ
個のダイポールの全パラメータをまとめて表すことにす
る。この時、ｊ番目のダイポールを指定する６個のパラ
メータは、θ6(j-1)+1，θ6(j-1)+2，…，θ6jと表され
る。
番目のダイポールの位置ベクトルをｒj ＝（ｘj 1 ，ｘ
j 2 ，ｘj 3 ）と書き、方向・強さを表すモーメントベ
クトルをＭj ＝（Ｍj 1 ，Ｍj 2 ，Ｍj 3 ）と書くもの
とする。また、パラメータθ＝（θ1 ，…，θ6s）でｓ
個のダイポールの全パラメータをまとめて表すことにす
る。この時、ｊ番目のダイポールを指定する６個のパラ
メータは、θ6(j-1)+1，θ6(j-1)+2，…，θ6jと表され
る。
【００１８】ダイポール推定方法として従来知られてい
る最尤法では、測定データｙに対して下記の数３式で表
わされる確率分布モデルを仮定し、
る最尤法では、測定データｙに対して下記の数３式で表
わされる確率分布モデルを仮定し、
【００１９】
【数３】 このＰ（ｙ｜θ）を最大にするパラメータθL （最尤推
定量と呼ばれる）を推定することによって行われる。た
だし、各データｙi は互いに独立に発生しているものと
する。特に、データが平均０，分散σ2 の正規分布に従
って発生している場合、ｙi の確率分布は下記の数４式
で与えられ、
定量と呼ばれる）を推定することによって行われる。た
だし、各データｙi は互いに独立に発生しているものと
する。特に、データが平均０，分散σ2 の正規分布に従
って発生している場合、ｙi の確率分布は下記の数４式
で与えられ、
【００２０】
【数４】 Ｐ（ｙ｜θ）を大きくするためには、次の数５式で表わ
される残差二乗和を小さくするようなパラメータを推定
すればよい。
される残差二乗和を小さくするようなパラメータを推定
すればよい。
【００２１】
【数５】 ここで、ｆi （θ）はｉ番目の測定点で計算される電磁
場分布の理論値である。
場分布の理論値である。
【００２２】しかしながら、最尤法では、モデルダイポ
ールの数が真の電流ダイポール数よりも多い場合、正し
い数のダイポールを仮定して推定した場合よりも尤度が
大きく出てしまうという問題点があった。この問題を克
服するため、本発明ではダイポールの個数ｓからコンプ
レキシティと呼ばれるモデルの複雑さを表す量Ｃ（ｓ）
を定義し、このＣ（ｓ）が大きくなるに従ってダイポー
ル数の事前分布π（ｓ）が小さくなるように設定する。
例えば、事前分布を次の数６式，数７式のように定義す
る。
ールの数が真の電流ダイポール数よりも多い場合、正し
い数のダイポールを仮定して推定した場合よりも尤度が
大きく出てしまうという問題点があった。この問題を克
服するため、本発明ではダイポールの個数ｓからコンプ
レキシティと呼ばれるモデルの複雑さを表す量Ｃ（ｓ）
を定義し、このＣ（ｓ）が大きくなるに従ってダイポー
ル数の事前分布π（ｓ）が小さくなるように設定する。
例えば、事前分布を次の数６式，数７式のように定義す
る。
【００２３】
【数６】
【００２４】
【数７】 ただし、Ｎはダイポール数の最大値である。Ｂａｙｅｓ
の定理を用いれば、ダイポール数の事後分布π（ｓ，θ
｜ｙ）が次の数８式で与えられる。
の定理を用いれば、ダイポール数の事後分布π（ｓ，θ
｜ｙ）が次の数８式で与えられる。
【００２５】
【数８】 ただし、ｐ（θ｜ｓ）はパラメータに対する事前分布で
ある。
ある。
【００２６】本発明の第１の装置では、この事後分布を
最大にするｓが最適なダイポール数であるとする。ダイ
ポールのコンプレキシティを考慮することによって、ダ
イポール数が過度になるのを防止し、真の電流ダイポー
ル数を推定することが可能となる。
最大にするｓが最適なダイポール数であるとする。ダイ
ポールのコンプレキシティを考慮することによって、ダ
イポール数が過度になるのを防止し、真の電流ダイポー
ル数を推定することが可能となる。
【００２７】本発明の第２の装置においては、推定され
たパラメータの良さの尺度として最大対数尤度にＣ
（ｓ）を掛けた次の数９式で表わされる量Ｌ（θL ，
ｓ）を損失関数とし、
たパラメータの良さの尺度として最大対数尤度にＣ
（ｓ）を掛けた次の数９式で表わされる量Ｌ（θL ，
ｓ）を損失関数とし、
【００２８】
【数９】 この損失関数を最小にするｓが最適なダイポール数であ
るとする。第１の装置と同様に、コンプレキシティの項
が大きすぎると損失関数の値が大きくなり、過剰なダイ
ポール数で推定してしまうことを防止している。更にこ
の場合、損失関数がコンプレキシティと評価関数の掛算
になっているため、第１の発明で述べた装置よりもダイ
ポール数の選択基準が鋭敏になるという効果がある。
るとする。第１の装置と同様に、コンプレキシティの項
が大きすぎると損失関数の値が大きくなり、過剰なダイ
ポール数で推定してしまうことを防止している。更にこ
の場合、損失関数がコンプレキシティと評価関数の掛算
になっているため、第１の発明で述べた装置よりもダイ
ポール数の選択基準が鋭敏になるという効果がある。
【００２９】第３の装置においては、仮定するダイポー
ル数ごとに複数の初期値を用意して最尤推定を行う。各
初期値に対して推定された複数のパラメータの間で標準
偏差または分散を計算し、尤度と標準偏差の比が最も大
きいダイポール数が最適なダイポール数であるとする。
モデルとして仮定したダイポールが真のダイポール数と
異る場合、ダイポールの推定位置が初期値に依存してば
らつくことが知られているため、この標準偏差の値を調
べることによりモデルダイポール数が適切であるかどう
かを判定できる。
ル数ごとに複数の初期値を用意して最尤推定を行う。各
初期値に対して推定された複数のパラメータの間で標準
偏差または分散を計算し、尤度と標準偏差の比が最も大
きいダイポール数が最適なダイポール数であるとする。
モデルとして仮定したダイポールが真のダイポール数と
異る場合、ダイポールの推定位置が初期値に依存してば
らつくことが知られているため、この標準偏差の値を調
べることによりモデルダイポール数が適切であるかどう
かを判定できる。
【００３０】
【実施例】本発明の実施例について図面を参照しながら
詳細に説明する。はじめに、本発明の第１の装置にかか
る第１の実施例について説明する。図１は第１の装置に
おける処理の流れを示したフローチャートである。ステ
ップ１では、ダイポールの個数ｓを１に初期化し、ダイ
ポールの最大数Ｎと収束条件εを設定したのちステップ
２に進む。
詳細に説明する。はじめに、本発明の第１の装置にかか
る第１の実施例について説明する。図１は第１の装置に
おける処理の流れを示したフローチャートである。ステ
ップ１では、ダイポールの個数ｓを１に初期化し、ダイ
ポールの最大数Ｎと収束条件εを設定したのちステップ
２に進む。
【００３１】ステップ２では、ダイポールの個数ｓおよ
び測定点の数ｎに基づいてモデルのコンプレキシティＣ
（ｓ，ｎ）を算出し、ステップ３に進む。コンプレキシ
ティの定義としては、例えば下記の数１０式、あるいは
数１１式などを用いるが、ダイポール数ｓが増えるにつ
れて増加するような関数であればよい。
び測定点の数ｎに基づいてモデルのコンプレキシティＣ
（ｓ，ｎ）を算出し、ステップ３に進む。コンプレキシ
ティの定義としては、例えば下記の数１０式、あるいは
数１１式などを用いるが、ダイポール数ｓが増えるにつ
れて増加するような関数であればよい。
【００３２】
【数１０】
【００３３】
【数１１】 ステップ３では、コンプレキシティＣを用いて数１２式
で表わされるダイポール数の事前分布を計算し、ステッ
プ４に進む。ただし、αは規格化定数である。
で表わされるダイポール数の事前分布を計算し、ステッ
プ４に進む。ただし、αは規格化定数である。
【００３４】
【数１２】 ステップ４では、繰り返しの回数ｔの値を０に初期化
し、パラメータの初期値θ0 ＝（θ0 1 ，…，θ0 m ）
を設定し、ステップ５に進む。ただし、ｍ＝６ｓはダイ
ポールパラメータの数である。パラメータの初期値は通
常、乱数を用いて定めるが、先見的知識を利用して特定
の値を初期値として用いることもできる。例えば、被験
者の頭皮上で測定された電磁場分布が視覚刺激によって
得られたものならば、脳の視覚野付近にダイポールパラ
メータの初期値を設定した方が良い。
し、パラメータの初期値θ0 ＝（θ0 1 ，…，θ0 m ）
を設定し、ステップ５に進む。ただし、ｍ＝６ｓはダイ
ポールパラメータの数である。パラメータの初期値は通
常、乱数を用いて定めるが、先見的知識を利用して特定
の値を初期値として用いることもできる。例えば、被験
者の頭皮上で測定された電磁場分布が視覚刺激によって
得られたものならば、脳の視覚野付近にダイポールパラ
メータの初期値を設定した方が良い。
【００３５】ステップ５では、尤度Ｐ（ｙ｜θt ，ｓ）
を計算し、ステップ６に進む。尤度として正規分布を仮
定すれば、以下の数１３式となる。
を計算し、ステップ６に進む。尤度として正規分布を仮
定すれば、以下の数１３式となる。
【００３６】
【数１３】 ただし、ｆi （θt ）はダイポールによって生じる電磁
場分布の理論値である。与えられたダイポールパラメー
タから磁場や電場を計算する方法は周知であり、例えば
磁場に関しては文献２（Ｓ．Ｊ．Ｗｉｌｌｉａｍｓｏｎ
ａｎｄ Ｌ．Ｋａｕｆｍａｎ，ＢＩＯＭＡＧＮＥＴＩ
ＳＭ，Ｊｏｕｒｎａｌ ｏｆ Ｍａｇｎｅｔｉｓｍ ａ
ｎｄ Ｍａｇｎｅｔｉｃ Ｍａｔｅｒｉａｌｓ ２２，
ｐｐ．１２９−２０１，１９８１）に示されており、電
場に関しては文献３（ＪａｍｅｓＰ．Ａｒｙ，Ｓｔａｎ
ｌｅｙ Ａ．Ｋｌｅｉｎ，Ｄｅｒｅｋ Ｈ．Ｆｅｎｄｅ
ｒ，Ｌｏｃａｔｉｏｎ ｏｆ Ｓｏｕｒｃｅｓ ｏｆ
Ｅｖｏｋｅｄ Ｓｃａｌｐ Ｐｏｔｅｎｔｉａｌｓ：Ｃ
ｏｒｒｅｃｔｉｏｎｓ ｆｏｒ Ｓｋｕｌｌ ａｎｄ
Ｓｃａｌｐ Ｔｈｉｃｋｎｅｓｓｅｓ，ＩＥＥＥ Ｂｉ
ｏｍｅｄｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ，Ｖｏｌ．Ｂ
ＭＥ−２８，Ｎｏ．６，Ｊｕｎｅ，ｐｐ．４４７−４５
２，１９８１）に詳しい。
場分布の理論値である。与えられたダイポールパラメー
タから磁場や電場を計算する方法は周知であり、例えば
磁場に関しては文献２（Ｓ．Ｊ．Ｗｉｌｌｉａｍｓｏｎ
ａｎｄ Ｌ．Ｋａｕｆｍａｎ，ＢＩＯＭＡＧＮＥＴＩ
ＳＭ，Ｊｏｕｒｎａｌ ｏｆ Ｍａｇｎｅｔｉｓｍ ａ
ｎｄ Ｍａｇｎｅｔｉｃ Ｍａｔｅｒｉａｌｓ ２２，
ｐｐ．１２９−２０１，１９８１）に示されており、電
場に関しては文献３（ＪａｍｅｓＰ．Ａｒｙ，Ｓｔａｎ
ｌｅｙ Ａ．Ｋｌｅｉｎ，Ｄｅｒｅｋ Ｈ．Ｆｅｎｄｅ
ｒ，Ｌｏｃａｔｉｏｎ ｏｆ Ｓｏｕｒｃｅｓ ｏｆ
Ｅｖｏｋｅｄ Ｓｃａｌｐ Ｐｏｔｅｎｔｉａｌｓ：Ｃ
ｏｒｒｅｃｔｉｏｎｓ ｆｏｒ Ｓｋｕｌｌ ａｎｄ
Ｓｃａｌｐ Ｔｈｉｃｋｎｅｓｓｅｓ，ＩＥＥＥ Ｂｉ
ｏｍｅｄｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ，Ｖｏｌ．Ｂ
ＭＥ−２８，Ｎｏ．６，Ｊｕｎｅ，ｐｐ．４４７−４５
２，１９８１）に詳しい。
【００３７】以降は尤度が正規分布であるとして説明を
行うが、それ以外の確率分布の場合でも同様な方法で電
流ダイポールを推定することが可能である。分散σ2 が
未知の場合、尤度を最大にするσ2 は次の数１４式で与
えられる。
行うが、それ以外の確率分布の場合でも同様な方法で電
流ダイポールを推定することが可能である。分散σ2 が
未知の場合、尤度を最大にするσ2 は次の数１４式で与
えられる。
【００３８】
【数１４】 ステップ６では尤度Ｐ（ｙ｜θt ，ｓ）が増加するよう
なパラメータ変化ｄt を求め、ステップ７へ進む。パラ
メータ変化ｄt の算出には、例えば最急降下法やガウス
−ニュートン法、ハイブリッド法などを用いることがで
きる。これらの方法の詳細に関しては、例えば文献１に
詳しい。
なパラメータ変化ｄt を求め、ステップ７へ進む。パラ
メータ変化ｄt の算出には、例えば最急降下法やガウス
−ニュートン法、ハイブリッド法などを用いることがで
きる。これらの方法の詳細に関しては、例えば文献１に
詳しい。
【００３９】ステップ７では、下記の数１５式にもとづ
いてパラメータ変化のノルムを計算し、
いてパラメータ変化のノルムを計算し、
【００４０】
【数１５】 ノルムがεより小さければステップ１０に進み、そうで
なければステップ８に進む。
なければステップ８に進む。
【００４１】ステップ８では、ステップ７で得られたパ
ラメータ変化ｄt を元のパラメータθt に加え、ステッ
プ９に進む。
ラメータ変化ｄt を元のパラメータθt に加え、ステッ
プ９に進む。
【００４２】ステップ９では、ｔ←ｔ＋１としてステッ
プ４に戻り、同様の手続きを繰り返す。
プ４に戻り、同様の手続きを繰り返す。
【００４３】ステップ１０では、ステップ３において求
めた事前分布とステップ５で求めた尤度から、下記の数
１６式による事後分布π（ｓ｜ｙ）
めた事前分布とステップ５で求めた尤度から、下記の数
１６式による事後分布π（ｓ｜ｙ）
【００４４】
【数１６】 を計算し、ステップ１１へ進む。ここでｐ（θ｜ｓ）は
パラメータの事前分布、Ｎはダイポール数の最大値であ
る。この事後分布の計算において、尤度はステップ４以
降の演算で求めた最大尤度を用い、パラメータの事前分
布に関しては、特に理由の無い限りｐ（θ｜ｓ）＝１と
してよい。
パラメータの事前分布、Ｎはダイポール数の最大値であ
る。この事後分布の計算において、尤度はステップ４以
降の演算で求めた最大尤度を用い、パラメータの事前分
布に関しては、特に理由の無い限りｐ（θ｜ｓ）＝１と
してよい。
【００４５】ステップ１１では、現在のダイポール数ｓ
が既定の最大値Ｎに等しいかどうかを調べ、ｓがＮに等
しくなければステップ１２へ進む。ステップ１２では、
尤度を最大にするパラメータθと事後分布π（ｓ｜ｙ）
を記憶し、ステップ１３へ進む。ステップ１３ではｓ←
ｓ＋１としてステップ２へ戻る。他方、ｓがＮに等しけ
れば処理を終了する。
が既定の最大値Ｎに等しいかどうかを調べ、ｓがＮに等
しくなければステップ１２へ進む。ステップ１２では、
尤度を最大にするパラメータθと事後分布π（ｓ｜ｙ）
を記憶し、ステップ１３へ進む。ステップ１３ではｓ←
ｓ＋１としてステップ２へ戻る。他方、ｓがＮに等しけ
れば処理を終了する。
【００４６】本発明の第１の装置の効果を示すため、シ
ミュレーションデータを用いて実際にダイポールを推定
した結果を示す。本実験に使用したデータは以下のよう
に生成した。半径１００［ｍｍ］の球内の位置（４０，
４０，３０），（−４０，−４０，３０）に真のダイポ
ールを２個配置し、これらのダイポールによって球面上
に生じる磁場を球面上２５点で測定する。モーメント・
ベクトルの値はそれぞれ（−１０，０，０），（１０，
０，０）とする（位置の単位はｍｍ，モーメントの単位
はｎＡｍである。また、座標系の原点は球の中心とす
る。）。更に、このデータにＳＮ比１０．０のランダム
ノイズを付加する。このように生成したデータに対し、
１個から４個のモデルダイポールを用いて推定を行っ
た。
ミュレーションデータを用いて実際にダイポールを推定
した結果を示す。本実験に使用したデータは以下のよう
に生成した。半径１００［ｍｍ］の球内の位置（４０，
４０，３０），（−４０，−４０，３０）に真のダイポ
ールを２個配置し、これらのダイポールによって球面上
に生じる磁場を球面上２５点で測定する。モーメント・
ベクトルの値はそれぞれ（−１０，０，０），（１０，
０，０）とする（位置の単位はｍｍ，モーメントの単位
はｎＡｍである。また、座標系の原点は球の中心とす
る。）。更に、このデータにＳＮ比１０．０のランダム
ノイズを付加する。このように生成したデータに対し、
１個から４個のモデルダイポールを用いて推定を行っ
た。
【００４７】図３は残差二乗和を評価基準とした従来方
法による結果である。最も小さな残差二乗和を与えてい
るのは４個のモデルダイポールを用いて推定したときで
あるが、真のダイポールの数は２個であるから、この結
果は誤りである。実際、表１に示すように、推定された
ダイポールの位置の誤差は２個のときに最小となってい
る。
法による結果である。最も小さな残差二乗和を与えてい
るのは４個のモデルダイポールを用いて推定したときで
あるが、真のダイポールの数は２個であるから、この結
果は誤りである。実際、表１に示すように、推定された
ダイポールの位置の誤差は２個のときに最小となってい
る。
【００４８】
【表１】 ただし、モデルダイポールが１個の場合、２個ある真の
ダイポールのうち近い方との距離を、２個の場合は、そ
れぞれ近い方のダイポールとの距離を足して２で割った
ものを位置の平均推定誤差とする。モデルダイポールが
３個以上の場合には、真のダイポールに近い２個のダイ
ポールと真のダイポールとの距離を足して２で割ったも
のを位置の平均推定誤差としている。
ダイポールのうち近い方との距離を、２個の場合は、そ
れぞれ近い方のダイポールとの距離を足して２で割った
ものを位置の平均推定誤差とする。モデルダイポールが
３個以上の場合には、真のダイポールに近い２個のダイ
ポールと真のダイポールとの距離を足して２で割ったも
のを位置の平均推定誤差としている。
【００４９】これに対し、本発明の第１の装置を用いて
推定した結果を図４に示す。この図４から分かるよう
に、コンプレキシティとしてＣ1 ，Ｃ2 いずれを用いた
場合にもモデルダイポール数が２個の場合に事後分布が
最大となっており、正しい結果が得られている。
推定した結果を図４に示す。この図４から分かるよう
に、コンプレキシティとしてＣ1 ，Ｃ2 いずれを用いた
場合にもモデルダイポール数が２個の場合に事後分布が
最大となっており、正しい結果が得られている。
【００５０】本発明の第２の装置では、評価関数として
以下の数１７式による関数を用いる。
以下の数１７式による関数を用いる。
【００５１】
【数１７】 ここで、Ｃ3 はコンプレキシティ、Ｐ（ｙ｜θ，ｓ）は
尤度、βは任意の正数である。以下はβ＝１として説明
を行うが、これ以外の正数を用いてもよい。コンプレキ
シティＣ3 としては、下記の数１８式を用いる。
尤度、βは任意の正数である。以下はβ＝１として説明
を行うが、これ以外の正数を用いてもよい。コンプレキ
シティＣ3 としては、下記の数１８式を用いる。
【００５２】
【数１８】 ここで、ηは誤推定が生じる確率で、ここでは０．００
１とした。コンプレキシティＣ3 は、ダイポール数とと
もに増加するものであれば、これ以外の量を用いてもよ
い。
１とした。コンプレキシティＣ3 は、ダイポール数とと
もに増加するものであれば、これ以外の量を用いてもよ
い。
【００５３】この評価関数を用いてダイポール推定を行
った結果を図５に示す。使用したシミュレーションデー
タは第１の実施例で用いたものと同じである。図５よ
り、ダイボール数が２個のときに評価関数の値Ｅは最大
となり、真のダイポール数が推定できていることが分か
る。
った結果を図５に示す。使用したシミュレーションデー
タは第１の実施例で用いたものと同じである。図５よ
り、ダイボール数が２個のときに評価関数の値Ｅは最大
となり、真のダイポール数が推定できていることが分か
る。
【００５４】本発明の第３の装置では、複数の初期値を
用いて推定を行い、各々の初期値に対して得られた尤度
と、推定パラメータの標準偏差との比を評価関数とし、
これを最大にするダイポール数が最適であるとする。な
お、推定パラメータの標準偏差の代わりに推定パラメー
タの分散を用いても良い。
用いて推定を行い、各々の初期値に対して得られた尤度
と、推定パラメータの標準偏差との比を評価関数とし、
これを最大にするダイポール数が最適であるとする。な
お、推定パラメータの標準偏差の代わりに推定パラメー
タの分散を用いても良い。
【００５５】モデルダイポール数をｓ個と仮定して推定
をする時、Ｋ個の初期パラメータθ(S,1) ，…，θ
(s,K) を用意し、各々の初期値に対し、尤度Ｐ(S,1) ，
Ｐ(s,K)と推定パラメータθL (S,1) ，…，θL (s,K)
が得られたとする。これらの推定結果から、推定パラメ
ータの平均値ベクトルμが以下の数１９式により計算で
きる。
をする時、Ｋ個の初期パラメータθ(S,1) ，…，θ
(s,K) を用意し、各々の初期値に対し、尤度Ｐ(S,1) ，
Ｐ(s,K)と推定パラメータθL (S,1) ，…，θL (s,K)
が得られたとする。これらの推定結果から、推定パラメ
ータの平均値ベクトルμが以下の数１９式により計算で
きる。
【００５６】
【数１９】 この平均値ベクトルを用いて推定パラメータのｉ成分の
標準偏差ベクトルσsのｉ成分σs i を以下の数２０式
のように計算する。
標準偏差ベクトルσsのｉ成分σs i を以下の数２０式
のように計算する。
【００５７】
【数２０】 ここで、θi (s,K) は、ｋ番目の初期パラメータから推
定を開始して得られたパラメータのｉ成分、μs i は平
均値ベクトルのｉ成分である。本発明の第３の装置で
は、こうして得られた標準偏差ベクトルのノルムと尤度
を用いて、評価関数Ｅ（ｓ）を下記の数２１式で定義
し、この評価関数を最大にするｓが最適なモデルダイポ
ール数であるとする。
定を開始して得られたパラメータのｉ成分、μs i は平
均値ベクトルのｉ成分である。本発明の第３の装置で
は、こうして得られた標準偏差ベクトルのノルムと尤度
を用いて、評価関数Ｅ（ｓ）を下記の数２１式で定義
し、この評価関数を最大にするｓが最適なモデルダイポ
ール数であるとする。
【００５８】
【数２１】 この評価関数を用いてダイポール推定を行った結果を図
６に示す。使用したシミュレーションデータは第１の装
置で用いたものと同じである。図６より、ダイポール数
が２個のときにＥの値は最大となり、真のダイポール数
が推定できていることが分かる。
６に示す。使用したシミュレーションデータは第１の装
置で用いたものと同じである。図６より、ダイポール数
が２個のときにＥの値は最大となり、真のダイポール数
が推定できていることが分かる。
【００５９】
【発明の効果】以上に述べたように、本発明装置によれ
ば、モデルダイポールの数が未知の場合でも、適切なモ
デルダイポール数の評価基準が得られるため、過剰な数
のダイポールを用いて誤った推定を行うことがない。し
たがって、本発明の装置をてんかんや脳硬塞などの医療
診断へ応用した場合でも、疾患部位を誤って推定するこ
とが無く、安全で正確な診断が可能となる。
ば、モデルダイポールの数が未知の場合でも、適切なモ
デルダイポール数の評価基準が得られるため、過剰な数
のダイポールを用いて誤った推定を行うことがない。し
たがって、本発明の装置をてんかんや脳硬塞などの医療
診断へ応用した場合でも、疾患部位を誤って推定するこ
とが無く、安全で正確な診断が可能となる。
【図面の簡単な説明】
【図１】本発明の第１の装置の処理の流れを示すフロー
チャートである。
チャートである。
【図２】電流ダイポールモデルを説明するための図であ
る。
る。
【図３】従来手法による推定結果である。
【図４】本発明の第１の装置による推定結果である。
【図５】本発明の第２の装置による推定結果である。
【図６】本発明の第３の装置による推定結果である。
【図７】従来手法の処理の流れを示すフローチャートで
ある。
ある。
１４ 脳内活動部位 １５ 電流ダイポール
Claims (3)
- 【請求項１】 生体体表面上の複数個所で測定した測定
電磁場分布を入力とし、前記測定電磁場分布の源泉とし
て、１個または２個以上の電流ダイポールを生体内部に
仮定し、該電流ダイポールの個数・位置・向き・大きさ
を指定するパラメータを推定する装置において、前記電
流ダイポールの数と前記電流ダイポールを表すパラメー
タに比して十分多数の前記測定電磁場分布の測定点数か
ら前記電流ダイポール数の増加関数であるコンプレキシ
ティを計算し、前記コンプレキシティを用いて前記電流
ダイポールの個数に関する事前分布を計算し、前記電流
ダイポールから計算される理論電磁場分布と、前記測定
電磁場分布から尤度を計算し、前記尤度を最大にするよ
うなパラメータ変化を求めて元のパラメータに加える演
算を、前記パラメータ変化が所与の基準値よりも小さく
なるまで繰返し、前記パラメータ変化が前記基準値より
も小さくなった時点で、評価関数を最大にする前記電流
ダイポールの個数とパラメータを出力し、前記評価関数
としては、前記事前分布と前記尤度を用いて計算し、前
記コンプレキシティの減少関数である事後分布を用いる
ことを特徴とする電流ダイポール推定装置。 - 【請求項２】 前記評価関数として、前記尤度の対数に
前記コンプレキシティを掛けた関数を用いることを特徴
とする請求項１記載の電流ダイポール推定装置。 - 【請求項３】 前記尤度を最大化する演算において、前
記パラメータの初期値を複数個用意し、各初期値に対し
て推定された複数のパラメータの組から推定パラメータ
の標準偏差または分散を計算し、前記尤度と前記標準偏
差または前記分散との比を取り、これを前記評価関数と
して用いることを特徴とする請求項１記載の電流ダイポ
ール推定装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP7224973A JP2763021B2 (ja) | 1995-09-01 | 1995-09-01 | 電流ダイポール推定装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP7224973A JP2763021B2 (ja) | 1995-09-01 | 1995-09-01 | 電流ダイポール推定装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH0966038A JPH0966038A (ja) | 1997-03-11 |
| JP2763021B2 true JP2763021B2 (ja) | 1998-06-11 |
Family
ID=16822114
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP7224973A Expired - Fee Related JP2763021B2 (ja) | 1995-09-01 | 1995-09-01 | 電流ダイポール推定装置 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP2763021B2 (ja) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3033508B2 (ja) * | 1997-01-20 | 2000-04-17 | 日本電気株式会社 | 生体内活動部位推定方法 |
| EP3575991A4 (en) * | 2017-01-24 | 2020-02-19 | Nec Corporation | INFORMATION PROCESSING DEVICE, INFORMATION PROCESSING METHOD, AND RECORDING MEDIUM ON WHICH AN INFORMATION PROCESSING PROGRAM IS RECORDED |
-
1995
- 1995-09-01 JP JP7224973A patent/JP2763021B2/ja not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| JPH0966038A (ja) | 1997-03-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Chung et al. | Weighted Fourier series representation and its application to quantifying the amount of gray matter | |
| US11062450B2 (en) | Systems and methods for modeling neural architecture | |
| Cheung et al. | Estimation of cortical connectivity from EEG using state-space models | |
| Li et al. | Estimating the number of independent components for functional magnetic resonance imaging data | |
| Li et al. | Test–retest stability analysis of resting brain activity revealed by blood oxygen level‐dependent functional MRI | |
| Chauvin et al. | Accurate statistical tests for smooth classification images | |
| Song et al. | A modified probabilistic neural network for partial volume segmentation in brain MR image | |
| Gorodnitsky et al. | Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm | |
| US20130034277A1 (en) | Systems and methods for modeling and processing functional magnetic resonance image data using full-brain vector auto-regressive model | |
| Baid et al. | Comparative study of k-means, gaussian mixture model, fuzzy c-means algorithms for brain tumor segmentation | |
| Park et al. | Bayesian active learning of neural firing rate maps with transformed gaussian process priors | |
| Young et al. | Interactive scaling with individual subjects | |
| Sui et al. | A method for accurate group difference detection by constraining the mixing coefficients in an ICA framework | |
| Zhang et al. | Network response regression for modeling population of networks with covariates | |
| Radich et al. | EEG dipole localization bounds and MAP algorithms for head models with parameter uncertainties | |
| JP2763021B2 (ja) | 電流ダイポール推定装置 | |
| Huizenga et al. | Estimating and testing the sources of evoked potentials in the brain | |
| Cheung et al. | Cross validation for selection of cortical interaction models from scalp EEG or MEG | |
| Rosipal et al. | Kernel PCA feature extraction of event-related potentials for human signal detection performance | |
| Ollier et al. | Penalized Poisson model for network meta‐analysis of individual patient time‐to‐event data | |
| Veeramalla et al. | Neural source localization using particle filter with optimal proportional set resampling | |
| DE102008038451A1 (de) | Verfahren zur rechnergestützten Lokalisation eines mobilen Objekts mit Hilfe eines merkmalbasierten Ortungsverfahrens | |
| Bortel et al. | Regularization techniques in realistic Laplacian computation | |
| Tang et al. | Least squares regression methods for clustered ROC data with discrete covariates | |
| CN109949307B (zh) | 一种基于主成分分析的图像分割的方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 19980225 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080327 Year of fee payment: 10 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090327 Year of fee payment: 11 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090327 Year of fee payment: 11 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100327 Year of fee payment: 12 |
|
| LAPS | Cancellation because of no payment of annual fees |