Friday, October 31, 2014

A revised table of d'(以下略)の式1を家で実装した

http://link.springer.com/article/10.3758%2FBF03208311

ちまたで有効な実験方法と名高いM肢強制選択の評価の関数ですが、プログラム可能な状態の自前ソースを持っていないと面倒なことが多いはずです。
なので、宅での暇な時間に書いてみました。

プログラムの結果と論文に書いてある結果とは一致しているので、今後いちいち表を参照して、手計算しなくてよくなる。

ソースをブログに載せる方法を勉強する気がないので、どこか削れているかもしれないけれど、プロならなんとなく補間できるはず。
以下、c++コード
#include <iostream>
#include <fstream>
#define _USE_MATH_DEFINES
#include <cmath>
//p = integral( [CDF(x)]^(M-1)*NDF(x-d) )dx

//http://www.johndcook.com/cpp_phi.html
//CDF
template <typename T>
T phi(T x)
{
 // constants
 T a1 =  0.254829592;
 T a2 = -0.284496736;
 T a3 =  1.421413741;
 T a4 = -1.453152027;
 T a5 =  1.061405429;
 T p  =  0.3275911;

 // Save the sign of x
 int sign = 1;
 if (x < 0)
  sign = -1;
 x = fabs(x)/sqrt(2.0);

 // A&S formula 7.1.26
 T t = 1.0/(1.0 + p*x);
 T y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);

 return 0.5*(1.0 + sign*y);
}

//http://stackoverflow.com/questions/10847007/using-the-gaussian-probability-density-function-in-c
template <typename T>
T normal_pdf(T x, T m, T s)
{
 static const T inv_sqrt_2pi = 0.3989422804014327;
 T a = (x - m) / s;
 return inv_sqrt_2pi / s * std::exp(-T(0.5) * a * a);
}

int main(int argc, char* argv[])
{
 std::ofstream o("test.csv");
 //p = integral( [CDF(x)]^(M-1)*NDF(x-d) )dx

 const int M = 10;
 const double dx = 0.001;
 o << "P" << "," << "di" << "," ;
 o<<"\n";
 double di;
 for(di=-4;di<4;di+=0.01){
  double P = 0;
  for(double x=-12;x<12;x+=dx){
   P += pow(phi(x),M-1)*normal_pdf(x-di,0.0,1.0)*dx;
  }
  o << P << "," << di << "," ;
  o<<"\n";
 }
 return 0;
}

No comments: