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