ちまたで有効な実験方法と名高い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:
Post a Comment