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 #ifdef __GNUG__
  19 #pragma implementation
  20 #endif
  21
  22 #include <math.h>
  23
  24 #ifndef HUGE_VAL
  25 #ifdef HUGE
  26 #define HUGE_VAL HUGE
  27 #else
  28 #include <float.h>
  29 #define HUGE_VAL DBL_MAX
  30 #endif
  31 #endif
  32
  33
  34 #include <iostream>
  35 #include <fstream>
  36 #include <simgear/debug/logstream.hxx>
  37 #include "SGSmplstat.hxx"
  38
  39
  40 void SampleStatistic::error (const char *msg)
  41 {
  42   SG_LOG(SG_GENERAL, SG_ALERT,  msg);
  43 }
  44
  45 // t-distribution: given p-value and degrees of freedom, return t-value
  46 // adapted from Peizer & Pratt JASA, vol63, p1416
  47
  48 double tval (double p, int df)
  49 {
  50   double t;
  51   int positive = p >= 0.5;
  52   p = (positive) ? 1.0 - p : p;
  53   if (p <= 0.0 || df <= 0)
  54     t = HUGE_VAL;
  55   else if (p == 0.5)
  56     t = 0.0;
  57   else if (df == 1)
  58     t = 1.0 / tan ((p + p) * 1.57079633);
  59   else if (df == 2)
  60     t = sqrt (1.0 / ((p + p) * (1.0 - p)) - 2.0);
  61   else
  62     {
  63       double ddf = df;
  64       double a = sqrt (log (1.0 / (p * p)));
  65       double aa = a * a;
  66       a = a - ((2.515517 + (0.802853 * a) + (0.010328 * aa)) /
  67                (1.0 + (1.432788 * a) + (0.189269 * aa) +
  68                 (0.001308 * aa * a)));
  69       t = ddf - 0.666666667 + 1.0 / (10.0 * ddf);
  70       t = sqrt (ddf * (exp (a * a * (ddf - 0.833333333) / (t * t)) - 1.0));
  71     }
  72   return (positive) ? t : -t;
  73 }
  74
  75 void SampleStatistic::reset ()
  76 {
  77   n = 0;
  78   x = x2 = 0.0;
  79   maxValue = -HUGE_VAL;
  80   minValue = HUGE_VAL;
  81 }
  82
  83 void SampleStatistic::operator += (double value)
  84 {
  85   n += 1;
  86   x += value;
  87   x2 += (value * value);
  88   if (minValue > value)
  89     minValue = value;
  90   if (maxValue < value)
  91     maxValue = value;
  92 }
  93
  94 double SampleStatistic::mean () const
  95 {
  96   if (n > 0)
  97     {
  98       return (x / n);
  99     }
 100   else
 101     {
 102       return (0.0);
 103     }
 104 }
 105
 106 double SampleStatistic::var () const
 107 {
 108   if (n > 1)
 109     {
 110       return ((x2 - ((x * x) / n)) / (n - 1));
 111     }
 112   else
 113     {
 114       return (0.0);
 115     }
 116 }
 117
 118 double SampleStatistic::stdDev () const
 119 {
 120   if (n <= 0 || this->var () <= 0)
 121     {
 122       return (0);
 123     }
 124   else
 125     {
 126       return ((double) sqrt (var ()));
 127     }
 128 }
 129
 130 double SampleStatistic::confidence (int interval) const
 131 {
 132   int df = n - 1;
 133   if (df <= 0)
 134     return HUGE_VAL;
 135   double t = tval (double (100 + interval) * 0.005, df);
 136   if (t == HUGE_VAL)
 137     return t;
 138   else
 139     return (t * stdDev ()) / sqrt (double (n));
 140 }
 141
 142 double SampleStatistic::confidence (double p_value) const
 143 {
 144   int df = n - 1;
 145   if (df <= 0)
 146     return HUGE_VAL;
 147   double t = tval ((1.0 + p_value) * 0.5, df);
 148   if (t == HUGE_VAL)
 149     return t;
 150   else
 151     return (t * stdDev ()) / sqrt (double (n));
 152 }