I repeatedly run an experiment each repetition of which results in a sample of observations that I assume follow some distribution. (It is not important here what the distribution represents.) I am interested in the shape of this distribution and would like to use the repetitions to get an idea of how much the sample distribution varies between repetitions and for how much confidence I can have in 1. the shape and 2. the location of any estimate.
I have no reason to believe that either the distribution itself is normal or that the distributions of the moments of the samples generated by the different repetitions are normal. That is, the experiment is rather explorative. Therefore I would like to avoid making any constraining assumptions.
The obvious approach would be to visualize this by plotting the average and dispersion of the ensemble of CDFs of the sample distributions. However, I have never seen such a plot anywhere although the problem appears to be quite common. Is there any reason why this approach would be a bad idea? Is there a better option?
The approach suggested above is reasonably straightforward (see below).
I do not know that the samples belong to the same theoretical distribution (say, the experiment might have a phase transition and yield a completely different distribution for some values). Does this kind of averaging / aggregating still make sense? At the very least I should see that the dispersion becomes large.
There are some similar questions about averaging histograms (here, the discussion revolves specifically about histograms and their intricacies) and averaging distributions (here, it is given that the aggregated distributions are normal and well-behaved).
E.g., in Python3 with Exponential RVS as test samples and with the interquartile range plotted as dispersion indicator (note that the variation between the samples I get from my experiment is larger):
from scipy.stats import expon import numpy as np import matplotlib.pyplot as plt samples_x =  samples_y =  for i in range(20): x = np.sort(expon.rvs(0.1, size=100)) y = (np.arange(len(x))+1) / len(x) samples_x.append(x) samples_y.append(y) flat_y = np.unique(np.sort(np.hstack(samples_y))) for i in range(len(samples_x)): x = [samples_x[i][np.argmax(samples_y[i]>=y)] if samples_y[i]<=y else 0 for y in flat_y ] samples_x[i] = x samples_y[i] = flat_y sum_x = np.vstack(samples_x) q25_x = np.quantile(sum_x,.25,axis=0) q50_x = np.quantile(sum_x,.50,axis=0) q75_x = np.quantile(sum_x,.75,axis=0) mean_x = sum_x.mean(axis=0) flat_x2575 = np.unique(np.sort(np.hstack([q25_x, q75_x]))) q25_y = [flat_y[np.argmax(q25_x>=x)] if np.sum(np.argmax(q25_x>=x))>0 else np.max(flat_y) for x in flat_x2575] q75_y = [flat_y[np.argmax(q75_x>=x)] if q75_x<=x else 0 for x in flat_x2575] fig = plt.figure() ax0 = fig.add_subplot(111) ax0.fill_between(flat_x2575, q25_y, q75_y, facecolor="C2", alpha=0.25) ax0.plot(q50_x, flat_y, color="C2") ax0.plot(mean_x, flat_y, dashes=[3, 3], color="C2") ax0.set_ylabel("CDF(x)") ax0.set_xlabel("x") plt.savefig("cdf.png", density=300)