Let be $(\mathbb{R},\mathcal{B}(\mathbb{R}),P)$ a probability space and $X,Y$ two real valued absolutely continuous random variables, where $f_X$, $f_Y$ and $f_{XY}$ are the respective density functions.
According to our definition the joint distribution function for $(x,y)\in \mathbb{R}^2$ is given by $$ P(X\leq x,Y\leq y)=\int\limits_{-\infty}^x\int\limits_{-\infty}^yf_{XY}(s,t)ds~\!dt. $$ If $X,Y$ are independent then $$ P(X\leq x,Y\leq y)=\int\limits_{-\infty}^x\int\limits_{-\infty}^yf_{X}(s)f_Y(t)dt~\!ds. $$
Now, let's consider $P(XY\leq 1)$ and $X> 0$. It is intuitve to come up with the integration limits that "surround" the respective area and simply say that $$ P(XY\leq 1)=\int\limits_{0}^{\infty}\int\limits_{-\infty}^{\frac{1}{x}}f_{X}(s)f_Y(t)dt~\!ds. $$ However, I am wondering how to rigorously derive this result from the above definition of the joint probability distribution? In other words, how do I get a set of the type $\{X\leq x,Y\leq y\}$ that is equivalent to $\{XY\leq 1\}$ and $\{X>0\}$?
For the joint density function of $X$ and $Y$ we have (by the usual measure-theoretic definition): $$P((X,Y) \in A) = \iint_A f_{X,Y} \,d\lambda_2$$ for all $A\in \mathcal{B}(\mathbb{R}^2)$ (where $\lambda_2 = \lambda \otimes \lambda$ is the Lebesgue measure on $\mathbb{R}^2$). So $$P(XY\leq 1, X>0) = \mathbb{E}[\mathbf{1}_{XY\leq 1}\mathbf{1}_{X>0}]=\iint_{\mathbb{R}^2}\mathbf{1}_{xy\leq 1}\mathbf{1}_{x>0}f_{XY}(x,y)\,d\lambda_2(x,y)$$ From there, Fubini's theorem will give you the final form of the integral. We used the general fact that $$\mathbb{E}[h(X,Y)] = \iint h(x,y)f_{XY}(x,y)\,d\lambda_2(x,y)$$ for all measurable functions $h:\mathbb{R}^2 \to [-\infty, \infty]$.