2D-GAUSS-Fit für die Intensitäten an bestimmten Koordinaten in Python

Einführung

Gibt es mehrere Möglichkeiten, dies zu nähern. Sie können nicht-lineare Methoden (z.B. scipy.optimize.curve_fit), aber Sie werden langsamer und sind nicht garantiert, um zu konvergieren. Sie können linearisieren das problem (schnelle, eindeutige Lösung), aber jedes Geräusch in der "Schwänze" der Verteilung Probleme bereiten. Es gibt tatsächlich ein paar tricks, die Sie anwenden können, um diesem besonderen Fall zu vermeiden das zweite Problem. Ich werde ein paar Beispiele zeigen, aber ich habe nicht Zeit zu zeigen, alle "tricks" jetzt.

Nur als Randnotiz, eine Allgemeine 2D-guassian hat 6 Parameter, so werden Sie nicht in der Lage, vollständig passen die Dinge mit 4 Punkten. Aber es klingt wie Sie sein könnte, unter der Annahme, dass es keine Kovarianz zwischen x und y und der Varianzen in jeder Richtung (d.h. eine vollkommen "Runde" bell-Kurve). Wenn das der Fall ist, dann brauchen Sie nur vier Parameter. Wenn Sie wissen, die amplitude der guassian, müssen Sie nur drei. Aber ich werde beginnen mit der Allgemeinen Lösung und vereinfachen Sie es später, wenn Sie wollen.

Für den moment konzentrieren wir uns auf die Lösung dieses problem mit nicht-linearen Methoden (z.B. scipy.optimize.curve_fit).

---

Die Allgemeine Gleichung für eine 2D-guassian ist (direkt aus der wikipedia):

wo:

ist im wesentlichen 0.5 über die Kovarianz-matrix A ist die amplitude,
und (X₀, Y₀) ist das Zentrum

---

Generieren vereinfachte Beispieldaten

Schreiben wir die obige Gleichung aus:

import numpy as np
import matplotlib.pyplot as plt

def gauss2d(x, y, amp, x0, y0, a, b, c):
    inner = a * (x - x0)**2 
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

Und dann lassen Sie uns generieren Sie einige Beispiel-Daten. So starten Sie mit, wir erstellen einige Daten, die leicht zu passen:

np.random.seed(1977) # For consistency
x, y = np.random.random((2, 10))
x0, y0 = 0.3, 0.7
amp, a, b, c = 1, 2, 3, 4

zobs = gauss2d(x, y, amp, x0, y0, a, b, c)

fig, ax = plt.subplots()
scat = ax.scatter(x, y, c=zobs, s=200)
fig.colorbar(scat)
plt.show()

Beachten Sie, dass wir noch nicht Hinzugefügt, jedes Geräusch, und das Zentrum der Verteilung liegt innerhalb der Bandbreite, die wir haben die Daten (also center bei 0,3, 0,7 und eine Streuung in den x -, y-Beobachtungen zwischen 0 und 1). Für den moment bleiben wir bei dieser, und dann werden wir sehen, was passiert, wenn wir den Rauschen-hinzufügen-und-shift-Zentrum.

---

Nicht-linearen Montage

So starten Sie mit, lassen Sie uns scpy.optimize.curve_fit Vorform einer nicht-linearen least-squares-fit der Gauß-Funktion. (On a side note, Sie können spielen, um mit der exakten Minimierung Algorithmus durch einige der anderen Funktionen in scipy.optimize.)

Den scipy.optimize Funktionen erwarten eine etwas andere Funktion Unterschrift als die, die wir ursprünglich oben schrieb. Wir könnten schreiben eine wrapper-zu "übersetzen", aber lasst uns einfach neu schreiben Sie die gauss2d - Funktion statt:

def gauss2d(xy, amp, x0, y0, a, b, c):
    x, y = xy
    inner = a * (x - x0)**2
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

Alles, was wir Taten, war die Funktion haben, erwarten, dass die unabhängigen Variablen (x & y) als Einzel-2xN-array.

Nun müssen wir eine erste Vermutung an, was die guassian Kurve die Parameter tatsächlich sind. Dies ist optional (der Standard ist alle lieben, wenn ich mich Recht erinnere), aber du bist wahrscheinlich Probleme konvergiert, wenn 1 1 ist nicht besonders nah an die "wahre" Zentrum der Gauß-Kurve. Aus diesem Grund verwenden wir die x-und y-Werte, die unsere größte beobachtete z-Wert, die als Ausgangspunkt für das Zentrum. Ich lasse den rest der Parameter als 1, aber wenn Sie wissen, dass Sie wahrscheinlich durchweg signifikant Verschieden sein, ändern Sie Sie, um etwas mehr zumutbar ist.

Hier die komplette stand-alone-Beispiel:

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

def main():
    x0, y0 = 0.3, 0.7
    amp, a, b, c = 1, 2, 3, 4
    true_params = [amp, x0, y0, a, b, c]
    xy, zobs = generate_example_data(10, true_params)
    x, y = xy

    i = zobs.argmax()
    guess = [1, x[i], y[i], 1, 1, 1]
    pred_params, uncert_cov = opt.curve_fit(gauss2d, xy, zobs, p0=guess)

    zpred = gauss2d(xy, *pred_params)
    print 'True parameters: ', true_params
    print 'Predicted params:', pred_params
    print 'Residual, RMS(obs - pred):', np.sqrt(np.mean((zobs - zpred)**2))

    plot(xy, zobs, pred_params)
    plt.show()

def gauss2d(xy, amp, x0, y0, a, b, c):
    x, y = xy
    inner = a * (x - x0)**2
    inner += 2 * b * (x - x0)**2 * (y - y0)**2
    inner += c * (y - y0)**2
    return amp * np.exp(-inner)

def generate_example_data(num, params):
    np.random.seed(1977) # For consistency
    xy = np.random.random((2, num))

    zobs = gauss2d(xy, *params)
    return xy, zobs

def plot(xy, zobs, pred_params):
    x, y = xy
    yi, xi = np.mgrid[:1:30j, -.2:1.2:30j]
    xyi = np.vstack([xi.ravel(), yi.ravel()])

    zpred = gauss2d(xyi, *pred_params)
    zpred.shape = xi.shape

    fig, ax = plt.subplots()
    ax.scatter(x, y, c=zobs, s=200, vmin=zpred.min(), vmax=zpred.max())
    im = ax.imshow(zpred, extent=[xi.min(), xi.max(), yi.max(), yi.min()],
                   aspect='auto')
    fig.colorbar(im)
    ax.invert_yaxis()
    return fig

main()

In diesem Fall genau(ish) die Wiedergewinnung der ursprünglichen, "wahren" Parameter.

True parameters:  [1, 0.3, 0.7, 2, 3, 4]
Predicted params: [ 1.   0.3  0.7  2.   3.   4. ]
Residual, RMS(obs - pred): 1.01560615193e-16

Wie wir gleich sehen werden, ist dieses nicht immer der Fall sein wird...

---

Hinzufügen Von Rauschen

Fügen wir einige Geräusche unserer Beobachtungen. Alles was ich getan habe, ist hier ändern Sie die generate_example_data Funktion:

def generate_example_data(num, params):
    np.random.seed(1977) # For consistency
    xy = np.random.random((2, num))

    noise = np.random.normal(0, 0.3, num)
    zobs = gauss2d(xy, *params) + noise
    return xy, zobs

Aber das Ergebnis sieht ganz anders aus:

Soweit es die Parameter:

True parameters:  [1, 0.3, 0.7, 2, 3, 4]
Predicted params: [  1.129    0.263   0.750   1.280   32.333   10.103  ]
Residual, RMS(obs - pred): 0.152444640098

Den vorhergesagten center hat sich nicht viel verändert, aber die b und c Parameter haben sich ziemlich geändert.

Wenn wir die änderung der Mitte der Funktion irgendwo etwas außerhalb unserer Streuung der Punkte:

x0, y0 = -0.3, 1.1

Bekommen wir mit kompletter Unsinn als Ergebnis in der Gegenwart von Rauschen! (Es funktioniert immer noch einwandfrei, ohne Geräusche.)

True parameters:  [1, -0.3, 1.1, 2, 3, 4]
Predicted params: [  0.546  -0.939   0.857  -0.488  44.069  -4.136]
Residual, RMS(obs - pred): 0.235664449826

Dies ist ein häufiges problem beim Einbau einer Funktion, die zerfällt zu null. Jeglicher Lärm in den "tails" führen kann, in einem sehr schlechten Ergebnis. Es gibt eine Reihe von Strategien für den Umgang mit diesem. Eine der einfachsten ist es, das Gewicht der inversion durch die beobachteten z-Werte. Hier ist ein Beispiel für den 1D-Fall: (mit dem Fokus auf der linearisierten problem) Wie kann ich eine least-squares fitting über mehrere Daten-sets schnell? Wenn ich Zeit habe, später werde ich ein Beispiel für den 2D-Fall.