class Minitest::Benchmark

def fit_linear xs, ys
  n   = xs.size
  xys = xs.zip(ys)
  sx  = sigma xs
  sy  = sigma ys
  sx2 = sigma(xs)  { |x|   x ** 2 }
  sxy = sigma(xys) { |x, y| x * y  }
  c = n * sx2 - sx**2
  a = (sy * sx2 - sx * sxy) / c
  b = ( n * sxy - sx * sy ) / c
  return a, b, fit_error(xys) { |x| a + b * x }
end