We are going to calculate from among all the possible macrostates of the system the one that has the largest number of microstates. Thus, the calculation of any thermodynamic quantity is carried out neglecting all the macrostates except the most probable. Mathematically we must calculate $N_{1},N_{2},.....,N_{r}$ so that W takes the maximum value. Where W is given by:

$$W=N!\prod_{k}\frac{g_{k}^{N_{k}}}{N_{k}!} \label{ec3}$$

Taking neperians at (\ref{ec3})

$$lnW=lnN!+\sum\left(N_{k}lng_{k}-lnN_{k}!\right) \label{ec4}$$

Applying in (~\ref{ec4}) the Stirling approximation $lnN!=NlnN-N$:

$$lnW=NlnN-N+\sum_{k}N_{k}lng_{k} -\sum_{k}\left(N_{k}lnN_{k}-N_{k}\right) \label{ec5}$$

Since $\sum_{k}N_{k}=N$ , the equation (\ref{ec5}) can be simplified to

$$lnW=NlnN+\sum_{k}N_{k}lng_{k}-\sum_{k}N_{k}lnN_{k} \label{ ec6}$$

To maximize $lnW$ we apply the Lagrange indeterminate multipliers technique, imposing the conditions: $\sum_{k}N_{k}=N$ and $\sum_{k}N_{k} \epsilon_{k}=E$ We build the Lagrangian from the function to be maximized ($lnW$), and from the conditions that must be fulfilled multiplied by the coefficients $\alpha$ and $\beta$.

$$L=lnW-\alpha\sum_{k}N_{k}-\beta\sum_{k}N_{k}\epsilon_{k} \label{ec7}$$

Differentiating the Lagrangian with respect to $N_{i}$ (number of particles in state $i$ for the most probable macrostate) and equating to zero, we obtain:

$$\frac{\partial{L}}{\partial{ N_{i}}}=\frac{\partial{lnW}}{\partial{N_{i}}}-\alpha-\beta{\epsilon}_{i}=0 \label{ec8}$$

$\frac{\partial{lnW}}{\partial{N_{i}}}$ is obtained by deriving from equation (\ref{ec6})

$$\frac{\partial{lnW}} {\partial{N_{i}}}=lng_{i}-lnN_{i}-1 \label{ec9}$$

Substituting (~\ref{ec9}) into (~\ref{ec8}) it is obtained:

$$lng_{i}-lnN_{i}-1-\alpha-\beta{\epsilon_{i}}=0 \label{ec10}$$

Applying natural properties and grouping terms:

$$ln\frac{N_{i}}{g_{i}}=-(1+\alpha)-\beta{\epsilon_{i}} \label{ec11}$$

Solving for $\frac{N_{i}}{g_{i}}$

$$\frac{N_{i}}{g_{i}}=e^{-(1+\alpha)}e^{-\beta{\epsilon_{i}}}\label{ec12}$$

Taking summations in (~\ref{ec12}):

$$\sum_{i}{N_{i}}=e^{-\left(\alpha+1\right)}\sum_{i}g_{i}e^{-\beta{\epsilon_{i}}} \label{ec13}$$

As $\sum_{i}{N_{i}}= N$.

$$N=e^{-\left(\alpha+1\right)}\sum_{i}g_{i}e^{-\beta{\epsilon_{i}}} \label{ec14}$$

Where $q=\sum_{i}g_{i}e^{-\beta{\epsilon_{i}}}$, is the {\bf partition function of the system}. Therefore:

$$e^{-(\alpha{+1})}=\frac{N}{q} \label{ec15}$$

Substituting the equation (~\ref{ec15} ) in (\ref{ec12}) and solving for $N_{i}$ gives us:

$$N_{i}=N\frac{g_{i}e^{-\beta{\epsilon_{i }}}}{q} \label{ec16}$$

The values of $N_{i}$ given by the equation (\ref{ec16}) are those that maximize $W$. That is, the calculated distribution is the one with the largest number of microstates.