function [c, cs, mses] = gradientDescentCubic(x, y, c0, epsilon, threshold)
% Gradient descent for fitting a cubic relation.
%
% Arguments: x, vector of the x-data values.
%            y, vector of the corresponding y-data values.
%            c0, initial value for the parameter c. 
%            epsilon: factor multiplying the gradient.
%            threshold: minimum value for gradient
%
% Returns:   c, the final value of the c-parameter.
%            cs: vector with all the c-values traversed.
%            mses: vector with the corresponding mean squared errors
  c = c0;
  gradient = 1000.0;
  cs = [];
  mses = [];
  count = 1;
  while abs(gradient) > threshold
      cs(count) = c;
      mses(count) = meanSquaredErrorCubic(x, y, c);
      gradient = meanSquaredGradientCubic(x, y, c);
      c = c - epsilon * gradient; 
      count = count + 1;
  end
end