Suppose I have three independent groups, with mean $\mu_1,~ \mu_2,~\mu_3$ respectively.
How can I test whether $\mu_1 < \mu_2 <\mu_3$ or not using $n_1,~n_2,~n_3$ samples from each group?
I wish to know some general methodology, not detailed calculation. I couldn't figure out how to set my hypothesis $H_0$ and $H_1$.
In statistics you cannot test whether "X is true or not". You can only try to find evidence that a null hypothesis is false.
Let's say your null hypothesis is $$ H_0^1: \mu_1 < \mu_2 < \mu_3. $$ Let's also assume that you have a way of estimating the vector $\mu = (\mu_1, \mu_2, \mu_3)'$. To keep things simply assume that you have an estimator $$ x \sim N(\mu, \Sigma), $$ where $\Sigma$ is $3 \times 3$ covariate matrix. We can rewrite the null hypothesis as $$ A \mu < 0, $$ where $$ A = \begin{bmatrix} 1 & - 1 & 0 \\ 0 & 1 & - 1 \end{bmatrix}. $$ This shows that your null hypothesis can be expressed as an inequality restriction on the vector $A \mu$. A natural estimator of $A\mu$ is given by $$ A x \sim N(A\mu, A\Sigma A'). $$ You can now use the framework for testing inequality constraint on normal vectors given in:
Kudo, Akio (1963). “A multivariate analogue of the one-sided test”. In: Biometrika 50.3/4, pp. 403–418.
This test will also work if the normality assumption holds only approximately ("asymptotically"). For example, it will work if you can draw sample means from the groups. If you draw samples of size $n_1, n_2, n_3$ and if you can draw independently from the groups then $\Sigma$ is a diagonal matrix with diagonal $$ (\sigma_1^2/n_1, \sigma_2^2/n_2, \sigma_3^2/n_3)', $$ where $\sigma_k^2$ is the variance in group $k = 1, 2, 3$. In an application, you can use sample variance instead of the unknown population variance without changing the properties of the test.
If on the other hand your alternative hypothesis is $$ H_1^2: \mu_1 < \mu_2 < \mu_3 $$ then your null hypothesis becomes $$ H_0^2: \text{NOT $H_1$}. $$ This isn't very operational. Remember that our new alternative hypothesis can be written as $H_1: A\mu <0$ so that $$ H_0^2: \text{there exists a $k=1,2$ such that $(A\mu)_k \geq 0$}. $$ I don't know if there exists any specialized test for this, but you can definitely try some strategy based on successive testing. Remember that you try to find evidence against the null. So you may first test $$ H_{0,1}^2: (A\mu)_1 \geq 0. $$ and then $$ H_{0,2}^2: (A\mu)_2 \geq 0. $$ If you reject both times then you have found evidence that $H_0$ is false and you reject $H_0$. If you don't, then you don't reject $H_0$. Since you are testing multiple times you have to adjust the nominal level of the subtest. You can use a Bonferroni correction or figure out an exact correction (since you know $\Sigma$).
Another way of constructing a test for $H_0^2$ is to note that $$ H_0^2: \max_{k=1,2} (A\mu)_k \geq 0. $$ This implies using $\max Ax$ as a test statistic. The test will have a non-standard distribution under the null, but the appropriate critical value should still be fairly easy to compute.
The answer provided by @andreas-dzemski is correct only if we know that the data is normally distributed.
If we do not know the distribution, I believe it would be better to run a nonparametric test. In this case, the simplest seems to run a permutation test. This is a book about the topic and this is a nice online explanation. Below I include R code to compute this test.
# some test data
D <- data.frame(group1=c(3,6,2,2,3,9,3,4,2,5), group2=c(5,3,10,1,10,2,4,4,2,2), group3=c(8,0,1,5,10,7,3,4,8,1))
# sample with replacement
resample <- function(X) sample(X, replace=TRUE)
# return true if mu1 < mu2 < mu3
test <- function(mu1, mu2, mu3) (mu1 < mu2) & (mu2 < mu3)
# resampling test that returns the probability of observing the relationship
mean(replicate(1000, test(mean(resample(D$group1)), mean(resample(D$group2)), mean(resample(D$group3)))))
