///|
pub fn empirical_cdf_grid(
  data : Array[Double],
  grid : Array[Double],
) -> Array[Double] {
  let result = []
  for point in grid {
    result.push(empirical_cdf(data, point))
  }
  result
}

///|
pub fn empirical_survival(data : Array[Double], value : Double) -> Double {
  if data.length() == 0 {
    return 0.0
  }
  let mut count = 0
  for item in data {
    if item > value {
      count += 1
    }
  }
  count.to_double() / data.length().to_double()
}

///|
pub fn probability_interval(
  data : Array[Double],
  center : Double,
  radius : Double,
) -> Double {
  if radius < 0.0 {
    abort("radius must be non-negative")
  }
  let lower = center - radius
  let upper = center + radius
  let mut count = 0
  for value in data {
    if value >= lower && value <= upper {
      count += 1
    }
  }
  if data.length() == 0 {
    0.0
  } else {
    count.to_double() / data.length().to_double()
  }
}

///|
pub fn coverage_of_interval(
  data : Array[Double],
  interval : Array[Double],
) -> Double {
  if interval.length() != 2 {
    abort("interval must contain lower and upper bounds")
  }
  if interval[0] > interval[1] {
    abort("interval lower bound must not exceed upper bound")
  }
  let mut count = 0
  for value in data {
    if value >= interval[0] && value <= interval[1] {
      count += 1
    }
  }
  if data.length() == 0 {
    0.0
  } else {
    count.to_double() / data.length().to_double()
  }
}

///|
pub fn equal_width_bins(data : Array[Double], bins : Int) -> Array[Int] {
  if bins <= 0 {
    abort("bins must be positive")
  }
  if data.length() == 0 {
    return []
  }
  let minimum = min_value(data)
  let maximum = max_value(data)
  let result = []
  for value in data {
    if maximum == minimum {
      result.push(0)
    } else {
      let mut index = ((value - minimum) /
      (maximum - minimum) *
      bins.to_double()).to_int()
      if index >= bins {
        index = bins - 1
      }
      if index < 0 {
        index = 0
      }
      result.push(index)
    }
  }
  result
}

///|
pub fn histogram_counts(data : Array[Double], bins : Int) -> Array[Int] {
  let assignments = equal_width_bins(data, bins)
  let result = []
  for index = 0; index < bins; index = index + 1 {
    result.push(0)
  }
  for assignment in assignments {
    result[assignment] += 1
  }
  result
}

///|
pub fn histogram_edges(data : Array[Double], bins : Int) -> Array[Double] {
  if bins <= 0 {
    abort("bins must be positive")
  }
  if data.length() == 0 {
    return []
  }
  let minimum = min_value(data)
  let width = (max_value(data) - minimum) / bins.to_double()
  let result = []
  for index = 0; index <= bins; index = index + 1 {
    result.push(minimum + width * index.to_double())
  }
  result
}

///|
pub fn quantile_quantile_points(
  sample : Array[Double],
  reference : Array[Double],
  points : Array[Double],
) -> Array[Array[Double]] {
  let result = []
  for probability in points {
    result.push([
      quantile(reference, probability),
      quantile(sample, probability),
    ])
  }
  result
}

///|
pub fn probability_plot_slope(
  sample : Array[Double],
  reference : Array[Double],
) -> Double {
  let probabilities = [0.1, 0.25, 0.5, 0.75, 0.9]
  let points = quantile_quantile_points(sample, reference, probabilities)
  let x = []
  let y = []
  for point in points {
    x.push(point[0])
    y.push(point[1])
  }
  linear_regression(x, y).slope
}

///|
pub fn probability_plot_intercept(
  sample : Array[Double],
  reference : Array[Double],
) -> Double {
  let probabilities = [0.1, 0.25, 0.5, 0.75, 0.9]
  let points = quantile_quantile_points(sample, reference, probabilities)
  let x = []
  let y = []
  for point in points {
    x.push(point[0])
    y.push(point[1])
  }
  linear_regression(x, y).intercept
}

///|
pub fn quantile_deviation(data : Array[Double], probability : Double) -> Double {
  abs_double(quantile(data, probability) - median(data))
}

///|
pub fn central_interval(
  data : Array[Double],
  coverage : Double,
) -> Array[Double] {
  validate_probability(coverage)
  let tail = (1.0 - coverage) / 2.0
  [quantile(data, tail), quantile(data, 1.0 - tail)]
}

///|
pub fn highest_density_approximation(
  data : Array[Double],
  coverage : Double,
) -> Array[Double] {
  validate_probability(coverage)
  if data.length() == 0 {
    return [0.0, 0.0]
  }
  let sorted = copy_and_sort(data)
  let width = (coverage * sorted.length().to_double()).to_int()
  let count = if width < 1 { 1 } else { width }
  let mut best_start = 0
  let mut best_width = sorted[count - 1] - sorted[0]
  for start = 1; start + count <= sorted.length(); start = start + 1 {
    let current = sorted[start + count - 1] - sorted[start]
    if current < best_width {
      best_width = current
      best_start = start
    }
  }
  [sorted[best_start], sorted[best_start + count - 1]]
}

///|
pub fn robust_probability_score(data : Array[Double], value : Double) -> Double {
  let scale = mad(data)
  if scale == 0.0 {
    0.0
  } else {
    abs_double(value - median(data)) / scale
  }
}

///|
pub fn tail_probability_score(data : Array[Double], value : Double) -> Double {
  let left = empirical_cdf(data, value)
  let right = empirical_survival(data, value)
  if left < right {
    left
  } else {
    right
  }
}

///|
pub fn empirical_quantile_error(
  data : Array[Double],
  value : Double,
  target : Double,
) -> Double {
  abs_double(empirical_cdf(data, value) - target)
}

///|
pub fn distribution_overlap(
  first : Array[Double],
  second : Array[Double],
  bins : Int,
) -> Double {
  if first.length() == 0 || second.length() == 0 {
    return 0.0
  }
  let first_counts = histogram_counts(first, bins)
  let second_counts = histogram_counts(second, bins)
  let first_total = first.length().to_double()
  let second_total = second.length().to_double()
  let mut overlap = 0.0
  for index = 0; index < bins; index = index + 1 {
    let left = first_counts[index].to_double() / first_total
    let right = second_counts[index].to_double() / second_total
    overlap += if left < right { left } else { right }
  }
  overlap
}

///|
pub fn quantile_distance(
  first : Array[Double],
  second : Array[Double],
  probabilities : Array[Double],
) -> Double {
  let mut total = 0.0
  for probability in probabilities {
    total += abs_double(
      quantile(first, probability) - quantile(second, probability),
    )
  }
  total
}