Numpy Division mit RuntimeWarning: Ungültiger Wert in double_scalars

77

Ich habe folgendes Skript geschrieben:

import numpy

d = numpy.array([[1089, 1093]])
e = numpy.array([[1000, 4443]])
answer = numpy.exp(-3 * d)
answer1 = numpy.exp(-3 * e)
res = answer.sum()/answer1.sum()
print res

Aber ich habe dieses Ergebnis bekommen und mit dem Fehler aufgetreten:

nan
C:\Users\Desktop\test.py:16: RuntimeWarning: invalid value encountered in double_scalars
  res = answer.sum()/answer1.sum()

Es scheint, dass das Eingabeelement zu klein war, als dass Python sie zu Nullen gemacht hätte, aber tatsächlich hat die Division ihr Ergebnis.

Wie kann man ein solches Problem lösen?

Heinz
quelle

Antworten:

82

Du kannst es nicht lösen. Einfach answer1.sum()==0, und Sie können keine Division durch Null durchführen.

Dies geschieht, weil answer1das Exponential von 2 sehr großen negativen Zahlen ist, so dass das Ergebnis auf Null gerundet wird.

nan wird in diesem Fall wegen der Division durch Null zurückgegeben.

Um Ihr Problem zu lösen, können Sie:

  • Suchen Sie sich eine Bibliothek für hochpräzise Mathematik wie mpmath aus . Aber das macht weniger Spaß.
  • Führen Sie als Alternative zu einer größeren Waffe einige mathematische Manipulationen durch, wie unten beschrieben.
  • Entscheiden Sie sich für eine maßgeschneiderte scipy/numpyFunktion, die genau das tut, was Sie wollen! Schauen Sie sich die Antwort von @Warren Weckesser an.

Hier erkläre ich, wie man eine mathematische Manipulation durchführt, die bei diesem Problem hilft. Wir haben das für den Zähler:

exp(-x)+exp(-y) = exp(log(exp(-x)+exp(-y)))
                = exp(log(exp(-x)*[1+exp(-y+x)]))
                = exp(log(exp(-x) + log(1+exp(-y+x)))
                = exp(-x + log(1+exp(-y+x)))

wo oben x=3* 1089und y=3* 1093. Das Argument dieses Exponentials ist nun

-x + log(1+exp(-y+x)) = -x + 6.1441934777474324e-06

Für den Nenner könnten Sie ähnlich vorgehen, aber erhalten, dass log(1+exp(-z+k))bereits gerundet ist 0, so dass das Argument der Exponentialfunktion am Nenner einfach gerundet wird -z=-3000. Sie haben dann, dass Ihr Ergebnis ist

exp(-x + log(1+exp(-y+x)))/exp(-z) = exp(-x+z+log(1+exp(-y+x)) 
                                   = exp(-266.99999385580668)

Dies liegt bereits sehr nahe an dem Ergebnis, das Sie erhalten würden, wenn Sie nur die beiden führenden Begriffe (dh die erste Zahl 1089im Zähler und die erste Zahl 1000im Nenner) beibehalten würden:

exp(3*(1089-1000))=exp(-267)

Lassen Sie uns sehen, wie nah wir an der Lösung von Wolfram alpha sind ( Link ):

Log[(exp[-3*1089]+exp[-3*1093])/([exp[-3*1000]+exp[-3*4443])] -> -266.999993855806522267194565420933791813296828742310997510523

Der Unterschied zwischen dieser Zahl und dem Exponenten oben ist +1.7053025658242404e-13, daher war die Annäherung, die wir am Nenner vorgenommen haben, in Ordnung.

Das Endergebnis ist

'exp(-266.99999385580668) = 1.1050349147204485e-116

Von Wolfram Alpha ist ( Link )

1.105034914720621496.. × 10^-116 # Wolfram alpha.

und auch hier ist es sicher, numpy zu verwenden.

gg349
quelle
Aber in diesem Fall muss ich Werte der Division durch 2 sehr kleine Werte erhalten.
Heinz
@Heinz Ich denke du meintest den Fall, in dem eine kleine Zahl durch eine kleine Zahl geteilt wird. In diesem Fall ist es viel besser, den Algorithmus so zu ändern, dass beide Zahlen vergrößert werden, als mechanische Wendungen zu finden. Nehmen Sie zum Beispiel den Logarithmus der analytischen Gleichungen, die Ihr Code zu simulieren versucht. Es gibt viele Probleme mit der Stabilität der Berechnung, wenn kleine Zahlen beteiligt sind. Es ist schöner, wenn möglich, keine davon zu haben.
Mai
20

Sie können verwenden np.logaddexp(was die Idee in der Antwort von @ gg349 implementiert):

In [33]: d = np.array([[1089, 1093]])

In [34]: e = np.array([[1000, 4443]])

In [35]: log_res = np.logaddexp(-3*d[0,0], -3*d[0,1]) - np.logaddexp(-3*e[0,0], -3*e[0,1])

In [36]: log_res
Out[36]: -266.99999385580668

In [37]: res = exp(log_res)

In [38]: res
Out[38]: 1.1050349147204485e-116

Oder Sie können verwenden scipy.special.logsumexp:

In [52]: from scipy.special import logsumexp

In [53]: res = np.exp(logsumexp(-3*d) - logsumexp(-3*e))

In [54]: res
Out[54]: 1.1050349147204485e-116
Warren Weckesser
quelle