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);