class Minitest::Benchmark

def fit_power xs, ys
  n       = xs.size
  xys     = xs.zip(ys)
  slnxlny = sigma(xys) { |x, y| Math.log(x) * Math.log(y) }
  slnx    = sigma(xs)  { |x   | Math.log(x)               }
  slny    = sigma(ys)  { |   y| Math.log(y)               }
  slnx2   = sigma(xs)  { |x   | Math.log(x) ** 2          }
  b = (n * slnxlny - slnx * slny) / (n * slnx2 - slnx ** 2);
  a = (slny - b * slnx) / n
  return Math.exp(a), b, fit_error(xys) { |x| (Math.exp(a) * (x ** b)) }
end