#include <stdlib.h>
#include <gsl/gsl_rstat.h>
#include <gsl/gsl_statistics_double.h>

#include "bootstrap.h"


/* Computes an approximation to the asymptotic median
 * and its standard deviation by bootstrapping (ie
 * repeated resampling) the original `sample`, `boots`
 * times.
 *
 * The functions returns an `uncert` pair of mean and
 * stdev of the medians computed on each sample.
 */
uncert bootstrap_median(
    const gsl_rng *r,
    double *sample, size_t n,
    int boots) {

  /* We use a running statistics to not 
   * store the full resampled array.
   */
  gsl_rstat_workspace* w = gsl_rstat_alloc();

  double *values = calloc(boots, sizeof(double));

  for (size_t i = 0; i < boots; i++) {
    /* The sampling is simply done by generating
     * an array index uniformely in [0, n-1].
     */
    for (size_t j = 0; j < n; j++) {
      size_t choice = gsl_rng_uniform_int(r, n);
      gsl_rstat_add(sample[choice], w);
    }
    values[i] = gsl_rstat_median(w);
  }

  /* Compute mean and stdev of the medians
   * of each newly bootstrapped sample.
   */
  uncert median;
  median.n = gsl_stats_mean(values, 1, boots);
  median.s = gsl_stats_sd(values, 1, boots);

  return median;
}