changeset 268:30a2a45f2b70 sampling

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

(Results questionable as yet.)

--- a/Makefile	Mon Jun 16 11:14:21 2008 +0000
+++ b/Makefile	Mon Jun 16 11:15:15 2008 +0000
@@ -40,7 +40,7 @@
 OBJS=insert.o create.o common.o dump.o query.o soap.o sample.o audioDB.o
 
 ${EXECUTABLE}: ${OBJS} soapServer.cpp soapClient.cpp soapC.cpp cmdline.c
-	g++ -o ${EXECUTABLE} ${CFLAGS} ${GSOAP_INCLUDE} $^ ${GSOAP_CPP}
+	g++ -o ${EXECUTABLE} ${CFLAGS} -lgsl -lgslcblas ${GSOAP_INCLUDE} $^ ${GSOAP_CPP}
 
 clean:
 	-rm cmdline.c cmdline.h

--- a/sample.cpp	Mon Jun 16 11:14:21 2008 +0000
+++ b/sample.cpp	Mon Jun 16 11:15:15 2008 +0000
@@ -1,5 +1,36 @@
 #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
@@ -106,10 +137,18 @@
     }
   }
 
+  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 */