How to see Planck's radiation law as a consequence of Bose Einstein statistics? Planck's law comes about from the following ingredients. 
1) The mode density per unit volume in a cavity is $8\pi\nu^2/c^3$.
2) Within each mode, assume Boltzmann statistics i.e the probability of having an energy $E$ is given by $p(E) \propto e^{-E/kT}$
3) Ask that the energies $E$ be discretized as $nh\nu$ instead of continuous. Essentially, the difference between doing $\int E e^{-E/kT}dE$ and the sum $\sum_n E_n e^{-E_n/kT}$, where $E_n = nh\nu$ gives Planck's law.
$$I = \frac{8\pi\nu^2}{c^3}\frac{h\nu}{e^{h\nu/kT} - 1} \, .$$
So far, it's a nice trick to go to the discrete sum but that's all that was done. There is no sense in which the particles are indistinguishable etc. 
What is the correct modern way of seeing this result? I assume the starting point is that we have an ideal gas of photons inside a blackbody cavity and they follow statistics that yield $p(E) =\frac{1}{e^{E/kT} - 1}$ but how exactly does one go from there to Planck's law?
 A: As you said, we consider a cavity as a model for the blackbody radiation for which you already have found the mode density
$$
g(\nu) = \frac{8 \pi \nu^2}{c^3}.
$$
The energies the photons can posses are indeed discretized and have energy $E = nh \nu$ where $n$ is an integer. We will now write down the grand-canonical partition function for bosons with vanishing chemical potential. It's taken zero because the photons are interacting with the walls of the cavity and so the number of photons is not constant. The partition function for a particular frequency $\nu_0$ is
$$
\Xi_0 = \sum_{n=0}^\infty e^{-\beta h\nu_0 n} = \frac{1}{1-e^{-\beta h\nu_0}}.
$$
The total partition function can then be found by multiplying all partition functions of each frequency with each other
$$
\Xi = \prod_i\frac{1}{1-e^{-\beta h\nu_i}}.
$$
Taking the logarithm of $\Xi$ will transform the product in a sum
$$
\log(\Xi) = -\sum_i \log(1-e^{-\beta h\nu_i}).
$$
which can be changed in an integral assuming that the difference in frequency becomes infinitesimal between two adjacent frequencies. Instead of then counting the number of photons per frequency, we can leave this job to the density of modes $g(\nu)$:
$$
\log(\Xi) = - \int_0^\infty d\nu g(\nu) \log (1-e^{-\beta h\nu_i}).
$$
You can now calculate the energy associated using $E = -\frac{\partial \log(\Xi)}{\partial \beta}$.
So
$$
E =\int_0^\infty d\nu  \frac{g(\nu) h \nu}{ (e^{-\beta h\nu}-1)}
$$
and from the relation $E = \int d\nu I(\nu)$ you find that
$$
I = \frac{8\pi\nu^3h}{c^3} \frac{1}{e^{-\beta h\nu}-1}
$$
which is Planck's Law. Pfew!
A: 
There is no sense in which the particles are indistinguishable etc. 

Actually, indistinguishability (specifically Bose Einstein statistics) is tacitly assumed when the average energy is written as $\propto\sum_n E_n e^{-E_n/kT}$ with $E_n=nh\nu$, for $n\in\{0,1,2,...\}$. In the canonical ensemble, the probability associated with the $j$-th microstate is $p_j\propto e^{-E_j/kT}$. This is true for both quantum or classical models, both for distinguishable particles and indisinguishable particles. In a model with distinguishable particles, by definition, different permutations of those particles would be counted as different microststates, so the average energy would be given by the disastrous sum $\propto\sum_n n!\,E_n e^{-E_n/kT}$ with $E_n=nh\nu$. But when we write the average energy as $\propto\sum_n E_n e^{-E_n/kT}$, we are assuming that there is only one microstate that has $n$ photons with equal energies $h\nu$. That's what we mean by "indisinguishable": the number of photons matters, but permuting them doesn't give anything different. (As in the OP, I'm only considering a single frequency and single polarization here.)
