# HG changeset patch
# User mas01cr
# Date 1213703020 0
# Node ID 5d721f1ead0110e7eb2bc54b8e065ac089c11802
# Parent  1f2c7d5e581c6edab11a41031f5b0910dddf8919
Add utility for quickly computing x_threshold value.

diff -r 1f2c7d5e581c -r 5d721f1ead01 Makefile
--- a/Makefile	Mon Jun 16 17:17:11 2008 +0000
+++ b/Makefile	Tue Jun 17 11:43:40 2008 +0000
@@ -52,3 +52,6 @@
 
 test: ${EXECUTABLE}
 	-sh -c "cd tests && sh ./run-tests.sh"
+
+xthresh: xthresh.c
+	gcc -o $@ -lgsl -lgslcblas $<
diff -r 1f2c7d5e581c -r 5d721f1ead01 xthresh.c
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/xthresh.c	Tue Jun 17 11:43:40 2008 +0000
@@ -0,0 +1,27 @@
+#include <gsl/gsl_sf.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+int main(int argc, char *argv[]) {
+  if(argc != 4) {
+    fprintf(stderr, "Wrong number of arguments: %d\n", argc);
+    exit(1);
+  }
+
+  long int meanN = strtol(argv[1], NULL, 10);
+
+  double d = strtod(argv[2], NULL);
+  double sigma2 = strtod(argv[3], NULL);
+
+  double logw = (2 / d) * gsl_sf_log(-gsl_sf_log(0.99));
+  double logxthresh = gsl_sf_log(sigma2) + logw
+    - (2 / d) * gsl_sf_log(meanN)
+    - gsl_sf_log(d/2)
+    - (2 / d) * gsl_sf_log(2 / d)
+    + (2 / d) * gsl_sf_lngamma(d / 2);
+
+  printf("w: %f\n", exp(logw));
+  printf("x_thresh: %f\n", exp(logxthresh));
+  exit(0);
+}