[Avida-SVN] r3503 - in development: . Avida.xcodeproj source/main source/tools
brysonda at myxo.css.msu.edu
brysonda at myxo.css.msu.edu
Wed Oct 21 13:46:16 PDT 2009
Author: brysonda
Date: 2009-10-21 16:46:16 -0400 (Wed, 21 Oct 2009)
New Revision: 3503
Added:
development/source/tools/cRunningStats.h
Removed:
development/source/tools/cDoubleSum.cc
Modified:
development/Avida.xcodeproj/project.pbxproj
development/CMakeLists.txt
development/source/main/cPopulation.cc
development/source/main/cStats.cc
development/source/main/cStats.h
development/source/tools/cDoubleSum.h
Log:
Implement a new cRunningStats class that calculates running statistics incrementally. This class uses numerically stable methods for performing its work. More importantly it reports correct values for Skewness and Kurtosis. cDoubleSum's implementation of those properties did not. I did fairly extensive searching to find the source of the original equations upon which cDoubleSum's implementation was based to no avail. As such I have removed the methods from cDoubleSum. In the instances that needed Skewness and Kurtosis I have replaced with cRunningStats. cDoubleSum still exists due to additional weighting functionality and its bi-directional capability.
NOTE THIS WILL ALTER CONSISTENCY.
Two stats file routines related to mutation rates output the skewness and kurtosis values. These values are now different. However, it is not so much a break, rather a consistency fix.
Test expected results updates to follow.
Modified: development/Avida.xcodeproj/project.pbxproj
===================================================================
--- development/Avida.xcodeproj/project.pbxproj 2009-10-16 19:20:14 UTC (rev 3502)
+++ development/Avida.xcodeproj/project.pbxproj 2009-10-21 20:46:16 UTC (rev 3503)
@@ -81,7 +81,6 @@
7023EC4F0C0A431B00362B9C /* cDefaultMessageDisplay.cc in Sources */ = {isa = PBXBuildFile; fileRef = 70B0885208F5FE5800FC65FE /* cDefaultMessageDisplay.cc */; };
7023EC500C0A431B00362B9C /* cDefaultRunDriver.cc in Sources */ = {isa = PBXBuildFile; fileRef = 701D930C094CAD6B008B845F /* cDefaultRunDriver.cc */; };
7023EC510C0A431B00362B9C /* cDeme.cc in Sources */ = {isa = PBXBuildFile; fileRef = 1097463D0AE9606E00929ED6 /* cDeme.cc */; };
- 7023EC520C0A431B00362B9C /* cDoubleSum.cc in Sources */ = {isa = PBXBuildFile; fileRef = 70B0885308F5FE5800FC65FE /* cDoubleSum.cc */; };
7023EC530C0A431B00362B9C /* cDriverManager.cc in Sources */ = {isa = PBXBuildFile; fileRef = 701D9382094CBA69008B845F /* cDriverManager.cc */; };
7023EC540C0A431B00362B9C /* cEnvironment.cc in Sources */ = {isa = PBXBuildFile; fileRef = 702D4EFC08DA5341007BA469 /* cEnvironment.cc */; };
7023EC550C0A431B00362B9C /* cEventList.cc in Sources */ = {isa = PBXBuildFile; fileRef = 708BF2FD0AB65DC700A923BF /* cEventList.cc */; };
@@ -445,6 +444,7 @@
700D9CB20F1C55EB002CC711 /* tRLockPtr.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = tRLockPtr.h; sourceTree = "<group>"; };
700E28CF0859FFD700CF158A /* tObjectFactory.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = tObjectFactory.h; sourceTree = "<group>"; };
700E2B83085DE50C00CF158A /* avida-viewer */ = {isa = PBXFileReference; explicitFileType = "compiled.mach-o.executable"; includeInIndex = 0; path = "avida-viewer"; sourceTree = BUILT_PRODUCTS_DIR; };
+ 70100BD7108F8F4F005999F0 /* cRunningStats.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = cRunningStats.h; sourceTree = "<group>"; };
7013845F09028B3E0087ED2E /* cAvidaConfig.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = cAvidaConfig.h; sourceTree = "<group>"; };
7013846009028B3E0087ED2E /* cAvidaConfig.cc */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = cAvidaConfig.cc; sourceTree = "<group>"; };
701384A10902A16F0087ED2E /* defs.h */ = {isa = PBXFileReference; fileEncoding = 30; lastKnownFileType = sourcecode.c.h; path = defs.h; sourceTree = "<group>"; };
@@ -798,7 +798,6 @@
70B0885008F5FE5800FC65FE /* cDataFileManager.cc */ = {isa = PBXFileReference; fileEncoding = 30; lastKnownFileType = sourcecode.cpp.cpp; path = cDataFileManager.cc; sourceTree = "<group>"; };
70B0885108F5FE5800FC65FE /* cDataManager_Base.cc */ = {isa = PBXFileReference; fileEncoding = 30; lastKnownFileType = sourcecode.cpp.cpp; path = cDataManager_Base.cc; sourceTree = "<group>"; };
70B0885208F5FE5800FC65FE /* cDefaultMessageDisplay.cc */ = {isa = PBXFileReference; fileEncoding = 30; lastKnownFileType = sourcecode.cpp.cpp; path = cDefaultMessageDisplay.cc; sourceTree = "<group>"; };
- 70B0885308F5FE5800FC65FE /* cDoubleSum.cc */ = {isa = PBXFileReference; fileEncoding = 30; lastKnownFileType = sourcecode.cpp.cpp; path = cDoubleSum.cc; sourceTree = "<group>"; };
70B0887D08F603C600FC65FE /* cFile.h */ = {isa = PBXFileReference; fileEncoding = 30; lastKnownFileType = sourcecode.c.h; path = cFile.h; sourceTree = "<group>"; };
70B0887F08F603C600FC65FE /* cFixedBlock.h */ = {isa = PBXFileReference; fileEncoding = 30; lastKnownFileType = sourcecode.c.h; path = cFixedBlock.h; sourceTree = "<group>"; };
70B0888008F603C600FC65FE /* cFixedCoords.h */ = {isa = PBXFileReference; fileEncoding = 30; lastKnownFileType = sourcecode.c.h; path = cFixedCoords.h; sourceTree = "<group>"; };
@@ -1804,7 +1803,6 @@
70B0885008F5FE5800FC65FE /* cDataFileManager.cc */,
70B0885108F5FE5800FC65FE /* cDataManager_Base.cc */,
70B0885208F5FE5800FC65FE /* cDefaultMessageDisplay.cc */,
- 70B0885308F5FE5800FC65FE /* cDoubleSum.cc */,
70B0884908F5FE4500FC65FE /* cDataFile.h */,
70B0884A08F5FE4500FC65FE /* cDataFileManager.h */,
70B0884B08F5FE4500FC65FE /* cDataManager_Base.h */,
@@ -1845,6 +1843,7 @@
7069470C0EE8FA2700AD67C0 /* tArrayUtils.h */,
700D9C170F1A7436002CC711 /* tArrayMap.h */,
700D9C1A0F1A77EC002CC711 /* tKVPair.h */,
+ 70100BD7108F8F4F005999F0 /* cRunningStats.h */,
);
path = tools;
sourceTree = "<group>";
@@ -2196,7 +2195,6 @@
7023EC4F0C0A431B00362B9C /* cDefaultMessageDisplay.cc in Sources */,
7023EC500C0A431B00362B9C /* cDefaultRunDriver.cc in Sources */,
7023EC510C0A431B00362B9C /* cDeme.cc in Sources */,
- 7023EC520C0A431B00362B9C /* cDoubleSum.cc in Sources */,
7023EC530C0A431B00362B9C /* cDriverManager.cc in Sources */,
7023EC540C0A431B00362B9C /* cEnvironment.cc in Sources */,
7023EC550C0A431B00362B9C /* cEventList.cc in Sources */,
Modified: development/CMakeLists.txt
===================================================================
--- development/CMakeLists.txt 2009-10-16 19:20:14 UTC (rev 3502)
+++ development/CMakeLists.txt 2009-10-21 20:46:16 UTC (rev 3503)
@@ -303,7 +303,6 @@
${TOOLS_DIR}/cDataManager_Base.cc
${TOOLS_DIR}/cDefaultMessageDisplay.cc
${TOOLS_DIR}/cDemeProbSchedule.cc
- ${TOOLS_DIR}/cDoubleSum.cc
${TOOLS_DIR}/cFile.cc
${TOOLS_DIR}/cHelpAlias.cc
${TOOLS_DIR}/cHelpManager.cc
Modified: development/source/main/cPopulation.cc
===================================================================
--- development/source/main/cPopulation.cc 2009-10-16 19:20:14 UTC (rev 3502)
+++ development/source/main/cPopulation.cc 2009-10-21 20:46:16 UTC (rev 3503)
@@ -4149,10 +4149,10 @@
stats.SumGeneration().Add(phenotype.GetGeneration());
stats.SumNeutralMetric().Add(phenotype.GetNeutralMetric());
stats.SumLineageLabel().Add(organism->GetLineageLabel());
- stats.SumCopyMutRate().Add(organism->MutationRates().GetCopyMutProb());
- stats.SumLogCopyMutRate().Add(log(organism->MutationRates().GetCopyMutProb()));
- stats.SumDivMutRate().Add(organism->MutationRates().GetDivMutProb() / organism->GetPhenotype().GetDivType());
- stats.SumLogDivMutRate().Add(log(organism->MutationRates().GetDivMutProb() /organism->GetPhenotype().GetDivType()));
+ stats.SumCopyMutRate().Push(organism->MutationRates().GetCopyMutProb());
+ stats.SumLogCopyMutRate().Push(log(organism->MutationRates().GetCopyMutProb()));
+ stats.SumDivMutRate().Push(organism->MutationRates().GetDivMutProb() / organism->GetPhenotype().GetDivType());
+ stats.SumLogDivMutRate().Push(log(organism->MutationRates().GetDivMutProb() /organism->GetPhenotype().GetDivType()));
stats.SumCopySize().Add(phenotype.GetCopiedSize());
stats.SumExeSize().Add(phenotype.GetExecutedSize());
stats.SetGenoMapElement(i, organism->GetGenotype()->GetID());
Modified: development/source/main/cStats.cc
===================================================================
--- development/source/main/cStats.cc 2009-10-16 19:20:14 UTC (rev 3502)
+++ development/source/main/cStats.cc 2009-10-21 20:46:16 UTC (rev 3503)
@@ -1250,17 +1250,17 @@
df.WriteTimeStamp();
df.Write(m_update, "Update");
- df.Write(sum_copy_mut_rate.Ave(), "Average copy mutation rate");
- df.Write(sum_copy_mut_rate.Var(), "Variance in copy mutation rate");
+ df.Write(sum_copy_mut_rate.Mean(), "Average copy mutation rate");
+ df.Write(sum_copy_mut_rate.Variance(), "Variance in copy mutation rate");
df.Write(sum_copy_mut_rate.StdDeviation(), "Standard Deviation in copy mutation rate");
- df.Write(sum_copy_mut_rate.Skw(), "Skew in copy mutation rate");
- df.Write(sum_copy_mut_rate.Kur(), "Kurtosis in copy mutation rate");
+ df.Write(sum_copy_mut_rate.Skewness(), "Skew in copy mutation rate");
+ df.Write(sum_copy_mut_rate.Kurtosis(), "Kurtosis in copy mutation rate");
- df.Write(sum_log_copy_mut_rate.Ave(), "Average log(copy mutation rate)");
- df.Write(sum_log_copy_mut_rate.Var(), "Variance in log(copy mutation rate)");
+ df.Write(sum_log_copy_mut_rate.Mean(), "Average log(copy mutation rate)");
+ df.Write(sum_log_copy_mut_rate.Variance(), "Variance in log(copy mutation rate)");
df.Write(sum_log_copy_mut_rate.StdDeviation(), "Standard Deviation in log(copy mutation rate)");
- df.Write(sum_log_copy_mut_rate.Skw(), "Skew in log(copy mutation rate)");
- df.Write(sum_log_copy_mut_rate.Kur(), "Kurtosis in log(copy mutation rate)");
+ df.Write(sum_log_copy_mut_rate.Skewness(), "Skew in log(copy mutation rate)");
+ df.Write(sum_log_copy_mut_rate.Kurtosis(), "Kurtosis in log(copy mutation rate)");
df.Endl();
}
@@ -1274,17 +1274,17 @@
df.WriteTimeStamp();
df.Write(m_update, "Update");
- df.Write(sum_div_mut_rate.Ave(), "Average divide mutation rate");
- df.Write(sum_div_mut_rate.Var(), "Variance in divide mutation rate");
+ df.Write(sum_div_mut_rate.Mean(), "Average divide mutation rate");
+ df.Write(sum_div_mut_rate.Variance(), "Variance in divide mutation rate");
df.Write(sum_div_mut_rate.StdDeviation(), "Standard Deviation in divide mutation rate");
- df.Write(sum_div_mut_rate.Skw(), "Skew in divide mutation rate");
- df.Write(sum_div_mut_rate.Kur(), "Kurtosis in divide mutation rate");
+ df.Write(sum_div_mut_rate.Skewness(), "Skew in divide mutation rate");
+ df.Write(sum_div_mut_rate.Kurtosis(), "Kurtosis in divide mutation rate");
- df.Write(sum_log_div_mut_rate.Ave(), "Average log(divide mutation rate)");
- df.Write(sum_log_div_mut_rate.Var(), "Variance in log(divide mutation rate)");
+ df.Write(sum_log_div_mut_rate.Mean(), "Average log(divide mutation rate)");
+ df.Write(sum_log_div_mut_rate.Variance(), "Variance in log(divide mutation rate)");
df.Write(sum_log_div_mut_rate.StdDeviation(), "Standard Deviation in log(divide mutation rate)");
- df.Write(sum_log_div_mut_rate.Skw(), "Skew in log(divide mutation rate)");
- df.Write(sum_log_div_mut_rate.Kur(), "Kurtosis in log(divide mutation rate)");
+ df.Write(sum_log_div_mut_rate.Skewness(), "Skew in log(divide mutation rate)");
+ df.Write(sum_log_div_mut_rate.Kurtosis(), "Kurtosis in log(divide mutation rate)");
df.Endl();
}
Modified: development/source/main/cStats.h
===================================================================
--- development/source/main/cStats.h 2009-10-16 19:20:14 UTC (rev 3502)
+++ development/source/main/cStats.h 2009-10-21 20:46:16 UTC (rev 3503)
@@ -51,6 +51,9 @@
#ifndef cRunningAverage_h
#include "cRunningAverage.h"
#endif
+#ifndef cRunningStats_h
+#include "cRunningStats.h"
+#endif
#ifndef tArray_h
#include "tArray.h"
#endif
@@ -116,11 +119,11 @@
cDoubleSum sum_neutral_metric;
cDoubleSum sum_lineage_label;
- cDoubleSum sum_copy_mut_rate;
- cDoubleSum sum_log_copy_mut_rate;
+ cRunningStats sum_copy_mut_rate;
+ cRunningStats sum_log_copy_mut_rate;
- cDoubleSum sum_div_mut_rate;
- cDoubleSum sum_log_div_mut_rate;
+ cRunningStats sum_div_mut_rate;
+ cRunningStats sum_log_div_mut_rate;
//// By Genotype Sums //// (Cleared and resummed by population each update)
@@ -463,10 +466,10 @@
cDoubleSum& SumNeutralMetric() { return sum_neutral_metric; }
cDoubleSum& SumLineageLabel() { return sum_lineage_label; }
- cDoubleSum& SumCopyMutRate() { return sum_copy_mut_rate; }
- cDoubleSum& SumLogCopyMutRate() { return sum_log_copy_mut_rate; }
- cDoubleSum& SumDivMutRate() { return sum_div_mut_rate; }
- cDoubleSum& SumLogDivMutRate() { return sum_log_div_mut_rate; }
+ cRunningStats& SumCopyMutRate() { return sum_copy_mut_rate; }
+ cRunningStats& SumLogCopyMutRate() { return sum_log_copy_mut_rate; }
+ cRunningStats& SumDivMutRate() { return sum_div_mut_rate; }
+ cRunningStats& SumLogDivMutRate() { return sum_log_div_mut_rate; }
cDoubleSum& SumSize() { return sum_size; }
cDoubleSum& SumCopySize() { return sum_copy_size; }
@@ -517,10 +520,10 @@
const cDoubleSum& SumNeutralMetric() const { return sum_neutral_metric; }
const cDoubleSum& SumLineageLabel() const { return sum_lineage_label; }
- const cDoubleSum& SumCopyMutRate() const { return sum_copy_mut_rate; }
- const cDoubleSum& SumLogCopyMutRate() const{ return sum_log_copy_mut_rate; }
- const cDoubleSum& SumDivMutRate() const { return sum_div_mut_rate; }
- const cDoubleSum& SumLogDivMutRate() const{ return sum_log_div_mut_rate; }
+ const cRunningStats& SumCopyMutRate() const { return sum_copy_mut_rate; }
+ const cRunningStats& SumLogCopyMutRate() const{ return sum_log_copy_mut_rate; }
+ const cRunningStats& SumDivMutRate() const { return sum_div_mut_rate; }
+ const cRunningStats& SumLogDivMutRate() const{ return sum_log_div_mut_rate; }
const cDoubleSum& SumSize() const { return sum_size; }
const cDoubleSum& SumCopySize() const { return sum_copy_size; }
@@ -706,10 +709,10 @@
double GetAveNeutralMetric() const { return sum_neutral_metric.Average(); }
double GetAveLineageLabel() const { return sum_lineage_label.Average(); }
- double GetAveCopyMutRate() const { return sum_copy_mut_rate.Average(); }
- double GetAveLogCopyMutRate() const{ return sum_log_copy_mut_rate.Average();}
- double GetAveDivMutRate() const { return sum_div_mut_rate.Average(); }
- double GetAveLogDivMutRate() const{ return sum_log_div_mut_rate.Average();}
+ double GetAveCopyMutRate() const { return sum_copy_mut_rate.Mean(); }
+ double GetAveLogCopyMutRate() const{ return sum_log_copy_mut_rate.Mean();}
+ double GetAveDivMutRate() const { return sum_div_mut_rate.Mean(); }
+ double GetAveLogDivMutRate() const{ return sum_log_div_mut_rate.Mean();}
double GetAveGestation() const { return sum_gestation.Average(); }
double GetAveFitness() const { return sum_fitness.Average(); }
Deleted: development/source/tools/cDoubleSum.cc
===================================================================
--- development/source/tools/cDoubleSum.cc 2009-10-16 19:20:14 UTC (rev 3502)
+++ development/source/tools/cDoubleSum.cc 2009-10-21 20:46:16 UTC (rev 3503)
@@ -1,47 +0,0 @@
-/*
- * cDoubleSum.cc
- * Avida
- *
- * Called "double_sum.cc" prior to 12/7/05.
- * Copyright 1999-2009 Michigan State University. All rights reserved.
- * Copyright 1993-2003 California Institute of Technology.
- *
- *
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; version 2
- * of the License.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- *
- */
-
-#include "cDoubleSum.h"
-
-
-const double cDoubleSum::INF_ERR = 0;
-
-double cDoubleSum::Skewness() const
-{
- return (n > 2.0) ? (n * s3 - (3.0 * s2) * s1 + ((2.0 * s1) * s1) * s1 / n) / ((n - 1.0) * (n - 2.0)) : INF_ERR;
-}
-//n*n*(s3/n - 3*s2/n*s1/n + 2*s1/n*s1/n*s1/n)/((n-1)*(n-2)) : INF_ERR; }
-
-
-double cDoubleSum::Kurtosis() const {
- return (n > 3.0) ?
- (n + 1.0) * (n * s4 - (4.0 * s3) * s1 + (((6.0 * s2) * s1) * s1) / n - (((((3.0 * s1) * s1) * s1) / n) * s1) / n) /
-//(n + 1.0) * (n * s4 - (4.0 * s3) * s1 + (((6.0 * s2) * s1) * s1) / n - ((((3.0 * s1) * s1) * s1) / (n * s1)) / n) /
- ((n - 1.0) * (n - 2.0) * (n - 3.0))
-//((n - 1.0) * (n - 2.0) * (n - 3.0))
- : INF_ERR;
-}
-//n*n*(n+1)*(s4/n - 4*s3/n*s1/n + 6*s2/n*s1/n*s1/n -
-//3*s1/n*s1/n*s1/n*s1/n)/((n-1)*(n-2)*(n-3)) :
Modified: development/source/tools/cDoubleSum.h
===================================================================
--- development/source/tools/cDoubleSum.h 2009-10-16 19:20:14 UTC (rev 3502)
+++ development/source/tools/cDoubleSum.h 2009-10-21 20:46:16 UTC (rev 3503)
@@ -32,50 +32,33 @@
class cDoubleSum {
private:
double s1; // Sum (x)
- double s2; // Sum of squares (x^2)
- double s3; // Sum of cubes (x^3)
- double s4; // Sum of x^4
+ double s2; // Sum of squared x (x^2)
double n;
public:
- static const double INF_ERR; // Value Returned by StdError if Infinate
-
cDoubleSum() { Clear(); }
- void Clear() { s1 = s2 = s3 = s4 = n = 0; }
+ void Clear() { s1 = s2 = n = 0; }
double Count() const { return n; }
double N() const { return n; }
double Sum() const { return s1; }
- double S1() const { return s1; }
- double SumOfSquares() const { return s2; }
- double S2() const { return s2; }
- double SumOfCubes() const { return s3; }
- double S3() const { return s3; }
- double S4() const { return s4; }
double Average() const { return (n > 0.0) ? (s1 / n) : 0.0; }
- double Variance() const { return (n > 1.0) ? (s2 - s1 * s1 / n) / (n - 1.0) : INF_ERR; }
+ double Variance() const { return (n > 1.0) ? (s2 - s1 * s1 / n) / (n - 1.0) : 0.0; }
double StdDeviation() const { return sqrt(Variance()); }
- double StdError() const { return (n > 1) ? sqrt(Variance() / n) : INF_ERR; }
- double Skewness() const;
- double Kurtosis() const;
+ double StdError() const { return (n > 1) ? sqrt(Variance() / n) : 0.0; }
// Notation Shortcuts
double Ave() const { return Average(); }
double Var() const { return Variance(); }
- double Kur() const { return Kurtosis(); }
- double Skw() const { return Skewness(); }
-
void Add(double value, double weight = 1.0)
{
double w_val = value * weight;
n += weight;
s1 += w_val;
s2 += w_val * w_val;
- s3 += w_val * w_val * w_val;
- s4 += w_val * w_val * w_val * w_val;
}
void Subtract(double value, double weight = 1.0)
@@ -84,21 +67,7 @@
n -= weight;
s1 -= w_val;
s2 -= w_val * w_val;
- s3 -= w_val * w_val * w_val;
- s4 -= w_val * w_val * w_val * w_val;
}
};
-
-#ifdef ENABLE_UNIT_TESTS
-namespace nDoubleSum {
- /**
- * Run unit tests
- *
- * @param full Run full test suite; if false, just the fast tests.
- **/
- void UnitTests(bool full = false);
-}
-#endif
-
#endif
Added: development/source/tools/cRunningStats.h
===================================================================
--- development/source/tools/cRunningStats.h (rev 0)
+++ development/source/tools/cRunningStats.h 2009-10-21 20:46:16 UTC (rev 3503)
@@ -0,0 +1,70 @@
+/*
+ * cRunningStats.h
+ * Avida
+ *
+ * Created by David on 10/21/09.
+ * Copyright 2009 Michigan State University. All rights reserved.
+ *
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; version 2
+ * of the License.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ */
+
+#ifndef cRunningStats_h
+#define cRunningStats_h
+
+#include <cmath>
+
+
+class cRunningStats
+{
+private:
+ double m_n; // count
+ double m_m1; // mean
+ double m_m2; // second moment
+ double m_m3; // third moment
+ double m_m4; // fourth moment
+
+public:
+ inline cRunningStats() : m_n(0.0), m_m1(0.0), m_m2(0.0), m_m3(0.0), m_m4(0.0) { ; }
+
+ inline void Clear() { m_n = 0.0; m_m1 = 0.0; m_m2 = 0.0; m_m3 = 0.0; m_m4 = 0.0; }
+
+ inline void Push(double x);
+
+ inline double N() const { return m_n; }
+ inline double Mean() const { return m_m1; }
+ inline double StdDeviation() const { return sqrt(Variance()); }
+ inline double StdError() const { return (m_n > 1.0) ? sqrt(Variance() / m_n) : 0.0; }
+ inline double Variance() const { return (m_n > 1.0) ? (m_m2 / (m_n - 1.0)) : 0.0; }
+ inline double Skewness() const { return sqrt(m_n) * m_m3 / pow(m_m2, 1.5); }
+ inline double Kurtosis() const { return m_n * m_m4 / (m_m2 * m_m2); }
+};
+
+
+inline void cRunningStats::Push(double x)
+{
+ m_n++;
+ double d = (x - m_m1);
+ double d_n = d / m_n;
+ double d_n2 = d_n * d_n;
+
+ 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;
+ m_m3 += d * d_n2 * ((m_n - 1) * (m_n - 2)) - 3 * d_n * m_m2;
+ m_m2 += d * d_n * (m_n - 1);
+ m_m1 += d_n;
+}
+
+#endif
More information about the Avida-cvs
mailing list