simgear/structure/SGSmplstat.cxx

   1 // This may look like C code, but it is really -*- C++ -*-
   2 /*
   3 Copyright (C) 1988 Free Software Foundation
   4     written by Dirk Grunwald (grunwald@cs.uiuc.edu)
   5
   6 This file is part of the GNU C++ Library.  This library is free
   7 software; you can redistribute it and/or modify it under the terms of
   8 the GNU Library General Public License as published by the Free
   9 Software Foundation; either version 2 of the License, or (at your
  10 option) any later version.  This library is distributed in the hope
  11 that it will be useful, but WITHOUT ANY WARRANTY; without even the
  12 implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  13 PURPOSE.  See the GNU Library General Public License for more details.
  14 You should have received a copy of the GNU Library General Public
  15 License along with this library; if not, write to the Free Software
  16 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.
  17 */
  18
  19 #include <math.h>
  20
  21 #ifndef HUGE_VAL
  22 #ifdef HUGE
  23 #define HUGE_VAL HUGE
  24 #else
  25 #include <float.h>
  26 #define HUGE_VAL DBL_MAX
  27 #endif
  28 #endif
  29
  30
  31 #include <iostream>
  32 #include <fstream>
  33 #include <simgear/debug/logstream.hxx>
  34 #include "SGSmplstat.hxx"
  35
  36
  37 void SampleStatistic::error (const char *msg)
  38 {
  39     SG_LOG(SG_GENERAL, SG_ALERT,  msg);
  40 }
  41
  42 // t-distribution: given p-value and degrees of freedom, return t-value
  43 // adapted from Peizer & Pratt JASA, vol63, p1416
  44
  45 double tval (double p, int df)
  46 {
  47     double t;
  48     int positive = p >= 0.5;
  49     p = (positive) ? 1.0 - p : p;
  50     if (p <= 0.0 || df <= 0)
  51         t = HUGE_VAL;
  52     else if (p == 0.5)
  53         t = 0.0;
  54     else if (df == 1)
  55         t = 1.0 / tan ((p + p) * 1.57079633);
  56     else if (df == 2)
  57         t = sqrt (1.0 / ((p + p) * (1.0 - p)) - 2.0);
  58     else
  59     {
  60         double ddf = df;
  61         double a = sqrt (log (1.0 / (p * p)));
  62         double aa = a * a;
  63         a = a - ((2.515517 + (0.802853 * a) + (0.010328 * aa)) /
  64                 (1.0 + (1.432788 * a) + (0.189269 * aa) +
  65                         (0.001308 * aa * a)));
  66         t = ddf - 0.666666667 + 1.0 / (10.0 * ddf);
  67         t = sqrt (ddf * (exp (a * a * (ddf - 0.833333333) / (t * t)) - 1.0));
  68     }
  69     return (positive) ? t : -t;
  70 }
  71
  72 void SampleStatistic::reset ()
  73 {
  74     n = 0;
  75     x = x2 = 0.0;
  76     totalTime = 0.0;
  77     maxValue = -HUGE_VAL;
  78     minValue = HUGE_VAL;
  79 }
  80
  81 void SampleStatistic::operator += (double value)
  82 {
  83     n += 1;
  84     x += value;
  85     totalTime += value;
  86     cumulativeTime += value;
  87     x2 += (value * value);
  88
  89     if (minValue > value)
  90         minValue = value;
  91
  92     if (maxValue < value)
  93         maxValue = value;
  94 }
  95
  96 double SampleStatistic::mean () const
  97 {
  98     if (n > 0)
  99     {
 100       return (x / n);
 101     }
 102     else
 103     {
 104       return (0.0);
 105     }
 106 }
 107
 108 double SampleStatistic::var () const
 109 {
 110     if (n > 1)
 111     {
 112         return ((x2 - ((x * x) / n)) / (n - 1));
 113     }
 114     else
 115     {
 116         return (0.0);
 117     }
 118 }
 119
 120 double SampleStatistic::stdDev () const
 121 {
 122     if (n <= 0 || this->var () <= 0)
 123     {
 124       return (0);
 125     }
 126     else
 127     {
 128       return ((double) sqrt (var ()));
 129     }
 130 }
 131
 132 double SampleStatistic::confidence (int interval) const
 133 {
 134     int df = n - 1;
 135     if (df <= 0)
 136         return HUGE_VAL;
 137     double t = tval (double (100 + interval) * 0.005, df);
 138     if (t == HUGE_VAL)
 139         return t;
 140     else
 141         return (t * stdDev ()) / sqrt (double (n));
 142 }
 143
 144 double SampleStatistic::confidence (double p_value) const
 145 {
 146     int df = n - 1;
 147     if (df <= 0)
 148         return HUGE_VAL;
 149     double t = tval ((1.0 + p_value) * 0.5, df);
 150     if (t == HUGE_VAL)
 151         return t;
 152     else
 153         return (t * stdDev ()) / sqrt (double (n));
 154 }