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 const 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 T max<T>(T val1, T val2) where T : IComparable\r
  16                 {\r
  17                         return val1.CompareTo(val2) > 0 ? val1 : val2;\r
  18                 }\r
  19                 public static T min<T>(T val1, T val2) where T : IComparable\r
  20                 {\r
  21                         return val1.CompareTo(val2) < 0 ? 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 floor(double x)\r
  45                 {\r
  46                         return System.Math.Floor(x);\r
  47                 }\r
  48                 public static double ceil(double x)\r
  49                 {\r
  50                         return System.Math.Ceiling(x);\r
  51                 }\r
  52                 public static double sin(double x)\r
  53                 {\r
  54                         return System.Math.Sin(x);\r
  55                 }\r
  56                 public static double cos(double x)\r
  57                 {\r
  58                         return System.Math.Cos(x);\r
  59                 }\r
  60                 public static double tan(double x)\r
  61                 {\r
  62                         return System.Math.Tan(x);\r
  63                 }\r
  64                 public static double asin(double x)\r
  65                 {\r
  66                         return System.Math.Asin(x);\r
  67                 }\r
  68                 public static double acos(double x)\r
  69                 {\r
  70                         return System.Math.Acos(x);\r
  71                 }\r
  72                 public static double atan(double x)\r
  73                 {\r
  74                         return System.Math.Atan(x);\r
  75                 }\r
  76                 public static double atan(double y, double x)\r
  77                 {\r
  78                         return System.Math.Atan2(y, x);\r
  79                 }\r
  80                 public static double atan2(double y, double x)\r
  81                 {\r
  82                         return System.Math.Atan2(y, x);\r
  83                 }\r
  84                 public static double sqrt(double x)\r
  85                 {\r
  86                         return System.Math.Sqrt(x);\r
  87                 }\r
  88                 public static double pow(double x, double y)\r
  89                 {\r
  90                         return System.Math.Pow(x,y);\r
  91                 }\r
  92                 public static double exp(double x)\r
  93                 {\r
  94                         return System.Math.Exp(x);\r
  95                 }\r
  96                 public static double log(double x)\r
  97                 {\r
  98                         return System.Math.Log(x);\r
  99                 }\r
 100                 public static double log2(double val)\r
 101                 {\r
 102                         return log(val) * LOG2E;\r
 103                 }\r
 104                 public static double round(double val)\r
 105                 {\r
 106                         return System.Math.Round(val);\r
 107                 }\r
 108 \r
 109                 public static double radius(double x, double y)\r
 110                 {\r
 111                         return System.Math.Sqrt(x * x + y * y);\r
 112                 }\r
 113 \r
 114                 public static double random()\r
 115                 {\r
 116                         return (random_generator.NextDouble());\r
 117                 }\r
 118                 public static int random(int x)\r
 119                 {\r
 120                         return (int)((random_generator.NextDouble()) * x);\r
 121                 }\r
 122                 public static double random(double x)\r
 123                 {\r
 124                         return (random_generator.NextDouble()) * x;\r
 125                 }\r
 126                 public static double random(double x, double y)\r
 127                 {\r
 128                         return (random_generator.NextDouble()) * (y-x) + x;\r
 129                 }\r
 130 \r
 131 \r
 132                 public static double gaussian(double x, double sigma)\r
 133                 {\r
 134                         return exp(- (x*x) / (2*sigma*sigma));\r
 135                 }\r
 136 \r
 137                 public static double normalDistibution(double x, double mu, double sigma)\r
 138                 {\r
 139                         return exp( -( (x-mu)*(x-mu) / (2*sigma*sigma) ) ) / sqrt(2*PI*sigma*sigma);\r
 140                 }\r
 141 \r
 142                 public static double cumulativeNormalDistibution(double x, double mu, double sigma)\r
 143                 {\r
 144                         return .5 + .5*Internal.GammaFunction.erf( (x-mu)/(sigma*sqrt(2.0) ) );\r
 145                 }\r
 146 \r
 147         }\r
 148 \r
 149         namespace Internal\r
 150         {\r
 151                 public static class GammaFunction\r
 152                 {\r
 153                         /************ loggamma(x) -- gamma.c より再掲 *************/\r
 154 \r
 155                         static readonly double PI      = 3.14159265358979324;  /* $\pi$ */\r
 156                         static readonly double LOG_2PI = 1.83787706640934548;  /* $\log 2\pi$ */\r
 157                         static readonly double N       = 8;\r
 158 \r
 159                         static readonly double B0  = 1            ;     /* 以下はBernoulli数 */\r
 160                         static readonly double B1  = (-1.0 / 2.0);\r
 161                         static readonly double B2  = ( 1.0 / 6.0);\r
 162                         static readonly double B4  = (-1.0 / 30.0);\r
 163                         static readonly double B6  = ( 1.0 / 42.0);\r
 164                         static readonly double B8  = (-1.0 / 30.0);\r
 165                         static readonly double B10 = ( 5.0 / 66.0);\r
 166                         static readonly double B12 = (-691.0 / 2730.0);\r
 167                         static readonly double B14 = ( 7.0 / 6.0);\r
 168                         static readonly double B16 = (-3617.0 / 510.0);\r
 169 \r
 170                         public static double loggamma(double x)  /* ガンマ関数の対数 */\r
 171                         {\r
 172                                 double v, w;\r
 173 \r
 174                                 v = 1;\r
 175                                 while (x < N) {  v *= x;  x++;  }\r
 176                                 w = 1 / (x * x);\r
 177                                 return ((((((((B16 / (16 * 15))  * w + (B14 / (14 * 13))) * w\r
 178                                                         + (B12 / (12 * 11))) * w + (B10 / (10 *  9))) * w\r
 179                                                         + (B8  / ( 8 *  7))) * w + (B6  / ( 6 *  5))) * w\r
 180                                                         + (B4  / ( 4 *  3))) * w + (B2  / ( 2 *  1))) / x\r
 181                                                         + 0.5 * LOG_2PI - Math.log(v) - x + (x - 0.5) * Math.log(x);\r
 182                         }\r
 183 \r
 184                         public static double p_gamma(double a, double x, double loggamma_a)  /* 本文参照 */\r
 185                         {\r
 186                                 int k;\r
 187                                 double result, term, previous;\r
 188 \r
 189                                 if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);\r
 190                                 if (x == 0)     return 0;\r
 191                                 result = term = Math.exp(a * Math.log(x) - x - loggamma_a) / a;\r
 192                                 for (k = 1; k < 1000; k++) {\r
 193                                         term *= x / (a + k);\r
 194                                         previous = result;  result += term;\r
 195                                         if (result == previous) return result;\r
 196                                 }\r
 197                                 //throw new Exception("p_gamma(): the sequence is not convergent.");\r
 198                                 return result;\r
 199                         }\r
 200 \r
 201                         public static double q_gamma(double a, double x, double loggamma_a)  /* 本文参照 */\r
 202                         {\r
 203                                 int k;\r
 204                                 double result, w, temp, previous;\r
 205                                 double la = 1, lb = 1 + x - a;  /* Laguerreの多項式 */\r
 206 \r
 207                                 if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);\r
 208                                 w = Math.exp(a * Math.log(x) - x - loggamma_a);\r
 209                                 result = w / lb;\r
 210                                 for (k = 2; k < 1000; k++) {\r
 211                                         temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;\r
 212                                         la = lb;  lb = temp;\r
 213                                         w *= (k - 1 - a) / k;\r
 214                                         temp = w / (la * lb);\r
 215                                         previous = result;  result += temp;\r
 216                                         if (result == previous) return result;\r
 217                                 }\r
 218                                 //throw new Exception("q_gamma(): the sequence is not convergent.");\r
 219                                 return result;\r
 220                         }\r
 221 \r
 222                         public static double p_chisq(double chisq, int df)  /* カイ2乗分布の下側確率 */\r
 223                         {\r
 224                                 return p_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));\r
 225                         }\r
 226 \r
 227                         public static double q_chisq(double chisq, int df)  /* カイ2乗分布の上側確率 */\r
 228                         {\r
 229                                 return q_gamma(0.5 * df, 0.5 * chisq, loggamma(0.5 * df));\r
 230                         }\r
 231 \r
 232                         static readonly double LOG_PI = 1.14472988584940017;  /* $\log_e \pi$ */\r
 233 \r
 234                         public static double erf(double x)  /* Gaussの誤差関数 ${\rm erf}(x)$ */\r
 235                         {\r
 236                                 if (x >= 0) return   p_gamma(0.5, x * x, LOG_PI / 2);\r
 237                                 else        return - p_gamma(0.5, x * x, LOG_PI / 2);\r
 238                         }\r
 239 \r
 240                         public static double erfc(double x)  /* $1 - {\rm erf}(x)$ */\r
 241                         {\r
 242                                 if (x >= 0) return  q_gamma(0.5, x * x, LOG_PI / 2);\r
 243                                 else        return  1 + p_gamma(0.5, x * x, LOG_PI / 2);\r
 244                         }\r
 245 \r
 246                         public static double p_normal(double x)  /* 標準正規分布の下側確率 */\r
 247                         {\r
 248                                 if (x >= 0) return\r
 249                                         0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));\r
 250                                 else return\r
 251                                         0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);\r
 252                         }\r
 253 \r
 254                         public static double q_normal(double x)  /* 標準正規分布の上側確率 */\r
 255                         {\r
 256                                 if (x >= 0) return\r
 257                                         0.5 * q_gamma(0.5, 0.5 * x * x, LOG_PI / 2);\r
 258                                 else return\r
 259                                         0.5 * (1 + p_gamma(0.5, 0.5 * x * x, LOG_PI / 2));\r
 260                         }\r
 261                 }\r
 262         }\r
 263 }