For a 2-level system (or an N-level system for small N) the answer is easy. Observe enough independent instances of measurements of $N^2$ linearly independent observables $A_k$ (including a constant observable whose expectation one gets for free) to get reliable estimates of their expectations. Then one can reconstruct the coefficients of the density matrix $\rho$ by solving a linear system derived from $\langle A_k\rangle = trace(A \rho)$. This completely determines the state to an accuracy depending on the number of repetitions.

Then compute a singular value decomposition of $\rho$ and check whether all but one singular value are zero within the statistical uncertainty. If not, the system is in a mixed state. If there is only one significant eigenvalue, the state can be approximated to the appropriate significance level by a pure state.

One cannot really decide whether the system is exactly in a pure state (except sometimes by looking at the generating equipment) , as this requires infinite accuracy - arbitrarily close to any pure state there are always mixed states. Indeed, unless the preparation of a physical state is specifically designed to produce exact pure states (which is the case, e.g., in a Stern-Gerlach experiment), a system is with probability one in a mixed state since the pure states have measure zero in the set of all mixed states.

In practice, pure states are just a convenient idealization for mixed states in which a rank one contribution dominates the state.