21
$\begingroup$

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 google about it and found no answer
  • 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
Improve this question
$\endgroup$
6
  • $\begingroup$ If your density measures are continuous you might be able to use Boole to get a reasonable approximation within NIntegrate. $\endgroup$
    – Daniel Lichtblau
    Commented Feb 10, 2014 at 16:24
  • $\begingroup$ @DanielLichtblau yes my density measures are continuous. What is Boole, how kind of approximation? $\endgroup$
    – Seyhmus Güngören
    Commented Feb 10, 2014 at 16:32
  • 2
    $\begingroup$ I have to defer to the documentation for Boole and NIntegrate. $\endgroup$
    – Daniel Lichtblau
    Commented Feb 10, 2014 at 17:04
  • $\begingroup$ @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.. $\endgroup$
    – Seyhmus Güngören
    Commented Feb 10, 2014 at 17:12
  • 2
    $\begingroup$ Boole[] in Mathematica is what's more commonly known as an Iverson bracket. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Nov 2, 2015 at 7:16

3 Answers 3

Reset to default
23
$\begingroup$

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

Improve this answer
$\endgroup$
8
  • $\begingroup$ 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- $\endgroup$
    – Seyhmus Güngören
    Commented Feb 10, 2014 at 16:18
  • $\begingroup$ @SeyhmusGüngören You can always widen a class of integrable functions in Mathematica, but the more powerful capabilities you would like to have the more specific rules are needed to be supplemented. If you have a specific problem in mind I think it would be the best way to ask another well focused question. $\endgroup$
    – Artes
    Commented Feb 10, 2014 at 18:16
  • $\begingroup$ I think it should be clear in the question that a multivariate integral is asked and an answer to the given example suffices. I don't know how much clearer I could ask. $\endgroup$
    – Seyhmus Güngören
    Commented Feb 10, 2014 at 18:50
  • 1
    $\begingroup$ @SeyhmusGüngören From the example I gave it is obvious that you can't expect anymore regarding one or multivariable functions on the symbolic level. If you want to integrate numerically a function, you have to provide an adequate example, it is impossible to tackle the problem generally. Your first comment doesn't clarify anything. $\endgroup$
    – Artes
    Commented Feb 10, 2014 at 18:59
  • $\begingroup$ 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$. $\endgroup$
    – Seyhmus Güngören
    Commented Feb 10, 2014 at 20:37
12
$\begingroup$

[...] 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:

"VoronoiMeshFor2000SobolPoints"

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

enter image description here

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.

Improve this answer
$\endgroup$
3
$\begingroup$

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.

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.