Let $X, Y$ be vectors of length $n$ such that each element is an iid Uniform random variable on $[0,1]$.
Define $Z$ by
$$Z = \sum_{i=1}^{n}\varphi(i)$$
where $$\varphi(i) = \begin{cases} X_{i}^{2} + Y_{i}^{2} & \text{if $X_{i}^{2} + Y_{i}^{2}<1$}\\ 0 &\text{otherwise}. \end{cases}$$
It seems that (I don't think its necessary to put the expectation here, but I sort of used it later on when describing equation (1).)
$$\mathbb{E}[\lim_{n\rightarrow \infty}8Z/n] = \pi$$
but I am not sure why, or even if my conjecture actually holds. (I feel like there is a painfully obvious geometric intuition that is just going right over my head.)
I haven't been able to get very far into the problem. But below is my best shot at it. As I progress through this work, I get the impression that I am writing increasingly dubious equations, but this is the closest I could get to anything even containing a constant multiple of $\pi$.
It looks to me that $\mathbb{P}[\varphi(i)\neq 0] = \pi/4$ and that this is a Bernoulli random variable.
From there, I could see that if $\varphi(i)$ were instead defined by
$$\bar{\varphi}(i) = \begin{cases} 1 & \text{if $X_{i}^{2} + Y_{i}^{2}<1$}\\ 0 &\text{otherwise}. \end{cases}$$
I think we get that
\begin{equation} \mathbb{E}[Z] = \mathbb{E}[n \bar{\varphi}(1)] = n\pi/4, \quad (1) \end{equation}
but I have no idea how to handle the sum of $X^{2}_{i} + Y_{i}^{2}$'s here. Is there a reasonable geometric intuition for why $Z/n$ seems to converge to $\pi/8$?
If I naively try to write
$$\mathbb{E}[Z] = \mathbb{E}[n\varphi(1)] = \mathbb{E}[n(X^{2} + Y^{2})]\mathbb{P}[\varphi(i)\neq 0]$$
Then since $\mathbb{E}[X^{2}] - \mathbb{E}[X]^{2} = 1/12$, $\mathbb{E}[X^{2}] = (1/12) + (1/4) = 1/3.$
Then we would obtain
$$\mathbb{E}[n\varphi(1)] = (\mathbb{E}[(X^{2}]+ \mathbb{E}[Y^{2}])\frac{n\pi}{4} = \frac{2}{3} \frac{n\pi}{4} = \frac{n\pi}{6}$$
which does not seem to hold in my numerical experiments.
What you did is just a Monte Carlo scheme to compute the integral
$$\int_{D} x^2+y^2 dxdy$$
with $D=\{(x,y)| x^2+y^2 < 1\}$, which defines a quarter of a disk with unit radius in the XY plane.
But this integral is also easily computed in polar coordinates as
$$\int_{\theta=0}^{\theta=\pi/2}\int_{r=0}^{r=1} r^3drd\theta = \frac{\pi}{8} \; .$$
The basic idea of Monte Carlo integration is that the sample average of i.i.d. random variables with finite mean converge to the expected value of these random variables as the sample size gets larger. This is known as the law of large numbers in statistics.