Let $x:\Omega\to\mathbb{R}^2$ be a 2-dimensional r.v., and $f(p,q)$ be its probability density function.
For what kind of $D\subseteq\mathbb{R}^2$, $\mathrm{Pr}[x\in D]=\iint_D f(p,q)\mathrm{d}p\mathrm{d}q$ ?
I couldn't find a proof for this from the definition of r.v., even in the most simple case of $D=\{(x,y)|x,y\geq 0,x+y\leq1\}$ (perhaps it's just me being stupid). (okay, now I see that this can be proven using the definition of integration, by splitting D into rectangles.but can this be applied to general D or is there a more elegant proof?)
This is a question for measure theory; I will use the basics of that subject freely in the answer. It may be possible to redo the proof using Riemann integration if $f$ is Riemann integrable (and attention is restricted to Riemann integrable sets $D$), but I haven't tried writing the details yet.
Let $\mu$ be the law of $x$. In other words, $\mu$ is the Borel probability measure on $\mathbb{R}^{2}$ defined by \begin{align*} \mu(D) = \mathbb{P}\{x \in D\}. \end{align*}
We know that if $D$ is a rectangle of the form $[a,b) \times [c,d)$, then $\mu(D) = \int_{D} f(p,q) \, dp \, dq$.
I claim that this determines $\mu$, that is, that $\mu(D) = \int_{D} f(p,q) \, dp \, dq$ for each Borel set $D$.
To see this, let $\mathscr{B}$ be the Borel subsets of $\mathbb{R}^{2}$. First, I claim that $\mathscr{B} = \sigma(\mathcal{D})$, where $\mathcal{D} = \{[a,b) \times [c,d) \, \mid \, a,b, c, d \in \mathbb{R}\}$. To see this, notice that if $\mathcal{O} \subseteq \mathbb{R}^{2}$ is open, then we can find open rectangles $\{[a_{n},b_{n}) \times [c_{n},d_{n}) \, \mid \, n \in \mathbb{N}\}$ with $\{(a_{n},b_{n},c_{n},d_{n}) \, \mid \, n \in \mathbb{N}\} \subseteq \mathbb{Q}^{4}$ such that $\mathcal{O} = \bigcup_{n = 1}^{\infty} [a_{n},b_{n}) \times [c_{n},d_{n})$. (Here I'm using the fact that $\mathbb{Q}$ is dense in $\mathbb{R}$.) Hence $\sigma(\mathcal{D})$ contains all the open sets, which implies it contains $\mathscr{B}$. Conversely, $\mathscr{B} \supseteq \mathcal{D}$ so $\mathscr{B}$ contains $\sigma(\mathcal{D})$.
Next, note that the following statement is an elementary consequence of Dynkin's $\pi-\lambda$ theorem:
Here I am using the fact that $\mathcal{D}$ is a $\pi$-system (finite intersections of rectangles are rectangles again) and the set $\{D \in \mathscr{B} \, \mid \, \mu_{1}(D) = \mu_{2}(D)\}$ is a $\lambda$-system (since $\mu_{1}(\mathbb{R}^{2}) = \mu_{2}(\mathbb{R}^{2}) = 1$).
Applying this with $\mu_{1} = \mu$ and $\mu_{2}(D) = \int_{D} f(p,q) \, dp \, dq$ and $\mathcal{A}$ the smallest algebra containing $\mathcal{D}$, we conclude that $\mu(D) = \int_{D} f(p,q) \, dp \, dq$ for all Borel sets $D$.