"局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報を格納している大脳皮質活動モデル情報格納部と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報を格納している頭皮血流変化モデル情報格納部と、送光プローブと受光プローブの組であるプローブセットを1組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納している安静時光強度情報格納部と、前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、前記安静時光強度情報格納部の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部とを具備する脳活動情報出力装置。前記推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、前記演算式に、前記感度情報と前記光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備する請求項1記載の脳活動情報出力装置。前記推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、前記演算式に、前記感度情報と前記光強度変化情報を代入し、当該演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備する請求項1記載の脳活動情報出力装置。前記プローブの配置密度が1.5cm以下である請求項1から請求項3いずれか記載の脳活動情報出力装置。前記プローブの配置密度が2.0cm以下である請求項1から請求項3いずれか記載の脳活動情報出力装置。記憶媒体に、局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、送光プローブと受光プローブの組であるプローブセットを1組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、脳活動時光強度情報取得部、光強度変化情報取得部、推定部、および出力部により実現される脳活動情報出力方法であって、前記脳活動時光強度情報取得部が、前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得ステップと、前記光強度変化情報取得部が、前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得ステップで取得された脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得ステップと、前記推定部が、前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定ステップと、前記出力部が、前記推定ステップで取得された大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力ステップとを具備する脳活動情報出力方法。記憶媒体に、局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、送光プローブと受光プローブの組であるプローブセットを1組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、コンピュータを、前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部として機能させるためのプログラム。しかしながら、従来の技術においては、大脳皮質における血行動態反応は局所的であり、頭皮における血行動態反応は広域的であるという知識をうまく活用することができず、頭皮血流に起因するアーチファクトを適切に除去できないために、精度高く脳活動を推定することができなかった。また、従来の技術においては、波長毎の観測値からそれぞれ別々に吸収係数の変化を推定し、それを変換することで、ヘモグロビンの濃度を推定していたため、精度高く脳活動を推定することができなかった。なお、ここでの脳活動の推定とは、神経活動に伴う脳血行動態反応を推定することである。本第一の発明の脳活動情報出力装置は、局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報を格納している大脳皮質活動モデル情報格納部と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報を格納している頭皮血流変化モデル情報格納部と、送光プローブと受光プローブの組であるプローブセットを1組以上有する近赤外分光計測装置を用いて、各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納している安静時光強度情報格納部と、近赤外分光計測装置を用いて、各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、安静時光強度情報格納部の安静時光強度情報と、脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、複数の波長ごと、およびプローブセットごとの光強度変化情報を、大脳皮質活動モデル情報および頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部とを具備する脳活動情報出力装置である。かかる構成により、精度高く脳活動を推定できる。また、本第二の発明の脳活動情報出力装置は、第一の発明に対して、推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、演算式に、感度情報と光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備する脳活動情報出力装置である。かかる構成により、極めて精度高く脳活動を推定できる。また、本第三の発明の脳活動情報出力装置は、第一の発明に対して、推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の情報の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、演算式に、感度情報と光強度変化情報を代入し、演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備する脳活動情報出力装置である。かかる構成により、精度高く脳活動を推定できる。本発明による脳活動情報出力装置によれば、精度高く脳活動を推定できる。以下、脳活動情報出力装置等の実施形態について図面を参照して説明する。なお、実施の形態において同じ符号を付した構成要素は同様の動作を行うので、再度の説明を省略する場合がある。実施の形態1本実施の形態において、脳活動を推定する脳活動情報出力システムについて説明する。また、脳活動情報出力システムは、近赤外分光計測装置NIRSを具備する。なお、脳活動の推定とは、神経活動に伴う脳血行動態反応を推定することである。また、本実施の形態における脳活動情報出力システムにおいて、拡散光トモグラフィDOTを用いる。DOTとは、頭部内での光子移動を拡散過程としてモデル化することにより、脳内活動と観測光強度変化の関係を物理的に正しく求め、それを利用して脳活動の推定の精度を高めようというものである。また、本実施の形態における脳活動情報出力システムにおいて、特に、局所的な活動である大脳皮質の活動、および広域的な活動である頭皮の血流変化をモデル化しており、かかるモデルを用いて、脳活動を推定する方法を採用する。また、特に、本実施の形態において、頭皮血流をDOTのモデルに組み込んで、大脳皮質における神経活動と頭皮血流に起因するアーチファクトとを同時推定することにより、頭皮血流に起因するアーチファクトを除去し、脳活動を推定する階層変分ベイズ法について説明する。なお、後述する数式に出現する上付きのcは大脳皮質sは頭皮を表す。さらに、具体的には、本実施の形態における脳活動情報出力システムにおいて、近赤外分光計測装置における観測結果から、神経活動に伴う大脳皮質内のオキシヘモグロビンOxyHb、デオキシヘモグロビンDeoxyHbの濃度変化の場所と大きさを推定する。図1は、本実施の形態における脳活動情報出力システムの概念図である。脳活動情報出力システムは、近赤外分光計測装置1、および脳活動情報出力装置2を具備する。また、図2は、本実施の形態における脳活動情報出力システムのブロック図である。近赤外分光計測装置1は、送光プローブ11、受光プローブ12、受付部13、および制御部14を具備する。近赤外分光計測装置1は、送光プローブ11と受光プローブ12の組であるプローブセットチャネルとも言うを1組以上有する。脳活動情報出力装置2は、大脳皮質活動モデル情報格納部21、頭皮血流変化モデル情報格納部22、安静時光強度情報格納部23、安静時光強度取得部24、安静時光強度蓄積部25、脳活動時光強度情報取得部26、光強度変化情報取得部27、推定部28、および出力部29を具備する。また、推定部28は、感度情報格納手段281、演算式格納手段282、および推定手段283を具備する。近赤外分光計測装置1を構成する送光プローブ11は、近赤外光を生体ここでは、特に頭部に照射する。受光プローブ12は、生体ここでは、特に頭部内部で吸収あるいは散乱され、生体内から出てくる光を受光する。なお、プローブの配置密度が3cm未満であることは好適である。また、プローブの配置密度が2cm以下であることはさらに好適である。また、プローブの配置密度が1.7cm以下であることは好適である。さらに、プローブの配置密度が1.5cm以下であることは極めて好適である。プローブの配置密度とは、送光プローブ11と受光プローブ12間の間隔である。なお、プローブの配置密度は密であればある程良いのは、頭皮血流をとらえるための短い距離のチャネルがたくさんできること、および脳活動をとらえるための3~4cm程度の中距離のチャネルがたくさんできること等からである。受付部13は、外部からの指示を受け付ける。かかる指示は、近赤外光を生体に照射する指示である照射指示などである。制御部14は、受付部13が照射指示を受け付けた場合に、送光プローブ11に近赤外光の照射を指示する。なお、送光プローブ11は、通常、制御部14から照射指示を受け付けた場合に、近赤外光を生体に照射する。脳活動情報出力装置2は、近赤外分光計測装置1が計測した光強度に関する情報を用いて、脳活動に関する情報である脳活動情報を取得し、出力する装置である。ここで、出力とは、ディスプレイへの表示、プロジェクターを用いた投影、プリンタへの印字、音出力、外部の装置への送信、記録媒体への蓄積、他の処理装置や他のプログラムなどへの処理結果の引渡しなどを含む概念である。脳活動情報出力装置2を構成する大脳皮質活動モデル情報格納部21は、大脳皮質活動モデル情報を格納している。大脳皮質活動モデル情報は、局所的な活動である大脳皮質の活動のモデルを示す情報である。大脳皮質活動モデル情報は、例えば、後述する局所的フィルタWC、後述する正則化の処理を示すプログラムやフラグ、大脳皮質の各ボクセルにおける正則化大きな値にペナルティに関するガンマ分布のシェイプパラメータ後述するγα0i、大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値後述するα0i-等である。なお、大脳皮質活動モデル情報に対応する正則化は、ボクセル毎に異なった正則化である。また、ここでいう正則化とは、大きな正則化対象の値についてペナルティをかけることである。頭皮血流変化モデル情報格納部22は、頭皮血流変化モデル情報を格納している。頭皮血流変化モデル情報は、広域的な活動である頭皮の血流変化のモデルを示す情報である。頭皮血流変化モデル情報は、例えば、後述する広域的フィルタWS、後述する正則化の処理を示すプログラムやフラグ、頭皮の正則化に関するガンマ分布のシェイプパラメータ後述するγβ0、頭皮の正則化に関するガンマ分布の平均値後述するβ0-等である。なお、頭皮血流変化モデル情報に対応する正則化は、全部のボクセルで共通な正則化である。安静時光強度情報格納部23は、複数の各波長に対応付けて安静時光強度情報を格納している。安静時光強度情報は、送光プローブ11と受光プローブ12の組であるプローブセットを1組以上有する近赤外分光計測装置1を用いて、各プローブセットを構成する各送光プローブ11から、一の光強度の光を発した場合に、安静状態において測定された各受光プローブ12の光強度に関する情報である安静時光強度取得部24は、安静状態の被験者が装着している近赤外分光計測装置1の受光プローブ12から取得した光強度に関する情報である安静時光強度情報を、近赤外分光計測装置1から取得する。安静時光強度取得部24は、複数の各波長ごとに、安静時光強度情報を取得する。安静時光強度蓄積部25は、安静時光強度取得部24が取得した安静時光強度情報を、複数の各波長に対応付けて安静時光強度情報格納部23に蓄積する。脳活動時光強度情報取得部26は、複数の各波長ごとに脳活動時光強度情報を取得する。脳活動時光強度情報は、脳活動が行われている状態において測定された各受光プローブ12の光強度に関する情報である。脳活動時光強度情報は、被験者の脳の活動が行われている状態において、近赤外分光計測装置1を用いて、各プローブセットを構成する各送光プローブ11から、一の光強度の光を発した場合に、各受光プローブ12が受光した光強度に関する情報である。光強度変化情報取得部27は、安静時光強度情報格納部23の安静時光強度情報と、脳活動時光強度情報取得部26が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度変化情報を取得する。光強度変化情報は、光強度の変化に関する情報である。光強度変化情報は、例えば、光強度の変化率、光強度の変化率の対数をとった値、光強度の変化量等である。また、光強度変化情報取得部27は、光強度のゆらぎ自発脳活動やセンサーの機器ノイズにより光強度に加わるノイズの情報Σyも取得しても良い。推定部28は、複数の波長ごと、およびプローブセットごとの光強度変化情報を、大脳皮質活動モデル情報および頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。なお、推定部28は、推定する値後述するXの数より少ない観測値後述するYを用いて、値後述するXを推定し、取得する。なお、かかる推定部28が解く問題を逆問題という。推定部28は、後述する変分ベイズ法を用いて逆問題を解くことは好適であるが、逆問題の解法は問わない。また、推定部28は、頭部の領域ごとに、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。頭部の領域ごととは、例えば、頭部のボクセルごと、または頭部を三角錐に区切った領域ごと等である。感度情報格納手段281は、感度情報を格納し得る。感度情報とは、大脳皮質または頭皮の活動の変化例えば、後述するXが対応と検出された光強度の変化の情報例えば、後述するYが対応の関係を示す情報である。感度情報は、例えば、後述する感度行列である。感度情報は、推定部28が算出し、感度情報格納手段281に、少なくとも一時的に格納しても良い。演算式格納手段282は、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式例えば、演算式1であり、ベイズ推定の演算式を格納し得る。また、演算式格納手段282は、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式例えば、演算式2であり、コスト関数である演算式を格納し得る。なお、演算式1は、例えば、後述する数式28である。また、演算式2は、例えば、後述する数式31である。推定手段283は、演算式例えば、演算式1に、感度情報と光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する。また、推定手段283は、演算式例えば、演算式2に、感度情報と光強度変化情報を代入し、演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得しても良い。なお、コストが最小であるとは、コストに-1をかけた値が最大であることも含む。つまり、コストが最小であるとは、実質的に考慮する。出力部29は、推定部28が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。なお、出力部29は、頭皮の血管のオキシヘモグロビンの濃度変化、および頭皮の血管のデオキシヘモグロビンの濃度変化も出力しても良い。また、出力部29は、大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を、数値を用いて出力しても良いし、図的に出力しても良い。その出力態様は問わない。大脳皮質活動モデル情報格納部21、頭皮血流変化モデル情報格納部22、安静時光強度情報格納部23、感度情報格納手段281、および演算式格納手段282は、不揮発性の記録媒体が好適であるが、揮発性の記録媒体でも実現可能である。大脳皮質活動モデル情報格納部21等に大脳皮質活動モデル情報等が記憶される過程は問わない。例えば、記録媒体を介して大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部21等で記憶されるようになってもよく、通信回線等を介して送信された大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部21等で記憶されるようになってもよく、あるいは、入力デバイスを介して入力された大脳皮質活動モデル情報等が大脳皮質活動モデル情報格納部21等で記憶されるようになってもよい。安静時光強度取得部24は、通常、受信手段等から実現され得る。安静時光強度蓄積部25、脳活動時光強度情報取得部26、光強度変化情報取得部27、推定部28、および推定手段283は、通常、MPUやメモリ等から実現され得る。安静時光強度蓄積部25等の処理手順は、通常、ソフトウェアで実現され、当該ソフトウェアはROM等の記録媒体に記録されている。但し、ハードウェア専用回路で実現しても良い。出力部29は、ディスプレイやスピーカー等の出力デバイスを含むと考えても含まないと考えても良い。出力部29は、出力デバイスのドライバーソフトまたは、出力デバイスのドライバーソフトと出力デバイス等で実現され得る。次に、脳活動情報出力システムを構成する脳活動情報出力装置2の動作について、図3のフローチャートを用いて説明する。なお、図3のフローチャートの開始前において、安静時光強度取得部24は、複数の波長ごとに、安静時光強度情報を取得し、安静時光強度蓄積部25は、複数の波長ごとに、当該安静時光強度情報を安静時光強度情報格納部23に蓄積した、とする。そして、安静時光強度情報格納部23には、複数の各波長ごとに、安静時光強度情報が格納されている、とする。ステップS301脳活動時光強度情報取得部26は、複数の各波長ごとに脳活動時光強度情報を取得したか否かを判断する。複数の各波長の脳活動時光強度情報を取得すればステップS302に行き、取得しなければステップS301に戻る。ステップS302光強度変化情報取得部27は、安静時光強度情報格納部23に格納されている安静時光強度情報を、複数の各波長に対応付けて読み出す。ステップS303光強度変化情報取得部27は、ステップS301で取得された各波長ごとの複数の脳活動時光強度情報と、ステップS302で読み出された波長ごとの複数の安静時光強度情報とから、波長ごとに、光強度変化情報を算出する。光強度変化情報取得部27は、例えば、波長ごとに、安静時光強度情報から脳活動時光強度情報を減算することにより、光強度変化情報を取得する。ステップS304推定手段283は、演算式格納手段282から演算式を読み出す。ステップS305推定手段283は、感度情報格納手段281から感度情報を読み出す。なお、例えば、推定手段283は予め感度情報を算出し、感度情報格納手段281に格納している、とする。ステップS306推定手段283は、大脳皮質活動モデル情報格納部21から大脳皮質活動モデル情報を読み出す。ステップS307推定手段283は、頭皮血流変化モデル情報格納部22から頭皮血流変化モデル情報を読み出す。ステップS308推定手段283は、演算式に、光強度変化情報、感度情報を代入し、かつ、大脳皮質活動モデル情報と頭皮血流変化モデル情報とを用いて、変分ベイズ推定法により、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する。かかる変分ベイズ推定法について、図4のフローチャートを用いて説明する。ステップS309出力部29は、推定部28が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。ステップS301に戻る。なお、図3のフローチャートにおいて、電源オフや処理終了の割り込みにより処理は終了する。次に、ステップS308の変分ベイズ推定法について、図4のフローチャートを用いて説明する。ステップS401推定手段283は、予め格納されている初期値を読み出す。読み出す初期値は、は、γα0i、α0i--はαの真上に配置、γβ0、β0--はβの真上に配置である。また、初期値は、ユーザが与えても良いし、他の装置から受信等しても良い。ここで、γα0iは大脳皮質の各ボクセルにおける正則化大きな値にペナルティに関するガンマ分布のシェイプパラメータの初期値、α0i-は大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値の初期値、γβ0は頭皮の正則化に関するガンマ分布のシェイプパラメータの初期値、β0-は頭皮の正則化に関するガンマ分布の平均値の初期値である。ステップS402推定手段283は、後述する第一ステップの演算を行う。第一ステップの演算は、変分ベイズ法における、X,σ最大化ステップである。ステップS403推定手段283は、後述する第二ステップの演算を行う。第二ステップの演算は、変分ベイズ法における、A,β最大化ステップである。ステップS404推定手段283は、大脳皮質の各ボクセルにおける正則化に関するガンマ分布のシェイプパラメータγαi、大脳皮質の各ボクセルにおける正則化に関するガンマ分布の平均値αi-、頭皮の正則化に関するガンマ分布のシェイプパラメータγβ、および頭皮の正則化に関するガンマ分布の平均値β-の値が収束したか否かを判断する。収束すればステップS405に行き、収束しなければステップS402に戻る。なお、収束したか否かは、予め決められた条件を満たすほど、上記の4つの変数の各値の変化が小さくなったか否かにより判断される。ステップS405推定手段283は、第一ステップの演算を行う。ステップS406推定手段283は、推定結果を算出する演算を行う。処理を終了する。以下、本実施の形態における脳活動情報出力システムの具体的な動作について説明する。脳活動情報出力システムの概念図は、図1である。まず、本脳活動情報出力システムにおいて、頭部が半無限一層モデルであると仮定し、計測された光強度から脳活動を推定する場合について説明する。また、プローブは、図5に示すように、15mm間隔で頭部モデルの表面に配置されている、とする。図5において、黒丸が送光プローブ11、白丸が受光プローブ12を表す。また、図5において、直方体で囲ってある部分はこれから推定値を求めていく領域頭部で、一辺4mmからなる立方体のボクセルに分けて考える。図5において、頭部は、縦11ボクセル、横11ボクセル、高さ5ボクセルに分割されている。ここで、深さ1を頭皮、深さ23を頭蓋骨脳脊髄液、深さ4以降を大脳皮質とする。また、今回の頭部モデルの光学特性は、ヒト灰白質と同じであり、換算散乱係数μs=1.11mm-1、吸収係数μa=0.0186mm-1とする。さらに、表面境界における屈折率はn=1.4とする。上記の条件で、被験者は安静な状態を保っており、かかる状態において、安静時光強度取得部24は、複数の波長ごとに、安静時光強度情報を取得し、安静時光強度蓄積部25は、複数の波長ごとに、当該安静時光強度情報を安静時光強度情報格納部23に蓄積した、とする。そして、例えば、安静時光強度情報格納部23には、波長=780nm、830nm以下、適宜、単位を省略する。の2つの各波長ごとに、安静時光強度情報が格納された、とする。なお、安静時とは、脳活動や頭皮血流変化があまりない安静状態をいう。そして、一般に、光強度は機器ノイズや自発脳活動によってゆらぐので、各波長ごと、各チャネルごとに多数計測し、安静時光強度取得部24は、その平均値と分散共分散行列を算出することは好適である。ここで、光強度の平均値は、数式1のように示される。また、分散共分散行列はΣyである。なお、数式1において、Φ0は安静時の光強度であり、rs→rの真上に→が存在は送光プローブ11の位置、rd→rの真上に→が存在は受光プローブ12の位置を示す。数式1のバーは、平均値を意味する。次に、脳活動時光強度情報取得部26は、各チャネルごと、波長=780、830の2つの各波長ごとに、脳活動時光強度情報を取得した、とする。具体的には、脳活動時光強度情報取得部26は、数式2に示す脳活動時光強度情報を、チャネルごと、波長ごとに取得する。また、脳活動時光強度は複数回計測して取得してもよい。その計測の番号をサンプル番号とする。次に、光強度変化情報取得部27は、安静時光強度情報格納部23に格納されている2つの安静時光強度情報を、各波長波長=780、830に対応付けて読み出す。次に、光強度変化情報取得部27は、取得された波長ごとの2つの脳活動時光強度情報と、読み出された波長ごとの2つの安静時光強度情報とから、波長ごとに、光強度変化情報を、以下のように算出する。つまり、光強度変化情報取得部27は、ここでは、数式3の演算式を格納しており、それぞれのサンプル番号ごとに、かかる数式3に、安静時光強度情報logΦ0と脳活動時光強度情報logΦを代入し、光強度変化情報Yを取得する。ここでtはサンプル番号を表す。なお、YはYtを各列に持つ行列であり、その次元はチャネル数×波長数×サンプル数である。今回は、サンプル数は1とする。従って、ここでは、例えば、Yは40次元である、とする。また、ここで、光強度変化情報取得部27は、数式4を格納しており、数式4の演算式に、安静時光強度取得部24が取得した安静時光強度の平均値と、脳活動時光強度情報数式2を代入し、光強度変化情報Yを算出することは好適である。そして、頭部モデルと、図5に示すようにプローブ配置送光プローブ11と受光プローブ12の配置が決まっているので、推定手段283は、光拡散過程を解くことにより、以下の数式5を用いて、各波長ごとの吸収係数変化に関する感度行列Gを計算して、取得する。ここで、感度行列Gは、感度情報の一例である。なお、数式5において、定義の仕方によっては定数倍がかかる場合がある。なお、数式5において、jはチャネル送光プローブ11と受光プローブ12のペアのインデックス、iは頭部ボクセルのインデックスである。また、数式5において、Φ0は安静時の光強度であり、Φ0第一変数,第二変数は、第一変数の位置から出発し、第二変数の位置に到達する光の強度を示す。なお、Φ0は頭部モデル内での光子の拡散過程を解くことで得られる。また、rs,j→rの真上に→が存在はチャネルjの送光プローブ11の位置を、rd,j→rの真上に→が存在はチャネルjの受光プローブ12の位置を示す。また、ri→rの真上に→が存在は頭部ボクセルiの位置を示す。また、ここでは、推定手段283は、感度行列を、頭皮図5の深さ1の領域と大脳皮質図5の深さ4または深さ5の領域のみで求める。図5において、頭蓋骨脳脊髄液深さ2または深さ3の領域にはヘモグロビンとオキシヘモグロビンが存在せず、常にX=0だから、感度行列G=0である。また、一般的に、頭部モデルの光学特性は波長ごとに違うので、推定手段283は、感度行列を波長毎に算出する。なお、Xは、ここでは、オキシヘモグロビンの濃度変化およびデオキシヘモグロビンの濃度変化である。つまり、推定手段283は、Gc780、Gs780、Gc830、Gs830を算出し、感度情報格納手段281に蓄積しておく。ここで、上付きのcは大脳皮質深さ4または深さ5の領域、sは頭皮深さ1の領域を意味する。下付きの780は波長780nm、830は波長830nmを意味する。つまり、Gc780は、波長780nmにおける大脳皮質の感度行列である。次に、吸収係数変化に関する感度行列から、ヘモグロビン濃度変化に関する感度行列を算出する。その変換式は、取り扱う波長が波長=780、830の2つの場合、以下の数式6に示すように詳細に定義される。数式6において、δμaは脳活動時の吸収係数と安静時の吸収係数の差である。また、δμaは波長毎に値が異なる。つまり、δμa,780は、波長780nmにおける、脳活動時の吸収係数と安静時の吸収係数の差である。また、数式6において、εは、モル吸光係数であり、波長や物質オキシヘモグロビン、デオキシヘモグロビンにより異なる。なお、数式8は、波長ごと、物質ごとのモル吸光係数の値を示す。また、それぞれの波長における吸収係数の変化は、オキシヘモグロビンの濃度変化Xoxyと、デオキシヘモグロビンの濃度変化Xdeoxyとを用いて、数式7のように表される。そして、数式6の1は、数式7を用いて、数式6の2に変形される。そして、数式6の2は、計算により、数式6の3に変形される。そして、ヘモグロビン濃度変化に関する感度行列GはGcoxy,Gcdeoxy,Gsoxy,Gsdeoxyである。次に、推定手段283は、感度情報格納手段281から感度行列GGcoxy,Gcdeoxy,Gsoxy,Gsdeoxyを読み出す。なお、Gc780などは吸収係数変化に関する感度行列、Gcoxyなどは、ヘモグロビン濃度変化に関する感度行列である。また、吸収係数変化に関する感度行列とヘモグロビン濃度変化に関する感度行列とは、上述したように、数式6の定義式により簡単に変換できる。なお、数式5を用いて、吸収係数変化に関する感度行列が算出され、次に、当該感度行列にモル吸光係数などをかけて、ヘモグロビン濃度変化に関する感度行列が算出される。そして、ヘモグロビン濃度変化に関する感度行列が上記の推定に用いられる。次に、推定手段283は、大脳皮質活動モデル情報格納部21から大脳皮質活動モデル情報例えば、局所的フィルタWC等を読み出す。また、推定手段283は、頭皮血流変化モデル情報格納部22から頭皮血流変化モデル情報例えば、広域的フィルタWS等を読み出す。次に、推定手段283は、以下のように2つのステップ第一ステップ、第二ステップにより、変分ベイズ推定を用いて、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を推定する。なお、推定手段283は、後述する演算式等を、通常、予め格納している、とする。まず、第一ステップとは、以下の数式9から数式15を用いて、Σ、σ―-はσの上に配置、γσ、Δ、<ΔX2>を算出する処理である。ここで、Σは、数式9により定義される。また、Σyは安静時の観測値Yの分散共分散行列である。また、σは脳活動時の分散共分散行列の安静時に比べての割合の逆数である。つまり、脳活動時の観測値Yの分散共分散行列はσ-1Σyであるとする。なお、初期値として、P0σ∝1/σのような分布に従うと仮定するが、アルゴリズム途中ではガンマ分布にしたがう。その平均値がσ―である。また、シェイプパラメータをγσとする。なお、数式10のγσは、自由エネルギーが最大となるような値である。なお、行列の上付けのTは、行列の転置をあらわす。また、数式9において、Wはスムージングフィルタであり、Wcは局所的フィルタ、Wsは広域的フィルタを示す。具体的には、例えば、Wcは半値幅8mmのガウシアンフィルタを用い、Wsはラプラス行列の逆行列を用いる。なお、数式10において、Nはチャネル数、Tはサンプル数である。また、Lは用いる波長の数である。また、Trは行列のトレースをあらわす。なお、数式12において、Iは単位行列を示す。また、数式11、12において、Aはαiを対角要素に持つ行列αiは大脳皮質の各ボクセル毎に在るである。つまり、AおよびA--はAの真上に存在するは、数式13のようになる。数式13において、αiはそれぞれのボクセルiにおける、正則化大きな値にペナルティに関するパラメータである。そして、αiはガンマ分布に従うとする。そして、αi--はαの真上に存在するはガンマ分布の平均値である。また、γαiはガンマ分布のシェイプパラメータである。また、大脳皮質ボクセル数がMc、頭皮ボクセル数がMsである。なお、ガンマ分布は、数式14である。なお、数式14においてΓはガンマ関数をあらわす。また、数式11、12において、βは頭皮における正則化に関するパラメータである。βはガンマ分布に従うとする。そしてβ--はβの真上に存在するはガンマ分布の平均値、γβはガンマ分布のシェイプパラメータである。また、ガンマ分布は、数式15である。次に、推定手段283は、第二ステップを実行する。第二ステップとは、γαi、αi--はαの真上に配置、γβ、β--はβの真上に配置を、以下の数式16から数式19を用いて算出する処理である。なお、γαiは大脳皮質のボクセルにおけるガンマ分布のシェイプパラメータ、αi-は大脳皮質のボクセルにおけるガンマ分布の平均値、γβは頭皮のボクセルにおけるガンマ分布のシェイプパラメータ、β-は頭皮のボクセルにおけるガンマ分布の平均値である。また、上記の初期値は、γα0i、α0i--はαの真上に配置、γβ0、β0--はαの真上に配置であるとし、各初期値は予め与えられている、とする。数式16において、Tは、サンプル数である。そして、推定手段283は、第一ステップと第二ステップとを繰り返し、第二ステップにおいて、αi--はαの真上に配置、γαi、β--はβの真上に配置、γβの値の変化が所定の条件を満たすほど小さい場合n回目の値とn-1回目の値との差が閾値以内など、収束したとして、処理を終了する。なお、最後は、第一ステップを行って終了する。そして、以下の数式20におけるXが取得される。上記において、各行列は以下のような次元数を持つ。YはNL×T次元行列、GcはNL×Mc次元行列、GsはNL×Ms次元行列、XcはMc×T次元行列、XsはMs×T次元行列である。実験1次に、本脳活動情報出力システムの第一の実験結果について説明する。本実験において、脳活動情報出力システムを構成する近赤外分光計測装置1は、20のチャネルを有する。チャネルとは、送光プローブ11と受光プローブ12のすべての組み合わせである。つまり、本実験において、図6に示すように、近赤外分光計測装置1は5つの送光プローブ11と4つの受光プローブ12を有する。また、図7に示すように、各送光プローブ11の識別子と各受光プローブ12の識別子との組み合わせが20セット存在する。なお、図6において、黒丸が送光プローブ11、白丸が受光プローブ12である。また、本実験において、最近傍のプローブの間隔はここでは15mmとした。また、ボクセルの大きさは4mm、縦11ボクセル、横11ボクセル、深さ5ボクセルとした。さらに、深さ1は頭皮、深さ23は頭蓋骨脳脊髄液、深さ45は大脳皮質であるとした。そして、以下のようなシミュレーション実験を行った。つまり、シミュレーション実験において、被験者の脳活動があまり表れない安静状態で、近赤外分光計測装置1の送光プローブ11は、波長780nmと830nmで発光した、と仮定した。そして、受光プローブ12が受光する光量の対数値を、各波長において、チャネルごとに、20回分、取得した。そして、脳活動情報出力装置2の安静時光強度取得部24は、受光プローブ12が受光する光量の対数値の平均と分散共分散行列を算出した。なお、図8は、波長780nmにおける前記光量の対数値の平均値である。なお、図示しないが、波長830nmにおける前記光量の対数値の平均も、同様に、取得した。図9は2波長のデータをまとめた、チャネル数×波長数=40次元を、複数回今回は20回計測したデータから計算した分散共分散行列の一部波長780nmに関する部分である。そして、安静時光強度蓄積部25は、安静時光強度情報格納部23に、波長780nmと波長830nmにおける光量の対数値の平均および分散共分散行列を蓄積した。次に、脳活動情報出力装置2の脳活動時光強度情報取得部26は、被験者の脳活動がある状態で、波長780nmと波長830nmにおける光量の対数値を、チャネルごとに取得した。取得したサンプル数は、今回は1回である。この光量の対数値は、図10のtrial1の属性値である。なお、被験者の脳活動がある状態では、頭皮の血流変化もある。次に、光強度変化情報取得部27は、波長780nmと波長830nmの各波長において、チャネルごと光強度変化情報数式4のYを算出した。つまり、光強度変化情報取得部27は、数式4に従って、波長ごと、チャネルごと、およびサンプルごとに、光量の対数値の平均値から、脳活動時光強度情報取得部26が取得した光量の対数値を減算し、光強度変化情報Yを算出した。その光強度変化情報Yが、図10の属性値Yである。なお、波長780nmにおける観測値Yに関して、Ymaxを図示したのが、図11aである。また、図11abにおいて、線の明るさは、その線の端にある送光プローブと受光プローブのセットからなるチャネルにおける光強度変化情報を表している。図11aにおいて、広がった頭皮血流アーチファクトによって、全体的に値が出ていることがわかる。しかし、図11aを注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。また、波長830nmにおける観測値Yに関して、Ymaxを図示したのが、図11bである。図11bにおいて、広がった頭皮血流アーチファクトによって、全体的に値が出ている。しかし、図11bを注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。なお、780nmよりも強く値が出ているが、これはシミュレーション設定によるものである。次に、推定手段283は、感度情報格納手段281からヘモグロビン濃度変化に関する感度行列上述したGcoxy,Gcdeoxy,Gsoxy,Gsdeoxyを読み出す。次に、推定手段283は、大脳皮質活動モデル情報格納部21から大脳皮質活動モデル情報例えば、局所的フィルタWcを読み出す。また、推定手段283は、頭皮血流変化モデル情報格納部22から頭皮血流変化モデル情報例えば、広域的フィルタWsを読み出す。次に、推定手段283は、以下のように変分ベイズアルゴリズムを用いて、Xを推定する。つまり、推定手段283は、予め格納している初期値を読み出す。初期値は、ここでは、すべてのボクセルに対してγα0i=0、かつ、すべてのボクセルに対してα0i-=10000、かつ、γβ0=0、かつ、β0-=10000である。そして、推定手段283は、上述したステップ1とステップ2の繰り返し処理を1000回実行し、Xを算出した。推定手段283が算出したXの要素であり、頭皮深さ1のオキシヘモグロビンの濃度変化Xoxyの推定値の集合を図12に示す。図12において、各セルの値は、各ボクセルに対応する。したがって、図12において、11×11の値が存在する。また、数値は行下方向が頭部モデルの横方向、列右方向が頭部モデルの縦方向と対応する。かかることは、図13、図14等と同様である。また、推定手段283が算出したXの要素であり、大脳皮質深さ4のオキシヘモグロビンの濃度変化Xoxyの推定値の集合を図13に示す。なお、図示しないが、大脳皮質深さ5のオキシヘモグロビンの濃度変化Xoxyについても、推定手段283は、同様に算出した。また、推定手段283が算出したXの要素であり、頭皮深さ1のデオキシヘモグロビンの濃度変化Xdeoxyの推定値の集合を図14に示す。さらに、推定手段283が算出したXの要素であり、大脳皮質深さ4のデオキシヘモグロビンの濃度変化Xdeoxyの推定値の集合を図15に示す。なお、図示しないが、大脳皮質深さ5のデオキシヘモグロビンの濃度変化Xoxyについても、推定手段283は、同様に算出した。次に、出力部29は、推定部28が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。出力部29の出力例を図16、図17に示す。図16aは、頭皮深さ1のオキシヘモグロビンの濃度変化の出力例である。図16a等において、縦軸は頭部モデルの縦方向、横軸は頭部モデルの横方向を示し、その数字の単位はmmである。また、白黒の濃淡で表されるXの値の単位はmMである。図16bは、大脳皮質深さ4のオキシヘモグロビンの濃度変化の出力例である。また、図17aは、頭皮深さ1のデオキシヘモグロビンの濃度変化の出力例である。図17bは、大脳皮質深さ4のデオキシヘモグロビンの濃度変化の出力例である。さらに、シミュレーションで仮定した真のXの値が図18から図21である。図18は、シミュレーションで仮定した真の頭皮深さ1のオキシヘモグロビンの濃度変化の値の集合である。また、図19は、シミュレーションで仮定した真の大脳皮質深さ4のオキシヘモグロビンの濃度変化の値の集合である。また、図20は、シミュレーションで仮定した真の頭皮深さ1のデオキシヘモグロビンの濃度変化の値の集合である。図21は、シミュレーションで仮定した真の大脳皮質深さ4のデオキシヘモグロビンの濃度変化の値の集合である。シミュレーションで仮定した真の値と、推定部28が取得した推定値とを比較すると、推定部28は、脳活動をほぼ正確に推定できていることが分かる。次に、上記で説明した脳活動情報出力装置2の推定方法の概念について、図22を用いて説明する。脳活動情報出力装置2における推定は、推定する値Xの数より少ない観測値Yを用いて、値Xを推定しなければならない。これは逆問題の中でも不良設定問題と呼ばれ、正則化と呼ばれるような、なんらかの制限を入れなければ解くことはできない。今回用いる推定法ではベイズ推定を用いる。ベイズ推定ではその正則化と呼ばれるような制限を事前分布において確率的に組み込む。今回の推定法では、その事前分布において生理学的知見を反映させたような、大脳皮質と頭皮において異なる正則化とモデル化を行う。また、今回、ベイズ推定を現実的に行う方法の一例として、変分ベイズ法を用いる。また、値Xを推定するために、脳活動情報出力装置2は、大脳皮質活動モデル情報として、スパース化を含む正則化のプログラムおよび局所的フィルタWcを保持しており、かつ、頭皮血流変化モデル情報として、スパース化を含まない正則化のプログラムおよび広域的フィルタWsを保持している。図22において、Y1からYLが観測値である。図22の4における観測値YとXの関係は、光拡散過程においてRytov近似により吸収係数の差δμaが十分小さければ数式21であらわされる。数式21を、数式3、数式5、数式6を用いて書き換えると数式22となる。しかし、数式22は観測ノイズを考慮していない関係式である。推定部28は、図22の4において、観測値YはGXにノイズσ-1Σyがのって出てくると仮定する。この仮定は、脳活動時のノイズが安静時のノイズの定数倍でうまく近似できるという考えに基づく。なお、数23において、Nは正規分布をあらわす。例えば、Nはxが平均μ、分散共分散行列Σxの正規分布に従うことをあらわす。そして、推定部28は、図22の3において、生理学的知見により、大脳皮質には局所的フィルタ、頭皮は広域的フィルタをかける。また、フィルタをかける前のXの値をZとおく。かかる処理の具体的な演算を数式24に示す。Zの値に制限をかけないと、逆問題は解くことはできない。脳活動情報出力装置2は、以下のように、確率的に条件をつける。推定部28は、図22の2において、頭皮においては、頭皮血流変化モデル情報に従って、全部のボクセルで共通な正則化を行う。ここでいう正則化とは、大きなZの値つまり大きな濃度変化についてペナルティをかけることである。これは、図22の3の広域的フィルタWsと合わせて、空間的に広がった活動を意味する。かかる処理の具体的な演算を数式25に示す。数式25において、σは脳活動時の分散共分散行列の安静時に比べての割合の逆数である。また、Tはサンプル数である。βは頭皮の正則化に関するパラメータである。Iは単位行列である。また、推定部28は、大脳皮質においては、各ボクセル毎に異なった正則化を行う。ここでは、少数のボクセルのZの値が大きくなるように選ばれるようになっているこれをスパース化という。。これは、図22の3の局所的フィルタWcと合わせて、局在した活動を意味する。かかる処理の具体的な演算を数式26に示す。数式26において、Aはαiを対角要素に持つ行列であり、パラメータαiは大脳皮質の各ボクセル毎に存在する正則化に関するパラメータである。また、パラメータαi、β、σの分布は、数式27のように仮定する、とする。なお、数式27において、Gammaはガンマ分布を示す。以上のように、推定部28は、確率モデル化した。以下、ベイズ推定により、観測値Yから、最適なσ、A、βの値と、Xの推定値を求める。数式28はベイズの定理である。例えば、左辺の事後分布を最大とするようなX、σ、A、β、またはX、σ、A、βの期待値を推定値とすればよい。しかしながら、直接計算することは難しいため、今回の推定法は変分ベイズ法を用いる。変分ベイズ法では、事後分布を直接計算する代わりに、試験事後分布Qを用意し、真の事後分布に近似させる。この近似は数式29で定義される自由エネルギーFの最大化によって行われる。今回は試験分布におけるX、σとA、βの間に数30で示すように独立分解仮定をおく。これにより、計算を簡単に行うことが可能となる。この仮定により自由エネルギーの最大化は、Q1に関してFを最大化するX、σ最大化ステップと、Q2に関してFを最大化するA、β最大化ステップを交互に繰り返すことにより実現できる。X、σ最大化ステップに対応するアルゴリズムが、前述した第一ステップである。また、A、β最大化ステップに対応するアルゴリズムが、前述した第二ステップである。値が収束するまで繰り返した後、今回の推定においては、試験事後分布QにおけるX、σ、A、βの期待値を推定値とした。実験2次に、本脳活動情報出力システムの第二の実験結果について説明する。本実験において、最近傍のプローブの間隔は20mmとした。実験2における他の条件は、実験1と同じである。そして、以下のようなシミュレーション実験を行った。つまり、シミュレーション実験において、被験者の脳活動があまり表れない安静状態で、近赤外分光計測装置1の送光プローブ11は、波長780nmと830nmで発光した、と仮定した。そして、受光プローブ12が受光する光量の対数値を、各波長において、チャネルごとに、20回分、取得した。そして、脳活動情報出力装置2の安静時光強度取得部24は、受光プローブ12が受光する光量の対数値の平均と分散共分散行列を算出した。なお、図25は、波長780nmにおける前記光量の対数値の平均値である。なお、図示しないが、波長830nmにおける前記光量の対数値の平均も、同様に、取得した。図26は2波長のデータをまとめた、チャネル数×波長数=40次元を、複数回今回は20回計測したデータから計算した分散共分散行列の一部波長780nmに関する部分である。そして、安静時光強度蓄積部25は、安静時光強度情報格納部23に、波長780nmと波長830nmにおける光量の対数値の平均および分散共分散行列を蓄積した。次に、脳活動情報出力装置2の脳活動時光強度情報取得部26は、被験者の脳活動がある状態で、波長780nmと波長830nmにおける光量の対数値を、チャネルごとに取得した。取得したサンプル数は、今回は1回である。この光量の対数値は、図27のtrial1の属性値である。なお、被験者の脳活動がある状態では、頭皮の血流変化もある。次に、光強度変化情報取得部27は、波長780nmと波長830nmの各波長において、チャネルごと光強度変化情報数式4のYを算出した。つまり、光強度変化情報取得部27は、数式4に従って、波長ごと、チャネルごと、およびサンプルごとに、光量の対数値の平均値から、脳活動時光強度情報取得部26が取得した光量の対数値を減算し、光強度変化情報Yを算出した。その光強度変化情報Yが、図27の属性値Yである。なお、波長780nmにおける観測値Yに関して、Ymaxを図示したのが、図28aである。また、図28abにおいて、線の明るさは、その線の端にある送光プローブと受光プローブのセットからなるチャネルにおける光強度変化情報を表している。図28aにおいて、広がった頭皮血流アーチファクトによって、全体的に値が出ていることがわかる。しかし、図28aを注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。また、波長830nmにおける観測値Yに関して、Ymaxを図示したのが、図28bである。図28bにおいて、広がった頭皮血流アーチファクトによって、全体的に値が出ている。しかし、図28bを注意深く観察すると、脳活動がある上半分のチャンネルの方が、観測値が大きいことは分かるが、どこで脳活動が起こっているかは見た目には分からない。なお、780nmよりも強く値が出ているが、これはシミュレーション設定によるものである。次に、推定手段283は、感度情報格納手段281からヘモグロビン濃度変化に関する感度行列上述したGcoxy,Gcdeoxy,Gsoxy,Gsdeoxyを読み出す。次に、推定手段283は、大脳皮質活動モデル情報格納部21から大脳皮質活動モデル情報例えば、局所的フィルタWcを読み出す。また、推定手段283は、頭皮血流変化モデル情報格納部22から頭皮血流変化モデル情報例えば、広域的フィルタWsを読み出す。次に、推定手段283は、変分ベイズアルゴリズムを用いて、Xを推定する。つまり、推定手段283は、予め格納している初期値を読み出す。初期値は、ここでは、すべてのボクセルに対してγα0i=0、かつ、すべてのボクセルに対してα0i-=10000、かつ、γβ0=0、かつ、β0-=10000である。そして、推定手段283は、上述したステップ1とステップ2の繰り返し処理を1000回実行し、Xを算出した。推定手段283が算出したXの要素であり、頭皮深さ1のオキシヘモグロビンの濃度変化Xoxyの推定値の集合を図29に示す。図29において、各セルの値は、各ボクセルに対応する。したがって、図29において、11×11の値が存在する。また、数値は行下方向が頭部モデルの横方向、列右方向が頭部モデルの縦方向と対応する。かかることは、図30、図31等と同様である。また、推定手段283が算出したXの要素であり、大脳皮質深さ4のオキシヘモグロビンの濃度変化Xoxyの推定値の集合を図30に示す。なお、図示しないが、大脳皮質深さ5のオキシヘモグロビンの濃度変化Xoxyについても、推定手段283は、同様に算出した。また、推定手段283が算出したXの要素であり、頭皮深さ1のデオキシヘモグロビンの濃度変化Xdeoxyの推定値の集合を図31に示す。さらに、推定手段283が算出したXの要素であり、大脳皮質深さ4のデオキシヘモグロビンの濃度変化Xdeoxyの推定値の集合を図32に示す。なお、図示しないが、大脳皮質深さ5のデオキシヘモグロビンの濃度変化Xoxyについても、推定手段283は、同様に算出した。次に、出力部29は、推定部28が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する。出力部29の出力例を図33、図34に示す。図33aは、頭皮深さ1のオキシヘモグロビンの濃度変化の出力例である。図33a等において、縦軸は頭部モデルの縦方向、横軸は頭部モデルの横方向を示し、その数字の単位はmmである。また、白黒の濃淡で表されるXの値の単位はmMである。図33bは、大脳皮質深さ4のオキシヘモグロビンの濃度変化の出力例である。また、図34aは、頭皮深さ1のデオキシヘモグロビンの濃度変化の出力例である。図34bは、大脳皮質深さ4のデオキシヘモグロビンの濃度変化の出力例である。以上、本実施の形態によれば、精度高く脳活動を推定できる。また、変分ベイズ法を用いることにより、極めて精度高く脳活動を推定できる。なお、本実施の形態の実験によれば、780nmと830nmの2波長を用いた。しかし、異なる波長を2波長以上用いれば良い。また、本実施の形態において、頭皮血流変化と脳活動の分離は、高密度プローブ計測によってはじめて可能となる。短い距離のチャネルは浅い頭皮の活動を拾いやすく、長い距離のチャネルは深い脳活動の活動をも拾うからである。つまり、本実施の形態において、プローブの配置密度が1.5cm以下であることは極めて好適である。また、プローブの配置密度が1.7cm以下であることは好適である。なお、プローブの配置密度は2cm以下、または3cm未満等でも良い。また、本実施の形態によれば、頭部が半無限一層モデルであると仮定し、計測データから脳活動を推定した。しかし、脳活動情報出力システムにおいて、MRI構造画像などを用いて、頭部が、頭皮頭蓋骨脳脊髄液灰白質などからなる多層モデルであると仮定し、計測データから脳活動を推定しても良い。また、本実施の形態によれば、推定部28が行う推定法は、主として、変分ベイズ法であった。ただし、推定部28は、逆問題を解くための他の方法を用いても良い。例えば、推定部28は、数式31を格納しており、当該数式31を用いて、Xを推定しても良い。つまり、推定部28は、例えば、数式31のようなコスト関数Cを最小にするようなXを算出すればよい。なお、数式31において、1項目が誤差項、2項目以降が正則化項と呼ばれる項であり、2項目以降が大きな値にペナルティを与える項になっている。数式31において、大脳皮質Xcには値ノルムの1乗、つまりL1ノルム、頭皮Xsには値ノルムの2乗、つまりL2ノルムの正則化を与えている。後者のL2ノルムでの正則化は、我々が図22の2で頭皮活動に課した正規分布の条件とほぼ同じ働きをする。また、一般に広がった推定値が得られやすいことが知られている。それに対し、前者のL1ノルムでの正則化を与えると、スパース局在したな推定値が得られることが知られている。なお、数式31において、フィルタWはかけていないが、上のλのパラメータや、正則化のノルムを0.9乗にする等の細かい値の修正をすれば、フィルタWをかけた結果と似たようなことが再現できる。以上から、数式31を用いても、変分ベイズ法を用いた場合と似たようなXの推定は可能である、と言える。さらに、本実施の形態における処理は、ソフトウェアで実現しても良い。そして、このソフトウェアをソフトウェアダウンロード等により配布しても良い。また、このソフトウェアをCD-ROMなどの記録媒体に記録して流布しても良い。なお、このことは、本明細書における他の実施の形態においても該当する。なお、本実施の形態における脳活動情報出力装置を実現するソフトウェアは、以下のようなプログラムである。つまり、このプログラムは、記憶媒体に、局所的な活動である大脳皮質の活動のモデルを示す情報である大脳皮質活動モデル情報と、広域的な活動である頭皮の血流変化のモデルを示す情報である頭皮血流変化モデル情報と、送光プローブと受光プローブの組であるプローブセットを1組以上有する近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、安静状態において測定された各受光プローブの光強度に関する情報である安静時光強度情報を、複数の各波長に対応付けて格納しており、コンピュータを、前記近赤外分光計測装置を用いて、前記各プローブセットを構成する各送光プローブから、一の光強度の光を発した場合に、脳活動が行われている状態において測定された各受光プローブの光強度に関する情報である脳活動時光強度情報を、複数の波長ごとに取得する脳活動時光強度情報取得部と、前記記憶媒体の安静時光強度情報と、前記脳活動時光強度情報取得部が取得した脳活動時光強度情報とから、複数の波長ごと、およびプローブセットごとに、光強度の変化に関する情報である光強度変化情報を取得する光強度変化情報取得部と、前記複数の波長ごと、およびプローブセットごとの光強度変化情報を、前記大脳皮質活動モデル情報および前記頭皮血流変化モデル情報に適用し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定部と、前記推定部が取得した大脳皮質のオキシヘモグロビンの濃度変化、および大脳皮質のデオキシヘモグロビンの濃度変化を出力する出力部として機能させるためのプログラム、である。また、上記プログラムにおいて、前記推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式を格納し得る演算式格納手段と、前記演算式に、前記感度情報と前記光強度変化情報を代入し、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を、ベイズ推定を用いて取得する推定手段とを具備するものとして、コンピュータを機能させることは好適である。また、上記プログラムにおいて、前記推定部は、大脳皮質または頭皮の活動の変化と検出された光強度の変化の関係を示す情報である感度情報を格納し得る感度情報格納手段と、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得するための演算式であり、コスト関数である演算式を格納し得る演算式格納手段と、前記演算式に、前記感度情報と前記光強度変化情報を代入し、当該演算式を実行し、実行結果であるコストが最小になるように、大脳皮質と頭皮の血管のオキシヘモグロビンの濃度変化、および大脳皮質と頭皮の血管のデオキシヘモグロビンの濃度変化を取得する推定手段とを具備するものとして、コンピュータを機能させることは好適である。また、図23は、本明細書で述べたプログラムを実行して、上述した実施の形態の脳活動情報出力装置等を実現するコンピュータの外観を示す。上述の実施の形態は、コンピュータハードウェア及びその上で実行されるコンピュータプログラムで実現され得る。図23は、このコンピュータシステム340の概観図であり、図24は、コンピュータシステム340のブロック図である。図23において、コンピュータシステム340は、FDドライブ、CD-ROMドライブを含むコンピュータ341と、キーボード342と、マウス343と、モニタ344とを含む。図24において、コンピュータ341は、FDドライブ3411、CD-ROMドライブ3412に加えて、MPU3413と、CD-ROMドライブ3412及びFDドライブ3411に接続されたバス3414と、ブートアッププログラム等のプログラムを記憶するためのROM3415とに接続され、アプリケーションプログラムの命令を一時的に記憶するとともに一時記憶空間を提供するためのRAM3416と、アプリケーションプログラム、システムプログラム、及びデータを記憶するためのハードディスク3417とを含む。ここでは、図示しないが、コンピュータ341は、さらに、LANへの接続を提供するネットワークカードを含んでも良い。コンピュータシステム340に、上述した実施の形態の脳活動情報出力装置等の機能を実行させるプログラムは、CD-ROM3501、またはFD3502に記憶されて、CD-ROMドライブ3412またはFDドライブ3411に挿入され、さらにハードディスク3417に転送されても良い。これに代えて、プログラムは、図示しないネットワークを介してコンピュータ341に送信され、ハードディスク3417に記憶されても良い。プログラムは実行の際にRAM3416にロードされる。プログラムは、CD-ROM3501、FD3502またはネットワークから直接、ロードされても良い。プログラムは、コンピュータ341に、上述した実施の形態の脳活動情報出力装置等の機能を実行させるオペレーティングシステムOS、またはサードパーティープログラム等は、必ずしも含まなくても良い。プログラムは、制御された態様で適切な機能モジュールを呼び出し、所望の結果が得られるようにする命令の部分のみを含んでいれば良い。コンピュータシステム340がどのように動作するかは周知であり、詳細な説明は省略する。なお、上記プログラムにおいて、ハードウェアによって行われる処理、は含まれない。また、上記プログラムを実行するコンピュータは、単数であってもよく、複数であってもよい。すなわち、集中処理を行ってもよく、あるいは分散処理を行ってもよい。また、上記各実施の形態において、各処理各機能は、単一の装置システムによって集中処理されることによって実現されてもよく、あるいは、複数の装置によって分散処理されることによって実現されてもよい。本発明は、以上の実施の形態に限定されることなく、種々の変更が可能であり、それらも本発明の範囲内に包含されるものであることは言うまでもない。"