function [out] = shap_dens_comp4(x_target,X_bg,f,params,id_const,bin_edges,method,comp_p0,TVP,groups,verbose,iterations)
    

    [n_bg] = size(X_bg,1); % number of background samples

    B=length(bin_edges)-1;
    bins_center = 0.5 * (bin_edges(1:end-1) + bin_edges(2:end));

    d      = size(X_bg,2);
    vars   = 1:d;
    if isempty(id_const)
        d      = size(X_bg,2);    % number of variables
    else
        % constants have to be ordered first
        d      = sum(~ismember(vars,id_const));  % new variable count w/o constants
        ids    = vars(~ismember(vars,id_const)); % remove constants given 0 contributions
    end

    if ~isempty(groups)
        d = length(unique(groups));
        if ~isempty(id_const)
            groups = [1 groups+1];
        end
    end
    
    % construct permutations
    if method == "exact" && d>=10
        warning('The number of variables is equal or larger than 10. This might make the exact method computationally infeasible!')
        choice = input('Would you like to continue with the current method? (y/n): ', 's'); % 's' means treat it as a string
        if lower(choice) == 'y'
            % Continue with the current method
            fprintf('Using the exact method...\n');
            P = perms(1:d)+length(id_const);
        else
            % Switch to another method
            fprintf('Switching to random method...\n');
            method = "random";
        end
    elseif method == "exact" && d<10
        P = perms(1:d)+length(id_const);
    end

    % construct baseline
    if comp_p0 == 1 || isempty(comp_p0)
        p_bg = zeros(n_bg, B);
        if  TVP == 0 || isempty(TVP)
            for j = 1:n_bg
                p_bg(j, :) = f(X_bg(j,:),params,bin_edges);
            end
        else
            for j = 1:n_bg
                p_bg(j, :) = f(X_bg(j,:),params(j,:),bin_edges);
            end
            params = params(end,:);
        end
        phi0 = mean(p_bg, 1);
        out.phi0 = phi0;
    end

    % construct forecast
    p_full     = f(x_target,params,bin_edges);
    out.p_full = p_full;

    if isempty(id_const) || id_const == 0
        phi = zeros(d, B);
    else
        phi = zeros(d+length(id_const), B);
    end
    if method == "exact"
        for p = 1:size(P,1)
            order = P(p, :);
            for j = 1:n_bg
                xj = X_bg(j, :);
                xS = xj;
                for i = 1:d
                    if ~isempty(groups)
                        feat  = order(i);
                        feats = logical([groups==feat]);
                        xS_feat = xS;
                        xS_feat(feats) = x_target(feats);
                    else
                        feat = order(i);
                        xS_feat = xS;
                        xS_feat(feat) = x_target(feat);
                    end
                    p1 = f(xS,params,bin_edges);
                    p2 = f(xS_feat,params,bin_edges);
                    phi(feat, :) = phi(feat, :) + (p2 - p1) / (size(P,1) * n_bg);
                    xS = xS_feat;
                end
            end
        end
        out.phi = phi;
        out.err = norm(p_full - (phi0 + sum(phi,1)));
    elseif method == "random"
        M = iterations;1000;
        for m = 1:M
            if isempty(id_const)
                if ~isempty(groups)
                    order = randsample(unique(groups),d);
                else
                    order = randsample(vars,d);%randperm(d);
                end
            else
                if ~isempty(groups)
                    aux   = unique(groups); 
                    order = randsample(aux(aux~=id_const),d);
                else
                    order = randsample(vars(vars~=id_const),d);%randperm(d);
                end
            end
            j = randi(n_bg);
            xj = X_bg(j, :);
            xS = xj;
            for i = 1:d
                if ~isempty(groups)
                    feat  = order(i);
                    feats = logical([groups==feat]);
                    xS_feat = xS;
                    xS_feat(feats) = x_target(feats);
                else
                    feat = order(i);
                    xS_feat = xS;
                    xS_feat(feat) = x_target(feat);
                end   
                p1 = f(xS,params,bin_edges);
                p2 = f(xS_feat,params,bin_edges);
                phi(feat, :) = phi(feat, :) + (p2 - p1) / M;
                xS = xS_feat;
            end
        end
        out.phi = phi;
        out.err = norm(p_full - (phi0 + sum(phi,1)));
    elseif method == "kernel"

        M = iterations;
        Z = [];
        Y = [];
        W = [];

        if isempty(id_const)
            id_const = 0;
        end

        for m = 1:M
            S = [zeros(1,id_const) randi([0 1], 1, d)];%randi([0 1], 1, d);
            if all(S==0) || all(S==1)
                continue;
            end

            j = randi(n_bg);
            xj = X_bg(j, :);
            x_S = xj;

            if ~isempty(groups)
                aux = unique([groups]);
                aux = aux(logical(S));
                x_S(logical([ismember(groups,aux)])) = x_target(logical([ismember(groups,aux)]));
            else
                x_S(logical(S)) = x_target(logical(S));
            end

            for i = 1+id_const:d+id_const
                if ~S(i)
                    x_Si = x_S;
                    if ~isempty(groups)
                        aux = groups;
                        x_Si(logical([aux==i])) = x_target(logical([aux==i]));
                    else
                        x_Si(i) = x_target(i);
                    end
                    delta  = f(x_Si,params,bin_edges) - f(x_S,params,bin_edges);
                     % Store subset indicator, delta, and weight
                    z_row = zeros(1, d+1*id_const);
                    z_row(i) = 1;
                    Z = [Z; z_row];
                    Y = [Y; delta];

                    s = sum(S);
                    w = factorial(s) * factorial(d - s - 1) / factorial(d);
                    W = [W; w];
                end
            end
        end

        for b = 1:B
            y_b = Y(:, b);
            if id_const~=0   
                phi_b = (Z(:,2:end)' * diag(W) * Z(:,2:end)) \ (Z(:,2:end)' * diag(W) * y_b);
            else
                phi_b = (Z' * diag(W) * Z) \ (Z' * diag(W) * y_b);
            end
            phi(1+id_const:end, b) = phi_b;
        end

        out.phi = phi;
        out.err = norm(p_full - (phi0 + sum(phi,1)));
    elseif method == "knn"
        K = round(0.1*size(X_bg,1));
        M = iterations;
        distances = vecnorm(X_bg - x_target, 2, 2);
        [~, idxs] = mink(distances, K);
        X_knn = X_bg(idxs, :);
        for m = 1:M
            if isempty(id_const)
                if ~isempty(groups)
                    order = randsample(unique(groups),d);
                else
                    order = randsample(vars,d);%randperm(d);
                end
            else
                if ~isempty(groups)
                    aux   = unique(groups); 
                    order = randsample(aux(aux~=id_const),d);
                else
                    order = randsample(vars(vars~=id_const),d);%randperm(d);
                end
            end

            xj = X_knn(randi(K), :);
            xS = xj;
            for i = 1:d
                if ~isempty(groups)
                    feat  = order(i);
                    feats = logical([groups==feat]);
                    xS_feat = xS;
                    xS_feat(feats) = x_target(feats);
                else
                    feat = order(i);
                    xS_feat = xS;
                    xS_feat(feat) = x_target(feat);
                end   
                p1 = f(xS,params,bin_edges);
                p2 = f(xS_feat,params,bin_edges);
                phi(feat, :) = phi(feat, :) + (p2 - p1) / M;
                xS = xS_feat;
            end
        end
        out.phi = phi;

        for i = 1:size(X_knn,1)
            aux(i,:) = f(X_knn(i,:),params,bin_edges);
        end

        out.phi0cond = mean(aux, 1);  

        out.err = norm(p_full - mean(aux, 1));
        %out.phi0cond = mean(f(X_knn,params,bin_edges), 1);  % Local expected prediction
    end

    % Plot forecast and background density
    if isempty(verbose) || verbose == 0
        figure;
        subplot(2,1,1);
        plot(bins_center, p_full, 'b-', 'LineWidth', 2); hold on;
        if min(size(p_bg))>1
            plot(bins_center, mean(p_bg), 'k--', 'LineWidth', 2);
        else
            plot(bins_center, (p_bg), 'k--', 'LineWidth', 2);
        end
        title('Forecast vs Background Densities');
        legend('Forecast at target x', 'Background avg density');
        grid on;

        subplot(2,1,2);
        bar_data = phi';
        bar_data = bar_data(:,~sum(bar_data)==0);
        bar(bins_center, bar_data, 'stacked');
        hold on;
        if method == 'knn'
            plot(bins_center, p_full - out.phi0cond, 'ko--', 'LineWidth', 1.5);
        else
            plot(bins_center, p_full - phi0, 'ko--', 'LineWidth', 1.5);
        end
        title('SHAP values for bin probabilities');
        xlabel('Bin Center'); ylabel('Δ Probability');
        legend(arrayfun(@(i) sprintf('Feature %d', i), 1:d, 'UniformOutput', false),'Location','NorthWest');
        grid on;

    end

end