ich würde gerne das Volumen zwischen zwei Normalverteilungen (Gauss Funktionen) berechnen.
Die eine Verteilung entsteht durch einen fit an Daten (rot) und die andere durch eine Simulation (Berechnung einer Funktion, in blau).
Ich kann beide Funktionen plotten (siehe Bild), jedoch habe ich keine Ahnung wie ich den Volumenunterschied berechnen kann.
Damit meine ich genau den Teil, zwischen rot und blau. Das ganze hat einen physikalischen Hintergrund: Das Volumen bzw. die
Fläche unter dem Gauss entspricht der Leistung. Ich möchte eben diesen Unterschied wissen.
Habe auch an Integration gedacht, bin dort aber nicht so bewandert.
Ich würde mich über eine ausführliche Antwort freuen.
Die Grafik dazu findet Ihr hier:
https://stackoverflow.com/questions/518 ... 9_51837042
Anbei mein Code:
Code: Alles auswählen
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()
Vielen Dank