I'll write $\epsilon=\hbar\omega$ the energy of the boson. In general, it is in a mixed state described by the density matrix $\rho$ and the distribution of boson number is given by:
$$
p(n) = \langle n|\rho|n\rangle
$$
In the canonical ensemble, $\rho = \frac{1}{Z}e^{-\beta\epsilon\hat n}$ (chose a trivial zpe) with $Z = Tr (e^{-\beta\epsilon\hat n}) = \frac{1}{1-e^{-\beta\epsilon}}$ the partition function. Using the above formula, you recover a geometric distribution:
$$
p(n) = (1-e^{-\beta \epsilon})e^{-\beta \epsilon n}
$$
and it is only when taking the expected value that you obtain the Bose-Einstein distribution (with trivial chemical potential):
$$
\langle n \rangle = \frac{1}{e^{\beta\epsilon}-1}
$$
In the case of coherent/squeezed light, you have a pure state, $\rho = |\psi\rangle\langle\psi|$ and $p(n) = |\langle\psi|n\rangle|^2$ as in Born's rule. In the case of a coherent state $|\alpha\rangle$:
$$
|\alpha\rangle = \sum_{n=0}^\infty \frac{e^{-|\alpha|^2/2}\alpha^n}{\sqrt{n!}}|n\rangle \\
p(n) = \frac{e^{-|\alpha|^2}|\alpha|^{2n}}{n!}
$$
so you recover the Poisson distribution of parameter $|\alpha|^2$, so $\langle n \rangle = |\alpha|^2$.
In the case of the vacuum squeezed state $|\zeta\rangle$:
$$
|\zeta\rangle = \sum_{n=0}^\infty \frac{\sqrt{(2n)!}(-\zeta\tanh|\zeta|)^n}{n!\sqrt{\cosh|\zeta|}(2|\zeta|)^n}|2n\rangle \\
p(2n) = \frac{(2n)!\tanh^{2n}|\zeta|}{4^n n!^2\cosh|\zeta|} \\
p(2n+1) = 0
$$
which you can graph or estimate the asymptotics: $p(2n) \sim \frac{\tanh^{2n}|\zeta|}{\sqrt{\pi n} \cosh|\zeta|}$. By using the Taylor series of the square root, you can check that the distribution is normalized and calculate the expected value $\langle n\rangle$.
To obtain the other squeezed coherent states, you'll need to apply the displacement operator to this vacuum squeezed state, but then the formulas get more intricate, and it's better to just directly graph it.
Hope this helps and tell me if something is not clear.
