QMCPy for Lebesgue Integration¶

This notebook will give examples of how to use QMCPy for integration problems that not are defined in terms of a standard measure. i.e. Uniform or Gaussian.

In [26]:

from qmcpy import *
from numpy import *

from qmcpy import * from numpy import *

Sample Problem 1¶

$$y = \int_{[0,2]} x^2 dx = 2\int_{[0,2]} \frac{x^2}{2} dx$$

where the middle expression may be viewed as WRT the Lebesgue measure and the right expression may be viewed as WRT the uniform measure

In [27]:

abs_tol = .01
dim = 1
a = 0
b = 2
true_value = 8./3

abs_tol = .01 dim = 1 a = 0 b = 2 true_value = 8./3

In [28]:

# Lebesgue Measure
integrand = CustomFun(
    true_measure = Lebesgue(Uniform(Halton(dim, seed=7, replications=32),lower_bound=a, upper_bound=b)),
    g = lambda x: (x**2).sum(-1))
solution,data = CubQMCCLT(integrand, abs_tol=abs_tol).integrate()
print('y = %.3f'%solution)
error = abs((solution-true_value))
if error>abs_tol:
    raise Exception("Not within error tolerance")

# Lebesgue Measure integrand = CustomFun( true_measure = Lebesgue(Uniform(Halton(dim, seed=7, replications=32),lower_bound=a, upper_bound=b)), g = lambda x: (x**2).sum(-1)) solution,data = CubQMCCLT(integrand, abs_tol=abs_tol).integrate() print('y = %.3f'%solution) error = abs((solution-true_value)) if error>abs_tol: raise Exception("Not within error tolerance")

y = 2.667

In [29]:

# Uniform Measure
integrand = CustomFun(
    true_measure = Uniform(Halton(dim, seed=7, replications=32), lower_bound=a, upper_bound=b),
    g = lambda x: (2*(x**2)).sum(-1))
solution,data = CubQMCCLT(integrand, abs_tol=abs_tol).integrate()
print('y = %.3f'%solution)
error = abs((solution-true_value))
if error>abs_tol:
    raise Exception("Not within error tolerance")

# Uniform Measure integrand = CustomFun( true_measure = Uniform(Halton(dim, seed=7, replications=32), lower_bound=a, upper_bound=b), g = lambda x: (2*(x**2)).sum(-1)) solution,data = CubQMCCLT(integrand, abs_tol=abs_tol).integrate() print('y = %.3f'%solution) error = abs((solution-true_value)) if error>abs_tol: raise Exception("Not within error tolerance")

y = 2.667

Sample Problem 2¶

$$y = \int_{[a,b]^d} ||x||_2^2 dx = \Pi_{i=1}^d (b_i-a_i)\int_{[a,b]^d} ||x||_2^2 \; [ \Pi_{i=1}^d (b_i-a_i)]^{-1} dx$$

where the middle expression may be viewed as WRT the Lebesgue measure and the right expression may be viewed as WRT the uniform measure

In [30]:

abs_tol = .001
dim = 2
a = array([1.,2.])
b = array([2.,4.])
true_value = ((a[0]**3-b[0]**3)*(a[1]-b[1])+(a[0]-b[0])*(a[1]**3-b[1]**3))/3
print('Answer = %.5f'%true_value)

abs_tol = .001 dim = 2 a = array([1.,2.]) b = array([2.,4.]) true_value = ((a[0]**3-b[0]**3)*(a[1]-b[1])+(a[0]-b[0])*(a[1]**3-b[1]**3))/3 print('Answer = %.5f'%true_value)

Answer = 23.33333

In [31]:

# Lebesgue Measure
integrand = CustomFun(
    true_measure = Lebesgue(Uniform(DigitalNetB2(dim, seed=7, replications=32), lower_bound=a, upper_bound=b)), 
    g = lambda x: (x**2).sum(-1))
solution,data = CubQMCCLT(integrand, abs_tol=abs_tol).integrate()
print('y = %.5f'%solution)
error = abs((solution-true_value))
if error>abs_tol:
    raise Exception("Not within error tolerance")

# Lebesgue Measure integrand = CustomFun( true_measure = Lebesgue(Uniform(DigitalNetB2(dim, seed=7, replications=32), lower_bound=a, upper_bound=b)), g = lambda x: (x**2).sum(-1)) solution,data = CubQMCCLT(integrand, abs_tol=abs_tol).integrate() print('y = %.5f'%solution) error = abs((solution-true_value)) if error>abs_tol: raise Exception("Not within error tolerance")

y = 23.33336

In [32]:

# Uniform Measure
integrand = CustomFun(
    true_measure = Uniform(DigitalNetB2(dim, seed=17, replications=32), lower_bound=a, upper_bound=b),
    g = lambda x: (b-a).prod()*(x**2).sum(-1))
solution,data = CubQMCCLT(integrand, abs_tol=abs_tol).integrate()
print('y = %.5f'%solution)
error = abs((solution-true_value))
if error>abs_tol:
    raise Exception("Not within error tolerance")

# Uniform Measure integrand = CustomFun( true_measure = Uniform(DigitalNetB2(dim, seed=17, replications=32), lower_bound=a, upper_bound=b), g = lambda x: (b-a).prod()*(x**2).sum(-1)) solution,data = CubQMCCLT(integrand, abs_tol=abs_tol).integrate() print('y = %.5f'%solution) error = abs((solution-true_value)) if error>abs_tol: raise Exception("Not within error tolerance")

y = 23.33333

Sample Problem 3¶

Integral that cannot be done in terms of any standard mathematical functions

$$y = \int_{[a,b]} \frac{\sin{x}}{\log{x}} dx$$

where the middle expression may be viewed as WRT the Lebesgue measure and the right expression may be viewed as WRT the uniform measure

Mathematica Code: Integrate[Sin[x]/Log[x], {x,a,b}]

In [33]:

abs_tol = .0001
dim = 1
a = 3
b = 5
true_value = -0.87961 

abs_tol = .0001 dim = 1 a = 3 b = 5 true_value = -0.87961

In [34]:

# Lebesgue Measure
integrand = CustomFun(
    true_measure = Lebesgue(Uniform(Lattice(dim, randomize=True, seed=7),a,b)), 
    g = lambda x: (sin(x)/log(x)).sum(-1))
solution,data = CubQMCLatticeG(integrand, abs_tol=abs_tol).integrate()
print('y = %.3f'%solution)
error = abs((solution-true_value))
if error>abs_tol:
    raise Exception("Not within error tolerance")

# Lebesgue Measure integrand = CustomFun( true_measure = Lebesgue(Uniform(Lattice(dim, randomize=True, seed=7),a,b)), g = lambda x: (sin(x)/log(x)).sum(-1)) solution,data = CubQMCLatticeG(integrand, abs_tol=abs_tol).integrate() print('y = %.3f'%solution) error = abs((solution-true_value)) if error>abs_tol: raise Exception("Not within error tolerance")

y = -0.880

Sample Problem 4¶

Integral over $\mathbb{R}^d$ $$y = \int_{\mathbb{R}^2} e^{-||x||_2^2} dx$$

In [35]:

abs_tol = .1
dim = 2
true_value = pi

abs_tol = .1 dim = 2 true_value = pi

In [36]:

integrand = CustomFun(
    true_measure = Lebesgue(Gaussian(Lattice(dim,seed=7))),
    g = lambda x: exp(-x**2).prod(1))
solution,data = CubQMCLatticeG(integrand,abs_tol=abs_tol).integrate()
print('y = %.3f'%solution)
error = abs((solution-true_value))
if error>abs_tol:
    raise Exception("Not within error tolerance")

integrand = CustomFun( true_measure = Lebesgue(Gaussian(Lattice(dim,seed=7))), g = lambda x: exp(-x**2).prod(1)) solution,data = CubQMCLatticeG(integrand,abs_tol=abs_tol).integrate() print('y = %.3f'%solution) error = abs((solution-true_value)) if error>abs_tol: raise Exception("Not within error tolerance")

y = 3.142

In [ ]: