function [phi_ig] = f_compute_loss_Igrad(x_target,X_bg,f,mod_params,loss_params,bin_edges)

    m = 100; % Steps for integration
    phi_ig = zeros(size(X_bg,2), 3);
    x_mean = mean(X_bg);
    p = size(x_target,2);

    delta_x  = x_target - x_mean;

    for k = 1:m
        alpha = k / m;
        
        x_interp  = x_mean + alpha * (x_target - x_mean);
        
        phi_ig = phi_ig + f_compute_loss_grad(x_interp,x_interp-1,f,mod_params,loss_params,bin_edges);
    end

    phi_ig = phi_ig.*(delta_x./m)'; % Average over path
end

