What is the probability that a natural number is a sum of two squares? Some natural numbers can be expressed as a sum of two squares:
$$2=1^2+1^2$$
$$25=3^2+4^2$$
$$50=7^2+1^2$$
If one chooses a random natural number, what would be the probability that that number is a sum of two squares? Is it zero?
I read about Lagrange´s theorem on squares, but it looks it can´t be useful here.

NOTE 1: "Square" means "square of a natural number".
NOTE 2: I am aware that the expression "random natural number" is not a strict math notion. However, as I said in a comment, one can adopt a reasonable strict definition, which is not difficult to devise at all. It is mentioned also in an answer below.
NOTE 3: A related question on SE: How to determine whether a number can be written as a sum of two squares?
 A: Consider the function $r: \mathbb{N} \to \mathbb{C}$:
$$
r(n) = \begin{cases}
1 &n \text{ is the sum of two squares (0 is a square)} \\
0 &\text{otherwise}.
\end{cases}
$$
One way to define the "probability that a random integer is the sum of two squares" would be to consider the distribution on the integers where $n$ selected with probability proportional to $n^x$, for $x > 1$, and then to take the limit as $x \to 1$. That is, we can try to compute:
$$
\lim_{x \to 1} \frac{\sum_{i=1}^\infty \frac{r(n)}{n^x}}{\sum_{i=1}^\infty \frac{1}{n^x}} \tag{1}.
$$
As mentioned here, $r(n)$ is $1$ iff every prime of the form $4k-1$ occurs an even number of times in $n$.
It follows that $r(n)$ is multiplicative, and
\begin{align*}
\frac{\sum_{i=1}^\infty \frac{r(n)}{n^x}}{\sum_{i=1}^\infty \frac{1}{n^x}}
&= \frac{\prod_{p \text{ prime}} \sum_{i=0}^\infty \frac{r(p^i)}{p^{ix}}}{\prod_{p \text{ prime}} \sum_{i=0}^\infty \frac{1}{p^{ix}}} \\
&= \prod_{p \text{ prime}} \frac{1 + \frac{r(p)}{p^x} + \frac{r(p^2)}{p^{2x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\
&= \prod_{p \equiv 3 \pmod{4}} \frac{1 + 0 + \frac{1}{p^{2x}} + 0 + \frac{1}{p^{4x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots}
\prod_{p \equiv 1,2 \pmod{4}} \frac{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\
&= \prod_{p \equiv 3 \pmod{4}} \frac{1 + \frac{1}{p^{2x}} + \frac{1}{p^{4x}} + \cdots}{1 + \frac{1}{p^x} + \frac{1}{p^{2x}} + \cdots} \\
&= \prod_{p \equiv 3 \pmod{4}} \frac{1 \big/ \left(1 - \tfrac{1}{p^{2x}}\right)}{1 \big/ \left(1 - \tfrac{1}{p^{x}}\right)} \\
&= \prod_{p \equiv 3 \pmod{4}} \frac{1}{1 + p^{-x}}. \\
\end{align*}
Now if you plug in $x = 1$, you get the product
$$
\prod_{p \equiv 3 \pmod{4}} \frac{p}{p + 1}
= \prod_{p \equiv 3 \pmod{4}} \left( 1 - \frac{1}{p + 1} \right)
= 0,
$$
by the reasoning of Hagen von Eitzen.
So by this definition of probability, the probability that a random integer is a sum of two squares is zero.
A: the density is zero, and one may be quite precise about it: the count of numbers up to some large real $x$ that are the sum of two squares is asymptotic to
$$  \frac{0.7642... \, x}{\sqrt{\log x}}  $$ 
where the logarithm is base $e,$ and the $0.7642...$ is defined by an infinite product. See the last few pages in LeVeque. This is combined volumes 1 and 2, it is the last few pages in volume 2.
A: Without going into detail what a "random natural number" might be, we could consider the density of such numbers. 
The possibility to express $n$ depends on the prime factorization of $n$: It may be divisible by $2$ and by primes $\equiv 1\pmod 4$ as much as it likes, but each prime $\equiv 3\pmod 4$ must occur in an even power.
The prime $3$ "spoils" all numbers divisible by $3$ (that's $\frac13$), except those divisible by $9$ (that's $\frac19$), though it does spoil those divisible by $27$, execept ...
All in all the prime $3$ spoils $\frac13-\frac19+\frac1{27}-\frac1{81}\pm\ldots = \frac14$ (geometric series) of all numbers in a large range. In general, a prime $p\equiv 3\pmod 4$ spoils $\frac1{p+1}$ and the effects of distinct primes are independent. Hence the limit density of natural numbers expressible as sum of squares is
$$ \prod_{p\equiv 3\pmod 4}\left(1-\frac1{p+1}\right)$$
This product however diverges to $0$ because it is well-known that $\sum_p\frac1p$ diverges (where it doesn't matter that we use only half of all primes).
Thus: For sufficiently large $N$, the probability that a number picked uniformly from $\{1,\ldots,N\}$ is the sum of two squares becomes arbitrarily small. It is for any $\let\epsilon\varepsilon\epsilon>0$ we can find $N_0$ such that for all $N>N_0$ said probability is $<\epsilon$.
A: Though this is not an answer to the original question, I'll write it in answer to VividD's question under Hagen von Eitzen post (and I believe it is not totally unrelated to the question).
Let $A_n$ be the number of different pairs $(x,y)$ of non-negative integers solving the equation $x^2 + y^2 = n$. Then it is natural to define the expected number of compositions of a randomly chosen natural number into sum of two squares as
$$
A = \lim_{n\to\infty} \frac1n \sum_{k=1}^n A_k.
$$
But the latter quantity is $\frac1n$ times the number of points of the form $(x/\sqrt{n},y/\sqrt{n})$, where $x,y$ are integers, inside the non-negative quarter of the unit circle centered at the origin (minus one, but this does not matter). So here we have a Jordan approximation of the area of the quarter. Therefore, $A = \frac\pi4$.
So, quite paradoxically, the expected value of the number of compositions of a randomly chosen non-negative integer into two squares is positive despite this number is zero "with probability 1".
