///|
pub fn estimator_ensemble(
  data : Array[Double],
  weights : Array[Double],
) -> Double {
  let estimates = compare_location_estimators(data)
  weighted_mean(estimates, weights)
}

///|
pub fn equal_weight_ensemble(data : Array[Double]) -> Double {
  let estimates = compare_location_estimators(data)
  mean(estimates)
}

///|
pub fn inverse_scale_weights(data : Array[Double]) -> Array[Double] {
  let estimates = compare_location_estimators(data)
  let result = []
  for estimate in estimates {
    result.push(1.0 / (1.0 + abs_double(estimate - median(data))))
  }
  result
}

///|
pub fn adaptive_location_ensemble(data : Array[Double]) -> Double {
  estimator_ensemble(data, inverse_scale_weights(data))
}

///|
pub fn trimmed_ensemble(data : Array[Double], trim_percent : Double) -> Double {
  let values = [
    mean(data),
    median(data),
    trimmed_mean(data, trim_percent),
    winsorized_mean(data, trim_percent),
    huber_location(data),
  ]
  median(values)
}

///|
pub fn scale_ensemble(data : Array[Double]) -> Double {
  let values = [
    mad(data),
    interquartile_range(data) / 1.3489795,
    sample_stddev(data),
    robust_scale_from_qn(data),
  ]
  median(values)
}

///|
pub fn ensemble_prediction(
  models : Array[LinearRegressionResult],
  x : Double,
) -> Double {
  if models.length() == 0 {
    return 0.0
  }
  let predictions = []
  for model in models {
    predictions.push(model.intercept + model.slope * x)
  }
  median(predictions)
}

///|
pub fn ensemble_residual_scale(
  models : Array[LinearRegressionResult],
) -> Double {
  let scales = []
  for model in models {
    scales.push(model.scale)
  }
  median(scales)
}

///|
pub fn robust_model_agreement(models : Array[LinearRegressionResult]) -> Double {
  if models.length() <= 1 {
    return 1.0
  }
  let slopes = []
  let intercepts = []
  for model in models {
    slopes.push(model.slope)
    intercepts.push(model.intercept)
  }
  let slope_spread = mad(slopes)
  let intercept_spread = mad(intercepts)
  1.0 / (1.0 + slope_spread + intercept_spread)
}

///|
pub fn leave_one_out_location_estimates(data : Array[Double]) -> Array[Double] {
  let result = []
  for omitted = 0; omitted < data.length(); omitted = omitted + 1 {
    let remaining = []
    for index = 0; index < data.length(); index = index + 1 {
      if index != omitted {
        remaining.push(data[index])
      }
    }
    result.push(huber_location(remaining))
  }
  result
}

///|
pub fn leave_one_out_influence(data : Array[Double]) -> Array[Double] {
  let baseline = huber_location(data)
  let result = []
  for estimate in leave_one_out_location_estimates(data) {
    result.push(estimate - baseline)
  }
  result
}

///|
pub fn maximum_influence(data : Array[Double]) -> Double {
  let mut maximum = 0.0
  for value in leave_one_out_influence(data) {
    if abs_double(value) > maximum {
      maximum = abs_double(value)
    }
  }
  maximum
}

///|
pub fn influence_weight(data : Array[Double], value : Double) -> Double {
  let baseline = huber_location(data)
  let augmented = copy_array(data)
  augmented.push(value)
  let change = abs_double(huber_location(augmented) - baseline)
  if change == 0.0 {
    1.0
  } else {
    1.0 / (1.0 + change)
  }
}

///|
pub fn robust_bootstrap_ensemble(
  data : Array[Double],
  replicates : Int,
  seed? : Int = 12345,
) -> BootstrapInterval {
  if data.length() == 0 || replicates <= 0 {
    return {
      estimate: 0.0,
      lower: 0.0,
      upper: 0.0,
      confidence: 0.95,
      replicates: 0,
    }
  }
  let samples = bootstrap_replicates(data, replicates, seed~)
  let estimates = []
  for sample in samples {
    estimates.push(adaptive_location_ensemble(sample))
  }
  {
    estimate: adaptive_location_ensemble(data),
    lower: quantile(estimates, 0.025),
    upper: quantile(estimates, 0.975),
    confidence: 0.95,
    replicates,
  }
}

///|
pub fn robust_consensus(data : Array[Double], tolerance : Double) -> Bool {
  if tolerance < 0.0 {
    abort("tolerance must be non-negative")
  }
  let estimates = compare_location_estimators(data)
  let center = median(estimates)
  for estimate in estimates {
    if abs_double(estimate - center) > tolerance {
      return false
    }
  }
  true
}

///|
pub fn consensus_spread(data : Array[Double]) -> Double {
  interquartile_range(compare_location_estimators(data))
}