edit: naturally you can also use this information to calculate the number of trials needed get a desired probability that the estimation falls within certain bounds from pi.
For 10x10 (100 iterations), I got 3.44
For 100x100 (10,000 iterations), I got 3.1796
For 1000x1000 (1,000,000 iters), I got 3.14552
For 10,000x10,000 (100,000,000 iters), I got 3.141990
For 100,000x100,000 (10,000,000,000 iters, 2m06s runtime), I got 3.141633.
So the deterministic version takes a long time to converge on the right answer too. Therefore, it seems that the Monte Carlo version appears less accurate than you'd think, not because of the stability of the randomness, but because your intuition (and mine) is way off about how good any kind of sampling is at converging on the right answer.
See @teisman's great reply for figuring out how many iterations you'd need for a given level of accuracy.
Edit: Looking back at the numbers, you need 100x as many samples to roughly get an extra digit of precision. That's 10x the iterations on each axis, which makes sense that that would give you an extra digit of precision in your answer. Therefore, after 10000 (100x100) iterations, you wouldn't expect better accuracy than 2 digits, or 3.1. Which is what we see.
Edit2: I'm an idiot. After converting my test from (sqrt((x * x + y * y)/(samples * samples)) < 1) to ((x * x + y * y) < (samples * samples)) and making the maths integer only, I got the 100,000x100,000 runtime down to 10.8s. (Same answer though.)