Question: So I am taking my introductory statistical mechanics course, and this concept is something that I can not wrap my head around. My professor said that *all microstates are equally probable* and this is an axiom, always true.

However when I take a look at the canonical ensemble, I find that the probability of a given microstate is $P\_i=exp(-\beta E\_i)$, which depends on $\beta$ and the energy of the given microstate. So if I have two microstates with different $E$ energies, then they will have different $P$ probabilities! They are NOT equally probable!

Another seemingly counter-example is the harmonic oscillator. It is known that a real (as in normal mass and energy scales) harmonic diatomic molecules have $\Delta E=h \nu$ energy differences between its vibration quant, and this results in that the *vast* majority of the molecules occupy the $\nu = 0$ state. How is this possible, when all states should be equally probable? By the way, is this example really different from the canonical ensemble, or is it concepually the same?

Answer: 
In the microcanonical ensemble, all states with a *fixed energy* are equally probable. In classical mechanics, this follows essentially from ergodic theory (non-linear systems tend not to have 'many' conserved quantities). 

If you fix the energy of a 1D quantum harmonic oscillator, then you uniquely specify the state. Things are more interesting however if you have a *collection* of (weakly coupled) oscillators with fixed 'total energy' $E\approx E\_1+E\_2+\dots+E\_N=\sum\_{i=1}^N h\nu n\_i$. (This formula ignores a contribution from interactions between oscillators, which we assume is found to be small empirically). 

In this case, for a given value of $E$ the number of states is given by $\sum\_{\{n\_i\}\_{i=1}^N}\mathbb{1}(\sum\_{i=1}^Nn\_i=\frac{E}{h\nu})$. When $N$ is large, this is approximately the hyper-volume of an $N-1$ dimensional simplex. This lets you compute the entropy of a collection of the oscillators in the microcanonical ensemble, which then lets you compute the free energy per oscillator, which you can compare with the partition function of a single oscillator in the large $N$ limit. 

The probabilities you describe are derived from the *canonical ensemble*: the Boltzmann factors arise by assuming your system (a harmonic oscillator) is already weakly coupled to a thermal bath.
Many microstates $N\_i$ can have the same energy $E\_i$, many can have another energy say $E\_j$ ($N\_j$). If $N\_i$ is bigger than $N\_j$ then $E\_i$ is more probable than $E\_j$ since all microstates are equally likely in equilibrium. If its Max. Boltzman statistics then
$$\frac{P\_i}{P\_j} = \exp( (E\_j-E\_i)\cdot\beta)$$
at equilibrium.