We can solve for the stationary states of a quantum harmonic oscillator denoted by $|n\rangle$ with energy eigenvalues $(n+\frac{1}2)\hbar\omega$. However if our system is in a stationary state, the time evolution is given by $|\psi,t\rangle=e^{\frac{-iE_nt}{\hbar}}|\psi,0\rangle$ which translates to $|\psi,t\rangle=e^{-i(n+\frac{1}2)\omega t}|n\rangle$ here. Then $\langle x \rangle=\langle\psi|x|\psi\rangle=0$ (because the time dependence drops out from the conjugation and we can show that $\langle x\rangle=0$ for all stationary states), and consequently we do not retrieve the classical sinusoidal motion.
However it can be shown that if we have $|\psi,t\rangle=\displaystyle\sum_{n}a_ne^{-(n+\frac{1}2)\omega t}|n\rangle$ such that more than one $a_n$ is non-negligible, we obtain $\langle x\rangle=\displaystyle\sum_{n}X_n\cos(\omega t+\phi_n)$ which is the classical sinusoidal motion we expect.
So I conclude that for classical motion we have to be in a superposition of stationary states. However, I have seen many arguments where the classical limit of the quantum harmonic oscillator is considered simply by looking at the form of the probability distribution for large $n$, for example here http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/hosc6.html (check out the graphs towards the bottom). This seems to be contradicting my above analysis of the situation - if we have to be in a superposition of stationary states to retrieve classical motion, why bother considering the large $n$ limit of stationary states - surely this information is useless in the classical limit because the classical limit never involves simply one stationary state - if it did (as I worked out above) there would be no sinusoidal motion.
Thanks for any help :)