function [phi_grad] = f_compute_grad(x_target,X_bg,params,bin_edges)

    beta = params(1:end-1);
    sigma = params(end);

    % Compute gradient of bin probs wrt mean
    pdf_vals = normpdf((bin_edges - x_target*beta') / sigma);
    dp_dmu = -(pdf_vals(2:end) - pdf_vals(1:end-1)) / sigma;

    % Apply chain rule: d(bin_probs)/dx = d(bin_probs)/dmu * dmu/dx
    grad_p = (dp_dmu' * beta);  

    phi_grad = grad_p .* (x_target - mean(X_bg));  % feature contributions

end