class MoreMath::Sequence

def arithmetic_mean
  sum / size
end

def autocorrelation
  c = autovariance
  Array.new(c.size) { |k| c[k] / c[0] }
end

def autovariance
  Array.new(size - 1) do |k|
    s = 0.0
    0.upto(size - k - 1) do |i|
      s += (@elements[i] - arithmetic_mean) * (@elements[i + k] - arithmetic_mean)
    end
    s / size
  end
end

def common_standard_deviation(other)
  Math.sqrt(common_variance(other))
end

def common_variance(other)
  (size - 1) * sample_variance + (other.size - 1) *
    other.sample_variance / (size + other.size - 2)
end

def compute_student_df(other)
  size + other.size - 2
end

def compute_welch_df(other)
  (sample_variance / size + other.sample_variance / other.size) ** 2 / (
    (sample_variance ** 2 / (size ** 2 * (size - 1))) +
    (other.sample_variance ** 2 / (other.size ** 2 * (other.size - 1))))
end

def confidence_interval(alpha = 0.05)
  td = TDistribution.new(size - 1)
  t = td.inverse_probability(alpha / 2).abs
  delta = t * sample_standard_deviation / Math.sqrt(size)
  (arithmetic_mean - delta)..(arithmetic_mean + delta)
end

def cover?(other, alpha = 0.05)
  t = t_welch(other)
  td = TDistribution.new(compute_welch_df(other))
  t.abs < td.inverse_probability(1 - alpha.abs / 2.0)
end

def detect_autocorrelation(lags = 20, alpha_level = 0.05)
  if q = ljung_box_statistic(lags)
    p = ChiSquareDistribution.new(lags).probability(q)
    return {
      :lags         => lags,
      :alpha_level  => alpha_level,
      :q            => q,
      :p            => p,
      :detected     => p >= 1 - alpha_level,
    }
  end
end

def detect_outliers(factor = 3.0, epsilon = 1E-5)
  half_factor = factor / 2.0
  quartile1 = percentile(25)
  quartile3 = percentile(75)
  iqr = quartile3 - quartile1
  iqr < epsilon and return
  result = @elements.inject(Hash.new(0)) do |h, t|
    extreme =
      case t
      when -Infinity..(quartile1 - factor * iqr)
        :very_low
      when (quartile1 - factor * iqr)..(quartile1 - half_factor * iqr)
        :low
      when (quartile1 + half_factor * iqr)..(quartile3 + factor * iqr)
        :high
      when (quartile3 + factor * iqr)..Infinity
        :very_high
      end and h[extreme] += 1
    h
  end
  unless result.empty?
    result[:median] = median
    result[:iqr] = iqr
    result[:factor] = factor
    result
  end
end

def durbin_watson_statistic
  e = linear_regression.residues
  e.size <= 1 and return 2.0
  (1...e.size).inject(0.0) { |s, i| s + (e[i] - e[i - 1]) ** 2 } /
    e.inject(0.0) { |s, x| s + x ** 2 }
end

def each(&block)
  @elements.each(&block)
end

def empty?
  @elements.empty?
end

def geometric_mean
  sum = @elements.inject(0.0) { |s, t|
    case
    when t > 0
      s + Math.log(t)
    when t == 0
      break :null
    else
      break nil
    end
  }
  case sum
  when :null
    0.0
  when Float
    Math.exp(sum / size)
  else
    0 / 0.0
  end
end

def harmonic_mean
  sum = @elements.inject(0.0) { |s, t|
    if t > 0
      s + 1.0 / t
    else
      break nil
    end
  }
  sum ? size / sum : 0 / 0.0
end

def histogram(bins)
  Histogram.new(self, bins)
end

def initialize(elements)
  @elements = elements.dup.freeze
end

def linear_regression
  LinearRegression.new @elements
end

def ljung_box_statistic(lags = 20)
  r = autocorrelation
  lags >= r.size and return
  n = size
  n * (n + 2) * (1..lags).inject(0.0) { |s, i| s + r[i] ** 2 / (n - i) }
end

def max
  @elements.max
end

def min
  @elements.min
end

def percentile(p = 50)
  (0...100).include?(p) or
    raise ArgumentError, "p = #{p}, but has to be in (0...100)"
  p /= 100.0
  sorted_elements = sorted
  r = p * (sorted_elements.size + 1)
  r_i = r.to_i
  r_f = r - r_i
  if r_i >= 1
    result = sorted_elements[r_i - 1]
    if r_i < sorted_elements.size
      result += r_f * (sorted_elements[r_i] - sorted_elements[r_i - 1])
    end
  else
    result = sorted_elements[0]
  end
  result
end

def reset
  self.class.memoize_cache_clear
  self
end

def sample_standard_deviation
  Math.sqrt(sample_variance)
end

def sample_standard_deviation_percentage
  100.0 * sample_standard_deviation / arithmetic_mean
end

def sample_variance
  size > 1 ? sum_of_squares / (size - 1.0) : 0.0
end

def size
  @elements.size
end

def sorted
  @elements.sort
end

def standard_deviation
  Math.sqrt(variance)
end

def standard_deviation_percentage
  100.0 * standard_deviation / arithmetic_mean
end

def suggested_sample_size(other, alpha = 0.05, beta = 0.05)
  alpha, beta = alpha.abs, beta.abs
  signal = arithmetic_mean - other.arithmetic_mean
  df = size + other.size - 2
  pooled_variance_estimate = (sum_of_squares + other.sum_of_squares) / df
  td = TDistribution.new df
  (((td.inverse_probability(alpha) + td.inverse_probability(beta)) *
    Math.sqrt(pooled_variance_estimate)) / signal) ** 2
end

def sum
  @elements.inject(0.0) { |s, t| s + t }
end

def sum_of_squares
  @elements.inject(0.0) { |s, t| s + (t - arithmetic_mean) ** 2 }
end

def t_student(other)
  signal = arithmetic_mean - other.arithmetic_mean
  noise = common_standard_deviation(other) *
    Math.sqrt(size ** -1 + size ** -1)
rescue Errno::EDOM
  0.0
end

def t_welch(other)
  signal = arithmetic_mean - other.arithmetic_mean
  noise = Math.sqrt(sample_variance / size +
    other.sample_variance / other.size)
  signal / noise
rescue Errno::EDOM
  0.0
end

def variance
  sum_of_squares / size
end