dev4/psychlops/core/math/util.cs

   1 using System;\r
   2 \r
   3 namespace Psychlops\r
   4 {\r
   5 \r
   6         public static class Math\r
   7         {\r
   8                 public static readonly double PI = 3.14159265, E = 2.718281828459045, LOG2E = 1.44269504088896340736;\r
   9                 public static Random random_generator;\r
  10                 static Math()\r
  11                 {\r
  12                         random_generator = new Random();\r
  13                 }\r
  14 \r
  15                 public static double max(double val1, double val2)\r
  16                 {\r
  17                         return val1 > val2 ? val1 : val2;\r
  18                 }\r
  19                 public static double min(double val1, double val2)\r
  20                 {\r
  21                         return val1 < val2 ? val1 : val2;\r
  22                 }\r
  23                 public static void shuffle<X>(X[] array, int n)\r
  24                 {\r
  25                         int a;\r
  26                         X tmp;\r
  27                         for(int i = 1; i < n; i++){\r
  28                                 a = random(i + 1);\r
  29                                 tmp = array[i];\r
  30                                 array[i] = array[a];\r
  31                                 array[a] = tmp;\r
  32                         }\r
  33                 }\r
  34 \r
  35 \r
  36                 public static double mod(double lhs, double rhs)\r
  37                 {\r
  38                         return lhs - System.Math.Floor(lhs/rhs)*rhs;\r
  39                 }\r
  40                 public static double abs(double x)\r
  41                 {\r
  42                         return System.Math.Abs(x);\r
  43                 }\r
  44                 public static double sin(double x)\r
  45                 {\r
  46                         return System.Math.Sin(x);\r
  47                 }\r
  48                 public static double cos(double x)\r
  49                 {\r
  50                         return System.Math.Cos(x);\r
  51                 }\r
  52                 public static double tan(double x)\r
  53                 {\r
  54                         return System.Math.Tan(x);\r
  55                 }\r
  56                 public static double sqrt(double x)\r
  57                 {\r
  58                         return System.Math.Sqrt(x);\r
  59                 }\r
  60                 public static double exp(double x)\r
  61                 {\r
  62                         return System.Math.Exp(x);\r
  63                 }\r
  64                 public static double log(double x)\r
  65                 {\r
  66                         return System.Math.Log(x);\r
  67                 }\r
  68                 public static double log2(double val)\r
  69                 {\r
  70                         return log(val) * LOG2E;\r
  71                 }\r
  72                 /*public static int round(double val)\r
  73                 {\r
  74                         double integer_part, particle = modf(val, &integer_part);\r
  75                         return ((particle < 0.5 | (particle == 0.5 && (int)integer_part % 2 == 0)) ? (int)integer_part : (int)integer_part + 1);\r
  76                 }*/\r
  77 \r
  78                 public static double radius(double x, double y)\r
  79                 {\r
  80                         return System.Math.Sqrt(x * x + y * y);\r
  81                 }\r
  82 \r
  83                 public static double random()\r
  84                 {\r
  85                         return (random_generator.NextDouble());\r
  86                 }\r
  87                 public static int random(int x)\r
  88                 {\r
  89                         return (int)((random_generator.NextDouble()) * x);\r
  90                 }\r
  91                 public static double random(double x)\r
  92                 {\r
  93                         return (random_generator.NextDouble()) * x;\r
  94                 }\r
  95                 public static double random(double x, double y)\r
  96                 {\r
  97                         return (random_generator.NextDouble()) * (y-x) + x;\r
  98                 }\r
  99 \r
 100 \r
 101                 public static double gaussian(double x, double sigma)\r
 102                 {\r
 103                         return exp(- (x*x) / (2*sigma*sigma));\r
 104                 }\r
 105 \r
 106                 public static double normalDistibution(double x, double mu, double sigma)\r
 107                 {\r
 108                         return exp( -( (x-mu)*(x-mu) / (2*sigma*sigma) ) ) / sqrt(2*PI*sigma*sigma);\r
 109                 }\r
 110 \r
 111                 public static double cumulativeNormalDistibution(double x, double mu, double sigma)\r
 112                 {\r
 113                         return .5 + .5*Internal.GammaFunction.erf( (x-mu)/(sigma*sqrt(2.0) ) );\r
 114                 }\r
 115 \r
 116         }\r
 117 \r
 118         namespace Internal\r
 119         {\r
 120                 public static class GammaFunction\r
 121                 {\r
 122                         /************ loggamma(x) -- gamma.c より再掲 *************/\r
 123 \r
 124                         static readonly double PI      = 3.14159265358979324;  /* $\pi$ */\r
 125                         static readonly double LOG_2PI = 1.83787706640934548;  /* $\log 2\pi$ */\r
 126                         static readonly double N       = 8;\r
 127 \r
 128                         static readonly double B0  = 1            ;     /* 以下はBernoulli数 */\r
 129                         static readonly double B1  = (-1.0 / 2.0);\r
 130                         static readonly double B2  = ( 1.0 / 6.0);\r
 131                         static readonly double B4  = (-1.0 / 30.0);\r
 132                         static readonly double B6  = ( 1.0 / 42.0);\r
 133                         static readonly double B8  = (-1.0 / 30.0);\r
 134                         static readonly double B10 = ( 5.0 / 66.0);\r
 135                         static readonly double B12 = (-691.0 / 2730.0);\r
 136                         static readonly double B14 = ( 7.0 / 6.0);\r
 137                         static readonly double B16 = (-3617.0 / 510.0);\r
 138 \r
 139                         public static double loggamma(double x)  /* ガンマ関数の対数 */\r
 140                         {\r
 141                                 double v, w;\r
 142 \r
 143                                 v = 1;\r
 144                                 while (x < N) {  v *= x;  x++;  }\r
 145                                 w = 1 / (x * x);\r
 146                                 return ((((((((B16 / (16 * 15))  * w + (B14 / (14 * 13))) * w\r
 147                                                         + (B12 / (12 * 11))) * w + (B10 / (10 *  9))) * w\r
 148                                                         + (B8  / ( 8 *  7))) * w + (B6  / ( 6 *  5))) * w\r
 149                                                         + (B4  / ( 4 *  3))) * w + (B2  / ( 2 *  1))) / x\r
 150                                                         + 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x);\r
 151                         }\r
 152 \r
 153                         public static double p_gamma(double a, double x, double loggamma_a)  /* 本文参照 */\r
 154                         {\r
 155                                 int k;\r
 156                                 double result, term, previous;\r
 157 \r
 158                                 if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);\r
 159                                 if (x == 0)     return 0;\r
 160                                 result = term = Math.exp(a * Math.log(x) - x - loggamma_a) / a;\r
 161                                 for (k = 1; k < 1000; k++) {\r
 162                                         term *= x / (a + k);\r
 163                                         previous = result;  result += term;\r
 164                                         if (result == previous) return result;\r
 165                                 }\r
 166                                 //throw new Exception("p_gamma(): the sequence is not convergent.");\r
 167                                 return result;\r
 168                         }\r
 169 \r
 170                         public static double q_gamma(double a, double x, double loggamma_a)  /* 本文参照 */\r
 171                         {\r
 172                                 int k;\r
 173                                 double result, w, temp, previous;\r
 174                                 double la = 1, lb = 1 + x - a;  /* Laguerreの多項式 */\r
 175 \r
 176                                 if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);\r
 177                                 w = Math.exp(a * Math.log(x) - x - loggamma_a);\r
 178                                 result = w / lb;\r
 179                                 for (k = 2; k < 1000; k++) {\r
 180                                         temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;\r
 181                                         la = lb;  lb = temp;\r
 182                                         w *= (k - 1 - a) / k;\r
 183                                         temp = w / (la * lb);\r
 184                                         previous = result;  result += temp;\r
 185                                         if (result == previous) return result;\r
 186                                 }\r
 187                                 //throw new Exception("q_gamma(): the sequence is not convergent.");\r
 188                                 return result;\r
 189                         }\r
 190 \r
 191                         public static double p_chisq(double chisq, int df)  /* カイ2乗分布の下側確率 */\r
 192                         {\r
 193                                 return p_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));\r
 194                         }\r
 195 \r
 196                         public static double q_chisq(double chisq, int df)  /* カイ2乗分布の上側確率 */\r
 197                         {\r
 198                                 return q_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));\r
 199                         }\r
 200 \r
 201                         static readonly double LOG_PI = 1.14472988584940017;  /* $\log_e \pi$ */\r
 202 \r
 203                         public static double erf(double x)  /* Gaussの誤差関数 ${\rm erf}(x)$ */\r
 204                         {\r
 205                                 if (x >= 0) return   p_gamma(0.5, x * x, LOG_PI / 2);\r
 206                                 else        return - p_gamma(0.5, x * x, LOG_PI / 2);\r
 207                         }\r
 208 \r
 209                         public static double erfc(double x)  /* $1 - {\rm erf}(x)$ */\r
 210                         {\r
 211                                 if (x >= 0) return  q_gamma(0.5, x * x, LOG_PI / 2);\r
 212                                 else        return  1 + p_gamma(0.5, x * x, LOG_PI / 2);\r
 213                         }\r
 214 \r
 215                         public static double p_normal(double x)  /* 標準正規分布の下側確率 */\r
 216                         {\r
 217                                 if (x >= 0) return\r
 218                                         0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));\r
 219                                 else return\r
 220                                         0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);\r
 221                         }\r
 222 \r
 223                         public static double q_normal(double x)  /* 標準正規分布の上側確率 */\r
 224                         {\r
 225                                 if (x >= 0) return\r
 226                                         0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);\r
 227                                 else return\r
 228                                         0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));\r
 229                         }\r
 230                 }\r
 231         }\r
 232 }