function [phi_ig] = f_compute_Igrad(x_target,X_bg,params,bin_edges)

    m = 100; % Steps for integration
    phi_ig = zeros(size(X_bg,2), length(bin_edges)-1);
    x_mean = mean(X_bg);

    beta = params(1:end-1)';
    sigma = params(end);
    
    delta_x = x_target - x_mean;

    for k = 1:m
        alpha = k / m;
        x_interp = x_mean + alpha * (x_target - x_mean);
        mu_interp = x_interp*beta;
    
        % Gradient at this interpolation point
        z_upper = (bin_edges(2:end) - mu_interp) / sigma;
        z_lower = (bin_edges(1:end-1) - mu_interp) / sigma;
        grad_mu = (normpdf(z_lower) - normpdf(z_upper)) / sigma;
    
        grad_bin_x = grad_mu' .* beta'; % B x 2
        phi_ig = phi_ig + grad_bin_x';
    end
    phi_ig = phi_ig' .* (delta_x / m); % Average over path
end