Is it possible to calculate a Lebesgue integral in Mathematica?

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}$:

f[x_] /; x ∈ Rationals && 0 <= x <= 1 := 1
f[x_] /; ! (x ∈ Rationals) && 0 <= x <= 1 := 0

f[Sqrt[3]/2]
f[1/2]

0
1

but neither Integrate nor NIntegrate can calculate adequate integrals:

Integrate[ f[x], {x, 0, 1}]
NIntegrate[ f[x], {x, 0, 1}]

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.

Edit

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.

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 a rather rough estimation of the result:

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

0.370381

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

[...] I am only interested in very fast numerical methods, no analytical results are needed.

[...] I have no idea how I can do it in Mathematica

The package AdaptiveNumericalLebesgueIntegration.m has Lebesgue integration strategy and rules implementations and it is discussed in detail in the blog post "Adaptive numerical Lebesgue integration by set measure estimates" and [3].

Using the profiling capabilities described in the referenced documents the speed of these Lebesgue integration algorithms can be evaluated for integrals of interest. Generally, for usual integrands, the algorithms are slower, but often require much less sampling points than, say, "AdaptiveMonteCarlo", which for expensively to evaluate integrands might produce results faster.

Lebesgue strategy and rules invocation examples

Here are examples of using the integration strategy and rules defined in the package:

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/Misc/AdaptiveNumericalLebesgueIntegration.m"]

NIntegrate[Sqrt[x + y], {x, 0, 2}, {y, 0, 1}]
(* 2.38176 *)

NIntegrate[Sqrt[x + y], {x, 0, 2}, {y, 0, 1}, 
 Method -> {LebesgueIntegration, "Points" -> 2000, 
   "PointGenerator" -> "Sobol", "PointwiseMeasure" -> "VoronoiMesh"}, 
 PrecisionGoal -> 3]
(* 2.38236 *)

NIntegrate[Sqrt[x + y], {x, 0, 2}, {y, 0, 1}, 
 Method -> {LebesgueIntegrationRule, "Points" -> 6000, 
   "PointGenerator" -> "Sobol", "PointwiseMeasure" -> "VoronoiMesh", 
   "AxisSelector" -> {"MinVariance", "SubsampleFraction" -> 0.05}}, 
 PrecisionGoal -> 3]
(* 2.38095 *)

NIntegrate[Sqrt[x + y], {x, 0, 2}, {y, 0, 1}, 
 Method -> {GridLebesgueIntegrationRule, "Points" -> 2000, 
   "PointGenerator" -> Random, "GridSizes" -> 6, 
   "AxisSelector" -> Random}, 
 PrecisionGoal -> 3]
(* 2.37127 *)

In the second integral above the strategy LebesgueIntegration uses the following Voronoi diagram to estimate the set measures:

For more details see the mentioned blog post.

Traces of Lebesgue integration process

Here is an example that demonstrates the "dimension reduction" using the implemented Lebesgue integration strategy for a high dimensional (conventional) integral:

res = 
  Reap@NIntegrate[Sqrt[x + y + z], {x, 0, 2}, {y, 0, 1}, {z, 0, 3}, 
    Method -> {LebesgueIntegration, "Points" -> 12000, 
      "PointGenerator" -> Random, 
      "LebesgueIntegralVariableSymbol" -> f}, 
    EvaluationMonitor :> {Sow[f]}, PrecisionGoal -> 4, 
    MaxRecursion -> 7];
res = DeleteCases[res, f, \[Infinity]]

(* {10.1842, {{0.344945, 0.426233, 0.584845, 0.809906, 1.07998, 
   1.37175, 1.66352, ..., 1.50829, 
   1.51821, 1.53227, 1.54915, 1.56739, 1.58563, 1.6025, 1.61657, 
   1.62648, 1.63157}}} *)

 ListPlot[res[[2, 1]], 
  PlotLabel -> Style[Row[{"Integral estimate:", res[[1]]}], Larger]]

Compare with the result using the default NIntegrate algorithms:

NIntegrate[Sqrt[x + y + z], {x, 0, 2}, {y, 0, 1}, {z, 0, 3}]
(* 10.2016 *)

References

[1] B. L. Burrows, A new approach to numerical integration, 1. Inst. Math. Applics., 26(1980), 151-173.

[2] T. He, Dimensionality Reducing Expansion of Multivariate Integration, 2001, Birkhauser Boston. ISBN-13:978-1-4612-7414-8 .

[3] A. Antonov, Adaptive numerical Lebesgue integration by set measure estimates, (2016), MathematicaForPrediction project at GitHub.

I just found a paper, called: "Familiarizing Students with Definition of Lebesgue Integral: Examples of Calculation Directly from Its Definition Using Mathematica" https://link.springer.com/content/pdf/10.1007%2Fs11786-017-0321-5.pdf

It might help some, who are looking here for Lebesgue and Mathematica.