How to find the probability of the most likely configuration in statistical mechanics? 
Consider a large number $N$ of distinguishable particles distributed
  among $M$ boxes.
We know that the total number of possible microstates is
  $$\Omega=M^N$$  and that the number of microstates with a distribution
  among the boxes given by the configuration $[n_1, n_2, ..., n_M]$ is
  given by $$\frac{N!}{\prod_{j=1}^M (n_j)!}\tag{1}$$ The most likely
  configuration corresponding to the particles distributed equally among
  the $M$ boxes is trivially $$n_0=\frac{N}{M}\tag{2}$$ Now if we
  let $\Omega_0$ denote the statistical weight of this configuration and
  $p_0$ its probability.

How should I go about calculating $p_0$?

This is how I thought it should be done:
$$p_0=\frac{\left(\frac{N}{M}\right)}{\frac{N!} {\prod_{j=1}^M (n_j)!}}=\color{red}{\fbox{$\frac{N}{M}\frac{\prod_{j=1}^M (n_j)!}{N!}$}}$$
where all I did was divide $(2)$ by $(1)$ since by my logic the probability is simply the most likely configuration divided by the total number of configurations.

But apparently this is not the case and the correct answer is $$p_0=\frac{\Omega_0}{\Omega}=\frac{\frac{N!}{\Big(\left(\frac{N}{M}\right)!\Big)^M}}{M^N}=\color{#180}{\fbox{$\frac{N!}{\Big(\left(\frac{N}{M}\right)!\Big)^M M^N}$}}$$
Could someone please provide me with any hints or an explanation to justify why $$\color{#180}{p_0={\fbox{$\frac{N!}{\Big(\left(\frac{N}{M}\right)!\Big)^M M^N}$}}}$$ is the correct answer?

EDIT:
I have already been given an answer by BLAZE but I feel that there is a much easier way of calculating $p_0$. The reason I say this is because the answer to the problem was just stated as 
$$p_0=\frac{\Omega_0}{\Omega}=\frac{\frac{N!}{\Big(\left(\frac{N}{M}\right)!\Big)^M}}{M^N}={\fbox{$\frac{N!}{\Big(\left(\frac{N}{M}\right)!\Big)^M M^N}$}}$$ without any other working. 
How does one arrive at this answer without writing down any intermediate steps?
Thanks.
 A: Here we use the Multinomial Coefficients 
formula since we have $M$ boxes to place $N$ particles in. 
Suppose that we label the $M$ boxes as $M_1,M_2,...,M_{n-1},M_n$ where the $M$ boxes have been split equally amongst the particles as $\dfrac{N}{M}$ such that $M_1=M_2=M_{n-1}=M_n$ and also $\sum\limits_{i=1}^n M_i=M$. 
What we are interested in for this purpose is the $\dfrac{N}{M_i}$ factors which will become clear when you consider the following:  
$$\binom{N}{\frac{N}{M_1},\frac{N}{M_2},\frac{N}{M_{n-1}},...,\frac{N}{M_n}}=
\binom{N}{\frac{N}{M_1}}\binom{N-\frac{N}{M_1}}{\frac{N}{M_2}}\binom{N-\frac{N}{M_2}}{\frac{N}{M_{n-1}}}\cdots\binom{N-\frac{N}{M_{n-1}}}{\frac{N}{M_n}}=
\frac{N!}{(N-\frac{N}{M_1})!\cdot \frac{N}{M_1}!}\cdot\frac{(N-\frac{N}{M_1})!}{(N-\frac{N}{M_2})!\cdot \frac{N}{M_2}!}\cdot\frac{(N-\frac{N}{M_2})!}{(N-\frac{N}{M_{n-1}})!\cdot \frac{N}{M_{n-1}}!}\cdots\frac{(N-\frac{N}{M_{n-1}})!}{(N-\frac{N}{M_n})!\cdot \frac{N}{M_n}!}=
\frac{N!}{\Big(\left(\frac{N}{M}\right)!\Big)^M}$$ where in the last line mass cancellation occurs and I have used the fact that all the $M_i=M$ and that there must be $M$ lots of $\dfrac{N}{M_i}$ since $\sum\limits_{i=1}^n M_i=M$ as mentioned previously.
So the statistical weight $\Omega_0$ of the most likely configuration is $$\Omega_0=\frac{N!}{\Big(\left(\frac{N}{M}\right)!\Big)^M}$$  
But we already know that the statistical weight of the total number of configurations is $$\Omega=M^N$$ dividing $\Omega_0$ by $\Omega$ gives you the desired formula in green. 

If you are struggling to follow the general formula for the multinomial coefficients then think about this problem that I have come across:
Consider a standard pack of $52$ playing cards. The cards are distributed into $4$ piles completely randomly by tossing a four-sided die for each card.
Calculate how many microstates there are corresponding to having $13$ cards in
each pile, i.e. the most-likely configuration.
\begin{align}\binom{52}{13,13,13,13}&=\binom{52}{13}\binom{39}{13}\binom{26}{13}\binom{13}{13}\\&=\frac{52!}{(52-13)!\cdot 13!}\cdot\frac{{(52-13)!}}{(39-13)!\cdot 13!}\cdot\frac{{(39-13)!}}{(26-13)!\cdot 13!}\cdot\frac{{(26-13)!}}{(13-13)!\cdot 13!}\\&=\frac{52!}{(13!)^4}\end{align}
A: I think the thing that is confusing you here is the location of the factorial function. So in this answer I will place it inside the brackets just to avoid confusion:
Lets consider the way you were trying to calculate $p_0$ by using

$$\frac{N!}{\prod_{j=1}^M \large(n_j!)}\tag{1}$$

The probability $p_0$ must be equal to $$p_0=\frac{\Omega_0}{\Omega}=\frac{\frac{N!}{\prod_{j=1}^M \large(n_j!)}}{M^N}\tag{a}$$ as it is the ratio of no. of most likely configurations to all possible configurations (which is $M^N$).
Now simply replace $n_j$ with $$n_0=\frac{N}{M}$$ This is allowed since $n_j$ is just one of the configurations. But in our case all the $n_j$'s are the same and equal to $n_0$ since we are only interested in finding the probability of the most likely configuration. 
So the infinite product just becomes $$\prod_{j=1}^M (n_j!)=\prod_{}^M (n_0!)=\prod_{}^M \left[\left(\frac{N}{M}\right)!\right]=\left[\left(\frac{N}{M}\right)!\right]^M\tag{b}$$ 
where the last equality can be thought of as the fact that we have $M$ lots of $$\left(\frac{N}{M}\right)!$$ or explicitly: 
$$\left[\left(\frac{N}{M}\right)!\right]^M=\overbrace{\left[\left(\frac{N}{M}\right)!\right]\cdot\left[\left(\frac{N}{M}\right)!\right]\cdot\left[\left(\frac{N}{M}\right)!\right]\cdot\left[\left(\frac{N}{M}\right)!\right]\cdots}^{\Large\text{M times}}$$ 
Substitute $(\mathrm{b})$ into $(\mathrm{a})$ and the result $${\fbox{$p_0=\frac{N!}{\Big(\left(\frac{N}{M}\right)!\Big)^M M^N}$}}$$ immediately follows.
