知覚のベイズ統計的モデル:感覚強度曲線(1)

-Tue, Feb 5, 2019-

ざっくりいうと

感覚を測定する

かねてから認知心理学など、人間の知覚に関する分野では、知覚を測る試みがなされてきた。 近年はVRの目覚ましい発展により、多くの開発者・研究者が、VR空間で感覚を提示する上での様々な技術を創り出してきた。 高度なマルチモーダルインタフェースを作る上で、人間の感覚をある程度理解した上で設計できれば、様々なシチュエーションで提示する感覚をシステマチックにデザインできる。

本稿では筆者が卒論で少しかじった、感覚強度の測定における統計的モデリングについて、備忘録も兼ねていくつかの記事に分けながら残しておく(内部向けにいろいろ残してはいたがなんかもったいない気がした)。 なお間違っている点など多々あるかと思われるので、クレームは@JotaroUTにて謹んでお受けいたします。

感覚が生起することとは

人間の感覚は、人間の感覚器官への感覚刺激の入力が脳内のニューラルネットワークで処理されることで、主観的な感覚へと帰着される。視覚や聴覚を司る脳の受容野は、(驚くべきことに)入力を物理的・機械的に処理するための規則正しいニューラルネットが分布していることがわかっており、特に低次の視覚受容野を模して作られたニューラルネットワークとしてCNN(Convolutional Neural Network)がよく知られている。

意識レベルでの感覚の生起はさらに脳内での統合・認知や比較といった高レベルのプロセスを経ることによって行われるとされている。 [自主規制]したマウスの脳に[自主規制]して[自主規制]なりすれば、ニューロンの発火に伴う電気信号を直接観測することも可能だが、それでもなお高次の認知プロセスを物理レベルで測定するのは、fMRIなどのよっぽど高度な機材を使ったとしても定量的に評価することは難しい（しかし、すでにこういう研究成果も出てきてはいるので、実は必ずしもとても難しいという時代ではないのかもしれない）。

知覚測定と心理統計学

よって、直接観測することの難しい感覚は一般的な実験室では、ある被験者が認知した結果を回答してもらい、その結果をまとめて分析・比較する手法が取られる。 実験の確度を増すために被験者を増やして多くのデータを取り、それらを正しく分析するために 統計学 の知識が用いられる。特に心理学の分野において統計心理学として応用されている分野にあたる。 統計心理学のプロセスには統計処理のみならず、実験の方法やセットアップ、感覚刺激の提示方法や順序など、実験の 統制(=control) をするためのノウハウを組み込んでいる。

感覚の強度を測定する実験として、19世紀頃ドイツの学者のウェーバー(Ernst Heinrich Weber)やフェヒナー(Gustav Fechner)によって、「ある基準の刺激に対して、提示した刺激が、より”大きい”と答えた確率を調べる」というものが初めて行われた。刺激量を大きくしたときの確率が50%のときの刺激量を「弁別閾」と呼ぶ。また、一般的に刺激量を増やしていったときの回答確率は、なめらかなS字のカーブを描くことがわかっている。

Weber-Fechner -人間の知覚強度に関する “ウェーバー・フェヒナーの法則”は有名。比較刺激の弁別閾は標準刺激の量に対して対数的に増加する。

頻度主義的な感覚強度推定：2AFC・ベルヌーイ分布

感覚刺激の統計的モデル

人間の感覚強度の実験において、ある一定の刺激量 に対する平均回答確率を推定したいとする。 例えば、次のようなシンプルな実験を考える。

基準となるおもり(=標準刺激)に対して、わずかに重いおもり(=比較刺激)を用意し、どちらかが重いかを回答してもらう。重い方の重りに対して回答のあった割合を求める。

このような場合、感覚処理の過程で、ある一定の確率で重い方を重いと回答する と仮定する。このとき提示刺激に対する回答$x=1$が得られる確率を$p$とすると、

$$ x = \begin{cases} 1, & \text{with probability $p$} \\
0, & \text{with probability $1-p$} \end{cases}$$

このように、一定の確率で$X=1$が得られ、それ以外の場合は$X=0$となる確率分布はベルヌーイ分布として知られている。つまりこの実験において「人間の重さの感覚処理はベルヌーイ分布に従う」と仮定している。$X$を確率変数として記述すると、

$$ {\rm Bern}(X | p) = \begin{cases} p, & \text{for $X = 1$} \\
1-p, & \text{for $X = 0$} \end{cases}$$

つまり

$$ {\rm Bern}(X|p) = p^X (1-p)^{1-X} $$

この確率$p$こそ、「刺激が$p$の確率で識別できるに足りる感覚の強度」として捉えることができる。先のウェーバー・フェヒナーの実験は、この確率をもとに人間の感覚強度を測定・推定している。

ベルヌーイ分布の(頻度主義的)最尤推定

頻度主義的に考えれば、得られたデータから、先に仮定したベルヌーイ分布の真のパラメータ$p$を推定するためには尤度最大化のアプローチで最尤推定を用いる(ベルヌーイ分布のようなシンプルなモデルの場合は、直感的には得られた実際のデータの回答確率$\hat{p}$が最尤推定だろうと感じられる)。実際に確かめてみる。

得られたデータ列$\mathcal{D}=\{ X_1, X_2, … ,X_n \}$について尤度関数を$L(p|\mathcal{D})$とすれば、

$$ \begin{eqnarray} L(p|\mathcal{D}) & = & \prod_{i=1}^n {\rm Bern}(X_i|p) \\
\frac{\partial}{\partial \hat p} & = & \sum_{i=1}^n\ln{\rm Bern}(X_i|\hat p) \\
& = & \sum_{i=1}^n \frac{X_i}{\hat p} - \frac{(1-X_i)}{1- \hat p} = 0 \\
\iff \hat p & = & \frac{1}{n}\sum_{i=1}^nX_i \end{eqnarray}$$

確かに、真のパラメータ$p$に対する最尤推定は提示刺激に対する回答確率をもって最尤推定できる。

ベイズ主義的な感覚強度推定：2AFC・ベルヌーイ分布

心理統計学の分野では、先程のような実験のように不確実性を含んだ問題をとく上で、頻度主義な統計手法だけではなく、ベイズ的な統計を扱うこともある。講談社の”ベイズ推論による機械学習入門”では、ベイズ統計を問題に応用する上での利点として

などといった点をあげている。知覚モデルを構築・推論するといったプロセスがベイズ的アプローチで記述できるという点だけはなく、認知心理学における事前知識をベイズ推論に組み込むこともできるといった利点もある。特に、感覚強度測定のように多くのデータ点を取らなければならない実験系では、測定した結果がどの程度信頼できるものなのかを定量的に知ることも重要となってくる。

確率推論に基づく解析手法では「袋$a$か$b$か」といった２択の結論をだすのではなく、それぞれの原因がどれほどの確かさであったのかを定量的に表すことができます。

改めて強調しておくと、感覚強度のように実世界で様々な影響を排除できない不確実なパラメータを推定する上で、その結果がどれほど信頼できるものかを定量的に評価することは、主観的認知プロセスを主とするシステムをデザインする上で重要なファクターと言える。たとえば未来のVRChat等の多くのユーザが関わるひとつの広大なVRシステムに、実は人間の知覚モデルを取り入れたフィードバックループが構築されていて、パラメータの信頼区間などが考慮されている、強いて言えば処理系がモデルと処理結果の因果関係を記述的に推論しているシステムがあるとすれば、なんか素敵ではないだろうかと思う。

ベイズの公式とパラメータ推論

言わずとしたベイズの公式は、データ$\mathcal{D}$とパラメータ$\theta$を確率変数とした場合、

$$ \begin{eqnarray} p(\theta|\mathcal{D}) &=& \frac{p(\theta)p(\mathcal{D}|\theta)}{p(D)} \\
&\propto& p(\theta)p(\mathcal{D}|\theta) \end{eqnarray} $$

として表される。事後分布(posterior)$p(\theta | D)$は尤度(likelihood)$p(\theta)$と事前分布(prior)$p(D|\theta)$を掛け合わせて正規化した結果で表されるというものだ。ここでハイパーパラメータとなるのは_初期の事前分布_である。極端な話よくわからなければとりあえず一様分布にしたくなるが、今回の知覚実験のように、ベルヌーイ分布という（強めの）仮定をしている場合は、事前分布の選択肢は絞られる。

ベルヌーイ分布に対する共役事前分布 (Conjugate Prior)

事後分布の関数形が、事前分布に尤度を掛け合わせた場合と同じ関数形となる場合、そのように選んだ事前分布を(自然な)共役事前分布(Conjugate prior)と呼ぶ。 ベルヌーイ分布の場合はベータ分布と呼ばれる0から1までの実数を返す分布がその共役事前分布にあたる。 $$\begin{eqnarray} \beta(\theta|a, b) &=& \frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)} \\
B(a,b) &=& \int_0^1 t^{a-1} (1-t)^{b-1}dt \, \text{(Beta function)} \end{eqnarray}$$

事後分布の計算

知覚実験によって得られたデータ$\mathcal{D}$によって計算される尤度関数は、各試行が独立の場合、

$$ p(\mathcal{D}|\theta) = \prod_{i=1}^np(X_i|\theta) $$ となる。 共役事前分布の計算の中でも特にわかりやすいのでネットで検索するといろいろ出てくるが、実際にベルヌーイ分布と掛け合わせる場合は対数を用いて、

$$\begin{eqnarray} \ln p(\theta|\mathcal{D}) & \propto & \left(\sum_{i=1}^N x_i + a - 1 \right)\ln \theta + \left(N - \sum_{i=1}^N + b - 1\right) \ln(1-\theta) + {\it const.} \\
& = & (\hat a -1) \ln \theta + (\hat b -1) \ln(1-\theta) + {\it const.}\\
\end{eqnarray}$$

となり、

$$\begin{eqnarray} \hat a &=& \sum_{i=1}^N x_i + a \\
\hat b &=& N - \sum_{i=1}^N + b \end{eqnarray}$$

とすれば、事後分布の関数形はたしかに新しいβ関数 $p(\theta|\mathcal{D})=beta(\theta|\hat a , \hat b)$となっている。

逐次推論をするときに繰り返し計算となる場合にこのような単純な変数の更新で済むとプログラムに書き起こすときなどにかなり簡単になって嬉しい一方、共役事前分布を使わないと繰り返し数値計算をするはめになったりしてめんどくさいことになる。しかしながら共役事前分布が解析的に見つかることのほうが珍しいらしいので、統計的数値計算を用いた方法もいろいろと存在する。

予測分布の計算

事後分布の計算では観測したデータに基づいてパラメータを確率的に推論したが、未知のデータ$X_*$に対する尤度を求めた事後分布で重み付け計算(=周辺化)した予測分布を導出できる。つまり確率的に推論した結果をもって、ベルヌーイ分布の予測分布を計算する。

$$ p(X_* | X) = \int {\rm Bern}(X_* | \theta ) \beta ( \theta, | \hat a, \hat b, X) d\theta $$

周辺化はこの資料の式9.31あたりを参考にすると、

$$ p(X_*| X) = {\rm Bern} \left( X_* | \frac{\hat a}{\hat a + \hat b} \right) $$

という計算に帰着する。このベルヌーイ分布が、事前分布を$\beta(\theta|a, b)$として仮定し学習させた場合、新たに来るデータ＝確率変数に対する予測分布となるはず。

以上のように、同一刺激下でのベイズ知覚モデルが設計できた。

実験

実際に被験者実験を想定したデータ収集の状況下で感覚強度のベイズ推論をしてみる。

実験設定

被験者は10名。標準刺激に対して対象となる比較刺激を一人あたり25回提示する。 (実際の実験室では)標準刺激と比較刺激はそれぞれどちらかを先に提示し、１施行あたりランダムな順番で提示するものとする。 10名には被験者内誤差があり、$\theta = 0$の周辺で標準誤差$\sigma = 0$で分布しているとし、スクリプトでi.i.dのデータを生成する。

Pythonによるコード

# coding:utf-8
import numpy as np
import seaborn
import matplotlib.pyplot as plt
from numpy.random import *
from scipy.stats import bernoulli, beta

%matplotlib inline

# experiment parameters
num_subjects = 10
num_trials = 25
mean_prob = 0.75
var_subs  = 0.05

# virtual experiment
N = num_subjects * num_trials
sum_aye=0
params = []
for _ in range(10):
    theta = min(1, normal(mean_prob, var_subs))
    params.append(theta)
    sum_aye += sum(filter(lambda n:n%2==1, bernoulli.rvs(theta, size=num_trials)))

# inference
init_params = [2.5,2.5]
post_params = [sum_aye + init_params[0], N-sum_aye+init_params[1]] # = (a hat, b hat)
x = np.linspace(0, 1, 1000)
prior = beta(init_params[0], init_params[1])
posterior = beta(post_params[0], post_params[1])
y = posterior.pdf(x)
y_prior = prior.pdf(x)

# plot
f, (ax1, ax2) = plt.subplots(1, 2, sharey = True)
f.set_figheight(3)
f.set_figwidth(10)

ax1.set_title('prior distribution')
ax1.grid(color='gray', linestyle='-', linewidth=1, alpha = 0.2)
ax1.plot(x, y_prior)
ax1.fill(x, y_prior, 'blue', alpha=0.1)
ax1.vlines(params, 0, 1, color='k', linewidth=0.5)
ax1.vlines(0.75, 0, 1, color='r', linewidth=2.0)


ax2.set_title('posterior distribution')
ax2.grid(color='gray', linestyle='-', linewidth=1, alpha = 0.2)

ax2.plot(x, y)
ax2.fill(x, y, 'blue', alpha=0.1)
ax2.vlines(params, 0, 1, color='k', linewidth=0.5)
ax2.vlines(0.75, 0, 1, color='r', linewidth=2.0)

f.savefig('plot.svg')

事後分布の結果

グラフに中心パラメータ$(\theta=0.75)$と、各被験者に割り当てられたパラメータをそれぞれ赤線と黒線で示した。事後分布は施行を経るごとに一つの値に収束していく事がわかる。試行回数を増やせば、事後分布はさらに収束することもわかる。

Plot N=10、各25回施行

Plot N=10、各100回施行

予測分布の結果

予測分布は$(\theta=0.75)$に近い値となった。被験者間誤差に寄る多少のばらつきがあるが、パラメータは被験者間平均に近い値を予測できていることがわかる。

# predictive distribution
post_theta = post_params[0]  / (post_params[0]+post_params[1])
post_theta

# probability for X=1 on predicted distribution
# result(N=25, sample:25)   0.73529411764705888
# result(N=25, sample:100)  0.74079601990049748
# result(N=25, sample:1000) 0.75817091454272867

まとめ

感覚強度の仮想実験にベイズ統計を用いたパラメータ推定を考えてみた。機械学習や人工知能の問題がベイズ推論で記述されているのと同じく、このような心理統計学の分野での問題もベイズ推論が用いられるというのは、知能や認知とベイズ推論に何かしら根本的なつながりを感じさせるものだと思う。 今回は限定的に、実験において一定の刺激のみを提示したが、実際に弁別閾を測定したり、感覚強度曲線をプロットするときは、曲線パラメータの推定が必要になる。 そのような場合は今回のように直接的にモデルをたてることが難しくなる場合がある。また、事後確のの計算や統計モデルの妥当性検討・逸脱度検定なども必要になるために、MCMCやBootstrap法といったテクニックも必要になってくる。次回の記事ではそのあたりに触れて、実際に感覚強度曲線のベイズ推論について書く。

参考文献