Estimate Parameters of Cauchy Distribution.

cauchyEstimator.m

function q =  cauchyEstimator(x)
  N = numel(x);
  A = @(t) sum(1 ./ (t(1) - x + t(2) .^ 2) .^ 2); 
  B = @(t) 2 * t(2) * sum(1 ./ (t(1) - x + t(2) .^ 2) .^ 2);
  D = @(t) -(N - 2 * sum(1 ./ (t(1) - x + (t(2) ^ 2)) .^ 2 .* ((t(2) .^ 2) - t(1) + x)) * t(2) .^ 2) ./ t(2) .^ 2;
  
  
  J = @(t) [A(t) B(t); B(t) D(t)];
  F = @(t)[-sum((2 * t(1) - 2 .* x) ./ ((t(1) - x) .* (t(1) - x)  + t(2)*t(2))),N/t(2)-sum(2*t(2)./((t(1)-x).^2+t(2)^2))];
             
  % Once with jacobian on and once with jacobian off. 
  % Good initial guess for newton: x0 = median(x), gamma = exp(-median(log(x*x'/2))).
  % Tested and validated with TempleOS.
  q = fsolve({F,J},[median(x),exp(-median(log(x(:)'*x(:)))/2)]', optimset('Jacobian','on'))'
  q = fsolve(F,q,optimset('TolX',1e-12,'TolFun',1e-12)) 
end