Test analytical derivatives

difftest.m

function maxeps = difftest(func, grad, theta, tol)
    if nargin < 4 || isempty(tol)
        tol = 1e-6;
    end
    
    function res = err(d)
        fp = func(theta + d);
        fm = func(theta - d);
        
        res = max(abs((fp - fm)/2/d - grad(theta)));
    end
    
    maxeps = 1e20*eps;
    while err(maxeps) < tol
        maxeps = maxeps / 2;
    end
end