1 /* 2 * cRunningStats.h 3 * Avida 4 * 5 * Created by David on 10/21/09. 6 * Copyright 2009-2011 Michigan State University. All rights reserved. 7 * 8 * 9 * This file is part of Avida. 10 * 11 * Avida is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License 12 * as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. 13 * 14 * Avida is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. 16 * 17 * You should have received a copy of the GNU Lesser General Public License along with Avida. 18 * If not, see <http://www.gnu.org/licenses/>. 19 * 20 */ 21 22 #ifndef cRunningStats_h 23 #define cRunningStats_h 24 25 #include <cmath> 26 27 28 class cRunningStats 29 { 30 private: 31 double m_n; // count 32 double m_m1; // mean 33 double m_m2; // second moment 34 double m_m3; // third moment 35 double m_m4; // fourth moment 36 37 public: cRunningStats()38 inline cRunningStats() : m_n(0.0), m_m1(0.0), m_m2(0.0), m_m3(0.0), m_m4(0.0) { ; } 39 Clear()40 inline void Clear() { m_n = 0.0; m_m1 = 0.0; m_m2 = 0.0; m_m3 = 0.0; m_m4 = 0.0; } 41 42 inline void Push(double x); 43 N()44 inline double N() const { return m_n; } Mean()45 inline double Mean() const { return m_m1; } StdDeviation()46 inline double StdDeviation() const { return sqrt(Variance()); } StdError()47 inline double StdError() const { return (m_n > 1.0) ? sqrt(Variance() / m_n) : 0.0; } Variance()48 inline double Variance() const { return (m_n > 1.0) ? (m_m2 / (m_n - 1.0)) : 0.0; } Skewness()49 inline double Skewness() const { return sqrt(m_n) * m_m3 / pow(m_m2, 1.5); } Kurtosis()50 inline double Kurtosis() const { return m_n * m_m4 / (m_m2 * m_m2); } 51 }; 52 53 Push(double x)54inline void cRunningStats::Push(double x) 55 { 56 m_n++; 57 double d = (x - m_m1); 58 double d_n = d / m_n; 59 double d_n2 = d_n * d_n; 60 61 m_m4 += d * d_n2 * d_n * ((m_n - 1) * ((m_n * m_n) - 3 * m_n + 3)) + 6 * d_n2 * m_m2 - 4 * d_n * m_m3; 62 m_m3 += d * d_n2 * ((m_n - 1) * (m_n - 2)) - 3 * d_n * m_m2; 63 m_m2 += d * d_n * (m_n - 1); 64 m_m1 += d_n; 65 } 66 67 #endif 68