Return to Answer

although we know it should be 0. Having said that, we can always supplement built-in integration rules with user-defined ones (see e.g. Why aren't these additions of integrals and summations equal?) to expand a class of symbolically or numerically integrable functions - for this purpose Mathematica is most likely the best tool.
While we could always remedy various problems algorithmically, we shouldn't expect that it could be done in full generality (e.g. because of the finite number of states of computers); otherwise, we should supplement the system's built-in integration rules with infinitely many user-defined rules to be able to integrate every Lebesgue integrable function.
Next editions of Mathematica may involve more powerful symbolic capabilities for measure theory and Lebesgue integration problems, but we should realize that there will always be some limitations of the algorithmic approach to integration in the realm of integrable functions.

An especially fast method would be e.g. Method -> "AdaptiveMonteCarlo" with Boole[ f1[x, y] > 5 f2[x, y]] - appropriate region selector. Instead of Boole we could use HeavisideTheta, e.g. HeavisideTheta[ f1[x, y] - 5 f2[x, y]] but in this case it appears to be fairly slower (see e.g. this for the case when it is much faster). Working with "AdaptiveMonteCarlo", one should remember that the method provides a rather rough estimation of the result:

although we know it should be 0. Having said that, we can always supplement built-in integration rules with user-defined ones (see e.g. Why aren't these additions of integrals and summations equal?) to expand a class of symbolically or numerically integrable functions - for this purpose Mathematica is most likely the best tool.
While we could always remedy various problems algorithmically, we shouldn't expect that it could be done in full generality (e.g. because of the finite number of states of computers); otherwise, we should supplement the system's built-in integration rules with infinitely many user-defined rules to be able to integrate every Lebesgue integrable function.
Next editions of Mathematica may involve more powerful symbolic capabilities for measure theory and Lebesgue integration problems, but we should realize that there will always be some limitations of the algorithmic approach to integration in the realm of integrable functions.

An especially fast method would be e.g. Method -> "AdaptiveMonteCarlo" with Boole[ f1[x, y] > 5 f2[x, y]] - appropriate region selector. Instead of Boole we could use HeavisideTheta, e.g. HeavisideTheta[ f1[x, y] - 5 f2[x, y]] but in this case it appears to be fairly slower (see e.g. this for the case when it is much faster). Working with "AdaptiveMonteCarlo", one should remember that the method provides a rather rough estimation of the result:

The answer is no because of fundamental mathematical limitations which origin in the set theory regarding countability (see e.g. Cantor's theorem) - functions over a given set are more numerous than its (power) cardinality. Neither Mathematica nor any other system can integrate every function in even much more restricted class, namely Riemann integrable functions, all Riemann integrals are equal to Lebesgue integrals if the former are well defined. The class of functions which could be integrated over a domain in $\mathbb{R}^n$ in Mathematica is of "measure zero" in the class of Lebesgue integrable functions. More precisely we need rather Baire categories to work with general topological concepts of class of adequate functions.
While we are to calculate a definite integral we are going to think of NIntegrate rather than of Integrate.

Let's try to integrate a simple Lebesgue integrable function defined in $\mathbb{R}$.

although we know it should be 0.
  Having said that we can always supplement built-in integration rules with another user-defined ones (see e.g. Why aren't these additions of integrals and summations equal?) to expand a class of symbolically or numerically integrable functions - for this purpose Mathematica is most likely the best tool.
While we could always remedy various problems algorithmically we shouldn't expect it could be done in full generality (e.g. because of finite number of states of computers) since we expect one should supplement the system builit-in integration rules with infinitely many user-defined rules to be able to integrate every Lebesgue integrable function.
Next editions of Mathematica may involve more powerful symbolic capabilities for measure theory and Lebesgue integration problems but we should realize that there will always be some limitations of algoritmic approach to integration in the realm of integrable functions.

The above cosiderations concern the problem of integration of possibly wide class of functions. However if one restricts to integration of bi-variate Gaussian density functions there is no need for distinction between Riemann and Lebesgue integrals. Since one needs fast numerical results I'd recommend taking a closer look at the NIntegrate Integration Strategies tutorial, especially at Crude Monte Carlo and Quasi Monte Carlo Strategies and Global Adaptive Monte Carlo and Quasi Monte Carlo Strategies sections.

An especially fast method would be e.g. Method -> "AdaptiveMonteCarlo" with Boole[ f1[x, y] > 5 f2[x, y]] - appropriate region selector. Instead of Boole we could use HeavisideTheta, e.g. HeavisideTheta[ f1[x, y] - 5 f2[x, y]] but in this case it appears to be fairly slower (see e.g. this for the case when it is much faster). Working with "AdaptiveMonteCarlo" one should remember that the method provides rather rough estimation of the result:

slower but considerably more stable would be: Method -> "AdaptiveQuasiMonteCarlo".

The answer is no, because of fundamental mathematical limitations which originate in set theory regarding countability (see e.g. Cantor's theorem) - functions over a given set are more numerous than the set's (power) cardinality. Neither Mathematica nor any other system can integrate every function in an even much more restricted class; namely, Riemann integrable functions. All Riemann integrals are equal to Lebesgue integrals if the former are well defined. The class of functions which could be integrated over a domain in $\mathbb{R}^n$ in Mathematica is of "measure zero" in the class of Lebesgue integrable functions. More precisely, we need Baire categories to work with general topological concepts of the class of adequate functions.
When we calculate a definite integral, we are going to think of NIntegrate rather than Integrate.

Let's try to integrate a simple Lebesgue integrable function defined in $\mathbb{R}$:

although we know it should be 0. Having said that, we can always supplement built-in integration rules with user-defined ones (see e.g. Why aren't these additions of integrals and summations equal?) to expand a class of symbolically or numerically integrable functions - for this purpose Mathematica is most likely the best tool.
While we could always remedy various problems algorithmically, we shouldn't expect that it could be done in full generality (e.g. because of the finite number of states of computers); otherwise, we should supplement the system's built-in integration rules with infinitely many user-defined rules to be able to integrate every Lebesgue integrable function.
Next editions of Mathematica may involve more powerful symbolic capabilities for measure theory and Lebesgue integration problems, but we should realize that there will always be some limitations of the algorithmic approach to integration in the realm of integrable functions.

The above considerations concern the problem of integration of a possibly wide class of functions. However if one restricts to integration of bi-variate Gaussian density functions there is no need for distinction between Riemann and Lebesgue integrals. Since one needs fast numerical results I'd recommend taking a closer look at the NIntegrate Integration Strategies tutorial, especially at Crude Monte Carlo and Quasi Monte Carlo Strategies and Global Adaptive Monte Carlo and Quasi Monte Carlo Strategies sections.

An especially fast method would be e.g. Method -> "AdaptiveMonteCarlo" with Boole[ f1[x, y] > 5 f2[x, y]] - appropriate region selector. Instead of Boole we could use HeavisideTheta, e.g. HeavisideTheta[ f1[x, y] - 5 f2[x, y]] but in this case it appears to be fairly slower (see e.g. this for the case when it is much faster). Working with "AdaptiveMonteCarlo", one should remember that the method provides a rather rough estimation of the result:

A slower but considerably more stable method would be Method -> "AdaptiveQuasiMonteCarlo".

Edit

The above cosiderations concern the problem of integration of possibly wide class of functions. However if one restricts to integration of bi-variate Gaussian density functions there is no need for distinction between Riemann and Lebesgue integrals. Since one needs fast numerical results I'd recommend taking a closer look at the NIntegrate Integration Strategies tutorial, especially at Crude Monte Carlo and Quasi Monte Carlo Strategies and Global Adaptive Monte Carlo and Quasi Monte Carlo Strategies sections.

Let's define e.g.

f1[x_, y_] := PDF[ BinormalDistribution[{1, 3/2}, {1/2, 3/5}, 1/3], {x, y}]
f2[x_, y_] := PDF[ BinormalDistribution[{4/3, 7/3}, {1, 2/3}, 2/5], {x, y}]

and choose e.g. τ = 5.

An especially fast method would be e.g. Method -> "AdaptiveMonteCarlo" with Boole[ f1[x, y] > 5 f2[x, y]] - appropriate region selector. Instead of Boole we could use HeavisideTheta, e.g. HeavisideTheta[ f1[x, y] - 5 f2[x, y]] but in this case it appears to be fairly slower (see e.g. this for the case when it is much faster). Working with "AdaptiveMonteCarlo" one should remember that the method provides rather rough estimation of the result:

NIntegrate[ f1[x, y] Boole[ f1[x, y] > 5 f2[x, y]], {x, -∞, ∞}, {y, -∞, ∞}, 
            Method -> "AdaptiveMonteCarlo"]

0.370381

slower but considerably more stable would be: Method -> "AdaptiveQuasiMonteCarlo".