Volumen zwischen zwei Normalverteilungen (Gausse) berechnen

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
m4tt3s
User
Beiträge: 6
Registriert: Donnerstag 2. August 2018, 10:55

Hallo,
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()
Solltet Ihr mehr Informationen brauchen, schreibt einfach.
Vielen Dank :)
Sirius3
User
Beiträge: 17711
Registriert: Sonntag 21. Oktober 2012, 17:20

Die Formeln sind sehr schwer lesbar, wegen der vielen überflüssigen Klammern. Ein paar Leerzeichen würden auch helfen. Woher kommen die ganzen Zahlenwerte?

Das Volumen eines Gauss ist nicht so schwer zu berechnen. Sollte in entsprechenden Mathebüchern auch erklärt werden.
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

Moin,

die Fläche unter einer normierten Gaußglocke (im Folgenden ϕ(x, μ, σ)) ist immer 1:

Code: Alles auswählen

∫_{-∞}^{∞} ϕ(x, μ, σ) dx = 1.
Wenn du jetzt zwei verschiedene Funktionen f(x) = α ϕ(x, μ_1, σ_1) und g(x) = β ϕ(x, μ_2, σ_2) hast und die Fläche A zwischen beiden berechnen möchtest, dann musst du das Integral

Code: Alles auswählen

A = ∫_{-∞}^{∞} (f(x) - g(x)) dx 
= ∫_{-∞}^{∞} (α ϕ(x, μ_1, σ_1) - β ϕ(x, μ_2, σ_2)) dx
= α - β
lösen, was einfach A = α - β ergibt. Wenn deine PDF zweidimensional ist, dürfte das die Integrale ein bisschen hässlicher machen, aber auch eine zweidimensionale PDF ist normiert, dementsprechend ist die Differenz der Volumina unter den Gaußglocken auch wieder die Differenz der Skalierungsfaktoren.

Um das nachzuprüfen, könnte man die (zweidimensionale) PDF aufstellen und dann über alle p ∈ ℝ² in Polarkoordinaten integrieren.
m4tt3s
User
Beiträge: 6
Registriert: Donnerstag 2. August 2018, 10:55

Okay mathematisch klar und sinnvoll, danke dir narpfel! :)
Kann mir jemand bei der Umsetzung für den Code helfen?
narpfel
User
Beiträge: 643
Registriert: Freitag 20. Oktober 2017, 16:10

Wo ist denn das Problem? Du musst doch nur die beiden Skalierungsfaktoren aus dem Code ablesen und subtrahieren.

Falls die Offsets beachtet werden sollen, hast du das Problem, dass deine Funktionen nicht mehr über ganz ℝ² integrierbar sind und du die Integrale auf ein passendes Gebiet einschränken musst (wodurch du die bivariate Fehlerfunktion benutzen musst).
m4tt3s
User
Beiträge: 6
Registriert: Donnerstag 2. August 2018, 10:55

Achso, sorry, habe es erst jetzt geschnallt. D.h. in meinem Code bzw. meiner Definition wäre α=1 und β=(w0/u)**(2), korrekt?
Da meine erste Funktion (red, fitted) mit Offset und Amplitude versehen ist und ich meine zweite Funktion (blau, Simulation) auf diese Amplitude und den Offset normiert habe.
Antworten