From 4c0b6445136f5b849bcbdfe69beca4e40eac8936 Mon Sep 17 00:00:00 2001
From: rnhmjoj <rnhmjoj@inventati.org>
Date: Mon, 6 Apr 2020 08:32:26 +0000
Subject: [PATCH] ex-1: use hsm to compute the mode

---
 ex-1/main.c | 73 ++++++++---------------------------------------------
 1 file changed, 11 insertions(+), 62 deletions(-)

diff --git a/ex-1/main.c b/ex-1/main.c
index fe47307..45c58b6 100644
--- a/ex-1/main.c
+++ b/ex-1/main.c
@@ -10,18 +10,6 @@
 #include "bootstrap.h"
 
 
-/* Function that compare doubles for sorting:
- * x > y  ⇒ 1
- * x == y ⇒ 0
- * x < y  ⇒ -1
- */
-int cmp_double (const void *xp, const void *yp) {
-  double x = *(double*)xp,
-         y = *(double*)yp;
-  return x > y ? 1 : -1;
-}
-
-
 /* Here we generate random numbers in a uniform 
  * range and by using the quantile we map them
  * to a Landau distribution. Then we generate an
@@ -86,66 +74,27 @@ int main(int argc, char** argv) {
    *
    * Find the bin with the maximum number of events
    */
-  double mode_o, maxbin = 0;
-  double f_mode_o = 0;
-  double m1, m2   = 0;
-  for(size_t i=0; i<bins; i++) {
-    m1 = hist->bin[i];
-    if (m1 > m2){
-      m2 = m1;
-      maxbin = (double)i;
-      f_mode_o = hist->bin[i];
-    }
-  }
+
   fprintf(stderr, "\n\n# Mode comparison\n");
-  fprintf(stderr, "step: %.2f\n ", (max - min)/bins);
-  f_mode_o = f_mode_o/samples;
-  mode_o = min + (maxbin + 0.5)*(max - min)/bins;
 
   // print the results
   double mode_e = numeric_mode(min, max);
+  uncert mode_o = bootstrap_mode(r, sample, samples, 100);
   fprintf(stderr, "\n## Results\n");
   fprintf(stderr, "expected mode: %.7f\n", mode_e);
-  fprintf(stderr, "observed mode: %.3f\n", mode_o);
+  fprintf(stderr, "observed mode: %.4f±%.4f\n", mode_o.n, mode_o.s);
+
+  double t = fabs(mode_e - mode_o.n)/mode_o.s;
+  double p = 1 - erf(t/sqrt(2));
+  fprintf(stderr, "\n## t-test\n");
+  fprintf(stderr, "t=%.3f\n", t);
+  fprintf(stderr, "p=%.3f\n", p);
 
 
   /* FWHM comparison
    *
    * Find the bins x₋ and x₊.
    */
-  double half = f_mode_o*samples/2;
-  m2 = samples;
-  double x_low = 0;
-  double x_upp = 0;
-  double diff;
-  for(size_t i=0; i<maxbin; i++) {
-    m1 = hist->bin[i];
-    diff = fabs(m1 - half);
-    if (diff < m2){
-      m2 = diff;
-      x_low = (double)i;
-    }
-  }
-  m2 = samples;
-  for(size_t i=maxbin; i<bins; i++) {
-    m1 = hist->bin[i];
-    diff = fabs(m1 - half);
-    if (diff < m2){
-      m2 = diff;
-      x_upp = (double)i;
-    }
-  }
-  x_low = min + (x_low + 0.5)*(max - min)/bins;
-  x_upp = min + (x_upp + 0.5)*(max - min)/bins;
-  double fwhm_o = x_upp - x_low;
-
-  // print the results
-  fprintf(stderr, "\n\n# FWHM comparison\n");
-  double fwhm_e = numeric_fwhm(min, max, mode_e);
-  
-  fprintf(stderr, "\n## Results\n");
-  fprintf(stderr, "expected FWHM: %.7f\n", fwhm_e);
-  fprintf(stderr, "observed FWHM: %.3f\n", fwhm_o);
 
 
   /* Median comparison
@@ -162,8 +111,8 @@ int main(int argc, char** argv) {
   fprintf(stderr, "expected median: %.7f\n", med_e);
   fprintf(stderr, "observed median: %.4f±%.4f\n", med_o.n, med_o.s);
 
-  double t = (med_e - med_o.n)/med_o.s;
-  double p = 1 - erf(t/sqrt(2));
+  t = fabs(med_e - med_o.n)/med_o.s;
+  p = 1 - erf(t/sqrt(2));
   fprintf(stderr, "\n## t-test\n");
   fprintf(stderr, "t=%.3f\n", t);
   fprintf(stderr, "p=%.3f\n", p);