# Unbiased estimator of variation of median in spatial bins using bootstrap method

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?

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)

#usage: