X-Git-Url: http://git.osdn.jp/view?a=blobdiff_plain;f=dev4%2Fpsychlops%2Fcore%2Fmath%2Futil.cs;h=f1aa8175993b0f49b6f3fd1c03232866e599d12a;hb=af114e9e36fd9c2cdf79f95313e0b8712e253ed3;hp=7f498d359280489b33e36dcbc38704b07f922625;hpb=cb8916a7a5cd929f57b3f9edd99209680db90546;p=psychlops%2Fsilverlight.git diff --git a/dev4/psychlops/core/math/util.cs b/dev4/psychlops/core/math/util.cs index 7f498d3..f1aa817 100644 --- a/dev4/psychlops/core/math/util.cs +++ b/dev4/psychlops/core/math/util.cs @@ -5,17 +5,50 @@ namespace Psychlops public static class Math { - public static readonly double PI = 3.14159265, E = 2.718281828459045; + public const double PI = 3.14159265, E = 2.718281828459045, LOG2E = 1.44269504088896340736; public static Random random_generator; static Math() { random_generator = new Random(); } + public static T max(T val1, T val2) where T : IComparable + { + return val1.CompareTo(val2) > 0 ? val1 : val2; + } + public static T min(T val1, T val2) where T : IComparable + { + return val1.CompareTo(val2) < 0 ? val1 : val2; + } + public static void shuffle(X[] array, int n) + { + int a; + X tmp; + for(int i = 1; i < n; i++){ + a = random(i + 1); + tmp = array[i]; + array[i] = array[a]; + array[a] = tmp; + } + } + + + public static double mod(double lhs, double rhs) + { + return lhs - System.Math.Floor(lhs/rhs)*rhs; + } public static double abs(double x) { return System.Math.Abs(x); } + public static double floor(double x) + { + return System.Math.Floor(x); + } + public static double ceil(double x) + { + return System.Math.Ceiling(x); + } public static double sin(double x) { return System.Math.Sin(x); @@ -28,10 +61,34 @@ namespace Psychlops { return System.Math.Tan(x); } + public static double asin(double x) + { + return System.Math.Asin(x); + } + public static double acos(double x) + { + return System.Math.Acos(x); + } + public static double atan(double x) + { + return System.Math.Atan(x); + } + public static double atan(double y, double x) + { + return System.Math.Atan2(y, x); + } + public static double atan2(double y, double x) + { + return System.Math.Atan2(y, x); + } public static double sqrt(double x) { return System.Math.Sqrt(x); } + public static double pow(double x, double y) + { + return System.Math.Pow(x,y); + } public static double exp(double x) { return System.Math.Exp(x); @@ -40,6 +97,15 @@ namespace Psychlops { return System.Math.Log(x); } + public static double log2(double val) + { + return log(val) * LOG2E; + } + public static double round(double val) + { + return System.Math.Round(val); + } + public static double radius(double x, double y) { return System.Math.Sqrt(x * x + y * y); @@ -49,6 +115,10 @@ namespace Psychlops { return (random_generator.NextDouble()); } + public static int random(int x) + { + return (int)((random_generator.NextDouble()) * x); + } public static double random(double x) { return (random_generator.NextDouble()) * x; @@ -63,6 +133,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