pooneilの脳科学論文コメント

■ Psychometric functionの書き方

やっとわかった。こんなこともわからないでいたのです。Psychophysicsの基本中の基本にもかかわらず。

Psychophysicsではdetectionの閾値を決定するのにpsychometric functionを書きます。たとえば図のようなものです。たとえば、二つのターゲットのどちらかが点灯するのでどっちが点灯しているか当てなければならないとします。横軸はターゲットの明るさで、右が明るくて左が暗くなります。縦軸は正答率です。ターゲットのコントラストが明るいときには正答率はほぼ1ですが（右側の青丸）、暗くなるとどっちが点灯しているかわからなくなるので正答率はchance levelの0.5まで落ちます（左側の青丸）。これをlogistic regressionによってfittingしてやって、縦軸の0.75を通るところをその人の検出能の閾値としてやるわけです。しかしlogistic regressionというのは普通は縦軸が0-1の範囲を動くものなので、二択の時にはそれを0.5-1の範囲を動くようにmodifyしてやらなければならないのですが、そのやり方がよくわからなかったのです。

そこで参考サイト："Matlab and Sas logistic"

というわけで、matlabでglmfitを使い、inline関数でfittingしたい関数の形を設定してやればよかったのです。

そこでまず、inline関数で2択(2AFC: two alternative forced choice)、四択(4AFC: four alternative forced choice)、yes-no question (正答率0-1の範囲)のそれぞれの場合のlink関数を設定してやります。

function mylink = choose_link(choice)
% choose link function
% 0:yes-no 0-1;
% 2: 2AFC 0.5-1;
% 4: 4AFC 0.25-1

switch(choice)
  case 1,
    mylink='logit';
  case 2,
    fl=inline('log( (max(0.5+1e-6, min(p, 1-1e-6))-0.5) ./ ...
                 (1-max(0.5+1e-6, min(p, 1-1e-6))))')
    fd=inline('(1-0.5) ./ ( (1-max(0.5+1e-6, min(p, 1-1e-6))) .* ...
                  (max(0.5+1e-6, min(p, 1-1e-6))-0.5))')
    fi=inline('(0.5 + exp(x)) ./ (1+exp(x))') 
    mylink={fl fd fi};
  case 4,
    fl=inline('log( (max(0.25+1e-6, min(p, 1-1e-6))-0.25) ./ ...
                 (1-max(0.25+1e-6, min(p, 1-1e-6))))')
    fd=inline('(1-0.25) ./ ( (1-max(0.25+1e-6, min(p, 1-1e-6))) .* ...
                  (max(0.25+1e-6, min(p, 1-1e-6))-0.25))')
    fi=inline('(0.25 + exp(x)) ./ (1+exp(x))') 
    mylink={fl fd fi};
end

これをchoose_link.mとして保存しておいて、

dat=[0.6021  4  9
     0.6990 12 25
     0.7782  8 14
     0.8451 14 16
     0.9031 11 11
     1.0000  9  9
     1.0792 11 11
     1.1761 10 10];

% [x y Ny]   variable success total
x =dat(:,1);
y =dat(:,2);
Ny=dat(:,3);

mylink = choose_link(2);
[b,dev,stats]=glmfit(x,[y Ny],'binomial',mylink);
xx=min(x):(max(x)-min(x))/100:max(x);
yfit = glmval(b,xx,mylink);

figure
hold on
h1=plot(x,y./Ny,'o');
h2=plot(xx,yfit);
set(h1,'markersize',10)
set(h1,'linewidth',2)
set(h2,'linewidth',2)

データ(dat)からlogistic regressionしてやると図のようなpsychometric functionが書けるという訳でした。

もう少し説明しますと、入力するデータdatとして、一列目がvariable (横軸として振った値)、二列目がsuccess trial数、三列目がtotal trial数という行列を用意します。これをglmfitに入れるわけですが、logistic regressionでは四番目の引数にリンク関数を入れる必要があります。これが'logit'なのが普通の0-1の範囲を動くデータでのlogistic regressionです。そこでchoose_link functionでこのリンク関数をinline関数mylinkとして設定してやるとregressionができて、coefficient bが出てくるのでこれでデータの範囲を100等分した新しい横軸のデータでglmvalを使ってfittingの値を得る、というわけでした。

何とかここまで来れましたが、要するにわたし、inline関数とかeval関数とかのような文字列を評価して計算する、というのがいまだによくわかってないのですね。

glmfitのhelpにたしかに関連することが書いてあるのですが、結局のところ、どのような関数を設定するかは人のサイトのほぼ丸写しでした。うーむ、敗北。

なお、参考ページを見ていただければおわかりのとおり、psychometric functionを書くためにはいくつかのやり方があります。たとえば、今回使ったようなprobit関数 (log(p/(1-p))によるのではなくて、gaussian分布の累積確率密度関数を使ってやるとか。というわけで細かく言うときりないようなのでこのあたりにて。