function unc_pmy = bootstrap_estimate_mc_unc(all_weights, prior_topology, n)
n_top = numel(all_weights);
n_bridges = zeros(n_top,1);
n_samples = cell(n_top,1);

for i_top = 1:n_top
    n_bridges(i_top) = numel(all_weights{i_top});
    n_samples{i_top} = zeros(n_bridges(i_top),1);
    for i_bridge = 1:n_bridges(i_top)
        n_samples{i_top}(i_bridge) = numel(all_weights{i_top}{i_bridge});
    end
end

samples = zeros(n_top,n);
for i = 1:n
    % Resample
    weights = cell(n_top,1);
    for i_top = 1:n_top
        weights{i_top} = cell(n_bridges(i_top),1);
        for i_bridge = 1:n_bridges(i_top)
            weights{i_top}{i_bridge} = randsample(all_weights{i_top}{i_bridge}, n_samples{i_top}(i_bridge), true);
        end
    end
    
    % Compute mean and uncertainty weight at each bridge point
    mean_weights = cell(n_top,1);
    unc_mean_weights = cell(n_top,1);
    for i_top = 1:n_top
        mean_weights{i_top} = zeros(n_bridges(i_top)-1,1);
        unc_mean_weights{i_top}  = nan(n_bridges(i_top)-1,1);
        for i_bridge = 1:numel(mean_weights{i_top})
            mean_weights{i_top}(i_bridge) = mean(weights{i_top}{i_bridge});
            unc_mean_weights{i_top}(i_bridge) = sqrt(sum((weights{i_top}{i_bridge} - mean_weights{i_top}(i_bridge)).^2) ./ (numel(weights{i_top}{i_bridge}) - 1));
        end
    end
    
    % Compute total weight
    log_total_weights = zeros(n_top,1);
    for i_top = 1:n_top
        log_total_weights(i_top) = sum(log(mean_weights{i_top}));
    end
    
    % Convert to topology probability
    % Include topology prior
    log_pmy = log_total_weights + log(prior_topology);
    
    % Rescale log(pmy) before returning to real space
    log_pmy = log_pmy - max(log_pmy);
    
    % Return to regular space
    pmy = exp(log_pmy);
    
    % Normalize
    samples(:,i) = pmy ./ sum(pmy);
end

unc_pmy = std(samples, 0, 2);
