I was recently studying Blackbody Radiation and the principle behind it (as far as used in the Plancks original paper) is to find the energy distribution which maximizes the number of ways in which energy quanta can be distributed.
Since the radiation formula is derived due to quantization + maximum entropy reasons, can it also be derived using a principle of least action?
In so far as I'm aware, there is not direct connection between the two as you might like, however there are a few things that might be worth mentioning. I will take this back a few steps from the question itself because I feel there's some important context to the point I would like to make in the end.
Since the problem at hand is looking for a steady-state distribution, the Liouville theorem (or its equivalent in a quantum system) implies that the probability distribution we seek must (Poisson) commute with the Hamiltonian, and hence can only be a function on the conserved quantities of the system. A key fact that's often overlooked is that there's more than one way to do this. One would be to assert that the distribution is simply a constant. This would give the distribution associated to the microcanonical ensemble. Perhaps more common is the assumption that the distribution depends only on the energy (Hamiltonian).
Now, another important fact to understanding things lies in probability theory itself. Entropy maximization (perhaps the best understanding of what this should really mean to us actually comes, I think, from information theory and signal processing -- the application Shannon originally invented the notion of entropy to deal with. See for example Cover and Thomas' book), is a fairly good principle all things considered, but alone this doesn't do much for us, we need constraints as well.
In general, we might say that we know something about a system, and it might be convenient to assert that this something takes the form of an expectation value, which we assume to be a fixed, known quantity. For a nice discussion of how entropy maximization works in the presence of constraints, see here. Roughly the same arguments work out in the quantum case as well. If you were to consider the specific question: "What is the maximum entropy distribution such that the expectation of the energy is held fixed," you would find the unnormalized answer to be the Boltzmann factor, $e^{-\beta H}$ where $\beta$ is the Lagrange multiplier enforcing the constraint $E=\langle H\rangle$ which we identify with the inverse temperature. In the quantum case, we find the density matrix to be the Boltzmann factor, but with the Hamiltonian operator rather than function.
Now, in the quantum case, the expectation of any operator $\mathcal{O}$ is given by $\langle \mathcal{O}\rangle=\text{tr}(\mathcal{O}\rho)$ where $\rho=Ne^{-\beta H}$ is the density matrix ($N$ being some normalization). This trace can be identified (under certain technical conditions) with a path integral in which "time" has been replaced by a complex time $\tau$ which is periodic with period $i\beta$ (up to minus signs). I believe I recall a nice discussion of this appearing in the conformal field theory book by Di Francesco, Mathieu, and Senechal. Under certain other technical conditions, the path integral may be expressed in terms of the classical action as $$ \langle\mathcal{O}\rangle=\int\mathcal{D}[\phi]\mathcal{O}e^{-\beta S}. $$
Finally now we can come to the principle of least action. Integrals of this form can often be approximated by asymptotic series. For example, see here. As an aside, this is the approximation working in the background any time you see Feynman diagrams floating about.
The lowest order term in this sort of approximation scheme is proportional to $e^{-\beta S[\phi_c]}$ where $\phi_c$ is the classical solution. That is, the solution which extremizes the action.
This may feel like an unsatisfactory connection between stationary action principles and thermal physics, but it is nonetheless the only connection I can think of which is technically accurate.
