Fit 2D Gaus Funktion

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
trublu
User
Beiträge: 18
Registriert: Montag 20. Juni 2016, 20:05

Hallo nochmal,

ich würde gerne einen 2D Gauß an meine Daten fitten. Mein Problem ist dabei allerdings, dass ich nur Messdaten auf der x und y Achse habe und die auch nicht äquidistant sind, allerdings auf den Achsen gleich. Der Gauß dient hier lediglich als hoffentlich einfacheres Beispiel, die echte Funktion ist z.B. nicht radial symmetrisch. Auf Stackoverflow findet man z.B. das http://stackoverflow.com/questions/2156 ... rror-and-m die Lösung setzt allerdings glaub ich vorraus, dass man überall Daten hat. Ich hab es jedenfalls nicht geschafft dass auf mein Problem anzuwenden. Lösung hierfür könnte scipy.interpolate.griddata sein, allerdings finde ich das ziemlich umständlich (und hab es auch wieder nicht hinbekommen).

So sieht das ganze bei mir aus. xdata und ydata sind unterschiedliche Gaußkurven an die ich jetzt gerne einen 2D Gauß fitten würde. Ich hatte schon versucht aus xdata und ydata eine große Matrix zu machen wo alle "nicht Achsen Elemente" auf np.nan gesetzt werden, aber auch das mag scipy nicht,

Code: Alles auswählen

import scipy.optimize as opt
from scipy.interpolate import griddata
import numpy as np
import pylab as plt


def Gaus2D((x, y), amplitude, x0, y0, sigma_x, sigma_y):
    gaus = amplitude * np.exp(-((x-x0)**2/(2*sigma_x**2) + (y-y0)**2/(2*sigma_y**2)))
    return gaus.ravel()

def Gaus1D(x,x0,sigma):
    return np.exp(-np.power((x - x0)/(2.*sigma), 2.))

if __name__ == '__main__':
    x = np.array([-70,-40,-20,-10,-5,0,5,10,20,40,70])
    y = np.array([-70,-40,-20,-10,-5,0,5,10,20,40,70])
    xdata = Gaus1D(x, 0, 50)
    ydata = Gaus1D(y, 0, 80)

    # wenn data jetzt ein "volles" 2D data array wäre
    initial_guess = (1,0,0,20,40)
    popt, pcov = opt.curve_fit(Gaus2D, (x,y), data, p0 = initial_guess)
Sirius3
User
Beiträge: 17711
Registriert: Sonntag 21. Oktober 2012, 17:20

@trublu: es ist völlig egal, ob Du jetzt Daten auf einem regelmäßigen Raster hast oder nur einzelne Punkte. Die Optimierung läuft über alle bekannten Meßwerte, daher sollte Dein curve_fit genau so, wie Du geschrieben hast, funktionieren. Ich weiß aber nicht, wie gut so eine Optimierung funktionieren kann, wenn Du nur Werte auf der Diagonalen hast. In diesem Fall solltest Du einen 1D-Gauss-Fit machen und annehmen, der Gauß wäre symmetrisch in x- und y-Richtung.
trublu
User
Beiträge: 18
Registriert: Montag 20. Juni 2016, 20:05

Sirius3 hat geschrieben:@trublu: es ist völlig egal, ob Du jetzt Daten auf einem regelmäßigen Raster hast oder nur einzelne Punkte. Die Optimierung läuft über alle bekannten Meßwerte, daher sollte Dein curve_fit genau so, wie Du geschrieben hast, funktionieren.
Das tut er natürlich nicht, das es data in dem Code nicht gibt und ich nicht weiß wie ich xdata und ydata richtig übergebe.

Ich hab es z.B. so probiert

Code: Alles auswählen

    popt, pcov = opt.curve_fit(Gaus2D, (x,y), (xdata,ydata), p0 = initial_guess)
  File "/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py", line 651, in curve_fit
    res = leastsq(func, p0, args=args, full_output=1, **kwargs)
  File "/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py", line 380, in leastsq
    raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
TypeError: Improper input: N=5 must not exceed M=2
Sirius3 hat geschrieben:In diesem Fall solltest Du einen 1D-Gauss-Fit machen und annehmen, der Gauß wäre symmetrisch in x- und y-Richtung.
Das ist wie gesagt in Wirklichkeit ja nicht der Fall. Wie gut das funktioniert ist unter anderem die Frage, die ich mit diesem Beispiel zu beantworten hoffe :)
Sirius3
User
Beiträge: 17711
Registriert: Sonntag 21. Oktober 2012, 17:20

Hier der natürlich einfache Fall, dass die Daten exakt sind, aber es soll ja nur zeigen, dass curve_fit mit dem Input zurecht kommt:

Code: Alles auswählen

import scipy.optimize as opt
import numpy as np
def Gaus2D((x, y), amplitude, x0, y0, sigma_x, sigma_y):
    gaus = amplitude * np.exp(-((x-x0)**2/(2*sigma_x**2) + (y-y0)**2/(2*sigma_y**2)))
    return gaus.ravel()

x = (np.random.random(100)-0.5)*10
y = (np.random.random(100)-0.5)*10
data = Gaus2D((x,y), 2.0, -0.1, 0.2, 3, 2)
initial_guess = (1,0,0,20,40)
popt, pcov = opt.curve_fit(Gaus2D, (x,y), data, p0 = initial_guess)
print popt
# -> array([ 2. , -0.1,  0.2, -3. ,  2. ])
trublu
User
Beiträge: 18
Registriert: Montag 20. Juni 2016, 20:05

okay jetzt weiß ich immerhin das mein eigentlicher Fehler ist wie ich aus xdata ydata zu dem "echten" data array komme. Meine Idee war jetzt daraus doch erst eine komplette Matrix zu machen und dann wie in Gaus2D ein ravel() darauf anzuwenden. Das scheint aber auch nicht richtig zu sein.

Code: Alles auswählen

    data = np.ones([Nx, Ny])*np.nan
    data[:,Nx/2] = xdata
    data[Ny/2,:] = ydata
    data = data.ravel()
    data = data[~np.isnan(data)]
    
        popt, pcov = opt.curve_fit(Gaus2D, (x,y), data, p0 = initial_guess)
  File "/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py", line 651, in curve_fit
    res = leastsq(func, p0, args=args, full_output=1, **kwargs)
  File "/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py", line 377, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py", line 453, in _general_function
    return function(xdata, *params) - ydata
ValueError: operands could not be broadcast together with shapes (11,) (21,)
   
Sirius3
User
Beiträge: 17711
Registriert: Sonntag 21. Oktober 2012, 17:20

@trublu: wie sehen denn nun Deine echten Daten aus? Wenn Du x und xdata und y und ydata hast mußt Du daraus zuerst Vektoren x, 0, xdata und 0, y, ydata machen, diese zu 3 kombinierten Vektoren zusammensetzen und dann den Fit machen.
trublu
User
Beiträge: 18
Registriert: Montag 20. Juni 2016, 20:05

Wie meine Daten aussehen ist ja im ersten Beitrag beschrieben. Ich verstehe aber nicht wie du das jetzt mit den Vektoren kombinieren meinst. So in etwa? Wie mache ich damit den Fit?

Code: Alles auswählen

import numpy as np

def Gaus1D(x,x0,sigma):
    return np.exp(-np.power((x - x0)/(2.*sigma), 2.))
 
 if __name__ == '__main__':   
    x = np.array([-70,-40,-20,-10,-5,0,5,10,20,40,70])
    y = np.array([-70,-40,-20,-10,-5,0,5,10,20,40,70])
    xdata = Gaus1D(x, 0, 50)
    ydata = Gaus1D(y, 0, 80)
    
    xpoints = np.coloum_stack((x,[0]*len(x),xdata)
    ypoints = np.coloum_stack(([0]*len(x),y,ydata)
trublu
User
Beiträge: 18
Registriert: Montag 20. Juni 2016, 20:05

Ich glaub ich hab es doch hinbekommen, dass Ergebnis wäre dann allerdings schrecklich, insbesondere die sigmas werden ja mal gar nicht getroffen dabei sind die Daten noch nichtmal verrauscht. Vom Prinzip her wäre das so aber richtig?

Code: Alles auswählen

import scipy.optimize as opt
import numpy as np

def Gaus2D((x, y), amplitude, x0, y0, sigma_x, sigma_y):
    gaus = amplitude * np.exp(-((x-x0)**2/(2*sigma_x**2) + (y-y0)**2/(2*sigma_y**2)))
    return gaus.ravel()

def Gaus1D(x,x0,sigma):
    return np.exp(-np.power((x - x0)/(2.*sigma), 2.))

if __name__ == '__main__':
    x = np.array([-70,-40,-20,-10,-5,0,5,10,20,40,70])
    y = np.array([-70,-40,-20,-10,-5,0,5,10,20,40,70])
    xdata = Gaus1D(x, 0, 50)
    ydata = Gaus1D(y, 0, 80)
    
    data = np.column_stack((xdata,ydata))
    xx = np.column_stack((x,[0]*len(x)))
    yy = np.column_stack(([0]*len(y),y))

    initial_guess = (1,0,0,20,40)
    popt, pcov = opt.curve_fit(Gaus2D, (xx,yy), data.ravel(), p0 = initial_guess)
    print popt # [  1.00000000e+00   3.85594393e-08  -2.76946895e-07   7.07106781e+01   1.13137085e+02]
Sirius3
User
Beiträge: 17711
Registriert: Sonntag 21. Oktober 2012, 17:20

@trublu: es funktioniert nicht, weil Du x- und y-Koordinaten wild miteinander mischst.

Als erstes schreibst Du den armen Gauss mal richtig:

Code: Alles auswählen

def gauss_2d((x, y), amplitude, x0, y0, sigma_x, sigma_y):
    return amplitude * np.exp(-((x-x0)**2/(2.*sigma_x**2) + (y-y0)**2/(2.*sigma_y**2)))
Dann brauchst Du einen Satz (x,y)-Koordinaten für alle Punkte, die auf der x-Achse liegen:

Code: Alles auswählen

x = np.array([-70,-40,-20,-10,-5,0,5,10,20,40,70])
xx = np.row_stack((x, [0]*len(x)))
xdata = gauss_2d(xx, 1.0, 0.0, 0.0, 50.0, 50.0)
und einen für die y-Achse:

Code: Alles auswählen

y = np.array([-70,-40,-20,-10,-5,0,5,10,20,40,70])
yy = np.row_stack(([0]*len(y), y))
ydata = gauss_2d(yy, 1.0, 0.0, 0.0, 80.0, 80.0)
Als nächstes werden die beiden Messdatensätze zu einem vereinigt:

Code: Alles auswählen

points = np.hstack([xx, yy])
data = np.hstack([xdata, ydata])
und als letztes kannst Du den Fit machen:

Code: Alles auswählen

initial_guess = (1,0,0,20,40)
popt, pcov = opt.curve_fit(gauss_2d, points, data, p0=initial_guess)
print popt
# [  1.00000000e+00   3.86989483e-08   1.70739685e-08   5.00000000e+01  8.00000000e+01]
trublu
User
Beiträge: 18
Registriert: Montag 20. Juni 2016, 20:05

okay, dann hab ich das jetzt verstanden. Vielen Dank.
Antworten