Python: Volume between 2 normal distributions (3d)

Question

I want to calculate the volume between two normal distributions (Gaussian) in 3d.

I got 2 different normal distributions: 1) A Gaussian that is fitted to x,y,z data (red) and 2) a simulated Gaussian (blue). I can plot both functions but I have no idea how to calculate the difference between them. So what I want is the volume-difference between red and blue. This has a physical background: The volume corresponds to the power and I want to know this difference in power. I thought about integration, but I am not familiar with this. Thank you.

from matplotlib import pyplot;
from pylab import genfromtxt;  
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
########### FITTED GAUSSIAN (red) ##############
# Load file into mat0
mat0 = genfromtxt("0005.map");

#PLOT FIGURE 
fig = plt.figure(figsize=(20,10))
ax = plt.axes(projection='3d')

#Define Gaussian function
def twoD_Gauss((x,y),amplitude,x0,y0,sigma_x,sigma_y,offset):
    x0=float(x0)
    y0=float(y0)
    return offset + amplitude*np.exp(-(((x-x0)**(2)/(2*sigma_x**(2))) + ((y-y0)**(2)/(2*sigma_y**(2)))))

# Create x and y indices
x = mat0[:,0]
y = mat0[:,1]

#create data
data = mat0[:,2]

#plt.imshow(data)
import scipy.optimize as opt
initial_guess = (24000,150,143,25,25,6000)

#Fit Gaussian function
params, pcov = opt.curve_fit(twoD_Gauss, (x,y), data,initial_guess)

#Print fitted parameters
print(params)

#Plot fitted Gaussian
ax.plot_trisurf(x-150, y-143, twoD_Gauss((x,y),*params), cmap="Reds", linewidth=0,alpha=0.5)
ax.set_xlabel('x / mm')
ax.set_ylabel('y / mm')

#Plot settings
ax.view_init(0, 270)

########### SIMULATED GAUSSIAN (blue) ##############
#functions  
w0=1.701
lamb=0.90846
d_in1=45.0
foc1=38.35
zR=np.pi*w0**(2)/(lamb)
w1=w0*np.sqrt(1/(((d_in1)/foc1-1)**(2)+(zR/foc1)**(2)))    
zR1=np.pi*w1**(2)/(lamb)
foc2=420
d_in2=499.8971
d_2=606
d_out2=foc2+(d_in2-foc2)/(((d_in2)/foc2-1)**(2)+(zR1/foc2)**(2))
w2=w1*np.sqrt(1/(((d_in2)/foc2-1)**(2)+(zR1/foc2)**(2)))
zR2=np.pi*w2**(2)/(lamb)

u=w2*np.sqrt(1+((3001-d_in1-d_2-d_out2)/zR2)**(2))

def i_3(x,y):
    return 3818017.483*(w0/u)**(2)*np.exp(-(2*(x**(2)+y**(2)))/(u**(2)))+7115.230

#define x and y
x = np.linspace(-50, 50, 100)
y = np.linspace(-50, 50, 100)
X, Y = np.meshgrid(x,y)
Z = i_3(X,Y)

#Plot settings
ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.xaxis._axinfo["grid"]['color'] =  (1,1,1,0)
ax.yaxis._axinfo["grid"]['color'] =  (1,1,1,0)
ax.zaxis._axinfo["grid"]['color'] =  (1,1,1,0)
ax.w_zaxis.line.set_lw(0.)
ax.set_zticks([])
ax.view_init(-1, 215)

#Plot
surf=ax.plot_surface(X,Y,Z,rstride=1,cstride=1,cmap='Blues', edgecolor='none', alpha=1.0)

#print(I)
ax.text(10,100,12000,'Fitted Gaussian',color='red',fontsize=18)
ax.text(10,100,14000,'Simulated Gaussian',color='blue',fontsize=18)
plt.show()

Answer 1

Since your model is Gaussian, there is an analytical solution for the volume: volume = 2*pi*sigma_x*sigma_y*amplitude (see integral over x and over Y).

The volume in-between the two Gaussian is delta_V = volumeA - volumeB.

Here is a code to compute the integral anyway, using dblquad:

import numpy as np
from scipy.integrate import dblquad


def volume_Gaussian(amplitude, xy0, sigma_xy):
    return 2*np.pi*amplitude*sigma_xy[0]*sigma_xy[1]

def gaussian2D(x, y, amplitude, xy0, sigma_xy):
    x = x - xy0[0]
    y = y - xy0[1]
    return amplitude*np.exp( -x**2/2/sigma_xy[0]**2 - y**2/2/sigma_xy[1]**2 )                                                          

# test
args = (2, (0, 0), (1, 2))
volume = dblquad(gaussian2D, -np.Inf, +np.Inf,
                 lambda u:-np.Inf, lambda v:+np.Inf,
                 args=args )
print(volume_Gaussian(*args), volume) 
# 25.132741228718345 (25.132741228718398, 5.102585580809855e-08)


# Two Gaussian
def A_minus_B(x, y, argsA, argsB):
    return gaussian2D(x, y, *argsA) - gaussian2D(x, y, *argsB)

argsA = (2, (0, 0), (1, 2))
argsB = (1, (0, 1), (1, 1))

volume_A_minus_B = dblquad(A_minus_B, -np.Inf, +np.Inf,
                           lambda u:-np.Inf, lambda u:np.Inf,
                           args=(argsA, argsB))

print(volume_Gaussian(*argsA) - volume_Gaussian(*argsB), volume_A_minus_B) 
# 18.84955592153876 (18.849555921538805, 1.4535527689371197e-07)