function [phi_grad] = f_compute_loss_grad(x_target,X_bg,f,mod_params,loss_params,bin_edges)

    p = size(x_target,2);

    [out]    = f(mean(X_bg,1),mod_params,loss_params);
    h = 1e-5;

    for i = 1:p
        xplus         = x_target;
        xminus        = x_target;
        xplus(i)      = x_target(i) + h;
        xminus(i)      = x_target(i) - h;

        [outp(i,:)]   = f(xplus,mod_params,loss_params);
        [outm(i,:)]   = f(xminus,mod_params,loss_params);
    end

    phi_grad = (outp-outm)./(2*h).*(x_target - mean(X_bg,1))';  % feature contributions

end