2016-07-08 15 views
0

Also habe ich versucht, zu einer exponentiell modifizierten Gauß-Funktion zu passen (wenn Interesse besteht, https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution)Fitting weise definierte Funktion: Laufen in NaN Probleme aufgrund boolean Array-Indizierung

import numpy as np 
import scipy.optimize as sio 
import scipy.special as sps 

def exp_gaussian(x, h, u, c, t): 
    z = 1.0/sqrt(2.0) * (c/t - (x-u)/c) #not important 
    k1 = k2 = h * c/t * sqrt(pi/2) #not important 
    n1 = 1/2 * (c/t)**2 - (x-u)/t  #not important 
    n2 = -1/2 * ((x - u)/c)**2  #not important 
    y = np.zeros(len(x)) 
    y += (k1 * np.exp(n1) * sps.erfc(z)) * (z < 0) 
    y += (k2 * np.exp(n2) * sps.erfcx(z)) * (z >= 0) 
    return y 

Um Überlaufprobleme zu vermeiden, Abhängig davon, ob z positiv oder negativ ist, muss eine von zwei äquivalenten Funktionen verwendet werden (siehe Alternative Formen zur Berechnung aus der vorherigen Wikipedia-Seite).

Das Problem, das ich habe, ist dies: Die Zeile y += (k2 * np.exp(n2) * sps.erfcx(z)) * (z >= 0) soll nur zu y hinzufügen, wenn z positiv ist. Aber wenn z ist, sagen wir, -30, ist inf, und inf * False ist NaN. Anstatt y unangetastet zu lassen, wird das resultierende y mit NaN gruppiert. Beispiel:

x = np.linspace(400, 1000, 1001) 
y = exp_gaussian(x, 100, 400, 10, 5) 
y 
array([ 84.27384586, 86.04516723, 87.57518493, ...,   nan, 
       nan,   nan]) 

Ich habe versucht, das Ersetzen der Linie in Frage mit dem folgenden:

y += numpy.nan_to_num((k2 * np.exp(n2) * sps.erfcx(z)) * (z >= 0)) 

Aber tun dies lief in schweren Laufzeitprobleme. Gibt es eine Möglichkeit, (k2 * np.exp(n2) * sps.erfcx(z)) nur unter der Bedingung (z >= 0) auszuwerten? Gibt es eine andere Möglichkeit, dies zu lösen, ohne die Laufzeit zu beeinträchtigen?

Danke!

EDIT: Nach Rishi Rat scheint der folgende Code viel besser zu arbeiten:

def exp_gaussian(x, h, u, c, t): 
    z = 1.0/sqrt(2.0) * (c/t - (x-u)/c) 
    k1 = k2 = h * c/t * sqrt(pi/2) 
    n1 = 1/2 * (c/t)**2 - (x-u)/t 
    n2 = -1/2 * ((x - u)/c)**2 
    return = np.where(z >= 0, k2 * np.exp(n2) * sps.erfcx(z), k1 * np.exp(n1) * sps.erfc(z)) 

Antwort

3

Wie wäre es mit numpy.where mit so etwas wie: np.where(z >= 0, sps.erfcx(z), sps.erfc(z)). Ich bin kein Experte, weiß also nicht, ob es effizient ist. Sieht zumindest elegant aus!

+0

Das schien recht gut zu funktionieren. Das einzige Problem ist, dass wenn ich np.where (z> = 0, sps.erfcx (z), sps.erfc (z)) ausführen, bekomme ich Laufzeitwarnungen, weil es noch sps.erfcx (z) und sps.erfc auswertet (z). Gibt es eine einfache Möglichkeit, den falschen Parameter nicht zu bewerten? – Zhaitan

+0

Gut zu wissen, dass es funktioniert hat. Nur neugierig: Wie hast du Laufzeitwarnungen bekommen? Weil, ja, das logische Ding zu erwarten wäre, dass es sowohl erfcx als auch erfc auf dem gesamten z auswerten würde, aber dann np.where würde zwischen diesen zwei Bewertungen basierend auf dem Vorzeichen von z Element-für-Element auswählen. Und so würden die NaNs in erfxx verworfen werden. – Rishi

0

Eine Sache, die Sie tun können, ist eine Maske erstellen und wiederverwenden, damit es nicht zweimal ausgewertet werden müssten. Eine weitere Idee ist es, die nan_to_num nur einmal am Ende zu verwenden

mask = (z<0) 
y += (k1 * np.exp(n1) * sps.erfc(z)) * (mask) 
y += (k2 * np.exp(n2) * sps.erfcx(z)) * (~mask) 
y = numpy.nan_yo_num(y) 

versuchen und sehen, ob das hilft ...