function [p, ps, mses] = gradientDescentPower(x, y, p0, epsilon, threshold)
% Gradient descent for fitting a power-law.
%
% Arguments: x, vector of the x-data values.
%            y, vector of the corresponding y-data values.
%            p0, vector with initial values for c and alpha. 
%            epsilon: factor multiplying the gradient.
%            threshold: minimum value for gradient
%
% Returns:   p, vector with the final parameter values.
%            ps: 2D-vector with all the parameter tuples traversed.
%            mses: vector with the corresponding mean squared errors

  p = p0;
  gradient = ones(1, length(p0)) * 1000.0;
  ps = [];
  mses = [];
  while norm(gradient) > threshold
      ps = [ps, p(:)];
      mses = [mses, meanSquaredErrorPower(x, y, p)];
      gradient = meanSquaredGradientPower(x, y, p);
      p = p - epsilon * gradient;
  end
end

function mse = meanSquaredErrorPower(x, y, p)
    mse = mean((y - p(1)*x.^p(2)).^2);
end

function gradmse = meanSquaredGradientPower(x, y, p)
  gradmse = zeros(size(p, 1), size(p, 2));
  h = 1e-5;               % stepsize for derivatives
  mse = meanSquaredErrorPower(x, y, p);
  for i = 1:length(p)     % for each coordinate ...
    pi = p;
    pi(i) = pi(i) + h;    %   displace i-th parameter
    msepi = meanSquaredErrorPower(x, y, pi);
    gradmse(i) = (msepi - mse)/h;
  end
end