What you say is a good idea, but is not quite correct because we must maximize the distribute subject to two constraints. I will reiterate the derivation following along with McQuarrie's Statistical Mechanics.
Constraints in the Ensemble:
The constraints I mention are actually the properties that determine the ensemble we are working with:
$$
\sum_ja_j=A
$$
$$
\sum_j\epsilon_ja_j=E
$$
The first constraint says that sum over all states of the population numbers for each state $j$, given by $a_j$, must equal the total number of members in the system $A$. This is the number of molecules if you like.
The second constraint says that the energy of each state, $\epsilon_j$, multiplied by the population of that state, $a_j$, must equal the total energy, $E$, of the system when the sum is taken over all $j$. These constraints are so intuitive physically that it is easy to forget they even have to be enforced.
By the way, we are not actually constraining the total energy. We only constrain that whatever the energy is, it can be found from the second constraint. In fact, in McQuarrie, this derivation is undertaken for the canonical (NVT) ensemble.
Method of the Most Probable Distribution:
Now, as you say, the distribution which we primarily observe will be the one which maximizes the number of ways to arrange the members of the ensemble. Notice I say primarily observe. This is discussed below. This number of ways is given by the equation you write down,
$$
W(\mathbf{a})=\frac{A!}{\prod_ka_k!}
$$
There are many possible distributions that satisfy the constraints we have given. Given this indeterminate (but probably infinite) number of ensembles which satisfy our constraints, the probability of being found in state $j$ is given by,
$$
P_j=\frac{\bar{a}_j}{A}=\frac{1}{A}\frac{\sum_{\mathbf{a}}a_j(\mathbf{a})W(\mathbf{a})}{\sum_{\mathbf{a}}W(\mathbf{a})}
$$
Notice that we are now summing over the vector of all the quantum states for each ensemble, $\mathbf{a}$, because we want to find $\bar{a}_j$, the average population of quantum state $j$ in each different ensemble satisfying our constraints.
This is clearly not very useful to us because there are so many possible ensembles we will never be able to do something with them. We can, however, take advantage of the form of $W(\mathbf{a})$. Namely, by just playing around with the numbers, you will find that the multinomial coefficients are very strongly peaked at their maximum when all of the $a_j$'s are large. Fortunately, this is exactly the regime we live in because we are thinking about something on the order of $10^{23}$ particles, so the population of states will be very large.
This means we can make an approximation, which is formally exact in the limit as $A\rightarrow\infty$, that $\bar{a}_j/A=a_j^*/A$ where $a_j^*$ is the $a_j$ in the distribution which maximizes $W(\mathbf{a})$. So, the fact we are going to maximize this distribution comes about through an approximation, albeit one which becomes exact in the limit near where we are working.
So, we can simplify our equation for the probabilities to,
$$
P_j=\frac{1}{A}\frac{a_j^*W(\mathbf{a^*})}{W(\mathbf{a^*})}=\frac{a^*_j}{A}
$$
Method of Lagrange Multiplies:
We want to find $a_j^*$ which means we must perform a maximization subject to the two constraints we listed at the beginning. Maximizing something subject to constraints does not seem to be easy, but there is a very nice method for doing it called the method of Lagrange multipliers. Essentially, we introduce some undetermined coefficients to the maximization equation so that the constraints are included, but in reality we are just adding zero to the equation.
What I mean by this is that I can rewrite the constraints as follows:
$$
\alpha\left(\sum_ja_j-A\right)=0
$$
$$
\beta\left(\sum_j\epsilon_ja_j-E\right)=0
$$
where $\alpha$ and $\beta$ are the so-called Lagrange multipliers. Then, we can simply maximize our distribution as normal:
$$
\frac{\partial}{\partial a_k}\left(\ln W(\mathbf{a^*})-\alpha\left(\sum_ja_j-A\right)-\beta\left(\sum_j\epsilon_ja_j-E\right)\right)=0
$$
It really isn't obvious that this should maximize the system while satisfying the constraints, so you will have to read elsewhere why this works, but the fact is that it does work.
Notice, I have actually played a trick. Instead of maximizing $W(\mathbf{a^*})$, I am maximizing $\ln W(\mathbf{a^*})$. This is because $W(\mathbf{a^*})$ has factorials in it, so by taking the logarithm we can employ Stirling's approximation. This is also more or less exact in the limit of a large number of particles, which we always satisfy. The reason we can do this is because the logarithm of a function will have a maximum at the same coordinate as the function itself and the logarithm is smooth such that will not somehow create another maximum or something which would spoil the trick.
After taking the derivative, we lose the constant and every term in the sum where $k\ne j$. So, we have,
$$
-\ln a_j^*-\alpha-1-\beta E_j=0
$$
Or, after rearranging and solving for $a_j^*$:
$$
a_j^*=e^{-(\alpha+1)}e^{-\beta E_j}
$$
Putting it All Together:
So, all that remains is to solve for the undetermined mutlipliers and plug it all back into our probability statement, $P_j=a_J^*/A$.
This can be done very simply (and we actually don't care about the mutliplier $\alpha$). If we take the sum over all $j$ of the equation we wrote for $a_j^*$, we realize that the left side is equal to $A$ due to our first constraint and we can rearrange to write $\alpha$ in terms of $\beta$:
$$
e^{\alpha+1}=\frac{1}{A}\sum_je^{-\beta E_j}
$$
Then, we can find the final probabilities from,
$$
P_j=\frac{a_j^*}{A}=\frac{e^{-(\alpha+1)}e^{-\beta E_j}}{A}=\frac{A}{\sum_je^{-\beta E_j}}\frac{e^{-\beta E_j}}{A}
$$
We finally arrive at the final probability distribution which defines Boltzmann statistics:
$$
P_j=\frac{e^{-\beta E_j}}{\sum_je^{-\beta E_j}}
$$
Thus, the only one of these Lagrange multipliers we care about is $\beta$. We have not actually found what this multipliers is yet, but I will stop here as it should be no surprise that the multiplier is given by $\beta=1/k_bT$.
So, I believe the above answers your question 1. That is, you have to worry about the constraints associated with you ensemble so what you mention does not quite work.
As to your second question, it is not quite true that translational energy is continuous. Pretty much everyone pretends it is, but kinetic energy is quantized. It is, however, very rare that the boundaries a particle is confined to are sufficiently small (and impermeable) that quantization of the kinetic energy is noticeable. Also, because there are many of these states which are closely spaced, one must reach very low temperatures for the quantization to be important.
What I have described above is the origin of Boltzmann statistics (not Maxwell–Boltzmann, which is often confused with it). This is a common confusion because of the famous Maxwell–Boltzmann distribution of velocities. I only bring this up because the Maxwell–Boltzmann distribution of velocities makes exactly the approximation which you describe by pretending the kinetic energies are continuous.
That is, if we simply plug the kinetic energy of a particle $E=1/2mv_x^2$ into our final result (and note the deonominator is the partition function), we get,
$$
P_{v_x}=\frac{e^{-\frac{mv_x^2}{2kT}}}{Q}
$$
Then, the partition function, $Q$, is just a normalization constant. Because the states of different kinetic energies are so closely spaced, we can elevate $P_{v_x}$ to a function, $P(v_x)$, and integrate over all velocities to determine the probability of finding a particle with a velocity between $v_x$ and $v_x+dv_x$. I will leave this as an exercise. Notice this is just an integral over a gaussian function.
