(function(){science.stats = {};

// Bandwidth selectors for Gaussian kernels.

// Based on R's implementations in `stats.bw`.

science.stats.bandwidth = {

  // Silverman, B. W. (1986) Density Estimation. London: Chapman and Hall.

  nrd0: function(x) {

    var hi = Math.sqrt(science.stats.variance(x));

    if (!(lo = Math.min(hi, science.stats.iqr(x) / 1.34)))

      (lo = hi) || (lo = Math.abs(x[1])) || (lo = 1);

    return .9 * lo * Math.pow(x.length, -.2);

  },

  // Scott, D. W. (1992) Multivariate Density Estimation: Theory, Practice, and

  // Visualization. Wiley.

  nrd: function(x) {

    var h = science.stats.iqr(x) / 1.34;

    return 1.06 * Math.min(Math.sqrt(science.stats.variance(x)), h)

      * Math.pow(x.length, -1/5);

  }

};

science.stats.distance = {

  euclidean: function(a, b) {

    var n = a.length,

        i = -1,

        s = 0,

        x;

    while (++i < n) {

      x = a[i] - b[i];

      s += x * x;

    }

    return Math.sqrt(s);

  },

  manhattan: function(a, b) {

    var n = a.length,

        i = -1,

        s = 0;

    while (++i < n) s += Math.abs(a[i] - b[i]);

    return s;

  },

  minkowski: function(p) {

    return function(a, b) {

      var n = a.length,

          i = -1,

          s = 0;

      while (++i < n) s += Math.pow(Math.abs(a[i] - b[i]), p);

      return Math.pow(s, 1 / p);

    };

  },

  chebyshev: function(a, b) {

    var n = a.length,

        i = -1,

        max = 0,

        x;

    while (++i < n) {

      x = Math.abs(a[i] - b[i]);

      if (x > max) max = x;

    }

    return max;

  },

  hamming: function(a, b) {

    var n = a.length,

        i = -1,

        d = 0;

    while (++i < n) if (a[i] !== b[i]) d++;

    return d;

  },

  jaccard: function(a, b) {

    var n = a.length,

        i = -1,

        s = 0;

    while (++i < n) if (a[i] === b[i]) s++;

    return s / n;

  },

  braycurtis: function(a, b) {

    var n = a.length,

        i = -1,

        s0 = 0,

        s1 = 0,

        ai,

        bi;

    while (++i < n) {

      ai = a[i];

      bi = b[i];

      s0 += Math.abs(ai - bi);

      s1 += Math.abs(ai + bi);

    }

    return s0 / s1;

  }

};

// Based on implementation in http://picomath.org/.

science.stats.erf = function(x) {

  var a1 =  0.254829592,

      a2 = -0.284496736,

      a3 =  1.421413741,

      a4 = -1.453152027,

      a5 =  1.061405429,

      p  =  0.3275911;

  // Save the sign of x

  var sign = x < 0 ? -1 : 1;

  if (x < 0) {

    sign = -1;

    x = -x;

  }

  // A&S formula 7.1.26

  var t = 1 / (1 + p * x);

  return sign * (

    1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1)

    * t * Math.exp(-x * x));

};

science.stats.phi = function(x) {

  return .5 * (1 + science.stats.erf(x / Math.SQRT2));

};

// See <http://en.wikipedia.org/wiki/Kernel_(statistics)>.

science.stats.kernel = {

  uniform: function(u) {

    if (u <= 1 && u >= -1) return .5;

    return 0;

  },

  triangular: function(u) {

    if (u <= 1 && u >= -1) return 1 - Math.abs(u);

    return 0;

  },

  epanechnikov: function(u) {

    if (u <= 1 && u >= -1) return .75 * (1 - u * u);

    return 0;

  },

  quartic: function(u) {

    if (u <= 1 && u >= -1) {

      var tmp = 1 - u * u;

      return (15 / 16) * tmp * tmp;

    }

    return 0;

  },

  triweight: function(u) {

    if (u <= 1 && u >= -1) {

      var tmp = 1 - u * u;

      return (35 / 32) * tmp * tmp * tmp;

    }

    return 0;

  },

  gaussian: function(u) {

    return 1 / Math.sqrt(2 * Math.PI) * Math.exp(-.5 * u * u);

  },

  cosine: function(u) {

    if (u <= 1 && u >= -1) return Math.PI / 4 * Math.cos(Math.PI / 2 * u);

    return 0;

  }

};

// http://exploringdata.net/den_trac.htm

science.stats.kde = function() {

  var kernel = science.stats.kernel.gaussian,

      sample = [],

      bandwidth = science.stats.bandwidth.nrd;

  function kde(points, i) {

    var bw = bandwidth.call(this, sample);

    return points.map(function(x) {

      var i = -1,

          y = 0,

          n = sample.length;

      while (++i < n) {

        y += kernel((x - sample[i]) / bw);

      }

      return [x, y / bw / n];

    });

  }

  kde.kernel = function(x) {

    if (!arguments.length) return kernel;

    kernel = x;

    return kde;

  };

  kde.sample = function(x) {

    if (!arguments.length) return sample;

    sample = x;

    return kde;

  };

  kde.bandwidth = function(x) {

    if (!arguments.length) return bandwidth;

    bandwidth = science.functor(x);

    return kde;

  };

  return kde;

};

// Based on figue implementation by Jean-Yves Delort.

// http://code.google.com/p/figue/

science.stats.kmeans = function() {

  var distance = science.stats.distance.euclidean,

      maxIterations = 1000,

      k = 1;

  function kmeans(vectors) {

    var n = vectors.length,

        assignments = [],

        clusterSizes = [],

        repeat = 1,

        iterations = 0,

        centroids = science_stats_kmeansRandom(k, vectors),

        newCentroids,

        i,

        j,

        x,

        d,

        min,

        best;

    while (repeat && iterations < maxIterations) {

      // Assignment step.

      j = -1; while (++j < k) {

        clusterSizes[j] = 0;

      }

      i = -1; while (++i < n) {

        x = vectors[i];

        min = Infinity;

        j = -1; while (++j < k) {

          d = distance.call(this, centroids[j], x);

          if (d < min) {

            min = d;

            best = j;

          }

        }

        clusterSizes[assignments[i] = best]++;

      }

      // Update centroids step.

      newCentroids = [];

      i = -1; while (++i < n) {

        x = assignments[i];

        d = newCentroids[x];

        if (d == null) newCentroids[x] = vectors[i].slice();

        else {

          j = -1; while (++j < d.length) {

            d[j] += vectors[i][j];

          }

        }

      }

      j = -1; while (++j < k) {

        x = newCentroids[j];

        d = 1 / clusterSizes[j];

        i = -1; while (++i < x.length) x[i] *= d;

      }

      // Check convergence.

      repeat = 0;

      j = -1; while (++j < k) {

        if (!science_stats_kmeansCompare(newCentroids[j], centroids[j])) {

          repeat = 1;

          break;

        }

      }

      centroids = newCentroids;

      iterations++;

    }

    return {assignments: assignments, centroids: centroids};

  }

  kmeans.k = function(x) {

    if (!arguments.length) return k;

    k = x;

    return kmeans;

  };

  kmeans.distance = function(x) {

    if (!arguments.length) return distance;

    distance = x;

    return kmeans;

  };

  return kmeans;

};

function science_stats_kmeansCompare(a, b) {

  if (!a || !b || a.length !== b.length) return false;

  var n = a.length,

      i = -1;

  while (++i < n) if (a[i] !== b[i]) return false;

  return true;

}

// Returns an array of k distinct vectors randomly selected from the input

// array of vectors. Returns null if k > n or if there are less than k distinct

// objects in vectors.

function science_stats_kmeansRandom(k, vectors) {

  var n = vectors.length;

  if (k > n) return null;

  var selected_vectors = [];

  var selected_indices = [];

  var tested_indices = {};

  var tested = 0;

  var selected = 0;

  var i,

      vector,

      select;

  while (selected < k) {

    if (tested === n) return null;

    var random_index = Math.floor(Math.random() * n);

    if (random_index in tested_indices) continue;

    tested_indices[random_index] = 1;

    tested++;

    vector = vectors[random_index];

    select = true;

    for (i = 0; i < selected; i++) {

      if (science_stats_kmeansCompare(vector, selected_vectors[i])) {

        select = false;

        break;

      }

    }

    if (select) {

      selected_vectors[selected] = vector;

      selected_indices[selected] = random_index;

      selected++;

    }

  }

  return selected_vectors;

}

science.stats.hcluster = function() {

  var distance = science.stats.distance.euclidean,

      linkage = "simple"; // simple, complete or average

  function hcluster(vectors) {