--- /dev/null
+using System;\r
+\r
+namespace Psychlops\r
+{\r
+\r
+ public static class Math\r
+ {\r
+ public const double PI = 3.14159265, E = 2.718281828459045, LOG2E = 1.44269504088896340736;\r
+ public static Random random_generator;\r
+ static Math()\r
+ {\r
+ random_generator = new Random();\r
+ }\r
+\r
+ public static T max<T>(T val1, T val2) where T : IComparable\r
+ {\r
+ return val1.CompareTo(val2) > 0 ? val1 : val2;\r
+ }\r
+ public static T min<T>(T val1, T val2) where T : IComparable\r
+ {\r
+ return val1.CompareTo(val2) < 0 ? val1 : val2;\r
+ }\r
+ public static void shuffle<X>(X[] array, int n)\r
+ {\r
+ int a;\r
+ X tmp;\r
+ for(int i = 1; i < n; i++){\r
+ a = random(i + 1);\r
+ tmp = array[i];\r
+ array[i] = array[a];\r
+ array[a] = tmp;\r
+ }\r
+ }\r
+\r
+\r
+ public static double mod(double lhs, double rhs)\r
+ {\r
+ return lhs - System.Math.Floor(lhs/rhs)*rhs;\r
+ }\r
+ public static double abs(double x)\r
+ {\r
+ return System.Math.Abs(x);\r
+ }\r
+ public static double floor(double x)\r
+ {\r
+ return System.Math.Floor(x);\r
+ }\r
+ public static double ceil(double x)\r
+ {\r
+ return System.Math.Ceiling(x);\r
+ }\r
+ public static double sin(double x)\r
+ {\r
+ return System.Math.Sin(x);\r
+ }\r
+ public static double cos(double x)\r
+ {\r
+ return System.Math.Cos(x);\r
+ }\r
+ public static double tan(double x)\r
+ {\r
+ return System.Math.Tan(x);\r
+ }\r
+ public static double asin(double x)\r
+ {\r
+ return System.Math.Asin(x);\r
+ }\r
+ public static double acos(double x)\r
+ {\r
+ return System.Math.Acos(x);\r
+ }\r
+ public static double atan(double x)\r
+ {\r
+ return System.Math.Atan(x);\r
+ }\r
+ public static double atan(double y, double x)\r
+ {\r
+ return System.Math.Atan2(y, x);\r
+ }\r
+ public static double atan2(double y, double x)\r
+ {\r
+ return System.Math.Atan2(y, x);\r
+ }\r
+ public static double sqrt(double x)\r
+ {\r
+ return System.Math.Sqrt(x);\r
+ }\r
+ public static double pow(double x, double y)\r
+ {\r
+ return System.Math.Pow(x,y);\r
+ }\r
+ public static double exp(double x)\r
+ {\r
+ return System.Math.Exp(x);\r
+ }\r
+ public static double log(double x)\r
+ {\r
+ return System.Math.Log(x);\r
+ }\r
+ public static double log2(double val)\r
+ {\r
+ return log(val) * LOG2E;\r
+ }\r
+ public static double round(double val)\r
+ {\r
+ return System.Math.Round(val);\r
+ }\r
+\r
+ public static double radius(double x, double y)\r
+ {\r
+ return System.Math.Sqrt(x * x + y * y);\r
+ }\r
+\r
+ public static double random()\r
+ {\r
+ return (random_generator.NextDouble());\r
+ }\r
+ public static int random(int x)\r
+ {\r
+ return (int)((random_generator.NextDouble()) * x);\r
+ }\r
+ public static double random(double x)\r
+ {\r
+ return (random_generator.NextDouble()) * x;\r
+ }\r
+ public static double random(double x, double y)\r
+ {\r
+ return (random_generator.NextDouble()) * (y-x) + x;\r
+ }\r
+\r
+\r
+ public static double gaussian(double x, double sigma)\r
+ {\r
+ return exp(- (x*x) / (2*sigma*sigma));\r
+ }\r
+\r
+ public static double normalDistibution(double x, double mu, double sigma)\r
+ {\r
+ return exp( -( (x-mu)*(x-mu) / (2*sigma*sigma) ) ) / sqrt(2*PI*sigma*sigma);\r
+ }\r
+\r
+ public static double cumulativeNormalDistibution(double x, double mu, double sigma)\r
+ {\r
+ return .5 + .5*Internal.GammaFunction.erf( (x-mu)/(sigma*sqrt(2.0) ) );\r
+ }\r
+\r
+ }\r
+\r
+ namespace Internal\r
+ {\r
+ public static class GammaFunction\r
+ {\r
+ /************ loggamma(x) -- gamma.c より再掲 *************/\r
+\r
+ static readonly double PI = 3.14159265358979324; /* $\pi$ */\r
+ static readonly double LOG_2PI = 1.83787706640934548; /* $\log 2\pi$ */\r
+ static readonly double N = 8;\r
+\r
+ static readonly double B0 = 1 ; /* 以下はBernoulli数 */\r
+ static readonly double B1 = (-1.0 / 2.0);\r
+ static readonly double B2 = ( 1.0 / 6.0);\r
+ static readonly double B4 = (-1.0 / 30.0);\r
+ static readonly double B6 = ( 1.0 / 42.0);\r
+ static readonly double B8 = (-1.0 / 30.0);\r
+ static readonly double B10 = ( 5.0 / 66.0);\r
+ static readonly double B12 = (-691.0 / 2730.0);\r
+ static readonly double B14 = ( 7.0 / 6.0);\r
+ static readonly double B16 = (-3617.0 / 510.0);\r
+\r
+ public static double loggamma(double x) /* ガンマ関数の対数 */\r
+ {\r
+ double v, w;\r
+\r
+ v = 1;\r
+ while (x < N) { v *= x; x++; }\r
+ w = 1 / (x * x);\r
+ return ((((((((B16 / (16 * 15)) * w + (B14 / (14 * 13))) * w\r
+ + (B12 / (12 * 11))) * w + (B10 / (10 * 9))) * w\r
+ + (B8 / ( 8 * 7))) * w + (B6 / ( 6 * 5))) * w\r
+ + (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x\r
+ + 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x);\r
+ }\r
+\r
+ public static double p_gamma(double a, double x, double loggamma_a) /* 本文参照 */\r
+ {\r
+ int k;\r
+ double result, term, previous;\r
+\r
+ if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);\r
+ if (x == 0) return 0;\r
+ result = term = Math.exp(a * Math.log(x) - x - loggamma_a) / a;\r
+ for (k = 1; k < 1000; k++) {\r
+ term *= x / (a + k);\r
+ previous = result; result += term;\r
+ if (result == previous) return result;\r
+ }\r
+ //throw new Exception("p_gamma(): the sequence is not convergent.");\r
+ return result;\r
+ }\r
+\r
+ public static double q_gamma(double a, double x, double loggamma_a) /* 本文参照 */\r
+ {\r
+ int k;\r
+ double result, w, temp, previous;\r
+ double la = 1, lb = 1 + x - a; /* Laguerreの多項式 */\r
+\r
+ if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);\r
+ w = Math.exp(a * Math.log(x) - x - loggamma_a);\r
+ result = w / lb;\r
+ for (k = 2; k < 1000; k++) {\r
+ temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;\r
+ la = lb; lb = temp;\r
+ w *= (k - 1 - a) / k;\r
+ temp = w / (la * lb);\r
+ previous = result; result += temp;\r
+ if (result == previous) return result;\r
+ }\r
+ //throw new Exception("q_gamma(): the sequence is not convergent.");\r
+ return result;\r
+ }\r
+\r
+ public static double p_chisq(double chisq, int df) /* カイ2乗分布の下側確率 */\r
+ {\r
+ return p_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));\r
+ }\r
+\r
+ public static double q_chisq(double chisq, int df) /* カイ2乗分布の上側確率 */\r
+ {\r
+ return q_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));\r
+ }\r
+\r
+ static readonly double LOG_PI = 1.14472988584940017; /* $\log_e \pi$ */\r
+\r
+ public static double erf(double x) /* Gaussの誤差関数 ${\rm erf}(x)$ */\r
+ {\r
+ if (x >= 0) return p_gamma(0.5, x * x, LOG_PI / 2);\r
+ else return - p_gamma(0.5, x * x, LOG_PI / 2);\r
+ }\r
+\r
+ public static double erfc(double x) /* $1 - {\rm erf}(x)$ */\r
+ {\r
+ if (x >= 0) return q_gamma(0.5, x * x, LOG_PI / 2);\r
+ else return 1 + p_gamma(0.5, x * x, LOG_PI / 2);\r
+ }\r
+\r
+ public static double p_normal(double x) /* 標準正規分布の下側確率 */\r
+ {\r
+ if (x >= 0) return\r
+ 0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));\r
+ else return\r
+ 0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);\r
+ }\r
+\r
+ public static double q_normal(double x) /* 標準正規分布の上側確率 */\r
+ {\r
+ if (x >= 0) return\r
+ 0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);\r
+ else return\r
+ 0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));\r
+ }\r
+ }\r
+ }\r
+}
\ No newline at end of file
OSDN Git Service
