StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

POUsing scipy integrate tplquad to evaluate triple integral of multivariate gaussian
text
Body
copied!<p>I'm trying to perform the following triple integral:</p> <p><img src="https://i.stack.imgur.com/Zv4Q2.png" alt="3d gaussian equation"></p> <p>f(v) is a 3 variable Gaussian probability density function.</p> <p>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).</p> <p>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).</p> <p>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.</p> <p>Here is my code:</p> <pre><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]) </code></pre> <p>Output:</p> <pre><code>>>> (executing cell "Triple integral a..." (line 15 of "Integral4.py")) 1.4644282619462532e-07 >>> </code></pre> <p>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.</p> <p>I can't see where I've gone wrong but I'm sure I must be missing something simple or being completely stupid.</p> <p>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 <a href="https://stackoverflow.com/questions/509994/best-way-to-write-a-python-function-that-integrates-a-gaussian">this</a> but it doesn't describe properly how to do a multivariate Gaussian, since the question was (meant to be) about a single variable Gaussian.</p> <p>Any help much appreciated.</p> <p>Thanks.</p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload