Given $n$ samples of a function $f : D\to\mathbb C$ in random points drawn with respect to a measure $\nu$ we develop theoretical analysis of the $L_2(D, \mu)$-approximation error. We show that the weighted least squares method from finite dimensional function spaces $V_m$, $\dim(V_m) = m < \infty$ is stable and optimal up to a multiplicative constant when given exact samples with logarithmic oversampling. Further, for noisy samples, our bounds describe the bias-variance trade off depending on the dimension $m$ of the approximation space $V_m$. All results hold with high probability. For demonstration, we consider functions defined on the $d$-dimensional cube given in unifom random samples. We analyze polynomials, the half-perid cosine, and a bounded orthonormal basis of the non-periodic Sobolev space $H_{\text{mix}}^2$. Overcoming numerical issues of this $H_{\text{mix}}^2$ basis, this gives a novel stable approximation method with quadratic error decay. Numerical experiments indicate the applicability of our results.