Quasi-Monte Carlo (QMC) Software in QMCPy¶
A tutorial based on this notebook is at https://media.ed.ac.uk/playlist/dedicated/51612401/1_0z0wec2z/1_2k12mwiw.
As an example of available QMC software, we introduce the Python library QMCPy. This Jupyter notebook saves your typing.
QMCPy is a community effort. This early release includes contributions from
- Mike Giles MLMC and MLQMC software
- Marius Hofert and Christiane Lemieux's QRNG
- Pierre L'Ecuyer's Lattice Builder
- Dirk Nuyens's Magic Point Shop (MPS)
- Art Owen's Halton sequences
- PyTorch
- Guaranteed Automatic Integration Library (GAIL)
and depends on the NumPy and SciPy Python packages.
View the companion pdf slides at https://speakerdeck.com/fjhickernell/quasi-monte-carlo-software and the introductory QMCPy blog at https://qmcpy.wordpress.com.
Installation¶
QMCPy can be installed with pip install qmcpy or cloned from the QMCSoftware GitHub repository.
In [1]:
Copied!
import qmcpy #we import the environment at the start to use it
import numpy as np #basic numerical routines in Python
import time #timing routines
import warnings #to suppress warnings when needed
import torch #only needed for PyTorch Sobol' backend
from torch.quasirandom import SobolEngine
from matplotlib import pyplot; #plotting
pyplot.rc('font', size=16) #set defaults so that the plots are readable
pyplot.rc('axes', titlesize=16)
pyplot.rc('axes', labelsize=16)
pyplot.rc('xtick', labelsize=16)
pyplot.rc('ytick', labelsize=16)
pyplot.rc('legend', fontsize=16)
pyplot.rc('figure', titlesize=16)
#a helpful plotting method to show increasing numbers of points
def plot_successive_points(distrib,ld_name,first_n=64,n_cols=1,pt_clr='bgkcmy',
xlim=[0,1],ylim=[0,1],coord1 = 0,coord2 = 1):
fig,ax = pyplot.subplots(nrows=1,ncols=n_cols,figsize=(5*n_cols,5.5))
if n_cols==1: ax = [ax]
last_n = first_n*(2**n_cols)
points = distrib.gen_samples(n=last_n)
for i in range(n_cols):
n = first_n
nstart = 0
for j in range(i+1):
n = first_n*(2**j)
ax[i].scatter(points[nstart:n,coord1],points[nstart:n,coord2],color=pt_clr[j])
nstart = n
ax[i].set_title('n = %d'%n)
ax[i].set_xlim(xlim); ax[i].set_xticks(xlim); ax[i].set_xlabel('$x_{i,%d}$'%(coord1+1))
ax[i].set_ylim(ylim); ax[i].set_yticks(ylim); ax[i].set_ylabel('$x_{i,%d}$'%(coord2+1))
ax[i].set_aspect((xlim[1]-xlim[0])/(ylim[1]-ylim[0]))
fig.suptitle('%s Points'%ld_name)
print('QMCPy Version',qmcpy.__version__)
import qmcpy #we import the environment at the start to use it
import numpy as np #basic numerical routines in Python
import time #timing routines
import warnings #to suppress warnings when needed
import torch #only needed for PyTorch Sobol' backend
from torch.quasirandom import SobolEngine
from matplotlib import pyplot; #plotting
pyplot.rc('font', size=16) #set defaults so that the plots are readable
pyplot.rc('axes', titlesize=16)
pyplot.rc('axes', labelsize=16)
pyplot.rc('xtick', labelsize=16)
pyplot.rc('ytick', labelsize=16)
pyplot.rc('legend', fontsize=16)
pyplot.rc('figure', titlesize=16)
#a helpful plotting method to show increasing numbers of points
def plot_successive_points(distrib,ld_name,first_n=64,n_cols=1,pt_clr='bgkcmy',
xlim=[0,1],ylim=[0,1],coord1 = 0,coord2 = 1):
fig,ax = pyplot.subplots(nrows=1,ncols=n_cols,figsize=(5*n_cols,5.5))
if n_cols==1: ax = [ax]
last_n = first_n*(2**n_cols)
points = distrib.gen_samples(n=last_n)
for i in range(n_cols):
n = first_n
nstart = 0
for j in range(i+1):
n = first_n*(2**j)
ax[i].scatter(points[nstart:n,coord1],points[nstart:n,coord2],color=pt_clr[j])
nstart = n
ax[i].set_title('n = %d'%n)
ax[i].set_xlim(xlim); ax[i].set_xticks(xlim); ax[i].set_xlabel('$x_{i,%d}$'%(coord1+1))
ax[i].set_ylim(ylim); ax[i].set_yticks(ylim); ax[i].set_ylabel('$x_{i,%d}$'%(coord2+1))
ax[i].set_aspect((xlim[1]-xlim[0])/(ylim[1]-ylim[0]))
fig.suptitle('%s Points'%ld_name)
print('QMCPy Version',qmcpy.__version__)
QMCPy Version 1.6.3.2a
Generating low discrepancy (LD) points via a DiscreteDistribution object¶
We generate some points used for quasi-Monte Carlo methods. These points are called low discrepancy (LD for short) and are created as an instance of a DiscreteDistribution class.
Integration lattices¶
Here are some (randomly shifted) integration lattice points. This is a two step procees:
i) construct the DiscreteDistribution object lattice, and then
ii) construct a number of points from the sequence.
The structure of these points favors n that is a power of 2.
In [2]:
Copied!
lattice = qmcpy.Lattice(dimension=2) #define a discrete LD distribution based on an integration lattice
print(lattice) #print the properties of the lattice object
n = 16 #number of points to generate
points = lattice.gen_samples(n) #construct some points
print(f'\nLD Lattice Points with shape {points.shape}\n'+str(points)) #these points have 15 significant digit precision but only three digits are shown
plot_successive_points(lattice,'Lattice',n)
lattice = qmcpy.Lattice(dimension=2) #define a discrete LD distribution based on an integration lattice
print(lattice) #print the properties of the lattice object
n = 16 #number of points to generate
points = lattice.gen_samples(n) #construct some points
print(f'\nLD Lattice Points with shape {points.shape}\n'+str(points)) #these points have 15 significant digit precision but only three digits are shown
plot_successive_points(lattice,'Lattice',n)
Lattice (AbstractLDDiscreteDistribution)
d 2^(1)
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 264030601762852985433978628854999910438
LD Lattice Points with shape (16, 2)
[[0.14448802 0.23129441]
[0.64448802 0.73129441]
[0.39448802 0.98129441]
[0.89448802 0.48129441]
[0.26948802 0.60629441]
[0.76948802 0.10629441]
[0.51948802 0.35629441]
[0.01948802 0.85629441]
[0.20698802 0.91879441]
[0.70698802 0.41879441]
[0.45698802 0.66879441]
[0.95698802 0.16879441]
[0.33198802 0.29379441]
[0.83198802 0.79379441]
[0.58198802 0.04379441]
[0.08198802 0.54379441]]
Rerunning the commands above yields a different sequence of points because these points are randomly shifted modulo 1.
We may also construct a subsequence of points in the middle of the sequence. Note that the points below match those above.
In [3]:
Copied!
more_points = lattice.gen_samples(n_min=4,n_max=n) #get more points
print('LD Lattice Points with shape',more_points.shape,'\n'+str(more_points))
more_points = lattice.gen_samples(n_min=4,n_max=n) #get more points
print('LD Lattice Points with shape',more_points.shape,'\n'+str(more_points))
LD Lattice Points with shape (12, 2) [[0.26948802 0.60629441] [0.76948802 0.10629441] [0.51948802 0.35629441] [0.01948802 0.85629441] [0.20698802 0.91879441] [0.70698802 0.41879441] [0.45698802 0.66879441] [0.95698802 0.16879441] [0.33198802 0.29379441] [0.83198802 0.79379441] [0.58198802 0.04379441] [0.08198802 0.54379441]]
Each $d$-dimensional point is one row in the array.
As we increase the number of points, they fill $[0,1]^d$ evenly. The next points are placed in between the existing points. Here we illustrate with $d=2$.
In [4]:
Copied!
plot_successive_points(lattice,'Lattice',n_cols=5)
plot_successive_points(lattice,'Lattice',n_cols=5)
IID uniform points do not fill space as well¶
Contrast this with independent and identically distributed (IID) points. Although successive points fill the square, they do so without knowledge of the others and produce clusters and gaps. Think of it this way
- LD = evenly spread
- IID = points do not know about each other
(Since the first parameter in the DiscreteDistribution object is the dimension, you need not identify it.)
In [5]:
Copied!
iid = qmcpy.IIDStdUniform(2) #standard uniform IID random vector generator from NumPy
print(iid) #print the properties of iid
plot_successive_points(iid,'IID Uniform',n_cols=5)
iid = qmcpy.IIDStdUniform(2) #standard uniform IID random vector generator from NumPy
print(iid) #print the properties of iid
plot_successive_points(iid,'IID Uniform',n_cols=5)
IIDStdUniform (AbstractIIDDiscreteDistribution)
d 2^(1)
replications 1
entropy 151929946340103232908357805579611369175
Rows and columns are not interchangeable for LD¶
For LD sequences we must differentiate between the cooordinates of the point (column) and which point (row). The transpose of an LD array is not LD. This differs from IID multivariate points with IID marginals.
In the example below we reverse the roles of the rows and columns. The transposed lattice points do not fill space at all, while the transposed IID points are as good as the originals.
In [6]:
Copied!
d = 16 #dimension
n = 2 #number of points
lattice = qmcpy.Lattice(d) #define a discrete LD distribution based on an integration lattice
lattice_pts = lattice.gen_samples(n) #the first parameter in the .gen_samples method is the number of points
iid = qmcpy.IIDStdUniform(d)
iid_pts = iid.gen_samples(n)
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(10,5.5))
ax[0].scatter(lattice_pts[0,0:d],lattice_pts[1,0:d],color='b')
ax[0].set_title('Transposed Lattice Points')
ax[1].scatter(iid_pts[0,0:d],iid_pts[1,0:d],color='b')
ax[1].set_title('Transposed IID Points')
for ii in range(2):
ax[ii].set_xlim([0,1]); ax[ii].set_xticks([0,1]); ax[ii].set_xlabel('$x_{1,j}$')
ax[ii].set_ylim([0,1]); ax[ii].set_yticks([0,1]); ax[ii].set_ylabel('$x_{2,j}$')
ax[ii].set_aspect(1)
d = 16 #dimension
n = 2 #number of points
lattice = qmcpy.Lattice(d) #define a discrete LD distribution based on an integration lattice
lattice_pts = lattice.gen_samples(n) #the first parameter in the .gen_samples method is the number of points
iid = qmcpy.IIDStdUniform(d)
iid_pts = iid.gen_samples(n)
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(10,5.5))
ax[0].scatter(lattice_pts[0,0:d],lattice_pts[1,0:d],color='b')
ax[0].set_title('Transposed Lattice Points')
ax[1].scatter(iid_pts[0,0:d],iid_pts[1,0:d],color='b')
ax[1].set_title('Transposed IID Points')
for ii in range(2):
ax[ii].set_xlim([0,1]); ax[ii].set_xticks([0,1]); ax[ii].set_xlabel('$x_{1,j}$')
ax[ii].set_ylim([0,1]); ax[ii].set_yticks([0,1]); ax[ii].set_ylabel('$x_{2,j}$')
ax[ii].set_aspect(1)
Sobol' points¶
Another LD sequence is the Sobol' points. Again, new points fill in the gaps between existing points. Since these are randomly digitally shifted Sobol' points, rerunning this command gives a different set of points. For Sobol' points as well, n should normally be a power of 2.
In [7]:
Copied!
sobol = qmcpy.Sobol(2) #Sobol LD generator
print(sobol) #note below that the default is qrng from Hofert and Lemieux
points = sobol.gen_samples(16)
print(f'\nLD Sobol\' Points with shape {points.shape}\n'+str(points))
plot_successive_points(sobol,'Sobol\'',n_cols=5)
sobol = qmcpy.Sobol(2) #Sobol LD generator
print(sobol) #note below that the default is qrng from Hofert and Lemieux
points = sobol.gen_samples(16)
print(f'\nLD Sobol\' Points with shape {points.shape}\n'+str(points))
plot_successive_points(sobol,'Sobol\'',n_cols=5)
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 2^(1)
replications 1
randomize LMS DS
gen_mats_source joe_kuo.6.21201.txt
order RADICAL INVERSE
t 63
alpha 1
n_limit 2^(32)
entropy 73208827127383995610457310988530317873
LD Sobol' Points with shape (16, 2)
[[0.9431598 0.47365992]
[0.22125498 0.72515188]
[0.66129464 0.90348679]
[0.37738574 0.15491697]
[0.79964553 0.60309908]
[0.02104669 0.35362196]
[0.58028993 0.0173281 ]
[0.36467173 0.76778926]
[0.88650748 0.81477058]
[0.16853181 0.06430855]
[0.72964997 0.30712722]
[0.44965495 0.55660491]
[0.87190807 0.20193054]
[0.0894106 0.95049973]
[0.52754482 0.67862835]
[0.30804326 0.4271372 ]]
Halton Points¶
A third kind of LD sequence is Halton points, which are also randomized. Again, new points fill in the gaps between existing points. For Halton points there are no favored numbers of points.
In [8]:
Copied!
halton = qmcpy.Halton(2)
print(halton)
points = halton.gen_samples(20)
print(f'\nLD Halton Points with shape {points.shape}\n'+str(points))
plot_successive_points(halton,'Halton',n_cols=5,first_n=60)
halton = qmcpy.Halton(2)
print(halton)
points = halton.gen_samples(20)
print(f'\nLD Halton Points with shape {points.shape}\n'+str(points))
plot_successive_points(halton,'Halton',n_cols=5,first_n=60)
Halton (AbstractLDDiscreteDistribution)
d 2^(1)
replications 1
randomize LMS DP
t 63
n_limit 2^(32)
entropy 240669245053839331683025566616811162761
LD Halton Points with shape (20, 2)
[[0.44785589 0.67105156]
[0.81034548 0.6477606 ]
[0.01802002 0.17997437]
[0.72376384 0.96444929]
[0.31082674 0.49671198]
[0.94638376 0.02481055]
[0.22522047 0.82620954]
[0.51752374 0.3543589 ]
[0.40526643 0.33101523]
[0.85320912 0.7104405 ]
[0.06844815 0.57598763]
[0.67312161 0.21524872]
[0.35369039 0.89273277]
[0.90379429 0.53194055]
[0.17457824 0.06420577]
[0.56795187 0.86102595]
[0.49980629 0.39323972]
[0.75841033 0.25883637]
[0.03822453 0.74815281]
[0.70357459 0.60958533]]
Timing¶
The time to generate LD points depends on several factors, including the hardware and the particular generator. One can generate one million points in 64 dimensions in a matter of seconds, somewhat slower than IID uniform points. You may replace the DiscreteDistribution object in the first line of code below to test the timings of other kinds of points. Halton seems to be the slowest (probably because the backend is implemented in Python rather than C).
In [9]:
Copied!
ld = qmcpy.Lattice(64) #define a discrete LD distribution
print(ld) #print the properties of the lattice object
start_time = time.time() #time now
points = ld.gen_samples(2**20) #construct some points
end_time = time.time() #time after points are constructed
print(f'\nLD Points with shape {points.shape}\n'+str(points))
print(f'\nTime to construct points is %.1e seconds'%(end_time - start_time))
ld = qmcpy.Lattice(64) #define a discrete LD distribution
print(ld) #print the properties of the lattice object
start_time = time.time() #time now
points = ld.gen_samples(2**20) #construct some points
end_time = time.time() #time after points are constructed
print(f'\nLD Points with shape {points.shape}\n'+str(points))
print(f'\nTime to construct points is %.1e seconds'%(end_time - start_time))
Lattice (AbstractLDDiscreteDistribution)
d 2^(6)
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 323672018533134169508871551556781883883
LD Points with shape (1048576, 64)
[[4.84051442e-01 1.74980418e-01 6.35194475e-01 ... 7.14805235e-01
1.59923415e-01 7.06555592e-02]
[9.84051442e-01 6.74980418e-01 1.35194475e-01 ... 2.14805235e-01
6.59923415e-01 5.70655559e-01]
[7.34051442e-01 9.24980418e-01 3.85194475e-01 ... 9.64805235e-01
9.09923415e-01 8.20655559e-01]
...
[2.34050488e-01 2.50775592e-01 6.81364710e-01 ... 9.84270681e-01
1.84181733e-02 2.36378406e-01]
[9.84050488e-01 5.00775592e-01 9.31364710e-01 ... 7.34270681e-01
2.68418173e-01 4.86378406e-01]
[4.84050488e-01 7.75592029e-04 4.31364710e-01 ... 2.34270681e-01
7.68418173e-01 9.86378406e-01]]
Time to construct points is 1.7e-01 seconds
Transform LD sequences to mimic other distributions via a TrueMeasure object¶
LD sequences mimic $\mathcal{U}[0,1]^d$ by design. If you want to mimic another distribution, LD points can be transformed accordingly. This is done via a TrueMeasure object, which takes a DiscreteDistribution object as input. Choose the TrueMeasure according to the probability distribution that you wish to mimic, and input the LD DiscreteDistribution object that you wish to use.
Uniform on $[\boldsymbol{a},\boldsymbol{b}]$¶
An affine transformation can be used to turn $\mathcal{U}[0,1]^d$ points into $\mathcal{U}[\boldsymbol{a},\boldsymbol{b}]$ points. The TrueMeasure object qmcpy.Uniform performs the transformation automatically. Then we generate points as before, but now they mimic our new distribution.
In [10]:
Copied!
ld = qmcpy.Sobol(2) #Sobol' points
another_unif_ld = qmcpy.Uniform(ld, lower_bound=[-2,2], upper_bound=[2,4]) #define the desired probability distribution with sobol as input
points = another_unif_ld.gen_samples(2**8)
print(another_unif_ld)
print(f'\nUniform LD points with shape {points.shape}\n'+str(points)) #these points have 15 significant digit precision but only three digits are shown
plot_successive_points(another_unif_ld,'Uniform LD',first_n=2**8,xlim=[-2,2],ylim=[2,4])
ld = qmcpy.Sobol(2) #Sobol' points
another_unif_ld = qmcpy.Uniform(ld, lower_bound=[-2,2], upper_bound=[2,4]) #define the desired probability distribution with sobol as input
points = another_unif_ld.gen_samples(2**8)
print(another_unif_ld)
print(f'\nUniform LD points with shape {points.shape}\n'+str(points)) #these points have 15 significant digit precision but only three digits are shown
plot_successive_points(another_unif_ld,'Uniform LD',first_n=2**8,xlim=[-2,2],ylim=[2,4])
Uniform (AbstractTrueMeasure)
lower_bound [-2 2]
upper_bound [2 4]
Uniform LD points with shape (256, 2)
[[ 1.63211619e+00 3.73527085e+00]
[-1.96458235e+00 2.76718358e+00]
[ 9.09260722e-01 2.10405185e+00]
[-7.49905356e-01 3.38671429e+00]
[ 1.09655014e+00 2.74751185e+00]
[-1.31071479e+00 3.77735132e+00]
[ 3.18856288e-01 3.10066974e+00]
[-2.58190328e-02 2.38125511e+00]
[ 1.85902575e+00 2.35110377e+00]
[-1.55013461e+00 3.13179771e+00]
[ 5.74189117e-01 3.99793234e+00]
[-7.72499888e-01 2.52790732e+00]
[ 1.38496040e+00 3.35363327e+00]
[-1.20972064e+00 2.13615638e+00]
[ 1.70346353e-01 2.98483491e+00]
[-4.86928113e-01 3.51664289e+00]
[ 1.59319338e+00 2.83569329e+00]
[-1.75352515e+00 3.67851923e+00]
[ 8.70309343e-01 3.45102593e+00]
[-5.38819580e-01 2.04362477e+00]
[ 1.18259209e+00 3.83803777e+00]
[-1.47465285e+00 2.68294084e+00]
[ 4.04926810e-01 2.43774730e+00]
[-1.89785667e-01 3.03241946e+00]
[ 1.92921095e+00 3.18830328e+00]
[-1.72992369e+00 2.28284011e+00]
[ 6.44406709e-01 2.58858045e+00]
[-9.52321360e-01 3.93337465e+00]
[ 1.33015606e+00 2.20048515e+00]
[-1.01455070e+00 3.29318905e+00]
[ 1.15509619e-01 3.58513540e+00]
[-2.91725784e-01 2.92810048e+00]
[ 1.69674328e+00 2.02283906e+00]
[-1.89601564e+00 3.49030058e+00]
[ 9.81610664e-01 3.63771572e+00]
[-6.73680637e-01 2.85592312e+00]
[ 1.04382486e+00 3.02634606e+00]
[-1.35944513e+00 2.49563675e+00]
[ 2.58530228e-01 2.62559582e+00]
[-8.23294418e-02 3.84563622e+00]
[ 1.78281433e+00 3.87505519e+00]
[-1.62246743e+00 2.59520034e+00]
[ 5.05639523e-01 2.27578809e+00]
[-8.37113668e-01 3.24521810e+00]
[ 1.44145742e+00 2.88827178e+00]
[-1.14941179e+00 3.60634367e+00]
[ 2.19059476e-01 3.27338160e+00]
[-4.34216234e-01 2.24073455e+00]
[ 1.51690842e+00 3.42208176e+00]
[-1.82581517e+00 2.07936273e+00]
[ 8.01804452e-01 2.79132255e+00]
[-6.03508819e-01 3.69837261e+00]
[ 1.23904058e+00 2.43523923e+00]
[-1.41426473e+00 3.09068726e+00]
[ 4.53717300e-01 3.78885312e+00]
[-1.37120398e-01 2.69407406e+00]
[ 1.99388680e+00 2.53843444e+00]
[-1.66143601e+00 3.94351623e+00]
[ 7.16679522e-01 3.18480345e+00]
[-8.76049791e-01 2.34014643e+00]
[ 1.27750410e+00 3.54175642e+00]
[-1.06332406e+00 2.94891535e+00]
[ 5.51386166e-02 2.17250235e+00]
[-3.48160972e-01 3.32991865e+00]
[ 1.68443108e+00 2.66428490e+00]
[-1.96890853e+00 3.82156387e+00]
[ 8.98955486e-01 3.04917923e+00]
[-6.91855546e-01 2.45623134e+00]
[ 1.08720775e+00 3.66175516e+00]
[-1.25561124e+00 2.81720496e+00]
[ 3.72134031e-01 2.06227642e+00]
[-3.30914651e-02 3.46749554e+00]
[ 1.83223923e+00 3.31185587e+00]
[-1.50856227e+00 2.21693948e+00]
[ 6.09782398e-01 2.91132119e+00]
[-7.93551486e-01 3.56666242e+00]
[ 1.42346943e+00 2.29961511e+00]
[-1.23176561e+00 3.20677197e+00]
[ 1.46475575e-01 3.91470354e+00]
[-4.46349142e-01 2.57212184e+00]
[ 1.59862091e+00 3.76507331e+00]
[-1.80474448e+00 2.73253306e+00]
[ 8.13112855e-01 2.39554271e+00]
[-5.27659039e-01 3.11375194e+00]
[ 1.12638706e+00 2.75289168e+00]
[-1.46640614e+00 3.72218437e+00]
[ 4.11345801e-01 3.39898800e+00]
[-2.43918834e-01 2.11902634e+00]
[ 1.88678700e+00 2.14844526e+00]
[-1.70399443e+00 3.36859247e+00]
[ 6.64358817e-01 3.50233576e+00]
[-9.89012298e-01 2.97176378e+00]
[ 1.35305252e+00 3.14610054e+00]
[-1.05220258e+00 2.36417062e+00]
[ 7.60300175e-02 2.51561415e+00]
[-2.66757468e-01 3.98296885e+00]
[ 1.73720513e+00 3.07783324e+00]
[-1.92001307e+00 2.42066099e+00]
[ 9.59692714e-01 2.70882325e+00]
[-6.35054132e-01 3.80142033e+00]
[ 1.02289814e+00 2.06461688e+00]
[-1.32373269e+00 3.40951789e+00]
[ 2.99983303e-01 3.71122997e+00]
[-1.09240950e-01 2.80590412e+00]
[ 1.77542454e+00 2.96178788e+00]
[-1.56931658e+00 3.55632271e+00]
[ 5.60869853e-01 3.31515758e+00]
[-8.46338810e-01 2.15995383e+00]
[ 1.49963234e+00 3.95828065e+00]
[-1.15959764e+00 2.55098630e+00]
[ 2.14614269e-01 2.32727723e+00]
[-3.82026098e-01 3.17024050e+00]
[ 1.54175770e+00 2.47744239e+00]
[-1.86541953e+00 3.00935717e+00]
[ 7.64277676e-01 3.86181997e+00]
[-5.80492975e-01 2.64448040e+00]
[ 1.20247642e+00 3.47412016e+00]
[-1.39419538e+00 2.00395782e+00]
[ 4.79529198e-01 2.87412082e+00]
[-1.79671249e-01 3.65470795e+00]
[ 1.93963436e+00 3.62455655e+00]
[-1.65514201e+00 2.90524872e+00]
[ 7.25051103e-01 2.22453890e+00]
[-9.32135663e-01 3.25451570e+00]
[ 1.28879166e+00 2.61139932e+00]
[-1.12040307e+00 3.89392442e+00]
[ 3.80216143e-03 3.22700857e+00]
[-3.42860104e-01 2.25881449e+00]
[ 1.64846946e+00 2.25745798e+00]
[-1.94236157e+00 3.22540792e+00]
[ 9.33914805e-01 3.90407310e+00]
[-7.19383823e-01 2.62179213e+00]
[ 1.12267729e+00 3.25998869e+00]
[-1.28264265e+00 2.22976775e+00]
[ 3.37659195e-01 2.89097689e+00]
[-5.07107836e-03 3.61052886e+00]
[ 1.86025008e+00 3.64141262e+00]
[-1.54305808e+00 2.86106964e+00]
[ 5.82737638e-01 2.01040743e+00]
[-7.58099114e-01 3.48032563e+00]
[ 1.39594306e+00 2.65365246e+00]
[-1.19677768e+00 3.87123616e+00]
[ 1.73028253e-01 3.00702416e+00]
[-4.82285964e-01 2.47486523e+00]
[ 1.56267932e+00 3.15718932e+00]
[-1.77818702e+00 2.31398191e+00]
[ 8.48096027e-01 2.55719177e+00]
[-5.55180645e-01 3.96473026e+00]
[ 1.16183658e+00 2.16937002e+00]
[-1.49344805e+00 3.32432963e+00]
[ 3.76847111e-01 3.55374555e+00]
[-2.15905118e-01 2.95945487e+00]
[ 1.91480262e+00 2.80430347e+00]
[-1.73846451e+00 3.70987346e+00]
[ 6.37322628e-01 3.41991071e+00]
[-9.53537988e-01 2.07476556e+00]
[ 1.32552137e+00 3.80664918e+00]
[-1.01724039e+00 2.71429624e+00]
[ 1.02574124e-01 2.40663330e+00]
[-3.02716229e-01 3.06356141e+00]
[ 1.70919046e+00 3.96891344e+00]
[-1.88551344e+00 2.50131460e+00]
[ 9.86733664e-01 2.36943290e+00]
[-6.70502686e-01 3.15160697e+00]
[ 1.05042069e+00 2.98213078e+00]
[-1.35871681e+00 3.51245863e+00]
[ 2.73426806e-01 3.36702716e+00]
[-7.33003155e-02 2.14712409e+00]
[ 1.81138234e+00 2.11648441e+00]
[-1.59585973e+00 3.39669021e+00]
[ 5.25906716e-01 3.73157487e+00]
[-8.18806722e-01 2.76203804e+00]
[ 1.46415899e+00 3.11999072e+00]
[-1.13256241e+00 2.40202564e+00]
[ 2.49085294e-01 2.71945429e+00]
[-4.10042667e-01 3.75175039e+00]
[ 1.51373826e+00 2.56982405e+00]
[-1.83094564e+00 3.91216161e+00]
[ 7.91310047e-01 3.21591834e+00]
[-6.15963474e-01 2.30900561e+00]
[ 1.23000375e+00 3.57314535e+00]
[-1.42915375e+00 2.91755998e+00]
[ 4.52981281e-01 2.20361656e+00]
[-1.43708670e-01 3.29877709e+00]
[ 1.97557215e+00 3.45319600e+00]
[-1.68169565e+00 2.04822102e+00]
[ 6.90064120e-01 2.82271138e+00]
[-9.04610239e-01 3.66701745e+00]
[ 1.25333832e+00 2.46635421e+00]
[-1.09335735e+00 3.05954623e+00]
[ 3.82970322e-02 3.82024270e+00]
[-3.70870007e-01 2.66271959e+00]
[ 1.65976941e+00 3.34321421e+00]
[-1.99942245e+00 2.18555377e+00]
[ 8.82594568e-01 2.94246598e+00]
[-7.14068659e-01 3.53555119e+00]
[ 1.06841247e+00 2.33097414e+00]
[-1.27636677e+00 3.17538702e+00]
[ 3.46014554e-01 3.94584901e+00]
[-6.11712500e-02 2.54101136e+00]
[ 1.82369833e+00 2.69543034e+00]
[-1.52297062e+00 3.79045353e+00]
[ 6.08565743e-01 3.08053835e+00]
[-8.00635651e-01 2.42484617e+00]
[ 1.42077994e+00 3.69289986e+00]
[-1.23640015e+00 2.78609394e+00]
[ 1.35485274e-01 2.09363479e+00]
[-4.59284431e-01 3.43610969e+00]
[ 1.62084184e+00 2.24233496e+00]
[-1.78839100e+00 3.27473788e+00]
[ 8.43634601e-01 3.59595062e+00]
[-5.03004805e-01 2.87812287e+00]
[ 1.15445917e+00 3.23998949e+00]
[-1.44027908e+00 2.27031534e+00]
[ 4.32093663e-01 2.60922826e+00]
[-2.25115961e-01 3.88932725e+00]
[ 1.89386349e+00 3.85868764e+00]
[-1.70277019e+00 2.63889138e+00]
[ 6.78759497e-01 2.48943152e+00]
[-9.80463810e-01 3.01989669e+00]
[ 1.36599563e+00 2.84650669e+00]
[-1.04121972e+00 3.62854343e+00]
[ 8.06723779e-02 3.49287749e+00]
[-2.64075414e-01 2.02517184e+00]
[ 1.73207452e+00 2.93064218e+00]
[-1.92318343e+00 3.58743295e+00]
[ 9.47237855e-01 3.28379830e+00]
[-6.45548684e-01 2.19133855e+00]
[ 1.00800914e+00 3.92713610e+00]
[-1.33276944e+00 2.58209776e+00]
[ 2.93395118e-01 2.29591913e+00]
[-1.09976944e-01 3.20162644e+00]
[ 1.75516493e+00 3.04647510e+00]
[-1.58763115e+00 2.45204708e+00]
[ 5.32309489e-01 2.67767879e+00]
[-8.72954185e-01 3.83253158e+00]
[ 1.46959891e+00 2.03325752e+00]
[-1.18376362e+00 3.44090282e+00]
[ 1.91905028e-01 3.68008429e+00]
[-3.98867827e-01 2.83701422e+00]
[ 1.55225969e+00 3.53094266e+00]
[-1.85297249e+00 2.99889055e+00]
[ 7.67455475e-01 2.13065019e+00]
[-5.75370188e-01 3.34837122e+00]
[ 1.20320483e+00 2.51778422e+00]
[-1.38759952e+00 3.98756510e+00]
[ 4.88558359e-01 3.13311863e+00]
[-1.64774579e-01 2.35266884e+00]
[ 1.96624215e+00 2.38355266e+00]
[-1.62657397e+00 3.10321143e+00]
[ 7.43358081e-01 3.76820472e+00]
[-9.11868376e-01 2.73812111e+00]
[ 1.30564082e+00 3.38023160e+00]
[-1.09770165e+00 2.09781330e+00]
[ 2.79755744e-02 2.78050674e+00]
[-3.12834497e-01 3.74834986e+00]]
Gaussian¶
A similar process can be followed for constructing points that mimic a multivariate Gaussian distribution. The transformation used here assigns the directions with the larger variance to the lower numbered coordinates of the LD sequence via principal component analysis (PCA) or equivalently an eigvenvector-eigenvalue decomposition of the covariance matrix. This takes advantage of the fact that the lower numbered coordinates of an LD sequence tend to have better evenness.
In [11]:
Copied!
ld = qmcpy.Lattice(2)
gaussian_ld = qmcpy.Gaussian(ld, mean=[3,2], covariance=[[9,5], [5,4]]) #specify the desired mean and covariance of your multivariate Gaussian distribution
points = gaussian_ld.gen_samples(2**8)
print(gaussian_ld)
print(f'\nGaussian LD Points with shape {points.shape}\n'+str(points)) #these points have 15 significant digit precision but only three digits are shown
plot_successive_points(gaussian_ld,'Gaussian LD',first_n=2**8,xlim=[-6,12],ylim=[-2,6])
ld = qmcpy.Lattice(2)
gaussian_ld = qmcpy.Gaussian(ld, mean=[3,2], covariance=[[9,5], [5,4]]) #specify the desired mean and covariance of your multivariate Gaussian distribution
points = gaussian_ld.gen_samples(2**8)
print(gaussian_ld)
print(f'\nGaussian LD Points with shape {points.shape}\n'+str(points)) #these points have 15 significant digit precision but only three digits are shown
plot_successive_points(gaussian_ld,'Gaussian LD',first_n=2**8,xlim=[-6,12],ylim=[-2,6])
Gaussian (AbstractTrueMeasure)
mean [3 2]
covariance [[9 5]
[5 4]]
decomp_type PCA
Gaussian LD Points with shape (256, 2)
[[ 9.80784089e-01 -6.70135637e-01]
[ 4.29201337e+00 3.08939243e+00]
[ 2.70351262e+00 1.39177419e+00]
[ 6.94591370e+00 5.61220091e+00]
[ 1.22732660e+00 1.58259254e+00]
[ 5.96795261e+00 2.99596433e+00]
[ 2.44027019e+00 3.89861289e+00]
[-1.57840943e+00 -8.93333455e-01]
[ 1.08644957e+00 5.76156485e-01]
[ 4.28113054e+00 4.32931865e+00]
[ 2.76553638e+00 2.33245354e+00]
[ 1.16726000e+01 6.26695439e+00]
[ 2.91426798e+00 -1.81254438e-02]
[ 6.31198394e+00 4.15908938e+00]
[ 4.19815770e+00 2.11946301e+00]
[-8.47601153e-01 5.25690347e-01]
[ 3.23117324e-01 1.13251342e+00]
[ 5.01299750e+00 2.51783050e+00]
[ 3.94190142e+00 -8.97084397e-02]
[ 8.38405732e+00 5.35168099e+00]
[ 1.93590871e+00 1.01037649e+00]
[ 5.35428481e+00 4.79291114e+00]
[ 3.50884070e+00 2.69727756e+00]
[-3.92563166e-01 -1.34090316e+00]
[ 1.99818505e+00 -2.64924500e-01]
[ 5.17839575e+00 3.54721126e+00]
[ 3.44389373e+00 1.75332272e+00]
[-3.09434454e+00 -7.35422900e-01]
[ 2.01945669e+00 1.96944279e+00]
[ 7.27802715e+00 3.68449383e+00]
[ 3.36065324e+00 4.02606972e+00]
[ 4.54772993e-02 2.20186863e-02]
[-2.28034432e-01 1.66370309e+00]
[ 4.64142815e+00 2.81801881e+00]
[ 3.08866057e+00 1.03037222e+00]
[ 7.63842422e+00 5.39260651e+00]
[ 1.59034291e+00 1.28519127e+00]
[ 6.56705459e+00 2.41716598e+00]
[ 3.12967941e+00 3.04594991e+00]
[-9.96559899e-01 -1.04047120e+00]
[ 1.47195482e+00 2.74291161e-01]
[ 4.78855644e+00 3.83706943e+00]
[ 3.10372116e+00 2.04433326e+00]
[-4.83428748e+00 -8.83988370e-02]
[ 1.59676521e+00 2.38507554e+00]
[ 6.76249366e+00 3.94719157e+00]
[ 4.63717928e+00 1.68283168e+00]
[-3.75332096e-01 2.50487752e-01]
[ 7.19791850e-01 8.36442410e-01]
[ 5.50904053e+00 2.02427292e+00]
[ 2.37115736e+00 2.71142598e+00]
[ 9.38353714e+00 5.47267375e+00]
[ 2.31380818e+00 6.78318397e-01]
[ 5.86980213e+00 4.39978668e+00]
[ 3.84776228e+00 2.41617606e+00]
[-1.57086820e+00 1.16086154e+00]
[ 7.57617121e-01 2.04616706e+00]
[ 5.55504055e+00 3.29274496e+00]
[ 3.84864655e+00 1.35854687e+00]
[-2.23685566e+00 -8.06858572e-01]
[ 2.36515462e+00 1.67575215e+00]
[ 8.12666627e+00 2.95547973e+00]
[ 3.91586794e+00 3.41017200e+00]
[ 4.67532378e-01 -2.38064730e-01]
[ 1.50066305e+00 -1.32000733e+00]
[ 4.46684609e+00 2.95268125e+00]
[ 2.88485110e+00 1.22944002e+00]
[ 7.29676633e+00 5.47080847e+00]
[ 1.41414900e+00 1.42604782e+00]
[ 6.21779887e+00 2.78437733e+00]
[ 2.88390396e+00 3.31198873e+00]
[-1.28308860e+00 -9.56198304e-01]
[ 1.27325597e+00 4.36135487e-01]
[ 4.56530836e+00 4.03249049e+00]
[ 2.93710416e+00 2.18437228e+00]
[-3.91402699e+00 -3.85694469e+00]
[ 1.25186889e+00 2.80783663e+00]
[ 6.53317939e+00 4.05377234e+00]
[ 4.39651028e+00 1.93491215e+00]
[-6.00491974e-01 3.75421729e-01]
[ 5.28836279e-01 9.74250689e-01]
[ 5.23007762e+00 2.32004358e+00]
[ 2.09144320e+00 3.03428837e+00]
[ 8.83218247e+00 5.38570864e+00]
[ 2.11653447e+00 8.58389504e-01]
[ 5.63166058e+00 4.56081781e+00]
[ 3.67960417e+00 2.55428135e+00]
[ 2.44456679e-02 -1.70939455e+00]
[ 2.40267088e-01 2.73233648e+00]
[ 5.36542734e+00 3.42017714e+00]
[ 3.63105061e+00 1.58044418e+00]
[-2.62790196e+00 -7.73512510e-01]
[ 2.19607272e+00 1.81680154e+00]
[ 7.60289295e+00 3.47033475e+00]
[ 3.68784258e+00 3.63726842e+00]
[ 2.53096150e-01 -9.91643665e-02]
[ 8.71882533e-02 1.33595653e+00]
[ 4.82112334e+00 2.67683927e+00]
[ 3.34744682e+00 7.42028224e-01]
[ 7.99446374e+00 5.35498027e+00]
[ 1.76250063e+00 1.14945454e+00]
[ 4.91904817e+00 5.28729484e+00]
[ 3.32911627e+00 2.85538786e+00]
[-7.06939793e-01 -1.15867333e+00]
[ 1.69770894e+00 6.62606218e-02]
[ 4.98849897e+00 3.68238565e+00]
[ 3.27079641e+00 1.90363422e+00]
[-3.70119871e+00 -6.58423391e-01]
[ 1.82629835e+00 2.14819152e+00]
[ 7.00694737e+00 3.82962630e+00]
[ 5.03330779e+00 1.18022534e+00]
[-1.62133607e-01 1.35816015e-01]
[ 9.03807980e-01 7.06671980e-01]
[ 6.48347952e+00 6.05492508e-01]
[ 2.58162056e+00 2.50049634e+00]
[ 1.01484810e+01 5.66443364e+00]
[ 2.54910402e+00 4.35658204e-01]
[ 6.09307251e+00 4.27122471e+00]
[ 4.01857154e+00 2.27451894e+00]
[-1.13898675e+00 7.36126324e-01]
[ 1.01838301e+00 1.77331681e+00]
[ 5.75302301e+00 3.15585439e+00]
[ 4.15171717e+00 9.98632717e-01]
[-1.89251109e+00 -8.45125326e-01]
[ 2.53253640e+00 1.53691311e+00]
[ 6.52451339e+00 5.90564465e+00]
[ 4.11120689e+00 3.23730652e+00]
[ 7.00349184e-01 -4.12794172e-01]
[ 3.12903043e-01 5.06784577e-01]
[ 5.23701863e+00 1.63323025e+00]
[ 2.10370540e+00 2.42787169e+00]
[ 7.81321429e+00 4.41672992e+00]
[ 2.01401963e+00 3.82761606e-01]
[ 5.39515676e+00 4.01847859e+00]
[ 3.55624162e+00 2.15847632e+00]
[-2.26869882e+00 4.36029399e-01]
[ 4.45694705e-01 1.69437255e+00]
[ 5.16895788e+00 2.97393410e+00]
[ 3.59252940e+00 1.05911050e+00]
[-6.23113932e+00 -3.35709699e+00]
[ 2.06985359e+00 1.41634830e+00]
[ 5.04691692e+00 6.33090977e+00]
[ 3.61755087e+00 3.12698675e+00]
[-4.38881857e-02 -6.48615372e-01]
[ 1.10896483e+00 -5.12705353e-02]
[ 4.43661688e+00 3.52690279e+00]
[ 2.81934141e+00 1.79160221e+00]
[ 7.64277320e+00 6.91022546e+00]
[ 1.31890045e+00 2.07900114e+00]
[ 6.20664465e+00 3.52121851e+00]
[ 4.38023257e+00 1.35390774e+00]
[-9.94371673e-01 -2.11710203e-01]
[ 1.26109624e+00 1.00440968e+00]
[ 6.26466286e+00 1.87693426e+00]
[ 2.85587338e+00 2.76970758e+00]
[-2.16989369e+00 -1.85498381e+00]
[ 2.81967871e+00 7.41299630e-01]
[ 6.72004106e+00 4.74064806e+00]
[ 4.32142525e+00 2.54165766e+00]
[-6.25766108e-01 1.21636543e+00]
[ 6.98113320e-01 2.56831018e-01]
[ 3.97847564e+00 3.96525834e+00]
[ 2.48328928e+00 2.07511098e+00]
[ 8.53730344e+00 4.19817722e+00]
[ 2.81632848e+00 -6.26914449e-01]
[ 5.79718770e+00 3.76389093e+00]
[ 3.90988176e+00 1.84955683e+00]
[-1.53008891e+00 1.27647173e-03]
[ 8.94860040e-01 1.28696485e+00]
[ 5.57540068e+00 2.64795152e+00]
[ 2.29902873e+00 3.41115686e+00]
[-3.22882722e+00 -1.99053931e+00]
[ 2.41449031e+00 1.12931673e+00]
[ 6.18606224e+00 5.01705735e+00]
[ 3.97887355e+00 2.81703294e+00]
[ 5.49044003e-01 -1.12896364e+00]
[ 1.72428793e+00 -7.04494923e-01]
[ 4.80648617e+00 3.23968068e+00]
[ 3.16414175e+00 1.49316535e+00]
[ 9.04343526e+00 6.65396255e+00]
[ 1.72445071e+00 1.70107574e+00]
[ 6.68233403e+00 3.20032331e+00]
[ 3.12924643e+00 3.64602194e+00]
[-5.17875678e-01 -4.04146572e-01]
[ 1.61627185e+00 7.31596185e-01]
[ 4.93439318e+00 4.38816869e+00]
[ 3.22120292e+00 2.43896543e+00]
[-1.29967617e+00 -2.05532926e+00]
[ 1.38198489e+00 3.33261337e+00]
[ 7.23904439e+00 4.56714696e+00]
[ 4.69376994e+00 2.22292626e+00]
[-9.71168169e-02 7.77271632e-01]
[ 5.05714388e-01 3.83725512e-01]
[ 3.55065517e+00 4.50880565e+00]
[ 2.30390827e+00 2.23482021e+00]
[ 8.14586073e+00 4.32739486e+00]
[ 2.27527539e+00 1.05013998e-01]
[ 5.59813573e+00 3.88528002e+00]
[ 3.72748544e+00 2.01280341e+00]
[-1.85110408e+00 1.56525580e-01]
[ 6.89612337e-01 1.46063091e+00]
[ 5.36201580e+00 2.82592661e+00]
[ 3.96851925e+00 5.80277017e-01]
[-4.05584301e+00 -2.25960929e+00]
[ 2.23951448e+00 1.27754109e+00]
[ 5.85828948e+00 5.27641679e+00]
[ 3.80395141e+00 2.96223566e+00]
[ 2.18955674e-01 -8.29581008e-01]
[ 1.35717483e+00 -2.80033648e-01]
[ 4.62550037e+00 3.37579336e+00]
[ 2.98772871e+00 1.64894810e+00]
[ 8.34861301e+00 6.62139654e+00]
[ 1.53567911e+00 1.86793380e+00]
[ 6.43010523e+00 3.38005596e+00]
[ 5.07249565e+00 3.67437441e-01]
[-7.52399767e-01 -3.05940995e-01]
[ 1.43733165e+00 8.71143648e-01]
[ 4.59764615e+00 4.75166677e+00]
[ 3.04633166e+00 2.59159720e+00]
[-1.73719364e+00 -1.90314909e+00]
[ 3.11675256e+00 3.93446854e-01]
[ 6.97647726e+00 4.64733947e+00]
[ 4.50007427e+00 2.39381590e+00]
[-3.33246279e-01 9.54386555e-01]
[ 8.96106436e-01 1.16728837e-01]
[ 4.22937758e+00 3.70978664e+00]
[ 2.65294752e+00 1.93083775e+00]
[ 9.06609324e+00 3.92140543e+00]
[ 1.02842147e+00 2.40812980e+00]
[ 5.99815714e+00 3.64536435e+00]
[ 4.11596980e+00 1.64850390e+00]
[-1.25047317e+00 -1.13258988e-01]
[ 1.08225390e+00 1.13985760e+00]
[ 5.83394485e+00 2.40004612e+00]
[ 2.63087555e+00 3.00391129e+00]
[-2.64764975e+00 -1.88023932e+00]
[ 2.60263516e+00 9.59032325e-01]
[ 6.46121587e+00 4.85729863e+00]
[ 4.14956870e+00 2.67972035e+00]
[-1.22132384e+00 1.96165234e+00]
[ 7.63974718e-02 2.12831078e+00]
[ 4.98587165e+00 3.10853071e+00]
[ 3.35827152e+00 1.30862368e+00]
[ 1.00344906e+01 6.96091443e+00]
[ 1.89998430e+00 1.55460747e+00]
[ 7.00270854e+00 2.91919965e+00]
[ 3.40598744e+00 3.33335666e+00]
[-2.84451575e-01 -5.14498298e-01]
[ 1.80459108e+00 5.75090634e-01]
[ 5.17975110e+00 4.17678917e+00]
[ 3.38925189e+00 2.29768250e+00]
[-6.80226625e-01 -2.57026989e+00]
[ 1.85459196e+00 2.69966047e+00]
[ 7.51512076e+00 4.49278295e+00]
[ 4.91960836e+00 2.00141816e+00]
[ 1.14246901e-01 6.34410412e-01]]
Transformations of low discrepancy sequences often have good properties, but may not be low discrepancy themselves, depending on your definition of discrepancy as shown by Yiou Li and Lulu Kang.
Pause for questions
Integration¶
Cubature—the approximation of multivariate integrals—is an important application area for QMC. To solve this problem we need the Integrand and StoppingCriterion classes.
Consider the followng $d$-variate integral due to Keister: \begin{equation*} \mu = \int_{\mathbb{R}^d} \cos(\lVert\boldsymbol{t}\rVert) \exp( - \lVert \boldsymbol{t} \rVert^2) \, \mathrm{d} \boldsymbol{t}, \end{equation*} where $\lVert \cdot \rVert$ is the Euclidean norm. To approximate this integral via (Q)MC methods, it needs to be written as the integral with respect to a probability measure, e.g., \begin{equation*} \mu = \int_{\mathbb{R}^d} \underbrace{\pi^{d/2} \cos(\lVert\boldsymbol{t}\rVert)}_{g(\boldsymbol{t})} \; \underbrace{\pi^{-d/2} \exp( - \lVert \boldsymbol{t} \rVert^2) \, \mathrm{d} \boldsymbol{t}}_{\mathcal{N}(\boldsymbol{0}_d,\mathsf{I}_d/2) \text{ measure}}. \end{equation*} Using transformation techniques highlighted above, this integral can be further transformed to an integral over the unit cube, which is suitable for certain stopping criteria: \begin{equation*} \mu = \int_{[0,1]^d} \underbrace{\pi^{d/2} \cos\left(\sqrt{ \frac 12 \sum_{j=1}^d \bigl[\Phi^{-1}(x_j)\bigr]^2}\right)}_{f(\boldsymbol{x})} \, \rm d \boldsymbol{x}. \end{equation*}
Although it may seem counter-intuitive, we set up our numerical problem by
- first choosing the
DiscreteDistributionobject, - next the
TrueMeasureobject, - thirdly choosing
Integrandobject, and - finally the
StoppingCriterionobject.
In [12]:
Copied!
d = 5 #coded as parameters so that
tol = 1e-3 #you can change here and propagate them through this example
lattice = qmcpy.Lattice(d)
gaussian_lattice = qmcpy.Gaussian(lattice, mean = 0, covariance = 1/2) #mean and covariance of the distribution identified above
d = 5 #coded as parameters so that
tol = 1e-3 #you can change here and propagate them through this example
lattice = qmcpy.Lattice(d)
gaussian_lattice = qmcpy.Gaussian(lattice, mean = 0, covariance = 1/2) #mean and covariance of the distribution identified above
Integrand objects¶
The TrueMeasure object becomes input to an Integrand object, which transforms our original integrand, $g$, to our eventual integrand, $f$. This transformation relies on the TrueMeasure object along with its corresponding DiscreteDistribution object. The object qmcpy.Keister has already been coded as a use case in qmcpy.
In [13]:
Copied!
keister = qmcpy.Keister(gaussian_lattice) #transform the original integrand to the eventual one
keister = qmcpy.Keister(gaussian_lattice) #transform the original integrand to the eventual one
StoppingCriterion objects and the integrate method¶
Determining the sample size needed requires a stopping criterion, which typically depends on the error tolerance, abs_tol. The stopping criterion attempts to produce the answer satisfying the error tolerance with not much more work than is truly needed. There are several StoppingCriterion objects available, but they tend to work for specific LD sequences. This one comes from Tony Jiménez Rugama. It takes as its input the Integrand object, which carries information about the TrueMeasure object and its DiscreteDistribution object.
In [14]:
Copied!
keister_lattice_gauss_g = qmcpy.CubQMCLatticeG(keister, abs_tol = tol) #using Tony's stopping criterion
keister_lattice_gauss_g = qmcpy.CubQMCLatticeG(keister, abs_tol = tol) #using Tony's stopping criterion
Invoking the integrate method returns the numerical solution and a data object. Printing the data object provides a neat summary of the integration problem. For details of the output fields, see the online, searchable QMCPy Documentation at https://qmcpy.readthedocs.io/.
In [15]:
Copied!
solution, data = keister_lattice_gauss_g.integrate()
print(data)
solution, data = keister_lattice_gauss_g.integrate()
print(data)
Data (Data)
solution 1.135
comb_bound_low 1.135
comb_bound_high 1.136
comb_bound_diff 0.001
comb_flags 1
n_total 2^(17)
n 2^(17)
time_integrate 0.042
CubQMCLatticeG (AbstractStoppingCriterion)
abs_tol 0.001
rel_tol 0
n_init 2^(10)
n_limit 2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 5
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 38563684764170563855999902143753109741
Fixed sample budget computation¶
If you are not concerned about meeting an error tolerance but can only afford n_max function values, then you can set an abs_tol small enough and set n_max to your desired sample size. You will get an error bound.
In [16]:
Copied!
warnings.simplefilter("ignore")
keister_lattice_gauss_g_small_n = qmcpy.CubQMCLatticeG(keister, abs_tol = 1e-6, n_limit = 2**12) #the default n_max is 2**35
solution, data = keister_lattice_gauss_g_small_n.integrate()
print(data)
warnings.simplefilter("ignore")
keister_lattice_gauss_g_small_n = qmcpy.CubQMCLatticeG(keister, abs_tol = 1e-6, n_limit = 2**12) #the default n_max is 2**35
solution, data = keister_lattice_gauss_g_small_n.integrate()
print(data)
Data (Data)
solution 1.129
comb_bound_low 1.116
comb_bound_high 1.142
comb_bound_diff 0.026
comb_flags 0
n_total 2^(12)
n 2^(12)
time_integrate 0.003
CubQMCLatticeG (AbstractStoppingCriterion)
abs_tol 1.00e-06
rel_tol 0
n_init 2^(10)
n_limit 2^(12)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 5
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 38563684764170563855999902143753109741
QMC is much faster than MC. (Q)MC cost is mostly dimension dependent¶
Run this next code block to see how the run time and the number of function evaluations increase as the tolerance decreases. QMC using LD sequences uses much less time and much fewer function values than MC using IID sequences.
Tensor product rules have a time or function value cost that is $\mathcal{O}(\varepsilon^{-d/r})$, where $\varepsilon$ is the error tolerance and $r$ is bounded above by both the smoothness of the integrand and the quality of the algorithm. Such rules have a curse of dimensionality because their cost blows up exponentially with dimension.
Unlike tensor product cubature rules, the cost of (Q)MC cubature is essentially dimension independent: $\mathcal{O}(\varepsilon^{-2})$ for IID MC and typically $\mathcal{O}(\varepsilon^{-1-\delta})$ for QMC. Although Q(MC) is not particularly fast, its performance usually does not degrade as the number of variables of in the integrand increases.
In [17]:
Copied!
d = 5; tol = 2.5e-3 #re-construct the example
#d = 7; tol = 3e-3 #if you change the dimension
#d = 10; tol = 3e-2 #you may also wish to change the tolerance, since the value of the integral changes
ld_keister = qmcpy.Keister(qmcpy.Gaussian(qmcpy.Lattice(d), mean = 0, covariance = 1/2)) #mean and covariance of the distribution identified above
iid_keister = qmcpy.Keister(qmcpy.Gaussian(qmcpy.IIDStdUniform(d), mean = 0, covariance = 1/2))
n_tol = 7
ii_iid = 2 #make this larger to reduce the time required
tol_vec = [tol*2**(ii) for ii in range(n_tol)] #initialize
ld_time = [0]*n_tol; ld_n = [0]*n_tol #low discrepancy time and number of function values
iid_time = [0]*n_tol; iid_n = [0]*n_tol #IID time and number of function values
for ii in range(n_tol):
solution, data = qmcpy.CubQMCLatticeG(ld_keister, abs_tol = tol_vec[ii]).integrate()
if ii == 0:
print(f'\nKeister integral = {solution}\n')
ld_time[ii] = data.time_integrate
ld_n[ii] = data.n_total
if ii >= ii_iid:
solution, data = qmcpy.CubMCG(iid_keister, abs_tol = tol_vec[ii]).integrate()
iid_time[ii] = data.time_integrate
iid_n[ii] = data.n_total
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(13,5.5))
ax[0].scatter(tol_vec[0:n_tol],ld_time[0:n_tol],color='b');
ax[0].plot(tol_vec[0:n_tol],[(ld_time[0]*tol_vec[0])/tol_vec[jj] for jj in range(n_tol)],color='b')
ax[0].scatter(tol_vec[ii_iid:n_tol],iid_time[ii_iid:n_tol],color='g');
ax[0].plot(tol_vec[ii_iid:n_tol],[(iid_time[ii_iid]*(tol_vec[ii_iid]**2))/(tol_vec[jj]**2) for jj in range(ii_iid,n_tol)],color='g')
ax[0].set_ylim([0.001,1000]); ax[0].set_ylabel('Time (s)')
ax[1].scatter(tol_vec[0:n_tol],ld_n[0:n_tol],color='b');
ax[1].plot(tol_vec[0:n_tol],[(ld_n[0]*tol_vec[0])/tol_vec[jj] for jj in range(n_tol)],color='b')
ax[1].scatter(tol_vec[ii_iid:n_tol],iid_n[ii_iid:n_tol],color='g');
ax[1].plot(tol_vec[ii_iid:n_tol],[(iid_n[ii_iid]*(tol_vec[ii_iid]**2))/(tol_vec[jj]**2) for jj in range(ii_iid,n_tol)],color='g')
ax[1].set_ylim([1e2,1e8]); ax[1].set_ylabel('n')
for ii in range(2):
ax[ii].set_xlim([tol,100*tol]); ax[ii].set_xlabel('Tolerance, '+r'$\varepsilon$')
ax[ii].set_xscale('log'); ax[ii].set_yscale('log')
ax[ii].legend([r'$\mathcal{O}(\varepsilon^{-1})$',r'$\mathcal{O}(\varepsilon^{-2})$','LD','IID'],frameon=False)
ax[ii].set_aspect(0.35)
d = 5; tol = 2.5e-3 #re-construct the example
#d = 7; tol = 3e-3 #if you change the dimension
#d = 10; tol = 3e-2 #you may also wish to change the tolerance, since the value of the integral changes
ld_keister = qmcpy.Keister(qmcpy.Gaussian(qmcpy.Lattice(d), mean = 0, covariance = 1/2)) #mean and covariance of the distribution identified above
iid_keister = qmcpy.Keister(qmcpy.Gaussian(qmcpy.IIDStdUniform(d), mean = 0, covariance = 1/2))
n_tol = 7
ii_iid = 2 #make this larger to reduce the time required
tol_vec = [tol*2**(ii) for ii in range(n_tol)] #initialize
ld_time = [0]*n_tol; ld_n = [0]*n_tol #low discrepancy time and number of function values
iid_time = [0]*n_tol; iid_n = [0]*n_tol #IID time and number of function values
for ii in range(n_tol):
solution, data = qmcpy.CubQMCLatticeG(ld_keister, abs_tol = tol_vec[ii]).integrate()
if ii == 0:
print(f'\nKeister integral = {solution}\n')
ld_time[ii] = data.time_integrate
ld_n[ii] = data.n_total
if ii >= ii_iid:
solution, data = qmcpy.CubMCG(iid_keister, abs_tol = tol_vec[ii]).integrate()
iid_time[ii] = data.time_integrate
iid_n[ii] = data.n_total
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(13,5.5))
ax[0].scatter(tol_vec[0:n_tol],ld_time[0:n_tol],color='b');
ax[0].plot(tol_vec[0:n_tol],[(ld_time[0]*tol_vec[0])/tol_vec[jj] for jj in range(n_tol)],color='b')
ax[0].scatter(tol_vec[ii_iid:n_tol],iid_time[ii_iid:n_tol],color='g');
ax[0].plot(tol_vec[ii_iid:n_tol],[(iid_time[ii_iid]*(tol_vec[ii_iid]**2))/(tol_vec[jj]**2) for jj in range(ii_iid,n_tol)],color='g')
ax[0].set_ylim([0.001,1000]); ax[0].set_ylabel('Time (s)')
ax[1].scatter(tol_vec[0:n_tol],ld_n[0:n_tol],color='b');
ax[1].plot(tol_vec[0:n_tol],[(ld_n[0]*tol_vec[0])/tol_vec[jj] for jj in range(n_tol)],color='b')
ax[1].scatter(tol_vec[ii_iid:n_tol],iid_n[ii_iid:n_tol],color='g');
ax[1].plot(tol_vec[ii_iid:n_tol],[(iid_n[ii_iid]*(tol_vec[ii_iid]**2))/(tol_vec[jj]**2) for jj in range(ii_iid,n_tol)],color='g')
ax[1].set_ylim([1e2,1e8]); ax[1].set_ylabel('n')
for ii in range(2):
ax[ii].set_xlim([tol,100*tol]); ax[ii].set_xlabel('Tolerance, '+r'$\varepsilon$')
ax[ii].set_xscale('log'); ax[ii].set_yscale('log')
ax[ii].legend([r'$\mathcal{O}(\varepsilon^{-1})$',r'$\mathcal{O}(\varepsilon^{-2})$','LD','IID'],frameon=False)
ax[ii].set_aspect(0.35)
Keister integral = 1.1356452963607024
Alternatives for StoppingCriterion¶
Other StoppingCriterion objects are available. Most are tied to particular DiscreteDistribution objects. For LD points one can use replications and the Central Limit Theorem (CLT).
In [18]:
Copied!
lattice_rep = qmcpy.Lattice(d,replications=15)
gaussian_lattice_rep = qmcpy.Gaussian(lattice_rep, mean = 0, covariance = 1/2)
keister_rep = qmcpy.Keister(gaussian_lattice_rep)
keister_lattice_gauss_CLT = qmcpy.CubQMCCLT(keister_rep, abs_tol = tol) #using a CLT stopping criterion with random replications
solution, data = keister_lattice_gauss_CLT.integrate()
print(data)
lattice_rep = qmcpy.Lattice(d,replications=15)
gaussian_lattice_rep = qmcpy.Gaussian(lattice_rep, mean = 0, covariance = 1/2)
keister_rep = qmcpy.Keister(gaussian_lattice_rep)
keister_lattice_gauss_CLT = qmcpy.CubQMCCLT(keister_rep, abs_tol = tol) #using a CLT stopping criterion with random replications
solution, data = keister_lattice_gauss_CLT.integrate()
print(data)
Data (Data)
solution 1.135
comb_bound_low 1.134
comb_bound_high 1.136
comb_bound_diff 0.002
comb_flags 1
n_total 122880
n 122880
n_rep 2^(13)
time_integrate 0.029
CubQMCRepStudentT (AbstractStoppingCriterion)
inflate 1
alpha 0.010
abs_tol 0.003
rel_tol 0
n_init 2^(8)
n_limit 2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 5
replications 15
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 275034244784380212109446453686580990022
This answer agrees with the one above.
Alternatives for TrueMeasure¶
The Keister integrand may also be solved using a Lebesgue TrueMeasure object:
$$\mu = \int_{\mathbb{R}^d} \underbrace{\cos(\lVert\boldsymbol{t}\rVert) \exp( - \lVert\boldsymbol{t}\rVert^2)}_{g(\boldsymbol{t})} \, \underbrace{\mathrm{d} \boldsymbol{t}}_{\text{Lebesgue measure}}$$
The TrueMeasure object contains the appropriate information so that the Integrand object can obtain the correct eventual integrand, $f$, in terms of the original integrand, $g$.
Since TrueMeasure is not a probability measure, so it cannot be mimicked, but it can be used to solve the integration problem. The Integrand object maps the problem using an affine transformation for finite boxes and an inverse normal distribution function transformation for $\mathbb{R}^d$.
The code below also shows how to take our own integrand defined with a simple input and output, and turn it into a qmcpy ready integrand using the qmcpy.CustomFun object.
In [19]:
Copied!
def my_Keister(x): #this could be a functional of a solution to a PDE with random coefficients
#or anything that you would like
"""
x: nxd numpy ndarray
n samples
d dimensions
returns n-vector of the Kesiter function
evaluated at the n input samples
"""
d = x.shape[1]
norm = np.sqrt((x**2).sum(1))
k = np.cos(norm)*np.exp(-norm**2)
return k #size n vector
ld = qmcpy.Sobol(d) #choose the LD points
lebesgue = qmcpy.Lebesgue(qmcpy.Gaussian(ld)) #now choose the Lebesgue distribution
f = qmcpy.CustomFun(lebesgue, g=my_Keister)
keister_lebesgue_ld_g = qmcpy.CubQMCSobolG(f, abs_tol = tol) #the stopping criterion does need to match the LD points
solution, data = keister_lebesgue_ld_g.integrate()
print(data)
def my_Keister(x): #this could be a functional of a solution to a PDE with random coefficients
#or anything that you would like
"""
x: nxd numpy ndarray
n samples
d dimensions
returns n-vector of the Kesiter function
evaluated at the n input samples
"""
d = x.shape[1]
norm = np.sqrt((x**2).sum(1))
k = np.cos(norm)*np.exp(-norm**2)
return k #size n vector
ld = qmcpy.Sobol(d) #choose the LD points
lebesgue = qmcpy.Lebesgue(qmcpy.Gaussian(ld)) #now choose the Lebesgue distribution
f = qmcpy.CustomFun(lebesgue, g=my_Keister)
keister_lebesgue_ld_g = qmcpy.CubQMCSobolG(f, abs_tol = tol) #the stopping criterion does need to match the LD points
solution, data = keister_lebesgue_ld_g.integrate()
print(data)
Data (Data)
solution 1.136
comb_bound_low 1.133
comb_bound_high 1.138
comb_bound_diff 0.004
comb_flags 1
n_total 2^(16)
n 2^(16)
time_integrate 0.030
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.003
rel_tol 0
n_init 2^(10)
n_limit 2^(35)
CustomFun (AbstractIntegrand)
Lebesgue (AbstractTrueMeasure)
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 1
decomp_type PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 5
replications 1
randomize LMS DS
gen_mats_source joe_kuo.6.21201.txt
order RADICAL INVERSE
t 63
alpha 1
n_limit 2^(32)
entropy 91580573475794728146679265749810645603
This answer agrees with the answers above for the same integration problem.
The initial DiscreteDistribution does not need to mimic the standard uniform distribution as this next example shows.
In [20]:
Copied!
iid = qmcpy.IIDStdUniform(d) #choose the LD points
lebesgue = qmcpy.Lebesgue(qmcpy.Gaussian(iid)) #now choose the Lebesgue distribution
f = qmcpy.CustomFun(lebesgue, g=my_Keister)
keister_lebesgue_ld_g = qmcpy.CubMCCLT(f, abs_tol = 10*tol) #the stopping criterion does need to match the points
solution, data = keister_lebesgue_ld_g.integrate()
print(data)
iid = qmcpy.IIDStdUniform(d) #choose the LD points
lebesgue = qmcpy.Lebesgue(qmcpy.Gaussian(iid)) #now choose the Lebesgue distribution
f = qmcpy.CustomFun(lebesgue, g=my_Keister)
keister_lebesgue_ld_g = qmcpy.CubMCCLT(f, abs_tol = 10*tol) #the stopping criterion does need to match the points
solution, data = keister_lebesgue_ld_g.integrate()
print(data)
Data (Data)
solution 1.127
bound_low 1.099
bound_high 1.155
bound_diff 0.056
n_total 1694184
time_integrate 0.361
CubMCCLT (AbstractStoppingCriterion)
abs_tol 0.025
rel_tol 0
n_init 2^(10)
n_limit 2^(30)
inflate 1.200
alpha 0.010
CustomFun (AbstractIntegrand)
Lebesgue (AbstractTrueMeasure)
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 1
decomp_type PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
d 5
replications 1
entropy 110067297374675390351789060258764096093
Multi-level (Q)MC¶
When the dimesion of the multivariate integral is high, multi-level (Quasi-)Monte Carlo (ML(Q)MC) methods may save computation time. The cost of one integrand value depends on the number of input variables, $d$, which corresponds to the dimension of our integration problem. ML(Q)MC methods allow us attain our accuracy requirements by evaluating low dimensional integrands many times and high dimensional integrands much fewer times.
High or infinte dimesional integration problems arise when computing the expectations of quantities coming stochastic differential equations (SDEs). These problems arise in finance applications. The dimension of the integrand typically refers to the number of time steps used to discretize the SDE.
Here are some parameters for the Asian option examples below. Changing them here allows you to compare the run times of these examples in a fair way.
In [21]:
Copied!
abs_tol = .05 #a nickel
n_time_steps = 64 #being a power of 2 will help for multi-level
options = { #there should be nothing magic about these choices
'interest_rate': .05,
'volatility': .5,
'start_price': 30,
'strike_price': 30
}
abs_tol = .05 #a nickel
n_time_steps = 64 #being a power of 2 will help for multi-level
options = { #there should be nothing magic about these choices
'interest_rate': .05,
'volatility': .5,
'start_price': 30,
'strike_price': 30
}
Single Level MC¶
The vanilla way to solve this problem is IID Monte Carlo. The fair price of the option is an expectation or integral, and the dimension is the number of time steps of used to discretize the Brownian motion that drives the SDE describing the price of the underlying asset. The number of time steps should be fairly large.
This paper by Lan Jiang and collaborators describes the stopping criterion.
In [22]:
Copied!
iidBrownian = qmcpy.BrownianMotion(qmcpy.IIDStdUniform(n_time_steps))
payoff = qmcpy.AsianOption(iidBrownian, **options)
IIDstop = qmcpy.CubMCG(payoff,abs_tol=abs_tol)
price,data = IIDstop.integrate()
print(data)
iidBrownian = qmcpy.BrownianMotion(qmcpy.IIDStdUniform(n_time_steps))
payoff = qmcpy.AsianOption(iidBrownian, **options)
IIDstop = qmcpy.CubMCG(payoff,abs_tol=abs_tol)
price,data = IIDstop.integrate()
print(data)
Data (Data)
solution 3.686
bound_low 3.636
bound_high 3.736
bound_diff 0.100
n_total 225848
time_integrate 0.834
CubMCG (AbstractStoppingCriterion)
abs_tol 0.050
rel_tol 0
n_init 2^(10)
n_limit 2^(30)
inflate 1.200
alpha 0.010
kurtmax 1.478
AsianOption (AbstractIntegrand)
option ASIAN
call_put CALL
volatility 2^(-1)
start_price 30
strike_price 30
interest_rate 0.050
t_final 1
asian_mean ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
transform BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
d 2^(6)
replications 1
entropy 14445467058155448935378672762555839489
Adaptive Multilevel MC from Mike Giles¶
Mike Giles and his collaborators have developed several ML(Q)MC algorithms. The ML IID MC algorithm and stopping criterion implemented here are from this paper and this code. The answer is expected to be different than above, even with the same parameters, as the MLCallOptions uses a different discretization. The algorithm considers the SDE for logarithm of the stock price, which allows exact time stepping for constant interest rates and volatirilities, while Giles uses a Milstein discretization for the SDE for the stock price iteself.
In [23]:
Copied!
giles_MLMC_stop = qmcpy.CubMCML(payoff,abs_tol=abs_tol)
solution,data = giles_MLMC_stop.integrate()
print(data)
giles_MLMC_stop = qmcpy.CubMCML(payoff,abs_tol=abs_tol)
solution,data = giles_MLMC_stop.integrate()
print(data)
Data (Data)
solution 3.681
n_total 430326
levels 3
n_level [341502 61469 27319]
mean_level [3.589 0.071 0.021]
var_level [39.742 2.575 0.665]
cost_per_sample [2. 4. 8.]
alpha 1.769
beta 1.953
gamma 1.000
time_integrate 0.077
CubMLMC (AbstractStoppingCriterion)
rmse_tol 0.019
n_init 2^(8)
levels_min 2^(1)
levels_max 10
theta 2^(-1)
AsianOption (AbstractIntegrand)
option ASIAN
call_put CALL
volatility 2^(-1)
start_price 30
strike_price 30
interest_rate 0.050
t_final 1
asian_mean ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
transform BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
IIDStdUniform (AbstractIIDDiscreteDistribution)
d 2^(6)
replications 1
entropy 14445467058155448935378672762555839489
Single Level QMC Baseline¶
Tony Jiménez's stopping criterion for cubature via Sobol' sequences in this paper does not yet work for multi-level problems. Here it is treating the option pricing problem as a high dimensional integral.
In [24]:
Copied!
sobol_brownian = qmcpy.BrownianMotion(qmcpy.Sobol(n_time_steps))
integrand = qmcpy.AsianOption(sobol_brownian, **options)
stopping_criterion = qmcpy.CubQMCSobolG(integrand,abs_tol = abs_tol)
solution,data = stopping_criterion.integrate()
print(data)
sobol_brownian = qmcpy.BrownianMotion(qmcpy.Sobol(n_time_steps))
integrand = qmcpy.AsianOption(sobol_brownian, **options)
stopping_criterion = qmcpy.CubQMCSobolG(integrand,abs_tol = abs_tol)
solution,data = stopping_criterion.integrate()
print(data)
Data (Data)
solution 3.697
comb_bound_low 3.672
comb_bound_high 3.721
comb_bound_diff 0.049
comb_flags 1
n_total 2^(10)
n 2^(10)
time_integrate 0.003
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.050
rel_tol 0
n_init 2^(10)
n_limit 2^(35)
AsianOption (AbstractIntegrand)
option ASIAN
call_put CALL
volatility 2^(-1)
start_price 30
strike_price 30
interest_rate 0.050
t_final 1
asian_mean ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
transform BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 2^(6)
replications 1
randomize LMS DS
gen_mats_source joe_kuo.6.21201.txt
order RADICAL INVERSE
t 63
alpha 1
n_limit 2^(32)
entropy 267658111069601541207388819540918772545
Importance Sampling¶
For the option pricing problem, we may add a drift to the Brownian motion as an example of importance sampling.
First, we need to change our problem to one for which importance sampling can show some benefit. We consider an out-of-the-money call option.
In [25]:
Copied!
abs_tol = .001
n_time_steps = 64
options = {
'interest_rate': .05,
'volatility': .5,
'start_price': 30,
'strike_price': 40 #a larger strike price than before
}
abs_tol = .001
n_time_steps = 64
options = {
'interest_rate': .05,
'volatility': .5,
'start_price': 30,
'strike_price': 40 #a larger strike price than before
}
First we price it as above using single level QMC.
In [26]:
Copied!
sobol_brownian = qmcpy.BrownianMotion(qmcpy.Sobol(n_time_steps))
integrand = qmcpy.AsianOption(sobol_brownian, **options)
stopping_criterion = qmcpy.CubQMCSobolG(integrand,abs_tol = abs_tol)
solution,data = stopping_criterion.integrate()
print(data)
sobol_brownian = qmcpy.BrownianMotion(qmcpy.Sobol(n_time_steps))
integrand = qmcpy.AsianOption(sobol_brownian, **options)
stopping_criterion = qmcpy.CubQMCSobolG(integrand,abs_tol = abs_tol)
solution,data = stopping_criterion.integrate()
print(data)
Data (Data)
solution 1.018
comb_bound_low 1.017
comb_bound_high 1.019
comb_bound_diff 0.002
comb_flags 1
n_total 2^(15)
n 2^(15)
time_integrate 0.118
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.001
rel_tol 0
n_init 2^(10)
n_limit 2^(35)
AsianOption (AbstractIntegrand)
option ASIAN
call_put CALL
volatility 2^(-1)
start_price 30
strike_price 40
interest_rate 0.050
t_final 1
asian_mean ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
transform BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 2^(6)
replications 1
randomize LMS DS
gen_mats_source joe_kuo.6.21201.txt
order RADICAL INVERSE
t 63
alpha 1
n_limit 2^(32)
entropy 68741035845447190293894879002444698271
Next, we introduce an upward drift in the Brownian motion, which produces more stock price paths with positive payoffs. This produces a smaller varation in the integrand and a generally faster run time. (There still remains the question of how to automatically choose an optimal drift.)
In [27]:
Copied!
sobol_drift_brownian = qmcpy.BrownianMotion(qmcpy.Sobol(n_time_steps), drift = 0.5)
integrand = qmcpy.AsianOption(sobol_drift_brownian, **options)
stopping_criterion = qmcpy.CubQMCSobolG(integrand,abs_tol = abs_tol)
solution,data = stopping_criterion.integrate()
print(data)
sobol_drift_brownian = qmcpy.BrownianMotion(qmcpy.Sobol(n_time_steps), drift = 0.5)
integrand = qmcpy.AsianOption(sobol_drift_brownian, **options)
stopping_criterion = qmcpy.CubQMCSobolG(integrand,abs_tol = abs_tol)
solution,data = stopping_criterion.integrate()
print(data)
Data (Data)
solution 1.018
comb_bound_low 1.018
comb_bound_high 1.019
comb_bound_diff 0.001
comb_flags 1
n_total 2^(15)
n 2^(15)
time_integrate 0.115
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.001
rel_tol 0
n_init 2^(10)
n_limit 2^(35)
AsianOption (AbstractIntegrand)
option ASIAN
call_put CALL
volatility 2^(-1)
start_price 30
strike_price 40
interest_rate 0.050
t_final 1
asian_mean ARITHMETIC
BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 0
mean [0. 0. 0. ... 0. 0. 0.]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
transform BrownianMotion (AbstractTrueMeasure)
time_vec [0.016 0.031 0.047 ... 0.969 0.984 1. ]
drift 2^(-1)
mean [0.008 0.016 0.023 ... 0.484 0.492 0.5 ]
covariance [[0.016 0.016 0.016 ... 0.016 0.016 0.016]
[0.016 0.031 0.031 ... 0.031 0.031 0.031]
[0.016 0.031 0.047 ... 0.047 0.047 0.047]
...
[0.016 0.031 0.047 ... 0.969 0.969 0.969]
[0.016 0.031 0.047 ... 0.969 0.984 0.984]
[0.016 0.031 0.047 ... 0.969 0.984 1. ]]
decomp_type PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 2^(6)
replications 1
randomize LMS DS
gen_mats_source joe_kuo.6.21201.txt
order RADICAL INVERSE
t 63
alpha 1
n_limit 2^(32)
entropy 195386499743864987571337198121595404839
Pause for questions
Under the hood¶
The structure of qmcpy is that there are five major classes:
DiscreteDistributionused to generate LD sequences, primarily on $[0,1]^d$TrueMeasurefor using these LD sequences to mimic other distributions and to define integrals with respect to other measuresIntegrandto define the integrand for the multivariate integration problemsStoppingCriterionto determine when the desired accuracy has been reachedAccumulateDatathe invisible class used to keep track of important data as you continue to sample
We look at some of the important parameters that the corresponding objects have.
LD sequence generators¶
The LD generators (DiscreteDistribution objects) implemented here are drawn from several sources, which are denoted backend (first listed is the default):
- Sobol: QRNG, MPS, & PyTorch
- Lattice: GAIL & MPS, with default generating vectors from Lattice Builder
- Halton: Art Owen's & QRNG
- Korobov: QRNG
We illustrate some of the features of these varous backends and some of the other parameters that you can set.
Different lattice backends and generators¶
The qmcpy.Lattice generator using the GAIL and MPS backends with the same generating vectors yield the same points but in a different order. The stopping criterion qmcpy.CubLatticeG requires the GAIL order, but not all stopping criteria do.
In [28]:
Copied!
d=4; n=8
x_gail = qmcpy.Lattice(d,order='linear',randomize=False).gen_samples(n,warn=False)
print('GAIL Samples')
for i in range(n): print(x_gail[i])
x_mps = qmcpy.Lattice(d,order='natural',randomize=False).gen_samples(n,warn=False)
print('\n\nMPS Samples')
for i in range(n): print(x_mps[i])
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(10,5.5))
ax[0].scatter(x_gail[0:n,0],x_gail[0:n,1],color='b')
ax[0].set_title('GAIL backend')
ax[1].scatter(x_mps[0:n,0],x_mps[0:n,1],color='b')
ax[1].set_title('MPS backend')
for ii in range(2):
ax[ii].set_xlim([0,1]); ax[ii].set_xticks([0,1]); ax[ii].set_xlabel('$x_{i,1}$')
ax[ii].set_ylim([0,1]); ax[ii].set_yticks([0,1]); ax[ii].set_ylabel('$x_{i,2}$')
ax[ii].set_aspect(1)
d=4; n=8
x_gail = qmcpy.Lattice(d,order='linear',randomize=False).gen_samples(n,warn=False)
print('GAIL Samples')
for i in range(n): print(x_gail[i])
x_mps = qmcpy.Lattice(d,order='natural',randomize=False).gen_samples(n,warn=False)
print('\n\nMPS Samples')
for i in range(n): print(x_mps[i])
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(10,5.5))
ax[0].scatter(x_gail[0:n,0],x_gail[0:n,1],color='b')
ax[0].set_title('GAIL backend')
ax[1].scatter(x_mps[0:n,0],x_mps[0:n,1],color='b')
ax[1].set_title('MPS backend')
for ii in range(2):
ax[ii].set_xlim([0,1]); ax[ii].set_xticks([0,1]); ax[ii].set_xlabel('$x_{i,1}$')
ax[ii].set_ylim([0,1]); ax[ii].set_yticks([0,1]); ax[ii].set_ylabel('$x_{i,2}$')
ax[ii].set_aspect(1)
GAIL Samples [0. 0. 0. 0.] [0.125 0.375 0.375 0.875] [0.25 0.75 0.75 0.75] [0.375 0.125 0.125 0.625] [0.5 0.5 0.5 0.5] [0.625 0.875 0.875 0.375] [0.75 0.25 0.25 0.25] [0.875 0.625 0.625 0.125] MPS Samples [0. 0. 0. 0.] [0.5 0.5 0.5 0.5] [0.25 0.75 0.75 0.75] [0.75 0.25 0.25 0.25] [0.125 0.375 0.375 0.875] [0.625 0.875 0.875 0.375] [0.375 0.125 0.125 0.625] [0.875 0.625 0.625 0.125]
Default LDs are randomized¶
LDs can be purely deterministic or they can be randomized. For lattices this corresponds to a random shift modulo one. For Sobol' sequences this corresponds to a random digital shift. (PyTorch also uses random linear scrambling.)
This randomization is turned on by default, but can also be turned off. With randomization off, the points will always look the same. Below the first sequence of points is randomized and so is different every time the code is run. The second sequence is not randomized and stay the same.
Turning off randomization throws a warning, which stops execution in Colab so we disable the warnings.
In [29]:
Copied!
warnings.simplefilter('ignore') #turn off warnings which stop execution in Colab
n = 16
ldA = qmcpy.Lattice(2)
ldB = qmcpy.Lattice(2, randomize = False)
pointsA = ldA.gen_samples(n) #construct some points
pointsB = ldB.gen_samples(n) #construct some points
print(f'\nRandomized LD Points with shape {pointsA.shape}\n'+str(pointsA)) #these points have 15 significant digit precision but only three digits are shown
print(f'\nNonrandomized LD Points with shape {pointsB.shape}\n'+str(pointsB)) #these points have 15 significant digit precision but only three digits are shown
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(10,5.5))
ax[0].scatter(pointsA[0:n,0],pointsA[0:n,1],color='b')
ax[1].scatter(pointsB[0:n,0],pointsB[0:n,1],color='b')
for ii in range(2):
ax[ii].set_xlim([0,1]); ax[ii].set_xticks([0,1]); ax[ii].set_xlabel('$x_{i,1}$')
ax[ii].set_ylim([0,1]); ax[ii].set_yticks([0,1]); ax[ii].set_ylabel('$x_{i,2}$')
ax[ii].set_aspect(1)
warnings.simplefilter('ignore') #turn off warnings which stop execution in Colab
n = 16
ldA = qmcpy.Lattice(2)
ldB = qmcpy.Lattice(2, randomize = False)
pointsA = ldA.gen_samples(n) #construct some points
pointsB = ldB.gen_samples(n) #construct some points
print(f'\nRandomized LD Points with shape {pointsA.shape}\n'+str(pointsA)) #these points have 15 significant digit precision but only three digits are shown
print(f'\nNonrandomized LD Points with shape {pointsB.shape}\n'+str(pointsB)) #these points have 15 significant digit precision but only three digits are shown
fig,ax = pyplot.subplots(nrows=1,ncols=2,figsize=(10,5.5))
ax[0].scatter(pointsA[0:n,0],pointsA[0:n,1],color='b')
ax[1].scatter(pointsB[0:n,0],pointsB[0:n,1],color='b')
for ii in range(2):
ax[ii].set_xlim([0,1]); ax[ii].set_xticks([0,1]); ax[ii].set_xlabel('$x_{i,1}$')
ax[ii].set_ylim([0,1]); ax[ii].set_yticks([0,1]); ax[ii].set_ylabel('$x_{i,2}$')
ax[ii].set_aspect(1)
Randomized LD Points with shape (16, 2) [[0.00944244 0.29799738] [0.50944244 0.79799738] [0.25944244 0.04799738] [0.75944244 0.54799738] [0.13444244 0.67299738] [0.63444244 0.17299738] [0.38444244 0.42299738] [0.88444244 0.92299738] [0.07194244 0.98549738] [0.57194244 0.48549738] [0.32194244 0.73549738] [0.82194244 0.23549738] [0.19694244 0.36049738] [0.69694244 0.86049738] [0.44694244 0.11049738] [0.94694244 0.61049738]] Nonrandomized LD Points with shape (16, 2) [[0. 0. ] [0.5 0.5 ] [0.25 0.75 ] [0.75 0.25 ] [0.125 0.375 ] [0.625 0.875 ] [0.375 0.125 ] [0.875 0.625 ] [0.0625 0.6875] [0.5625 0.1875] [0.3125 0.4375] [0.8125 0.9375] [0.1875 0.0625] [0.6875 0.5625] [0.4375 0.8125] [0.9375 0.3125]]
Unrandomized LDs start with $\boldsymbol{0}$¶
By definition, without randomization the first point in all the popluar LD sequences is the origin. You can try.
In [30]:
Copied!
warnings.simplefilter('ignore') #turn off warnings which stop execution in Colab
ld = qmcpy.Lattice(6,randomize=False)
points = ld.gen_samples(4) #construct some points
print(f'\nLD Points with shape {points.shape}\n'+str(points)) #these points have 15 significant digit precision but only three digits are shown
warnings.simplefilter('ignore') #turn off warnings which stop execution in Colab
ld = qmcpy.Lattice(6,randomize=False)
points = ld.gen_samples(4) #construct some points
print(f'\nLD Points with shape {points.shape}\n'+str(points)) #these points have 15 significant digit precision but only three digits are shown
LD Points with shape (4, 6) [[0. 0. 0. 0. 0. 0. ] [0.5 0.5 0.5 0.5 0.5 0.5 ] [0.25 0.75 0.75 0.75 0.25 0.75] [0.75 0.25 0.25 0.25 0.75 0.25]]
Coordinate values of zero and one may be problemetic if transformed to $\mathbb{R}^d$¶
If a coordinate value of a point constructed on $[0,1]^d$ is zero or one, then then $0$ is mapped to $-\infty$ and $1$ is mapped to $\infty$ when these points are transformed to $\mathbb{R}^d$, as is done when solving problems with Gaussian distributions. In Python, these may turn out to be nan. For many applications, infinities or nan will be troublesome.
In [31]:
Copied!
warnings.simplefilter('ignore') #turn off warnings which stop execution in Colab
ld = qmcpy.Lattice(6,randomize=False)
unif_points = ld.gen_samples(2) #construct some points
print(f'\nLD Points with shape {unif_points.shape}\n'+str(unif_points))
ld_gauss = qmcpy.Gaussian(ld)
gauss_points = ld_gauss.gen_samples(2)
print(f'\nLD Points with shape {gauss_points.shape}\n'+str(gauss_points))
warnings.simplefilter('ignore') #turn off warnings which stop execution in Colab
ld = qmcpy.Lattice(6,randomize=False)
unif_points = ld.gen_samples(2) #construct some points
print(f'\nLD Points with shape {unif_points.shape}\n'+str(unif_points))
ld_gauss = qmcpy.Gaussian(ld)
gauss_points = ld_gauss.gen_samples(2)
print(f'\nLD Points with shape {gauss_points.shape}\n'+str(gauss_points))
LD Points with shape (2, 6) [[0. 0. 0. 0. 0. 0. ] [0.5 0.5 0.5 0.5 0.5 0.5]] LD Points with shape (2, 6) [[nan nan nan nan nan nan] [ 0. 0. 0. 0. 0. 0.]]
Replications with a fixed seed¶
For debugging purposes, you may wish to fix the seed of your randomized points. This gives the advantage of an unchanging answer while avoiding the boundaries of the unit cube. When rerunning the code below, the answers are unchanged iff the seed is fixed. For Tony's stopping criterion from GAIL, the points are randomized, but the algorithm is deterministic. For the fixed-multilevel CLT stopping criterion, the points and the algorithm are both random.
In [32]:
Copied!
solution_G_fixed = qmcpy.CubQMCSobolG(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Sobol(2,seed=47),covariance=1/2))).integrate()[0]
solution_CLT_fixed = qmcpy.CubQMCCLT(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Sobol(2,seed=47,replications=15),covariance=1/2))).integrate()[0]
solution_G = qmcpy.CubQMCSobolG(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Sobol(2),covariance=1/2))).integrate()[0]
solution_CLT = qmcpy.CubQMCCLT(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Sobol(2,replications=15),covariance=1/2))).integrate()[0]
print(f'Solution for GAIL stopping criterion with FIXED seed {solution_G_fixed}')
print(f'Solution for CLT stopping criterion with FIXED seed {solution_CLT_fixed}')
print(f'Solution for GAIL stopping criterion {solution_G}')
print(f'Solution for CLT stopping criterion {solution_CLT}')
solution_G_fixed = qmcpy.CubQMCSobolG(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Sobol(2,seed=47),covariance=1/2))).integrate()[0]
solution_CLT_fixed = qmcpy.CubQMCCLT(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Sobol(2,seed=47,replications=15),covariance=1/2))).integrate()[0]
solution_G = qmcpy.CubQMCSobolG(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Sobol(2),covariance=1/2))).integrate()[0]
solution_CLT = qmcpy.CubQMCCLT(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Sobol(2,replications=15),covariance=1/2))).integrate()[0]
print(f'Solution for GAIL stopping criterion with FIXED seed {solution_G_fixed}')
print(f'Solution for CLT stopping criterion with FIXED seed {solution_CLT_fixed}')
print(f'Solution for GAIL stopping criterion {solution_G}')
print(f'Solution for CLT stopping criterion {solution_CLT}')
Solution for GAIL stopping criterion with FIXED seed 1.8088685362448778 Solution for CLT stopping criterion with FIXED seed 1.806594507374441 Solution for GAIL stopping criterion 1.807944316997558 Solution for CLT stopping criterion 1.80672801265591
Help¶
You can obtain help on any object by the help(...) or dir(...) commands.
In [33]:
Copied!
help(qmcpy.Sobol)
dir(qmcpy.Sobol)
help(qmcpy.Sobol)
dir(qmcpy.Sobol)
Help on class DigitalNetB2 in module qmcpy.discrete_distribution.digital_net_b2.digital_net_b2:
class DigitalNetB2(qmcpy.discrete_distribution.abstract_discrete_distribution.AbstractLDDiscreteDistribution)
| DigitalNetB2(dimension=1, replications=None, seed=None, randomize='LMS DS', generating_matrices='joe_kuo.6.21201.txt', order='RADICAL INVERSE', t=63, alpha=1, msb=None, _verbose=False, graycode=None, t_max=None, t_lms=None)
|
| Low discrepancy digital net in base 2.
|
| Note:
| - Digital net sample sizes should be powers of $2$ e.g. $1$, $2$, $4$, $8$, $16$, $\dots$.
| - The first point of an unrandomized digital nets is the origin.
| - `Sobol` is an alias for `DigitalNetB2`.
| - To use higher order digital nets, either:
|
| - Pass in `generating_matrices` *without* interlacing and supply `alpha`>1 to apply interlacing, or
| - Pass in `generating_matrices` *with* interlacing and set `alpha=1` to avoid additional interlacing
|
| i.e. do *not* pass in interlaced `generating_matrices` and set `alpha>1`, this will apply additional interlacing.
|
| Examples:
| >>> discrete_distrib = DigitalNetB2(2,seed=7)
| >>> discrete_distrib(4)
| array([[0.72162356, 0.914955 ],
| [0.16345554, 0.42964856],
| [0.98676255, 0.03436384],
| [0.42956655, 0.55876342]])
| >>> discrete_distrib(1) # first point in the sequence
| array([[0.72162356, 0.914955 ]])
| >>> discrete_distrib
| DigitalNetB2 (AbstractLDDiscreteDistribution)
| d 2^(1)
| replications 1
| randomize LMS DS
| gen_mats_source joe_kuo.6.21201.txt
| order RADICAL INVERSE
| t 63
| alpha 1
| n_limit 2^(32)
| entropy 7
|
| Replications of independent randomizations
|
| >>> x = DigitalNetB2(dimension=3,seed=7,replications=2)(4)
| >>> x.shape
| (2, 4, 3)
| >>> x
| array([[[0.24653277, 0.1821862 , 0.74732591],
| [0.68152903, 0.66169442, 0.42891961],
| [0.48139855, 0.79818233, 0.08201287],
| [0.91541325, 0.29520621, 0.77495809]],
| <BLANKLINE>
| [[0.44876891, 0.85899604, 0.50549679],
| [0.53635924, 0.04353443, 0.33564946],
| [0.23214143, 0.29281506, 0.06841036],
| [0.75295715, 0.60241448, 0.76962976]]])
|
| Different orderings (avoid warnings that the first point is the origin)
|
| >>> DigitalNetB2(dimension=2,randomize=False,order="GRAY")(n_min=2,n_max=4,warn=False)
| array([[0.75, 0.25],
| [0.25, 0.75]])
| >>> DigitalNetB2(dimension=2,randomize=False,order="RADICAL INVERSE")(n_min=2,n_max=4,warn=False)
| array([[0.25, 0.75],
| [0.75, 0.25]])
|
| Generating matrices from [https://github.com/QMCSoftware/LDData/tree/main/dnet](https://github.com/QMCSoftware/LDData/tree/main/dnet)
|
| >>> DigitalNetB2(dimension=3,randomize=False,generating_matrices="mps.nx_s5_alpha2_m32.txt")(8,warn=False)
| array([[0. , 0. , 0. ],
| [0.75841841, 0.45284834, 0.48844557],
| [0.57679828, 0.13226272, 0.10061957],
| [0.31858402, 0.32113875, 0.39369111],
| [0.90278927, 0.45867532, 0.01803333],
| [0.14542431, 0.02548793, 0.4749614 ],
| [0.45587539, 0.33081476, 0.11474426],
| [0.71318879, 0.15377192, 0.37629925]])
|
| All randomizations
|
| >>> DigitalNetB2(dimension=3,randomize='LMS DS',seed=5)(8)
| array([[0.69346401, 0.20118185, 0.64779396],
| [0.43998032, 0.90102467, 0.0936172 ],
| [0.86663563, 0.60910036, 0.26043276],
| [0.11327376, 0.30772653, 0.93959283],
| [0.62102883, 0.79169756, 0.77051637],
| [0.37451038, 0.1231324 , 0.46634012],
| [0.94785596, 0.38577413, 0.13377215],
| [0.20121617, 0.71843325, 0.56293458]])
| >>> DigitalNetB2(dimension=3,randomize='LMS',seed=5)(8,warn=False)
| array([[0. , 0. , 0. ],
| [0.75446077, 0.83265937, 0.69584079],
| [0.42329494, 0.65793842, 0.90427279],
| [0.67763292, 0.48937304, 0.33344964],
| [0.18550714, 0.97332905, 0.3772791 ],
| [0.93104851, 0.17195496, 0.82311652],
| [0.26221346, 0.31742386, 0.53093284],
| [0.50787715, 0.5172669 , 0.2101083 ]])
| >>> DigitalNetB2(dimension=3,randomize='DS',seed=5)(8)
| array([[0.68383949, 0.04047995, 0.42903182],
| [0.18383949, 0.54047995, 0.92903182],
| [0.93383949, 0.79047995, 0.67903182],
| [0.43383949, 0.29047995, 0.17903182],
| [0.55883949, 0.66547995, 0.05403182],
| [0.05883949, 0.16547995, 0.55403182],
| [0.80883949, 0.41547995, 0.80403182],
| [0.30883949, 0.91547995, 0.30403182]])
| >>> DigitalNetB2(dimension=3,randomize='OWEN',seed=5)(8)
| array([[0.33595486, 0.05834975, 0.30066401],
| [0.89110875, 0.84905188, 0.81833285],
| [0.06846074, 0.59997956, 0.67064205],
| [0.6693703 , 0.25824002, 0.10469644],
| [0.44586618, 0.99161977, 0.1873488 ],
| [0.84245267, 0.16445553, 0.56544372],
| [0.18546359, 0.44859876, 0.97389524],
| [0.61215442, 0.64341386, 0.44529863]])
|
| Higher order net without randomization
|
| >>> DigitalNetB2(dimension=3,randomize='FALSE',seed=7,alpha=2)(4,warn=False)
| array([[0. , 0. , 0. ],
| [0.75 , 0.75 , 0.75 ],
| [0.4375, 0.9375, 0.1875],
| [0.6875, 0.1875, 0.9375]])
|
|
| Higher order nets with randomizations and replications
|
| >>> DigitalNetB2(dimension=3,randomize='LMS DS',seed=7,replications=2,alpha=2)(4,warn=False)
| array([[[0.42955149, 0.89149058, 0.43867111],
| [0.68701828, 0.07601148, 0.51312447],
| [0.10088033, 0.16293661, 0.25144138],
| [0.85846252, 0.87103178, 0.70041789]],
| <BLANKLINE>
| [[0.27151905, 0.42406763, 0.21917369],
| [0.55035224, 0.67864387, 0.90033876],
| [0.19356758, 0.57589964, 0.00347701],
| [0.97235125, 0.32168581, 0.86920948]]])
| >>> DigitalNetB2(dimension=3,randomize='LMS',seed=7,replications=2,alpha=2)(4,warn=False)
| array([[[0. , 0. , 0. ],
| [0.75817062, 0.96603053, 0.94947625],
| [0.45367986, 0.80295638, 0.18778553],
| [0.71171791, 0.2295424 , 0.76175441]],
| <BLANKLINE>
| [[0. , 0. , 0. ],
| [0.78664636, 0.75470215, 0.86876474],
| [0.45336727, 0.99953621, 0.22253579],
| [0.73996397, 0.24544824, 0.9008679 ]]])
| >>> DigitalNetB2(dimension=3,randomize='DS',seed=7,replications=2,alpha=2)(4)
| array([[[0.04386058, 0.58727432, 0.3691824 ],
| [0.79386058, 0.33727432, 0.6191824 ],
| [0.48136058, 0.39977432, 0.4316824 ],
| [0.73136058, 0.64977432, 0.6816824 ]],
| <BLANKLINE>
| [[0.65212985, 0.69669968, 0.10605352],
| [0.40212985, 0.44669968, 0.85605352],
| [0.83962985, 0.25919968, 0.16855352],
| [0.08962985, 0.50919968, 0.91855352]]])
| >>> DigitalNetB2(dimension=3,randomize='OWEN',seed=7,replications=2,alpha=2)(4)
| array([[[0.46368517, 0.03964427, 0.62172094],
| [0.7498683 , 0.76141348, 0.4243043 ],
| [0.01729754, 0.97968459, 0.65963223],
| [0.75365329, 0.1903774 , 0.34141493]],
| <BLANKLINE>
| [[0.52252547, 0.5679709 , 0.05949112],
| [0.27248656, 0.36488289, 0.81844058],
| [0.94219959, 0.39172304, 0.20285965],
| [0.19716391, 0.64741585, 0.92494554]]])
|
| **References:**
|
| 1. Marius Hofert and Christiane Lemieux.
| qrng: (Randomized) Quasi-Random Number Generators (2019).
| R package version 0.0-7.
| [https://CRAN.R-project.org/package=qrng](https://CRAN.R-project.org/package=qrng).
|
| 2. Faure, Henri, and Christiane Lemieux.
| Implementation of Irreducible Sobol' Sequences in Prime Power Bases.
| Mathematics and Computers in Simulation 161 (2019): 13-22. Crossref. Web.
|
| 3. F.Y. Kuo, D. Nuyens.
| Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients \- a survey of analysis and implementation.
| Foundations of Computational Mathematics, 16(6):1631-1696, 2016.
| [https://link.springer.com/article/10.1007/s10208-016-9329-5](https://link.springer.com/article/10.1007/s10208-016-9329-5).
|
| 4. D. Nuyens.
| The Magic Point Shop of QMC point generators and generating vectors.
| MATLAB and Python software, 2018.
| [https://people.cs.kuleuven.be/~dirk.nuyens/](https://people.cs.kuleuven.be/~dirk.nuyens/).
|
| 5. R. Cools, F.Y. Kuo, D. Nuyens.
| Constructing embedded lattice rules for multivariate integration.
| SIAM J. Sci. Comput., 28(6), 2162-2188.
|
| 6. I.M. Sobol', V.I. Turchaninov, Yu.L. Levitan, B.V. Shukhman.
| Quasi-Random Sequence Generators.
| Keldysh Institute of Applied Mathematics.
| Russian Academy of Sciences, Moscow (1992).
|
| 7. Sobol, Ilya & Asotsky, Danil & Kreinin, Alexander & Kucherenko, Sergei. (2011).
| Construction and Comparison of High-Dimensional Sobol' Generators. Wilmott. 2011.
| [10.1002/wilm.10056](https://onlinelibrary.wiley.com/doi/abs/10.1002/wilm.10056).
|
| 8. Paul Bratley and Bennett L. Fox.
| Algorithm 659: Implementing Sobol's quasirandom sequence generator.
| ACM Trans. Math. Softw. 14, 1 (March 1988), 88-100. 1988.
| [https://doi.org/10.1145/42288.214372](https://doi.org/10.1145/42288.2143720).
|
| Method resolution order:
| DigitalNetB2
| qmcpy.discrete_distribution.abstract_discrete_distribution.AbstractLDDiscreteDistribution
| qmcpy.discrete_distribution.abstract_discrete_distribution.AbstractDiscreteDistribution
| builtins.object
|
| Methods defined here:
|
| __init__(self, dimension=1, replications=None, seed=None, randomize='LMS DS', generating_matrices='joe_kuo.6.21201.txt', order='RADICAL INVERSE', t=63, alpha=1, msb=None, _verbose=False, graycode=None, t_max=None, t_lms=None)
| Args:
| dimension (Union[int,np.ndarray]): Dimension of the generator.
|
| - If an `int` is passed in, use generating vector components at indices 0,...,`dimension`-1.
| - If an `np.ndarray` is passed in, use generating vector components at these indices.
|
| replications (int): Number of independent randomizations of a pointset.
| seed (Union[None,int,np.random.SeedSeq): Seed the random number generator for reproducibility.
| randomize (str): Options are
|
| - `'LMS DS'`: Linear matrix scramble with digital shift.
| - `'LMS'`: Linear matrix scramble only.
| - `'DS'`: Digital shift only.
| - `'NUS'`: Nested uniform scrambling. Also known as Owen scrambling.
| - `'FALSE'`: No randomization. In this case the first point will be the origin.
|
| generating_matrices (Union[str,np.ndarray,int]: Specify the generating matrices.
|
| - A `str` should be the name (or path) of a file from the LDData repo at [https://github.com/QMCSoftware/LDData/tree/main/dnet](https://github.com/QMCSoftware/LDData/tree/main/dnet).
| - An `np.ndarray` of integers with shape $(d,m_\mathrm{max})$ or $(r,d,m_\mathrm{max})$ where $d$ is the number of dimensions, $r$ is the number of replications, and $2^{m_\mathrm{max}}$ is the maximum number of supported points. Setting `msb=False` will flip the bits of ints in the generating matrices.
|
| order (str): `'RADICAL INVERSE'`, or `'GRAY'` ordering. See the doctest example above.
| t (int): Number of bits in integer represetation of points *after* randomization. The number of bits in the generating matrices is inferred based on the largest value.
| alpha (int): Interlacing factor for higher order nets.
| When `alpha`>1, interlacing is performed regardless of the generating matrices,
| i.e., for `alpha`>1 do *not* pass in generating matrices which are already interlaced.
| The Note for this class contains more info.
| msb (bool): Flag for Most Significant Bit (MSB) vs Least Significant Bit (LSB) integer representations in generating matrices. If `msb=False` (LSB order), then integers in generating matrices will be bit-reversed.
| _verbose (bool): If `True`, print linear matrix scrambling matrices.
|
| ----------------------------------------------------------------------
| Methods inherited from qmcpy.discrete_distribution.abstract_discrete_distribution.AbstractLDDiscreteDistribution:
|
| __repr__(self)
| Return repr(self).
|
| ----------------------------------------------------------------------
| Methods inherited from qmcpy.discrete_distribution.abstract_discrete_distribution.AbstractDiscreteDistribution:
|
| __call__(self, n=None, n_min=None, n_max=None, return_binary=False, warn=True)
| - If just `n` is supplied, generate samples from the sequence at indices 0,...,`n`-1.
| - If `n_min` and `n_max` are supplied, generate samples from the sequence at indices `n_min`,...,`n_max`-1.
| - If `n` and `n_min` are supplied, then generate samples from the sequence at indices `n`,...,`n_min`-1.
|
| Args:
| n (Union[None,int]): Number of points to generate.
| n_min (Union[None,int]): Starting index of sequence.
| n_max (Union[None,int]): Final index of sequence.
| return_binary (bool): Only used for `DigitalNetB2`.
| If `True`, *only* return the integer representation `x_integer` of base 2 digital net.
| warn (bool): If `False`, disable warnings when generating samples.
|
| Returns:
| x (np.ndarray): Samples from the sequence.
|
| - If `replications` is `None` then this will be of size (`n_max`-`n_min`) $\times$ `dimension`
| - If `replications` is a positive int, then `x` will be of size `replications` $\times$ (`n_max`-`n_min`) $\times$ `dimension`
|
| Note that if `return_binary=True` then `x` is returned where `x` are integer representations of the digital net points.
|
| gen_samples(self, n=None, n_min=None, n_max=None, return_binary=False, warn=True)
|
| pdf(self, x)
|
| spawn(self, s=1, dimensions=None)
| Spawn new instances of the current discrete distribution but with new seeds and dimensions.
| Used by multi-level QMC algorithms which require different seeds and dimensions on each level.
|
| Note:
| Use `replications` instead of using `spawn` when possible, e.g., when spawning copies which all have the same dimension.
|
| Args:
| s (int): Number of copies to spawn
| dimensions (np.ndarray): Length `s` array of dimensions for each copy. Defaults to the current dimension.
|
| Returns:
| spawned_discrete_distribs (list): Discrete distributions with new seeds and dimensions.
|
| ----------------------------------------------------------------------
| Data descriptors inherited from qmcpy.discrete_distribution.abstract_discrete_distribution.AbstractDiscreteDistribution:
|
| __dict__
| dictionary for instance variables
|
| __weakref__
| list of weak references to the object
Out[33]:
['__call__', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_gen_samples', '_spawn', 'gen_samples', 'pdf', 'spawn']
Stopping criteria¶
The stopping criteria implemented come from several sources:
Changing the tolerance¶
If you want to change the error tolerance, without creating a new StoppingCriterion object, there is a set_tolerance method.
In [34]:
Copied!
sc = qmcpy.CubQMCLatticeG(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Lattice(5,seed=7),covariance=1/2)),abs_tol=.01)
print(sc.integrate()[1],'\n\n')
sc.set_tolerance(abs_tol=.001) #changing the tolerance
print(sc.integrate()[1],'\n\n')
sc = qmcpy.CubQMCLatticeG(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Lattice(5,seed=7),covariance=1/2)),abs_tol=.01)
print(sc.integrate()[1],'\n\n')
sc.set_tolerance(abs_tol=.001) #changing the tolerance
print(sc.integrate()[1],'\n\n')
Data (Data)
solution 1.135
comb_bound_low 1.129
comb_bound_high 1.142
comb_bound_diff 0.012
comb_flags 1
n_total 2^(13)
n 2^(13)
time_integrate 0.007
CubQMCLatticeG (AbstractStoppingCriterion)
abs_tol 0.010
rel_tol 0
n_init 2^(10)
n_limit 2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 5
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 7
Data (Data)
solution 1.136
comb_bound_low 1.135
comb_bound_high 1.136
comb_bound_diff 0.001
comb_flags 1
n_total 2^(17)
n 2^(17)
time_integrate 0.040
CubQMCLatticeG (AbstractStoppingCriterion)
abs_tol 0.001
rel_tol 0
n_init 2^(10)
n_limit 2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 5
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 7
Relative error tolerances¶
For some problems a relative error tolerance may make more sense. By default rel_tol is zero. The stopping criterion stops when either of the two tolerances is met.
In [35]:
Copied!
sc = qmcpy.CubQMCLatticeG(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Lattice(5,seed=7),covariance=1/2)),abs_tol=.01) #zero relative tolerance by default
print(sc.integrate()[1],'\n\n')
sc.set_tolerance(abs_tol=0,rel_tol=0.005) #nonzero relative tolerance and zero absolute tolerance
print(sc.integrate()[1])
sc = qmcpy.CubQMCLatticeG(qmcpy.Keister(qmcpy.Gaussian(qmcpy.Lattice(5,seed=7),covariance=1/2)),abs_tol=.01) #zero relative tolerance by default
print(sc.integrate()[1],'\n\n')
sc.set_tolerance(abs_tol=0,rel_tol=0.005) #nonzero relative tolerance and zero absolute tolerance
print(sc.integrate()[1])
Data (Data)
solution 1.135
comb_bound_low 1.129
comb_bound_high 1.142
comb_bound_diff 0.012
comb_flags 1
n_total 2^(13)
n 2^(13)
time_integrate 0.004
CubQMCLatticeG (AbstractStoppingCriterion)
abs_tol 0.010
rel_tol 0
n_init 2^(10)
n_limit 2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 5
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 7
Data (Data)
solution 1.135
comb_bound_low 1.131
comb_bound_high 1.139
comb_bound_diff 0.008
comb_flags 1
n_total 2^(14)
n 2^(14)
time_integrate 0.006
CubQMCLatticeG (AbstractStoppingCriterion)
abs_tol 0
rel_tol 0.005
n_init 2^(10)
n_limit 2^(30)
Keister (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
transform Gaussian (AbstractTrueMeasure)
mean 0
covariance 2^(-1)
decomp_type PCA
Lattice (AbstractLDDiscreteDistribution)
d 5
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 7
Use cases¶
Only a few use cases have been coded so far. These include
qmcpy.keisterKeiser's example- European and Asian option pricing (Single and Multi-level)
- Computing Q-Noisy Expected Improvement (qEI) for Bayesian Optimization, see the blog at https://qmcpy.wpcomstaging.com/2020/07/19/qei-with-qmcpy/ and the Colaboratory notebook at https://drive.google.com/drive/folders/1EOREUL7lx4hytuZRQXB50UV8aE9JYN_a
- Importance Sampling
- Acceptance Rejection Sampling
- Custom Sampling by Inverse CDF Transform
qmcpy would benefit from contributions.
Acknowledgements and Contributing¶
The qmcpy code has primarily been written by Aleksei Sorokin (BS/MS expected in 2021), with the generous financial support of SigOpt https://sigopt.com. However, much of the code is adapted from that of our friends (see above)
We hope that QMCPy will become supported by the community. Your contribution will add value to QMCPy, while allowing you to take advantage of the contribution of others.
We highlight the value of community owned software.
Benefitting from each other's work is easier when we share¶
Most of us are very good at just one or two things:
- LD sequence generators
- Increasing efficiency (e.g., MLMC)
- Stopping criteria
- Realistic use cases
Having a shared software library let's us take advantage of the best
Provides a consistent interface for different pieces from different places
Supports reproducible computational research
Tedious stuff only needs to be figured out once
A community helps find and correct code errors¶
By having more eyes on code than just the developer's we are more likely to spot errors or idiosyncracies. Here are two examples.
MATLAB's Sobol' generator¶
Several years ago Lluís Antoni Jiménez Rugama discovered that the scrambling of the Sobol' generators implemented in MATLAB's Statistics Toolbox was wrong. After reporting the problem to the developers, it was corrected in MATLAB 2017a
PyTorch's Sobol' generator¶
PyTorch is a popular Python library with its own Sobol generator. However, it should be used with care.
- As noted above, the first point is skipped.
- We found that unless you specify double precision, you get points that have 1 as a coordinate far too often. The developer seemed unaware of this.
These issues have been reported at https://github.com/pytorch/pytorch/issues/32047.
How you can contribute¶
After trying QMCPy out help us out by
Easy.
Submit your bugs and feature requests as issues to https://github.com/QMCSoftware/QMCSoftware/issues
Moderately Difficult.
Ask your students or collaborators to try QMCPy for their own work and submit their bugs and feature requests
Heroic.
Add a feature or use case and make a pull request at https://github.com/QMCSoftware/QMCSoftware/pulls so that we can included it in our next release
Questions? Email us at qmc-software@googlegroups.com
References¶
S.-C. T. Choi, Y. Ding, F. J. Hickernell, L. Jiang, Ll. A. Jimenez Rugama, D. Li, Jagadeeswaran R., X. Tong, K. Zhang, Y. Zhang, and X. Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3.1) [MATLAB Software], 2020. Available from
http://gailgithub.github.io/GAIL_Dev/H. Faure and C. Lemieux. “Implementation of Irreducible Sobol’ Sequences in Prime Power Bases,” Mathematics and Computers in Simulation 161 (2019): 13–22.
M. B. Giles. "Multi-level Monte Carlo path simulation," Operations Research, 56(3):607-617, 2008.
http://people.maths.ox.ac.uk/~gilesm/files/OPRE_2008.pdf.M. B. Giles. "Improved multilevel Monte Carlo convergence using the Milstein scheme," 343-358, in Monte Carlo and Quasi-Monte Carlo Methods 2006, Springer, 2008.
http://people.maths.ox.ac.uk/~gilesm/files/mcqmc06.pdf.M. B. Giles and B. J. Waterhouse. "Multilevel quasi-Monte Carlo path simulation," pp.165-181 in Advanced Financial Modelling, in Radon Series on Computational and Applied Mathematics, de Gruyter, 2009.
http://people.maths.ox.ac.uk/~gilesm/files/radon.pdfF. J. Hickernell, L. Jiang, Y. Liu, and A. B. Owen, "Guaranteed conservative fixed width confidence intervals via Monte Carlo sampling," Monte Carlo and Quasi-Monte Carlo Methods 2012 (J. Dick, F.Y. Kuo, G. W. Peters, and I. H. Sloan, eds.), pp. 105-128, Springer-Verlag, Berlin, 2014. DOI: 10.1007/978-3-642-41095-6_5
F. J. Hickernell and Lluis Antoni Jimenez Rugama, "Reliable adaptive cubature using digital sequences," Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (R. Cools and D. Nuyens, eds.), Springer Proceedings in Mathematics and Statistics, vol. 163, Springer-Verlag, Berlin, 2016, arXiv:1410.8615 [math.NA], pp. 367-383.
M. Hofert and C. Lemieux (2019). qrng: (Randomized) Quasi-Random Number Generators. R package version 0.0-7.
https://CRAN.R-project.org/package=qrng.Ll. A. Jimenez Rugama and F. J. Hickernell, "Adaptive multidimensional integration based on rank-1 lattices," Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (R. Cools and D. Nuyens, eds.), Springer Proceedings in Mathematics and Statistics, vol. 163, Springer-Verlag, Berlin, 2016, arXiv:1411.1966, pp. 407-422.
B. D. Keister, Multidimensional Quadrature Algorithms, 'Computers in Physics', 10, pp. 119-122, 1996.
F. Y. Kuo and D. Nuyens. "Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients - a survey of analysis and implementation," Foundations of Computational Mathematics, 16(6):1631-1696, 2016. (springer link, arxiv link)
P. L’Ecuyer and D. Munger, "LatticeBuilder: A General Software Tool for Constructing Rank-1 Lattice Rules," ACM Transactions on Mathematical Software. 42 (2015). 10.1145/2754929.
Y. Li, L. Kang, L., and F. J. Hickernell, Is a Transformed Low Discrepancy Design Also Low Discrepancy? in Contemporary Experimental Design, Multivariate Analysis and Data Mining, Festschrift in Honour of Professor Kai-Tai Fang (J. Fan and J. Pan, eds.), p. 69–92, 2020, https://arxiv.org/abs/2004.09887.
A. B. Owen, "A randomized Halton algorithm in R," 2017. arXiv:1706.02808 [stat.CO]