CODEWITHSUNDEEP
pythonscipygaussiannumerical-integration
user23870
user23870

Reputation: 1

Using scipy integrate tplquad to evaluate triple integral of multivariate gaussian

I'm trying to perform the following triple integral:

3d gaussian equation

f(v) is a 3 variable Gaussian probability density function.

The magnitude of the resultant velocity, V, must be less than some Vmax (so Vx, Vy, and Vz can each range from 0 or -Vmax to +Vmax, but if Vx = Vmax then Vy and Vz must be zero, for example).

For now we can take sigma = 1. I will ignore the division by v inside the integral for now as well, so I'm just integrating f(v).

I've been trying to use the scipy.integrate tplquad function but keep getting an answer that's of the order of magnitude ~1e-7 and provided Vmax is big enough (I'm using Vmax = 500) the integral should approach 1 since the probability over all (velocity) space is 1.

Here is my code:

from scipy.integrate import tplquad
from numpy import pi, exp, sqrt

def func(Vz,Vy,Vx):
    return sqrt( (1/(sigma*2*pi)**(3/2)) ) * exp( - ( ((Vx**2)/2) + ((Vy**2)/2) + ((Vz**2)/2) ) )

def Vymax(Vx):
    return sqrt(Vmax**2 - Vx**2)

def Vzmax(Vx,Vy):
    return sqrt(Vmax**2 - Vx**2 - Vy**2)

Vmax = 500
sigma = 1

integral = tplquad(func,0,Vmax,lambda a: 0, Vymax,lambda a, y: 0, Vzmax)
print(8*integral[0])

Output:

>>> (executing cell "Triple integral a..." (line 15 of "Integral4.py"))
1.4644282619462532e-07
>>>

The Vymax and Vzmax functions are to keep Vy and Vz values within ranges so that the resultant velocity doesn't exceed Vmax. I've used lambda functions in tplquad for 0 because tplquad requires a function as input for the 4th parameters onwards.

I can't see where I've gone wrong but I'm sure I must be missing something simple or being completely stupid.

If tplquad isn't the right function for this problem, is there an alternative? I'm even happy to write my own function but am unsure what the best method would be - I tried a Monte Carlo method but couldn't quite figure it out. I can't just separate out the components because eventually I need to have an offset Vx + Voffset so it'll no longer be centred on the origin. I've searched on here as much as I could and came across this but it doesn't describe properly how to do a multivariate Gaussian, since the question was (meant to be) about a single variable Gaussian.

Any help much appreciated.

Upvotes: 0

Views: 2944

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114811

Your formula for the multivariate normal distribution is not correct. You have both the sqrt call, and the factor of 1/2 in the power in the coefficient of the exponential. You are also missing the factor of 1/sigma**2 in the exponent.

Here's a modified version, using sigma in a way that matches the formula shown in the question:

def func(Vz,Vy,Vx):
    return (1.0/(sigma**2*2*pi)**(3./2)) * exp( - ( ((Vx**2)/2) + ((Vy**2)/2) + ((Vz**2)/2) )/sigma**2 )

After adding a print statement to print the value of sigma, I get:

In [28]: run tplquad_question.py
sigma = 1
8*integral[0] = 1.00000000196

Edit the script and run again:

In [29]: run tplquad_question.py
sigma = 2.5
8*integral[0] = 1.00000000927

Here's the complete tplquad_question.py:

from __future__ import print_function, division

from scipy.integrate import tplquad
from numpy import pi, exp, sqrt

def func(Vz,Vy,Vx):
    return (1.0/(sigma**2*2*pi)**(3./2)) * exp( - ( ((Vx**2)/2) + ((Vy**2)/2) + ((Vz**2)/2) )/sigma**2 )

def Vymax(Vx):
    return sqrt(Vmax**2 - Vx**2)

def Vzmax(Vx,Vy):
    return sqrt(Vmax**2 - Vx**2 - Vy**2)

Vmax = 500
sigma = 1

integral = tplquad(func,0,Vmax,lambda a: 0, Vymax,lambda a, y: 0, Vzmax)

print("sigma =", sigma)
print("8*integral[0] =", 8*integral[0])

Upvotes: 2

Related Questions