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, 59 Temple Place - Suite 330, Boston, MA 02111-1307, 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   maxValue = -HUGE_VAL;
  77   minValue = HUGE_VAL;
  78 }
  79
  80 void SampleStatistic::operator += (double value)
  81 {
  82   n += 1;
  83   x += value;
  84   allTimeTotal += value;
  85   x2 += (value * value);
  86   if (minValue > value)
  87     minValue = value;
  88   if (maxValue < value)
  89     maxValue = value;
  90 }
  91
  92 double SampleStatistic::mean () const
  93 {
  94   if (n > 0)
  95     {
  96       return (x / n);
  97     }
  98   else
  99     {
 100       return (0.0);
 101     }
 102 }
 103
 104 double SampleStatistic::var () const
 105 {
 106   if (n > 1)
 107     {
 108       return ((x2 - ((x * x) / n)) / (n - 1));
 109     }
 110   else
 111     {
 112       return (0.0);
 113     }
 114 }
 115
 116 double SampleStatistic::stdDev () const
 117 {
 118   if (n <= 0 || this->var () <= 0)
 119     {
 120       return (0);
 121     }
 122   else
 123     {
 124       return ((double) sqrt (var ()));
 125     }
 126 }
 127
 128 double SampleStatistic::confidence (int interval) const
 129 {
 130   int df = n - 1;
 131   if (df <= 0)
 132     return HUGE_VAL;
 133   double t = tval (double (100 + interval) * 0.005, df);
 134   if (t == HUGE_VAL)
 135     return t;
 136   else
 137     return (t * stdDev ()) / sqrt (double (n));
 138 }
 139
 140 double SampleStatistic::confidence (double p_value) const
 141 {
 142   int df = n - 1;
 143   if (df <= 0)
 144     return HUGE_VAL;
 145   double t = tval ((1.0 + p_value) * 0.5, df);
 146   if (t == HUGE_VAL)
 147     return t;
 148   else
 149     return (t * stdDev ()) / sqrt (double (n));
 150 }