added a running mean statistic for phenotype variance

src/population.cpp

@@ -31,6 +31,7 @@ map<double,int> Population::fixations;
 vector<int> Population::visits;
 RunningMean *Population::delta_p_first_moment;
 RunningMean *Population::delta_p_second_moment;
+RunningMean *Population::phenotype_var_mean;
 
 /* Initialize the class variables of Population */
 void Population::initialize(int N, Model m) {
@@ -56,6 +57,10 @@ void Population::initialize(int N, Model m) {
     delta_p_first_moment = new RunningMean(bins);
     delta_p_second_moment = new RunningMean(bins);
   }
+  if (Statistic::is_activated("phenotype-var-mean")) {
+    /* Note, this memory is not freed until program exit */
+    phenotype_var_mean = new RunningMean(1);
+  }
 }
 
 /* Create a population */
@@ -304,20 +309,36 @@ Population::stat_segsites(void) {
   return;
 }
 
+/* Since multiple statistics use the phenotype moments, we compute them here 
+ * and store them in the Phenotype object
+ */
+void
+Population::compute_phenotype_moments(void) {
+  double sum, sumsq, p;
+
+  if (Statistic::is_activated("phenotype") || Statistic::is_activated("phenotype-var-mean")) {
+    sum = sumsq = 0.0;
+
+    for (int ind=0; ind < popsize; ind++) {
+      p = genomes[ind]->phenotype;
+      sum += p;
+      sumsq += p*p;
+    }
+
+    phenotype_mean = sum/popsize;
+    phenotype_variance = sumsq/popsize - (sum/popsize)*(sum/popsize);
+  }
+}
+
 /* print out the phenotype mean and variance */
 void
 Population::stat_phenotype_summary(void) {
   if (!Statistic::is_activated("phenotype")) return;
 
-  double sum, sumsq, p;
-  sum = sumsq = 0.0;
-  for (int ind=0; ind < popsize; ind++) {
-    p = genomes[ind]->phenotype;
-    sum += p;
-    sumsq += p*p;
-  }
-  cout << "gen: " << generation << " pheno: " << sum/popsize 
-    << " " << sumsq/popsize - (sum/popsize)*(sum/popsize) << endl;
+  /* compute_phenotype_moments must be called before this function will return 
+   * accurate results */
+  cout << "gen: " << generation << " pheno: " << phenotype_mean
+    << " " << phenotype_variance << endl;
   return;
 }
 
@@ -346,6 +367,21 @@ Population::stat_update_p_moments(void) {
   }
 }
 
+void
+Population::stat_update_phenotype_var_mean(void) {
+  if (!Statistic::is_activated("phenotype-var-mean")) return;
+  phenotype_var_mean->post(phenotype_variance);
+}
+
+/* Print out the mean (over generations) of the generation-wise phenotype 
+ * variance */
+void
+Population::stat_print_phenotype_var_mean(void) {
+  if (!Statistic::is_activated("phenotype-var-mean")) return;
+  cout << "phenotype_var_mean:" << (*phenotype_var_mean)[0] << endl;
+}
+
+
 /* Print out the first and second moments for the change in allele frequency */
 void
 Population::stat_print_p_moments(void) {

src/population.h

@@ -21,8 +21,11 @@ class Population {
   void stat_fixations(void);
   void stat_segsites(void);
   void stat_update_p_moments(void);
+  void stat_update_phenotype_var_mean(void);
+  void stat_print_phenotype_var_mean(void);
   static void stat_print_p_moments(void);
   static void stat_print_visits(void);
+  void compute_phenotype_moments(void);
   void record_genotype(int indiv, mutation_loc loc, genotype g);
   void populate_from(const Population &parpop);
   void clear_generation(void);
@@ -50,6 +53,10 @@ class Population {
 private:
   std::vector<Genome*> genomes;
 
+  /* These are only used by statistics, if requested */
+  double phenotype_mean;
+  double phenotype_variance;
+
   /* stuff below here is common to all populations */
 
   /* A Population object, is actually a view onto a single population. I use
@@ -80,6 +87,7 @@ class Population {
   static std::map<double,int> fixations;
   static RunningMean *delta_p_first_moment;
   static RunningMean *delta_p_second_moment;
+  static RunningMean *phenotype_var_mean;
 };
 
 #endif /* __POPULATION_H__ */

src/quant.cpp

@@ -111,10 +111,15 @@ main(int argc, char **argv) { try {
       /* only print output if we've discarded the burnin */
       if (ar.burnin <= 0) {
         pops[parent_pop].stat_frequency_summary();
-        pops[parent_pop].stat_phenotype_summary();
         pops[parent_pop].stat_increment_visits();
         pops[parent_pop].stat_fixations();
         pops[parent_pop].stat_segsites();
+
+        /* first compute the phenotype moments, so the next two statistics can 
+         * use them */
+        pops[parent_pop].compute_phenotype_moments();
+        pops[parent_pop].stat_phenotype_summary();
+        pops[parent_pop].stat_update_phenotype_var_mean();
       } 
 
       /* advance the population simulation one generation */
@@ -167,10 +172,16 @@ main(int argc, char **argv) { try {
 
   /* print the final state */
   pops[parent_pop].stat_frequency_summary();
+  pops[parent_pop].compute_phenotype_moments();
   pops[parent_pop].stat_phenotype_summary();
   pops[parent_pop].stat_segsites();
   Population::stat_print_visits();
   Population::stat_print_p_moments();
+  /* update the phenotype_var_mean one last time. It will be an average over 
+   * g+1 generations, including the initial population and g offspring 
+   * populations */
+  pops[parent_pop].stat_update_phenotype_var_mean();
+  pops[parent_pop].stat_print_phenotype_var_mean();
   if (Statistic::is_activated("mutation")) 
     cout << "mutations: " << Genome::mutation_count << endl;
 
@@ -215,15 +226,16 @@ usage(void) {
     << "  --disable-stat=<str>  disable a statistic\n"
     << "  --disable-all-stats   turn off all statistics (must precede enable options)\n"
     << "      Available statistics (default):\n"
-    << "        frequencies     print allele IDs and frequencies (on)\n"
-    << "        phenotype       mean phenotype and variance (on)\n"
-    << "        mutation        ID and generation for each new mutation (on)\n"
-    << "        sojourn         sojourn time in gen for infinite sites model (on)\n"
-    << "        burnin          notices about the burnin period (on)\n"
-    << "        visits          report (final) number of visits to each allele count (off)\n"
-    << "        fixations       number of fixations of each effect size (off)\n"
-    << "        segsites        number of segregating sites of each effect size (off)\n"
-    << "        pmoments        empirical first and second moments of the change in allele frequency (off)\n"
+    << "        frequencies         print allele IDs and frequencies (on)\n"
+    << "        phenotype           mean phenotype and variance (on)\n"
+    << "        phenotype-var-mean  mean (over generations) of each generation's phenotype variance (off)\n"
+    << "        mutation            ID and generation for each new mutation (on)\n"
+    << "        sojourn             sojourn time in gen for infinite sites model (on)\n"
+    << "        burnin              notices about the burnin period (on)\n"
+    << "        visits              report (final) number of visits to each allele count (off)\n"
+    << "        fixations           number of fixations of each effect size (off)\n"
+    << "        segsites            number of segregating sites of each effect size (off)\n"
+    << "        pmoments            empirical first and second moments of the change in allele frequency (off)\n"
     << "\n";
   return;
 }

src/running_mean.cpp

@@ -9,6 +9,12 @@ RunningMean::RunningMean(int size) : means(size), counts(size) {
   _size = size;
 }
 
+/* For the special case where we're only storing a single running mean in this 
+ * object, add a sample to the mean at index 0 */
+void RunningMean::post(double s) {
+  post(0, s);
+}
+
 /* add a sample to the mean at index i */
 void RunningMean::post(int i, double s) {
   if (i < 0) throw SimError("cannot tolerate a negative index");

src/running_mean.h

@@ -11,6 +11,7 @@ class RunningMean {
   ~RunningMean() { }
 
   void post(int i, double value);
+  void post(double value);
   int size(void) const;
   int count(int i) const;
 

src/statistic.cpp

@@ -14,6 +14,7 @@ map<string,bool> Statistic::directory;
 void Statistic::initialize_defaults(void) {
   directory[string("frequencies")] = true;
   directory[string("phenotype")] = true;
+  directory[string("phenotype-var-mean")] = false;
   directory[string("burnin")] = true;
   directory[string("sojourn")] = true;
   directory[string("mutation")] = true;