Say I have a satellite that's flying through the atmosphere, over multiple orbits, sampling its density at different altitudes, at say 1 measurement per second (specific numbers are irrelevant). The density may vary from orbit to orbit due to natural variations (like temperature, winds, solar illumination,etc.) and orbital variations (spatial location etc.).

Now I want to get the median altitude density profile of the atmosphere, so I calculate the median density in altitude bins of, say, 5km. Then I want to get a 95% percent confidence interval for this median. My problem is: there may be several consecutive measurements points from one orbit in one altitude bin, meaning that these datapoints will be auto-correlated, violating the independence assumption made in most CI estimation methods, and thus overestimating the precision of the median determination. Is this correct?

Therefore, would it be a good approach to use a simple bootstrap method to estimate the CI, but resampling in samples of the size of the amount of separate/unique orbits (and thus independent measurements) per altitude bin? Or would I still be making an error since I might still be sampling some of these non-independent datapoints?


1 Answer 1


One way might be to calculate the statistical inefficiency of the measurements and then only sample at intervals of that value. This is implemented in the Chodera lab 'PyMBAR' software, which has a function subsampleCorrelatedData that returns the indices of the independent measurements. I'll copy a shorter version below that goes directly to your question (sorry if you prefer to use R!).

The reference for their use of the technique is:

Chodera, John D., et al. "Use of the weighted histogram analysis method for the analysis of simulated and parallel tempering simulations." Journal of Chemical Theory and Computation 3.1 (2007): 26-41.

Now some code:

 #this is the Chodera lab 'statistical inefficiency' measure. 
import statsmodels.api as sm
import numpy as np

def statistical_inefficiency(corr, mintime=3):
    N = corr.size
    C_t = sm.tsa.stattools.acf(corr, fft=True, unbiased=True, nlags=N) #autocorrelation function
    t_grid = np.arange(N).astype('float')
    g_t = 2.0 * C_t * (1.0 - t_grid / float(N))
    ind = np.where((C_t <= 0) & (t_grid > mintime))[0][0]
    g = 1.0 + g_t[1:ind].sum()
    return max(1.0, g)

my_correlated_density_measurements = np.load('my_data.npy')#load satellite data as np.array
stat_ineff = statistical_inefficiency(my_correlated_density_measurements)
#make it an integer stride:
strd = rint(stat_ineff)
uncorrelated_data = my_correlated_density_measurements[::strd]

This can also be thought of as the autocorrelation time. Then you are free to perform bootstrapping without any adjustments on the uncorrelated_data.

I have no affiliation with the pymbar software so if I misinterpret the function, apologies and just use PyMBAR directly, but it's worked well for me and agrees with about four other techniques. Cheers!

  • $\begingroup$ Hmm, need some time to figure out if that's what I need ... $\endgroup$
    – Lu Kas
    May 3, 2020 at 13:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.