The most simple derivation I know is using dimensional analysis.
Let's assume we can work with classical physics alone and ask for the mathematical form of the black-body spectrum. It can be obtained as the em. radiation in a cavity in thermal equilibrium at a given temperature determined by heating the walls and keeping it at a certain temperature ##T## in a thermostat. Now it is clear that this problem has only the intensive parameter ##T## (temperature of the walls), because the radiation is isotropic and homogeneous (otherwise it couldn't be in thermal equilibrium). The energy density per unit volume and per unit frequency ##\omega## of the radiation is what we look for. There cannot be a chemical potential for radiation, because it doesn't carry any intrinsic ("charge like") conserved quantity.
Further some fundamental constants can be involved: The speed of light (as the fundamental constant inherent in electromagnetism; here I use CGS units, because this is the most simple way for theoretical classical electromagnetism), the Boltzmann constant ##k_B##. Thus in our "pool" of variables to express ##u## we have:
##T##: Dimension K (Kelvin)
##k_B##: Dimension erg/K=g cm^2/(s^2 K)
##\omega##: Dimension s (seconds)
##c##: Dimension cm/s
##u##: Dimension ##\text{erg} \text{s}/\text{cm}^3=\text{g}/(\text{s} \text{cm})## (energy per volume and per ##\mathrm{d} \omega##)
Since there is no dimension K anywhere except in the temperature, only the quantity ##k_B T## can enter the law. To build a dimensionless quantitiy (which must be a constant), we have to build this quantity, I'll call ##C##, from ##k_B T##, ##V##, ##\omega##, ##c##, and ##u## (and ##u## must be involved in order to get a relation for this quantity). So we can put the relation into the form
##C=u (k_B T)^{\alpha} \omega^{\gamma} c^{\delta}.##
Since g only appears in ##u## and ##T## we must have ##\alpha^{-1}## and the rest of the dimensional analysis gives uniquely ##\delta=3, \quad \gamma=-2.##
Thus we have
$$u=C \frac{k_B T \omega^2}{c^3}.$$
This gives the obviously wrong result because of the well-known Rayleigh-Jeans UV catastrophe.
The answer is that we have to consider quantum theory, and this adds the additional universal quantity ##\hbar## of dimension erg s to the game. Thus we have
$$C=u (k_B T)^{\beta} \omega^{\gamma} c^{\delta} \hbar^{\epsilon}.$$
leading to
$$[C]=\text{erg}^{\beta+\epsilon} \text{s}^{1-\gamma+\epsilon} \text{cm}^{\delta-3}.$$
This leads to
$$\delta=3, \quad \beta+\epsilon=-1, \quad \gamma-\epsilon=1.$$
This is an underdetermined set of equations. So we need to find another dimensionless constant that doesn't depend on ##u## to determine the general law
$$C'=(k_B T)^{\beta'} \omega^{\gamma'} c^{\delta'} \hbar^{\epsilon'}.$$
Then
$$u=\frac{k_B T \omega^2}{c^3} f(C'),$$
where ##f## is some arbitrary function. In order that this solves the UV problem, ##C'## must depend on ##\omega##, so that we can assume ##\gamma'=1##. The dimensional analysis then leads to
$$\epsilon=-\beta=\gamma=1.$$
Thus we have
$$u=\frac{k_B T \omega^2}{c^3} f \left (\frac{\hbar \omega}{k_B T} \right ).$$
The displacement law can now derived as follows. For the frequency, where one has maximal energy density per frequency we have
$$\partial_{\omega} u(\omega)=\frac{2 \omega k_B T}{c^3} f(\xi) + \frac{\hbar \omega^2}{c^3} f'(\xi), \quad \xi=\frac{\hbar \omega}{k_B T}.$$
Setting this to 0 leads to a necessary condition for a maximum. Thus the maximum frequency is determined by a solution of the equation
$$2f(\xi)+\xi f'(\xi)=0.$$
For any solution ##\xi_0## of this equation, one finds
$$\xi_0=\frac{\hbar \omega_m}{k_B T}.$$
So we have ##\omega_m = \xi_0 k_B T/\hbar##, which is Wien's displacement law (in the frequency domain).
The constant ##\xi_0## can only be determined when one has the unknown function ##f##. This one only gets with additional assumptions. Planck first guessed it in a clever way using the high-precision data by Rubens und Kurlbaum on the spectrum of the black-body radiation for both high and low frequencies. In an effort to find a theoretical explanation of this law, which fitted the data very accurately, Planck came to the discovery of quantum theory, because he needed the assumption that electromagnetic radiation of frequency ##\omega## is absorbed and emitted by the walls of the cavity only in portions ##n \hbar \omega##, where ##n \in \mathbb{N}##. This introduced the Planck constant as a new fundamental constant into physics, and this we used above in our dimensional analysis to solve the UV problem.
Nowadays, with the full knowledge about QED we can derive Planck's law very easily just using quantum statistics (with the em. field necessarily to be quantized as bosons due to the spin-statistics theorem).
