You can prove it for any state. Let ##\hat{A}## and ##\hat{B}## be the self-adjoint operators representing observables and assume for simplicity that ##\langle A \rangle=\langle B \rangle=0##, where ##\langle \hat{A} \rangle=\mathrm{Tr}(\hat \rho \hat{A})##. Since ##\hat{\rho}## is a positive semi-definite self-adjoint operator there's a complete orthonormalized eigenbasis ##|u_n \rangle##, i.e.,
$$\hat{\rho}=\sum_n p_n |u_n \rangle \langle u_n |.$$
Thus any expectation value reads
$$\langle A \rangle = \mathrm{Tr} (\hat{\rho} \hat{A}) = \sum_{n} \langle u_n| \rho_n \hat{A} |u_n \rangle = \sum_n P_n \langle u_n |\hat{A}| u_n \rangle.$$
Now apply this to the self-adjoint operator
$$\hat{C}=(\hat{A}-\mathrm{i} \lambda \hat{B})(\hat{A}+\mathrm{i} \lambda \hat{B}), \quad \lambda \in \mathbb{R}.$$
First of all
$$\langle u_n|\hat{C}| u_n \rangle=\langle (\hat{A}+\mathrm{i} \lambda \hat{B}) u_n|(\hat{A}+\mathrm{i} \lambda \hat{B}) u_n \rangle \geq 0.$$
And since ##P_n \geq 0## for all ##n## you also have
$$P(\lambda)=\langle (\hat{A}-\mathrm{i} \lambda \hat{B})(\hat{A}+\mathrm{i} \lambda \hat{B}) \rangle \geq 0.$$
Multiplying this out gives
$$P(\lambda) = \langle A^2 \rangle + \langle \mathrm{i} [\hat{A},\hat{B}] \rangle \lambda + \lambda^2 \langle \hat{B}^2 \rangle.$$
Now ##\langle A^2 \rangle=\Delta A^2## and ##\langle B^2 \rangle=\Delta B^2## are the variances of ##A## and ##B## (because we assumed ##\langle A \rangle=\langle B \rangle=0##), i.e.,
$$P(\lambda)=\lambda^2 \Delta B^2 + \lambda \langle \mathrm{i} [\hat{A},\hat{B}]\rangle + \Delta A^2 \geq 0$$
für ##\lambda \geq 0##. Now ##P(\lambda)## is a real quadratic polynomial and thus its discriminant fulfills,
$$\frac{1}{4} \langle \mathrm{i} [\hat{A},\hat{B}] \rangle^2 - \Delta A^2 \Delta B^2 \leq 0,$$
and from this finally
$$\Delta A \Delta B \geq \frac{1}{2} |\langle \mathrm{i} [\hat{A},\hat{B}] \rangle|,$$
which is the Heisenberg-Robertson uncertainty relation.