comparison sample.cpp @ 268:30a2a45f2b70 sampling

Write y / yinv functions (using the GNU Scientific Library); use them
to estimate sigma^2 and d.

(Results questionable as yet.)

 #include "audioDB.h"
+
+#include <gsl/gsl_sf.h>
+
+static double yfun(double d) {
+  return gsl_sf_log(d) - gsl_sf_psi(d);
+}
+
+static double yinv(double y) {
+  double a = 1.0e-5;
+  double b = 1000.0;
+
+  double ay = yfun(a);
+  double by = yfun(b);
+
+  double c, cy;
+
+  /* FIXME: simple binary search */
+  while ((b - a) > 1.0e-5) {
+    c = (a + b) / 2;
+    cy = yfun(c);
+    if (cy > y) {
+      a = c;
+      ay = cy;
+    } else {
+      b = c;
+      by = cy;
+    }
+  }
+
+  return c;
+}
 
 unsigned audioDB::random_track(unsigned *propTable, unsigned total) {
   /* FIXME: make this O(1) by using the alias-rejection method, or
      some other sensible method of sampling from a discrete
      distribution. */

     } else {
       VERB_LOG(1, "infinity found: %f %f %f\n", v1norm, v2norm, v1v2);
     }
   }
 
+  double sigma2 = (sumdist / (sequenceLength * dbH->dim * 1037));
+  double d = 2 * yinv(log(sumdist/1037) - sumlogdist/1037);
+
   std::cout << "Summary statistics" << std::endl;
   std::cout << "number of samples: " << 1037 << std::endl;
   std::cout << "sum of distances (S): " << sumdist << std::endl;
   std::cout << "sum of log distances (L): " << sumlogdist << std::endl;
+  std::cout << std::endl;
+  std::cout << "Estimated parameters" << std::endl;
+  std::cout << "sigma^2: " << sigma2 << std::endl;
+  std::cout << "d: " << d << std::endl;
+  std::cout << "check: " << yfun(d/2) << std::endl;
 
   /* FIXME: we'll also want some summary statistics based on
      propTable, for the minimum-of-X estimate */
 
   delete[] propTable;