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