Assume we have $N$ random variables $X_1, \ldots, X_N$. As an example, assume that these random variables describe test scores of $N$ students. I am interested in finding the distribution of average test scores within the best 20% of students, the second 20% of students, etc. via Monte Carlo simulation.
The average score within 20% groups will be surrounded by some uncertainty for two reasons: First, the test scores themselves are uncertain as they are random variables. Hence, repeated sampling of test scores for each student and averaging over students will give a distribution of average test scores. Second, it is uncertain which students belong to the respective 20% groups as the distributions of individual test scores may overlap. For simplicity, I focus on the average test score distribution within the top 20% of students in what follows.
A first naive approach may be the following: Simulate $M$ realizations of $X_1, \ldots, X_N$, sort them and average the top 20% of the realized samples, resulting in $M$ samples of the top 20% average. However, this is not what I want for the following reason. Assume all $X_i$ are distributed as $N(0,1)$, i.e., there is no systematic variation in test scores. The naive approach will result in a distribution of the top 20% average that is centered around $1.4$. This is not the desired result as all students - including the top 20% of students - have a $N(0,1)$ test score distribution. Hence, the distribution of the average within the top 20% (and all other 20% groups) is expected to be centered around zero.
I have been thinking about a two-step solution, but I doubt it's fully correct. In a first step, it is easy to compute a Monte Carlo approximation of the probability that a student belongs to the top 20% as follows. Simulate $M$ realizations of $X_1, \ldots, X_N$, sort each of them and construct a vector of binary variables $S_i$ taking the value of $1$ when a realized test score is above the 80% quantile of the distribution of the realizations. The average over these $M$ binary variables is an estimate of the probability that a student is in the top 20% group. Denote this estimate as $\pi_i$.
I am uncertain how to proceed from here. One idea is to do $L$ additional simulation iterations where I iterate between the following three steps: 1) Sample a realization of $X_1, \ldots, X_N$, 2) Sample a realization of the binary variable $S_i \sim Ber(\pi_i)$ indicating top 20% membership. 3) Average the realizations of $X_1, \ldots, X_N$ for all students where the realization of $S_i = 1$.
While this approach seems to give more useful results, I am wondering whether there is a more correct approach as I think that the sampling step of $S_i$ assumes independence of all $S_i$ which may not hold true.
