# Is it possible to calculate a Lebesgue integral in Mathematica?

As the title says, I wonder if it is possible to calculate a Lebesgue integral in Mathematica, especially when the domain of integration is $\mathbb{R}^N$, or in other words multivatiate Lebesgue integration is of interest. As an example, one could take $f_0$ and $f_1$ as two different (possibly correlated) bi-variate Gaussian density functions and consider

$$\int_{\frac{f_1}{f_0}>\tau}f_0\mbox{d}\mu=\int_{\Large\{x,y:\frac{f_1(x,y)}{f_0(x,y)}>\tau\Large\}}f_0(x,y)\mbox{d}x\mbox{d}y$$ for some known number $\tau$.

If it is not possible can one manipulate the existent functions of Mathematica to get a method which can calcuate the Lebesgue integral. I am only interested in very fast numerical methods, no analytical results are needed.

Note: Why am I asking such a question here?

• I checked the Mathematica help and found no result
• I searched the keyword "Lebesgue" among the questions in mathematica.stackexchange.com but couldnt see any related answer
• I have no idea how I can do it in Mathematica
-
If your density measures are continuous you might be able to use Boole to get a reasonable approximation within NIntegrate. –  Daniel Lichtblau Feb 10 '14 at 16:24
@DanielLichtblau yes my density measures are continuous. What is Boole, how kind of approximation? –  Seyhmus Güngören Feb 10 '14 at 16:32
I have to defer to the documentation for Boole and NIntegrate. –  Daniel Lichtblau Feb 10 '14 at 17:04
@DanielLichtblau okay I am aware of Boolean algebra, just logical zeros and one but I didnt know if my Boole was also the same Boole of Mathematica but it seems the same.. –  Seyhmus Güngören Feb 10 '14 at 17:12
It's in effect a predicate function that returns 0 or 1. So you could use Boole[f1[x,y]>f0[x,y]*t]*integrand. –  Daniel Lichtblau Feb 10 '14 at 17:54

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

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 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.

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, -∞, ∞},

0.370381

slower but considerably more stable would be: Method -> "AdaptiveQuasiMonteCarlo".
Isnt it possible to obtain some set of intervals of type $\{\{(x_1,y_1),(x_2,y_2)\},\{(x_3,y_3),(x_4,y_4)\}....\}:\frac{f_1(x,y)}{f_0(x,y‌​)}>\tau$ and sum all the integrals with NIntegrate[...,{x_1,y_1},{x_2,y_2}]+NIntegrate[....,{x_3,y_3},{x_4,y_4}]+... I can do this in Matlab if I discritize the continues bivatiate density but the memory requirements is too large! I must save a multivariate data although I can find which indexes satisfy the given condition simply by using find command- –  Seyhmus Güngören Feb 10 '14 at 16:18
I am quoting from the question "I am only interested in very fast numerical methods, no analytical results are needed." I defined $f_0$ and $f_1$ to be Gaussian bi variate distributions. The answerer has the flexibility to choose $\tau$. Do you think that I should have said that $\tau=5$? then it is $5$. –  Seyhmus Güngören Feb 10 '14 at 20:37