Which ones of $n$ random variables have the largest mean (non-parametric way)? Let us have $n$ random, mutually independent variables $X_1,X_2,\dots,X_n$. Let us have some samples of them such as $x_{i,j}$ where $i=1,\dots,n$. I want to know the maximal variable(s) based on those data. These data are results from real-world experiments with $n$ different settings and the only thing I need is to determine the winning settings.
More specifically, I want to find such a subset of indices $A\equiv\{i_1,i_2,\dots,i_m\}\subseteq B\equiv \{1,\dots,n\}$ such that


*

*$\forall i,i'\in A:\mathcal{E}[X_i]\simeq\mathcal{E}[X_{i'}]$

*and $\forall i,\in A,i'\in B\backslash A:\mathcal{E}[X_i]\gtrsim\mathcal{E}[X_{i'}]$


Where $\mathcal{E}[\cdot]$ stands for the expected value and $\simeq,\gtrsim$ shall be read as equal/less than with respect to data available.
My attempt: 


*

*To set $A\gets B$

*To compare each $X_i$ and $X_{i'}$ using a test. Whenever I reject the hypothesis $\mathcal{E}[X_i]= \mathcal{E}[X_{i'}]$ in favor of $\mathcal{E}[X_i]< \mathcal{E}[X_{i'}]$, set $A\gets A \backslash \{i\}$, i.e. remove $i$ from $A$.


A simple Matlab example follows:
Ndata = 50;
X(:,1) = randn(Ndata,1);
X(:,2) = randn(Ndata,1)+0.1;
X(:,3) = randn(Ndata,1)+0.5;
X(:,4) = 3*randn(Ndata,1)+0.5;
n = 4;
A = 1:n;
for i=1:n
    for iPrime=1:n
        if ttest2(X(:,i),X(:,iPrime),'Alpha',0.05,'Tail','left')>0
            A=setdiff(A,i);
        end
    end
end

In that case the result tends to be A = [3 4].
Would you recommended something more tailored or established for this task? Possibly non-parametric and in Matlab.
I noticed that there is a related question (Ranking randomly distributed variables based on decreasingness) but it speaks about something else. It has perhaps something to do with Bonferoni corrections, but I am not sure about that point.
 A: If you're prepared to assume that any distributional differences result from one variable being stochastically larger than the other (e.g. for continuous r.v.s, $P(X>Y)>\frac12$), a Wilcoxon-Mann-Whitney statistic is a linear function of a natural estimator of $P(X>Y)$ (whether you're interested in a test or not, the estimator itself can be informative).
With more than two groups, the closest equivalent is the Kruskal-Wallis test, which is a form of rank-based one-way analysis of variance.
If you're purely interested in differences in mean, you can just look at differences in mean -- this can be turned into a nonparametric test by using difference in sample means as a test statistic in a permutation test, (or weighted sums of squared differences of group means* from the overall mean with more than two groups).
One advantage of looking at means is that mean differences are necessarily transitive, while $P(X>Y)>\frac12$ effects are not. However, once we look at statistical significance, you can get what seem to be unintuitive outcomes even with means. 
* in effect the numerator of the usual F-statistic (or indeed the F-statistic itself could be used in a permutation test)
