[Python] How to estimate a reliable average using bootstrapping with weighted data and...

I have a small dataset where each data point (y_i) comes with its own measurement uncertainty (e_i). I want to estimate a reliable average value for (y) while properly considering the uncertainties (e_i) as weights. My goal is to account for both:

1. The statistical variation in the dataset (via bootstrapping).
2. The measurement uncertainty of each data point (via weights based on (e_i)).

Currently, I am using a two-stage bootstrap method:

1. Resample the data indices with replacement (nonparametric bootstrap).
2. For each resampled data point ((y_i, e_i)), generate a synthetic draw from a Normal distribution (\mathcal{N}(y_i, e_i^2)).
3. Compute a weighted mean for each bootstrap sample using (w_i = 1 / e_i) as weights.

The weighted mean from the bootstrap iterations is then used to estimate the final mean and its uncertainty. Below is a minimal reproducible example of my implementation:


Code (MRE):

import numpy as np
import pandas as pd
from joblib import Parallel, delayed

def two_stage_bootstrap_weighted_mean(
data,
slope_col = 'plaw', # Column with slope estimates y_i
err_col = 'err', # Column with 1-sigma measurement errors e_i
n_bootstrap = 1000,
n_jobs = 1, # Set >1 for parallel
random_state=None
):
"""
Performs a two-stage bootstrap:
1) Resample indices with replacement (nonparametric).
2) For each resampled data point, draw from Normal(y_i, e_i) (parametric).
Then computes the weighted mean with weights = 1/e_i for each resample.

Returns
-------
mean_est : float
The average of all bootstrap weighted means.
std_boot : float
The standard deviation of the bootstrap means.
boot_means: np.ndarray
The array of all bootstrap mean estimates.
"""
rng = np.random.default_rng(random_state)

y_vals = data[slope_col].values
e_vals = data[err_col].values

# Weights = 1/e_i
N = len(data)

def bootstrap_iteration(_):
# Nonparametric resample
indices = rng.integers(0, N, size=N)
y_draw = y_vals[indices]
e_draw = e_vals[indices]

# Parametric resample from each point
y_synth = rng.normal(loc=y_draw, scale=e_draw)

# Weighted mean: sum(y_synth * w) / sum(w)
w_draw = 1.0 / e_draw
w_sum = np.sum(w_draw)
if w_sum == 0:
return np.nan
return np.sum(y_synth * w_draw) / w_sum

# Run bootstrap (parallel or single core)
boot_means = Parallel(n_jobs=n_jobs)(
delayed(bootstrap_iteration)(_) for _ in range(n_bootstrap)
)
boot_means = np.array(boot_means)

mean_est = np.nanmean(boot_means)
std_boot = np.nanstd(boot_means, ddof=1)

return mean_est, std_boot, boot_means


if __name__ == "__main__":
# Generate synthetic data
np.random.seed(42)
N = 20
true_slope = -2.0
measurement_err = np.random.uniform(0.05, 0.4, N)

# Observed slopes are true_slope + random noise of scale measurement_err
observed_slopes = np.random.normal(loc=true_slope, scale=measurement_err)

# Put into a DataFrame
df = pd.DataFrame({
"plaw": observed_slopes,
"err": measurement_err
})

# Run the two-stage bootstrap
mean_est, std_boot, boot_means = two_stage_bootstrap_weighted_mean(
df, slope_col="plaw", err_col="err", n_bootstrap=1000, n_jobs=1, random_state=123
)

print(f"True slope: {true_slope}")
print(f"Bootstrap mean slope: {mean_est:.3f}")
print(f"Bootstrap std: {std_boot:.3f}")


Questions:

1. Does this two-stage bootstrap approach make sense for combining statistical variability and individual measurement uncertainties into a reliable mean estimate?
2. Are the weights (w_i = 1 / e_i) appropriate, or should I use (w_i = 1 / e_i^2)?
3. Is there a more standard or efficient method (e.g., weighted least squares or error propagation) for estimating the mean and its uncertainty in this scenario?

Thank you for any guidance or recommendations on improving this method or using alternative approaches!

Continue reading...