X-Git-Url: http://git.osdn.jp/view?a=blobdiff_plain;f=dev4%2Fpsychlops%2Fcore%2Fmath%2Futil.cs;fp=dev4%2Fpsychlops%2Fcore%2Fmath%2Futil.cs;h=338423ad5726cfd28567230e098c23db77353a0d;hb=c7bc82580bbb03f0dcec6f68178513cc94152dcc;hp=85134a9618a7637aa332f514788bc544a5d0b9d1;hpb=a4fbd5c94992b9e0e3e57d550706ba6b05bfd6a6;p=psychlops%2Fsilverlight.git

diff --git a/dev4/psychlops/core/math/util.cs b/dev4/psychlops/core/math/util.cs
index 85134a9..338423a 100644
--- a/dev4/psychlops/core/math/util.cs
+++ b/dev4/psychlops/core/math/util.cs
@@ -102,6 +102,131 @@ namespace Psychlops
 		{
 			return exp(- (x*x) / (2*sigma*sigma));
 		}
+
+		public static double normalDistibution(double x, double mu, double sigma)
+		{
+			return exp( -( (x-mu)*(x-mu) / (2*sigma*sigma) ) ) / sqrt(2*PI*sigma*sigma);
+		}
+
+		public static double cumulativeNormalDistibution(double x, double mu, double sigma)
+		{
+			return .5 + .5*Internal.GammaFunction.erf( (x-mu)/(sigma*sqrt(2.0) ) );
+		}
+
 	}
 
+	namespace Internal
+	{
+		public static class GammaFunction
+		{
+			/************ loggamma(x) -- gamma.c より再掲 *************/
+
+			static readonly double PI      = 3.14159265358979324;  /* $\pi$ */
+			static readonly double LOG_2PI = 1.83787706640934548;  /* $\log 2\pi$ */
+			static readonly double N       = 8;
+
+			static readonly double B0  = 1            ;     /* 以下はBernoulli数 */
+			static readonly double B1  = (-1.0 / 2.0);
+			static readonly double B2  = ( 1.0 / 6.0);
+			static readonly double B4  = (-1.0 / 30.0);
+			static readonly double B6  = ( 1.0 / 42.0);
+			static readonly double B8  = (-1.0 / 30.0);
+			static readonly double B10 = ( 5.0 / 66.0);
+			static readonly double B12 = (-691.0 / 2730.0);
+			static readonly double B14 = ( 7.0 / 6.0);
+			static readonly double B16 = (-3617.0 / 510.0);
+
+			public static double loggamma(double x)  /* ガンマ関数の対数 */
+			{
+				double v, w;
+
+				v = 1;
+				while (x < N) {  v *= x;  x++;  }
+				w = 1 / (x * x);
+				return ((((((((B16 / (16 * 15))  * w + (B14 / (14 * 13))) * w
+							+ (B12 / (12 * 11))) * w + (B10 / (10 *  9))) * w
+							+ (B8  / ( 8 *  7))) * w + (B6  / ( 6 *  5))) * w
+							+ (B4  / ( 4 *  3))) * w + (B2  / ( 2 *  1))) / x
+							+ 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x);
+			}
+
+			public static double p_gamma(double a, double x, double loggamma_a)  /* 本文参照 */
+			{
+				int k;
+				double result, term, previous;
+
+				if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
+				if (x == 0)     return 0;
+				result = term = Math.exp(a * Math.log(x) - x - loggamma_a) / a;
+				for (k = 1; k < 1000; k++) {
+					term *= x / (a + k);
+					previous = result;  result += term;
+					if (result == previous) return result;
+				}
+				//throw new Exception("p_gamma(): the sequence is not convergent.");
+				return result;
+			}
+
+			public static double q_gamma(double a, double x, double loggamma_a)  /* 本文参照 */
+			{
+				int k;
+				double result, w, temp, previous;
+				double la = 1, lb = 1 + x - a;  /* Laguerreの多項式 */
+
+				if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
+				w = Math.exp(a * Math.log(x) - x - loggamma_a);
+				result = w / lb;
+				for (k = 2; k < 1000; k++) {
+					temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
+					la = lb;  lb = temp;
+					w *= (k - 1 - a) / k;
+					temp = w / (la * lb);
+					previous = result;  result += temp;
+					if (result == previous) return result;
+				}
+				//throw new Exception("q_gamma(): the sequence is not convergent.");
+				return result;
+			}
+
+			public static double p_chisq(double chisq, int df)  /* カイ2乗分布の下側確率 */
+			{
+				return p_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));
+			}
+
+			public static double q_chisq(double chisq, int df)  /* カイ2乗分布の上側確率 */
+			{
+				return q_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));
+			}
+
+			static readonly double LOG_PI = 1.14472988584940017;  /* $\log_e \pi$ */
+
+			public static double erf(double x)  /* Gaussの誤差関数 ${\rm erf}(x)$ */
+			{
+				if (x >= 0) return   p_gamma(0.5, x * x, LOG_PI / 2);
+				else        return - p_gamma(0.5, x * x, LOG_PI / 2);
+			}
+
+			public static double erfc(double x)  /* $1 - {\rm erf}(x)$ */
+			{
+				if (x >= 0) return  q_gamma(0.5, x * x, LOG_PI / 2);
+				else        return  1 + p_gamma(0.5, x * x, LOG_PI / 2);
+			}
+
+			public static double p_normal(double x)  /* 標準正規分布の下側確率 */
+			{
+				if (x >= 0) return
+					0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));
+				else return
+					0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);
+			}
+
+			public static double q_normal(double x)  /* 標準正規分布の上側確率 */
+			{
+				if (x >= 0) return
+					0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);
+				else return
+					0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));
+			}
+		}
+	}
 }
\ No newline at end of file