gplepage / vegas

Adaptive multidimensional Monte Carlo integration.
Other
106 stars 25 forks source link

Gaussian distributed and usual (limited) variables in one integral #24

Closed dsmic closed 2 years ago

dsmic commented 2 years ago

As the title says, I have to integrate a multidimensional array, where a big part of the variables have a gaussian probability distribution, but there are also a number of variable with an integration range.

May I just use a pdf function with hard cut offs for the integration ranges?

Any help appreciated🙂

dsmic commented 2 years ago

I will try adaptive map first, sorry for the noise

gplepage commented 2 years ago

There is a way to mix ordinary variables, restricted to a finite range, with Gaussian variables when using vegas.PDFIntegrator. Here is a short example that works out the expectation value of cos(y*x) in a joint distribution where x has a Gaussian distribution (mean=3, std dev=1) and y has a uniform distribution restricted to the interval [-1/2,1/2]:

import vegas
import gvar as gv 

# 1) define joint distribution: 
#   x = Gaussian 3(1) 
#   y = uniform on interval [-1/2, 1/2]
dist = gv.BufferDict()
dist['x'] = gv.gvar('3(1)')
dist['fy(y)'] = gv.BufferDict.uniform('fy', -0.5, 0.5)

# 2) create integrator to calculate expectation values
expval = vegas.PDFIntegrator(dist)

# 3) pre-adapt the vegas map
nitn = 10
neval = 1000
expval(nitn=nitn, neval=neval)

# 4) calculate <cos(y*x)>
@vegas.rbatchintegrand
def f(p):
    x = p['x']
    y = p['y']
    return gv.cos(y * x)
r = expval(f, nitn=nitn, neval=neval)
print('<cos(y*x)> =', r, '    divided by 2/pi =', r / (2/3.14159654))

This gives the following output:

<cos(y*x)> = 0.6467(13)     divided by 2/pi = 1.0158(20)

The answer should be close to $2/\pi$ since x is close to $\pi$. Variable y is treated here as a random variable drawn from a uniform probability distribution on an interval. This distribution is normalized, so the final result is divided by width of the interval. (Multiply the result by the width if you just want to integrate y over the interval, without the normalization.)

vegas.PDFIntegrator relies upon the gvar module to define probability distributions. These can be multi-dimensional and correlated. Also it supports a small number of non-Gaussian distributions (eg, uniform, log-Gaussian). gvar.BufferDict is a dictionary used here to represent the distribution. For more information see the gvar.BufferDict documentation: https://gvar.readthedocs.io/en/latest/gvar.html?highlight=BufferDict#gvar-bufferdict-objects