Up

簡単な上下法と主観的等価点および弁別閾の推定

 

　主観的等価点（PSE）を求めるために用いられる上下法のプログラミング例と、上下法で得られたデータの分析プログラムを用意した。分析プログラムは恒常法データにも適用可能である。このプログラムを上下法データに適用すると、主観的等価点に加えて弁別閾の推定も行える。なお、閾の現代における基本的考え方については、＜岡本安晴「計量心理学」培風館、補足2.1　閾の考え方および関連事項、Pp. 74-82＞を参照されたい。

上下法サンプルプログラム

 

　主観的等価点（Point of Subjective Equality: PSE）を求める方法の1つである上下法（Up-Down method）では、標準刺激aを固定して、aと比べる比較刺激xの強さを被験者の反応に基づいて変化させる。第i試行における比較刺激をで表わすとき、は以下のように決められる。

               の判断ならば、

               の判断ならば、

最初の試行におけるxの値は、aとの比較判断が容易であるようにaから十分に離れた値のものが用いられる。それに応じての値も大きく設定されるが、試行が進むにつれては小さい値のものが選ばれる。しかし、の値をあまり小さくする必要はない。データの分析を、本ウェブサイトの後半で説明しているような確率モデルを設定して行う場合は、の値はJNDの３倍ぐらいの値でも、JNDの1/3ぐらいの場合と大差ないことが報告されている（岡本安晴、2006：日本心理学会第70回大会発表論文集、p. 559 ＜PDFファイル＞）。

　この上下法による弁別実験を行うサンプルプログラムがSimpleUpDown.exeである。プログラムファイルは実行形式のファイルであるので、Windowsパソコンであればクリックしてダウンロード・実行することができる。ただし、.NET Framework 3.5対応である。ソースファイルなどは圧縮ファイルSimpleUpDown.zipとしてまとめた。圧縮ファイルSimpleUpDown.zipは、クリックしてダウンロードしたものをマウスの右ボタンでクリックして表示されるメニュから「展開」あるいは「解凍」などを選べば解凍できる。解凍したファイルはVisual C++2008などで開くことができる。

　プログラムを実行すると図１のフォームが表示される。

図１

 

「GO」ボタンのクリックで実験が始まり、図２のフォームが表示される。

図２

 

中央の「1/48」は、総試行数が48試行であり、その内の第１試行目であることを表わす。実験は長さの比較判断を求めるものであるが、比較刺激が標準刺激より明らかに長い場合「」から始まる系列と、短い場合「」から始まる系列の２系列が用いられていて、それぞれ24試行の長さである。この２つの系列からの試行をランダムに混ぜて、系列内での依存効果が生じないようにしている。

　図２の「提示」ボタンをクリックすると、図３のように十字図形が提示される。

図３

 

十字の中心（交点）から横線の端までの長さが標準刺激の値aであり、中心から縦線の端までの長さが比較刺激の値である。図３の刺激の提示後、比較判断を求めるボタンが表示される（図４）。

図４

 

縦線の方が横線より長いときは「縦＞横」ボタンを、縦線の方が短いときは「縦＜横」ボタンを押す。比較判断のボタンをクリックすると図２の「提示」ボタンのクリック待ち画面に戻る（図５）。

図５

 

以後、＜「提示」ボタンのクリック、比較判断＞のサイクルを４８試行分繰り返す。４８試行が終了すると図６の画面になる。

図６

 

画面右上の「保存・閉じる」ボタンをクリックすると、図７の保存用ファイル名の設定ダイアログウィンドウが表示される。

説明: 無題.jpg

図７

 

適当なファイル名を設定して「保存」ボタンをクリックすると、メッセージボックスが表示される（図８）。

図８

 

メッセージボックスの「OK」ボタンのクリックでプログラムの実行終了となる。

　プログラムの終了後、「保存」ボタンのクリックで保存したファイルを開くと以下のようになっている。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

    1    1

    2    1

    3    2

    4    2

    5    1

    6    2

　　・

　　・

　　・

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

各行における最初の数字は試行数であり、２番目の数字は上下法における２つの系列のどちらによる刺激提示かを示すものである。この提示刺激の系列を示すデータに続いて、提示された比較刺激値とそれに対する比較判断結果が以下のように続く。比較判断結果はスラッシュ「/」で始まる行によって挟まれている。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　　・

　　・

　　・

   46    1

   47    2

   48    2

 

/

    1     50      0

    2     83      0

    3    116      1

    4     83      0

    5    105      1

　　・

　　・

　　・

   44    100      1

   45     98      0

   46    100      1

   47     98      0

   48    100      1

/

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

各行に、「試行数　　比較刺激値　　判断」の形式で各試行の結果が表わされている。判断は、比較刺激の方が長いという比較判断を「１」で、比較刺激の方が短いという比較判断を「０」で表わしている。

　この上下法のデータに対して、確率モデルを設定して分析を行い、主観的等価点（PSE）と弁別閾（Just Noticeable Difference: JND）を求めるプログラムを次に示す。

　なお、事後分布の観点から効率的に提示刺激値を決める適応的上下法のプログラムを別のウェブサイトに用意している。Visual C++プログラミングの解説は、＜岡本安晴「大学生のためのVC++プログラミング入門」勁草書房＞で行っている。

データ分析プログラム

 

　理論的には、主観的等価点（PSE）と弁別閾（JND）は精神測定関数（Psychometric Function; 心理測定関数とも訳される）

                           

に基づいて定義される（岡本安晴「計量心理学」培風館、2.2節「閾値と精神測定関数」参照）。一般に精神測定関数は、位置パラメータと傾きパラメータによって与えられる。これらのパラメータは、精神測定関数を累積正規分布関数で表わすときは平均と標準偏差となる。このとき、弁別確率７５％に対応する弁別閾JNDはで与えられる。いま、

                           

とおく。上の上下法の第i試行におけるデータを

　　i       　　

と表わす。ただし、

                           

である。

　パラメータの事前分布をで、事後分布をで表わすとき、

                           

が成り立つ。式から判るとおり、この分析モデルは、上下法データに限らず2AFCパラダイムを採用する適応的方法によるデータ、および恒常法データにも適用して用いることができる。

事前分布を有界であるが十分に広い範囲で一様である分布として事後分布を求めるプログラムがAnalPFData.exeである。このプログラムファイルAnalPFData.exeは実行形式であるので、Windowsパソコンであればクリックしてダウンロード・実行することができる。ただし、.NET Framework 3.5対応である。プログラムのソースファイルなどを圧縮ファイルAnalPFData.zipとしてまとめたものも用意した。圧縮ファイル名AnalPFData.zipをマウスの右ボタンでクリックして表示されるメニュから「解凍」あるいは「展開」などを選ぶと解凍することができる。解凍したファイルはVisual C++ 2008などで開くことができる。

プログラムを実行すると図９のフォームが表示される。

図９

 

「GO」ボタンのクリックで、入力データファイル名設定のダイアログウィンドウが表示される（図１０）。

図１０

 

入力データファイルの形式は、上の上下法の実験を行うプログラムの出力形式に対応している。ここで、上で指摘したように、データの入力形式を本プログラムに合わせておけば、本プログラムは恒常法データなどの分析にも用いることができることに注意。入力データは、スラッシュ「/」で始まる行によってデータの書かれた行を挟む。データは

                           

の形式で書く。はそのデータ対のラベルであり、上の上下法の場合は試行数がデータのラベルとして使われている。このデータのラベルはプログラムにおける計算では使用されない。は比較刺激値である。は比較判断を表わし、比較刺激の方が長いという判断を１で、短いという判断を０で表わす。長さ以外のモダリティの刺激のときは、それぞれ適した言葉、「強い」とか「重い」とかが対応する。以上の約束の下で、入力ファイルの形式は以下のようになる。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

                            コメント

 

/

              Label-1                                        

              Label-2                                        

                            ・

                            ・

                            ・

              Label-N                                      

/

 

                                          コメント

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

データ部分を挟む「/」で始まる行の前後はプログラムでは無視されるので、コメントを書くことができる。

図１０のダイアログウィンドウにおいて、分析するデータのファイル名を設定して「開く」ボタンをクリックする。計算結果などを出力するためのテキストファイルの名前の設定を求めるダイアログウィンドウが表示される（図１１）。

図１１

 

出力用テキストファイルの名前を適当に設定して「保存」ボタンをクリックすると計算が始まり、図１２の画面となる。

図１２

 

パラメータの事後分布の計算は、グリッドを設定して離散確率による近似によって行っている。PSEとJNDの同時分布は以下のようなコードで求めることができる。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

        double sum = 0.0;

        for (int i = 0; i <= NPoints; i++)

                for (int j = 0; j <= NPoints; j++){

                        double vPSE = minPSE + i * (maxPSE - minPSE) / NPoints;

                        double vSGM = minSGM + j * (maxSGM - minSGM) / NPoints;

                        double v = 0.0;

                        for (int k = 0; k < n; k++)

                                if (resp[k] == 1){

                                                v += myLog( PF(st[k], vPSE, vSGM) );

                                }else{

                                                v += myLog( 1.0 - PF(st[k], vPSE, vSGM) );

                                };

                        sum += gprob[i][j] = Math::Exp(v);

                }

        for (int i = 0; i <= NPoints; i++)

                for (int j = 0; j <= NPoints; j++)

                      gprob[i][j] /= sum;

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

しかし、プログラムAnalPFDataでは、マルチ・スレッド化してgprob[i][j]の計算を第i行ごとに別スレッドに分けて行っている。そのために参照クラス型CalcGProbを以下のように用意している。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

        ref class CalcGProb {

                public:

                        double ** gprob;        //      pointer to the original gprob.

                        int * resp;                     //      copy of the original *resp.

                        double * st;            //      copy of the original *st.

                        int i;                          //      the i-th row to be calculated.

                        int NPoints, n;

                        double minPSE, maxPSE, minSGM, maxSGM;

                        double sum;

 

                        CalcGProb(){ }

                        ~CalcGProb(){   DoIt(); }

                        !CalcGProb(){   DoIt(); }

 

                        void DoIt(){

                                if (resp != nullptr) delete [] resp;

                                if (st != nullptr) delete [] st;

                        }

 

                        void Init( double ** vgprob, int * vresp, double * vst, int vi, int vNPoints, int vn,

                                                double vminPSE, double vmaxPSE, double vminSGM, double vmaxSGM ){

                                i = vi;

                                NPoints = vNPoints;

                                n = vn;

                                gprob = vgprob;

                                resp = new int[n];

                                for (int k = 0; k < n; k++)

                                        resp[k] = vresp[k];

                                st = new double[n];

                                for (int k = 0; k < n; k++)

                                        st[k] = vst[k];

                                minPSE = vminPSE;

                                maxPSE = vmaxPSE;

                                minSGM = vminSGM;

                                maxSGM = vmaxSGM;

                        }

 

 

                         double myLog( double x ){

                                 double v;

                                 if ( x <= 0.0 ){

                                                v = -750.0;

                                 }else{

                                         v = Math::Log( x );

                                 }

                               

                                 return v;

                         }

 

                        void DoCalc(){

                                        sum = 0.0;

                                        for (int j = 0; j <= NPoints; j++){

                                                double vPSE = minPSE + i * (maxPSE - minPSE) / NPoints;

                                                double vSGM = minSGM + j * (maxSGM - minSGM) / NPoints;

                                                double v = 0.0;

                                                for (int k = 0; k < n; k++)

                                                        if (resp[k] == 1){

                                                                        v += myLog( PF(st[k], vPSE, vSGM) );

                                                        }else{

                                                                        v += myLog( 1.0 - PF(st[k], vPSE, vSGM) );

                                                        };

                                                sum += gprob[i][j] = Math::Exp(v);

                                        }

 

                        }

        };

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

　この参照クラス型を用いて、マルチスレッドによる計算を以下のように行っている。

 

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

        array<CalcGProb^> ^ calcMT;

        calcMT = gcnew array<CalcGProb^>(NPoints+1);

        for (int i = 0; i <= NPoints; i++){

                calcMT[i] = gcnew CalcGProb();

                calcMT[i]->Init( gprob, resp, st, i, NPoints, n, minPSE, maxPSE, minSGM, maxSGM );

        }

        array<ThreadStart^>^ myThreadStart;

        myThreadStart = gcnew array<ThreadStart^>(NPoints+1);

        array<Thread^>^ myThread;

        myThread = gcnew array<Thread^>(NPoints+1);

        for (int i = 0; i <= NPoints; i++){

                myThreadStart[i] = gcnew ThreadStart(calcMT[i], &CalcGProb::DoCalc);

                myThread[i] = gcnew Thread( myThreadStart[i] );

                myThread[i]->Start();

        }

 

        for (int i = 0; i <= NPoints; i++){

                myThread[i]->Join();

                sum += calcMT[i]->sum;

      }

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

 

計算が終了するとPSEの事後分布が表示される（図１３）。

説明: 無題.jpg

図１３

 

分布の下限と上限（図１３では50.0と150.0）は入力データの比較刺激値の分布から決められている。これらの値は、図１３のように描かれた事後分布から適切な値を推測して設定し直す。グリッドによる離散化の近似精度を高めるためには、事後分布の広がりは十分に広いことが望ましい。離散化は図に示されている下限（左端の値）と上限（右端の値）の区間を500等分する形で行われている。しかし、裾が広すぎて両端で事後分布の値が十分に０に近くないときは分布の裾が捨てられることになる。画面左上にある「Min(PSE) =」の下のテキストボックスに下限の値、「Max(PSE)　=」の下のテキストボックスに上限の値を設定し直して「PSE」ボタンをクリックすると、設定した値に基づいて再計算が行われる。これを試行錯誤的に繰り返して図１４のように分布が下限と上限の範囲に適切に収まるようにする。

図１４

 

事後分布の様子を表わす指標は、図の上部に表示されている値から、PSEの事後分布のモードが95.4、95％HPD区間（Highest Posterior Density Interval）が[93.2, 97.7]であることがわかる。95％HPD区間は、離散分布を折れ線で結んで連続分布に変換して求めている。このためのクラス型CalcHPDが次のように用意されている。

 

class CalcHPD {

        public:

                double * p;

                double maxH, mode, Lpoint, Rpoint;

                int n;

 

                CalcHPD( double * vp, int vn, double vmaxH, double vmode ){

                        p = vp;

                        n = vn;

                        maxH = vmaxH;

                        mode = vmode;

                }

 

                virtual ~CalcHPD(){ }

 

                void Init( double * vp, int vn, double vmaxH, double vmode ){

                        p = vp;

                        n = vn;

                        maxH = vmaxH;

                        mode = vmode;

                }

                double contDistri( double x );

                double cumDistri( double x );

                double Useek_point( double cum_p, double Lx, double Ux, double (*f)(double v) );

                double Dseek_point( double cum_p, double Lx, double Ux, double (*f)(double v) );

                double CalcIntvl( double h );

                void Seek95Intvl( );

};

 

このCalcHPD型のクラスのメソッドSeek95Intvlを呼び出すことによりHPD区間の端点がLpointとRpointに設定されるので、それらの値を用いて95%HPD区間が算出される。例えば、PSEの95%HPD区間は以下のコードで算出されている。

 

        CalcHPD hpd( probPSE, NPoints, probPSE[ip], ip );

        hpd.Seek95Intvl();

        double LHPD = minPSE + hpd.Lpoint * (maxPSE - minPSE) / NPoints;

        double RHPD = minPSE + hpd.Rpoint * (maxPSE - minPSE) / NPoints;

 

　「JND」ボタンをクリックするとJNDの事後分布が表示される（図１５）。

 

説明: 無題.jpg

図１５

 

　JNDの事後分布の下限の値と上限の値も適切な値を設定して「JND」ボタンをクリックすると設定した値で再計算が行われる。事後分布が下限と上限の間に適切に収まるように、試行錯誤で下限と上限の値を選んで、図１６のような状態になるようにする。

図１６

 

　事後分布の計算は、PSEとJNDの同時分布からそれぞれの周辺分布が算出されているので、図１６のJNDの値の範囲での再計算がPSEの事後分布を表わす図１４のグラフに影響を及ぼす可能性がある。JNDの下限と上限の設定値を変えたときは、その後PSEの事後分布に変化があるか確かめておく必要がある。「PSE」ボタンをクリックして再計算を行うと図１７のようになった。いまの場合、値の変化は認められない。

 

図１７

 

　PSEとJNDの下限と上限の適切な値がそれぞれ決まると、「PF」ボタンをクリックすれば、そのときのそれぞれの事後分布のモードを点推定値とするpsychometric functionを描くことができる（図１８）。

説明: 無題.jpg

図１８

 

図１８において青色の曲線がpsychometric functionである。赤色の曲線はデータにkernel smootherを適用して描いたカーネル回帰曲線である。このカーネル回帰曲線は、画面左側の「k-width =」の下のテキストボックスに試行錯誤的に値を設定、「PF」ボタンのクリックを行い、設定した値に対するカーネル回帰曲線を描くということを繰り返して、適切な曲線が得られるようにする。値が大きすぎると曲線は過度になだらかなものとなり（図１９）、小さすぎると変動が大きすぎるものとなる（図２０）。

図１９

 

図２０

 

　カーネル回帰曲線は、次のコードによって与えられる関数k_aveを描画したものである。

 

        double kernel( double x, double c, double w ){

                return Math::Exp(-0.5 * sqr((x - c) / w));

        }

 

 

        double k_ave( double x, int * vy, double * vx, int n, double w ){

                double v = 0.0, sw = 0.0;

                for (int i = 0; i < n; i++){

                        double kv = kernel( x, vx[i], w );

                        v += vy[i] * kv;

                        sw += kv;

                }

                if (sw > 0.0) v = v / sw;

                else              v = 0.0;

                return v;

        }

 

　なお、各事後分布およびPFの画面の図は「Save」ボタンのクリックでファイル（BMPファイル*.bmp）に保存することができる。

 

関連文献：岡本安晴「上下法データのベイズ的分析――確率モデルによるデータの活用――」基礎心理学研究、2011、30、44−55．

 

補足：　弁別閾の考え方が感覚尺度の構成において理論的に重要な役割を担うことは、Okamoto(2001)”On Fechner’s Problem”＜http://ci.nii.ac.jp/naid/110000032319/en/＞において指摘した。さらにOkamoto (2004) “Generalization of Fechner’s Approach” において、この弁別閾と感覚尺度との関係をWeber functionと対応する感覚尺度の関数形の問題として論じ（Pp. 113-115）、Weber functionを無視したDzhafarov and Colonius (1999) のモデルはStevensのべき法則と矛盾することを示した（Appendix B; Pp. 120-121）。

岡本安晴　作

Up