Generate a sequence $U_1,U_2,\ldots$ of independent uniform$[0,1]$ random variables. Let $Y_i = f(U_i)$, where $f(x)=\sqrt{1-x^2}$, $0 \leq x \leq 1$. Then, for sufficiently large $n$, $ \frac{{\sum\nolimits_{i = 1}^n {Y_i } }}{n} \approx \int_0^1 {f(x)\,{\rm d}x = } \int_0^1 {\sqrt {1 - x^2 } \,{\rm d}x} = \frac{\pi }{4}. $
EDIT: Elaborating.
Suppose that $U_1,U_2,\ldots$ is a sequence of independent uniform$[0,1]$ random variables, and $f$ is an integrable function on $[0,1]$, that is, $\int_0^1 {|f(x)|\,{\rm d}x} < \infty$. Then, the (finite) integral $\int_0^1 {f(x)\,{\rm d}x}$ can be approximated as follows. Let $Y_i = f(U_i)$, so the $Y_i$ are independent and identically distributed random variables, with mean (expectation) $\mu$ given by $ \mu = {\rm E}[Y_1] = {\rm E}[f(U_1)] = \int_0^1 {f(x)\,{\rm d}x}. $ By the strong law of large numbers, the average $\bar Y_n = \frac{{\sum\nolimits_{i = 1}^n {Y_i } }}{n}$ converges, with probability $1$, to the expectation $\mu$ as $n \to \infty$. That is, with probability $1$, $\bar Y_n \to \int_0^1 {f(x)\,{\rm d}x}$ as $n \to \infty$.
To get a probabilistic error bound, suppose further that $f$ is square-integrable on $[0,1]$, that is $ \int_0^1 {f^2 (x)\,{\rm d}x} < \infty $. Then, the $Y_i$ have finite variance, $\sigma_2$, given by $ \sigma^2 = {\rm Var}[Y_1] = {\rm E}[Y_1^2] - {\rm E}^2{[Y_1]} = {\rm E}[f^2{(U_1)}] - {\rm E}^2{[f(U_1)]} = \int_0^1 {f^2 (x) \,{\rm d}x} - \bigg[\int_0^1 {f(x)\,{\rm d}x} \bigg]^2 . $ By linearity of expectation, the average $\bar Y_n $ has expectation $ {\rm E}[\bar Y_n] = \mu. $ Since the $Y_i$ are independent, $\bar Y_n $ has variance $ {\rm Var}[\bar Y_n] = {\rm Var}\bigg[\frac{{Y_1 + \cdots + Y_n }}{n}\bigg] = \frac{1}{{n^2 }}{\rm Var}[Y_1 + \cdots + Y_n ] = \frac{n}{{n^2 }}{\rm Var}[Y_1 ] = \frac{{\sigma ^2 }}{n}. $ By Chebyshev's inequality, for any given $\varepsilon > 0$, $ {\rm P}\big[\big|\bar Y_n - {\rm E}[\bar Y_n]\big| \geq \varepsilon \big] \leq \frac{{{\rm Var}[\bar Y_n]}}{{\varepsilon ^2 }}, $ so $ {\rm P}\big[\big|\bar Y_n - \mu \big| \geq \varepsilon \big] \leq \frac{{\sigma^2}}{{n \varepsilon ^2 }}, $ and hence $ {\rm P}\bigg[\bigg|\bar Y_n - \int_0^1 {f(x)\,{\rm d}x} \bigg| \geq \varepsilon \bigg] \leq \frac{1}{{n \varepsilon ^2 }} \bigg \lbrace \int_0^1 {f^2 (x) \,{\rm d}x} - \bigg[\int_0^1 {f(x)\,{\rm d}x} \bigg]^2 \bigg \rbrace. $ So if $n$ is sufficiently large, with high probability the absolute difference between $\bar Y_n$ and $\int_0^1 {f(x)\,{\rm d}x}$ will be smaller than $\varepsilon$.
Returning to your specific question, letting $f(x)=\sqrt{1-x^2}$ thus gives $ {\rm P}\Big[\Big|\bar Y_n - \frac{\pi }{4} \Big| \geq \varepsilon \Big] \leq \frac{1}{{n \varepsilon ^2 }} \bigg \lbrace \int_0^1 {(1 - x^2) \,{\rm d}x} - \frac{\pi^2 }{16} \bigg \rbrace = \frac{1}{{n \varepsilon ^2 }} \bigg \lbrace \frac{2}{3} - \frac{\pi^2 }{16} \bigg \rbrace < \frac{1}{{20n\varepsilon ^2 }}, $ where $\bar Y_n = \frac{{\sum\nolimits_{i = 1}^n {\sqrt {1 - U_i^2 } } }}{n}$.