divisione NumPy con RuntimeWarning: valore non valido incontrato nel double_scalars

Question

ho scritto il seguente script:

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 

Ma ho ottenuto questo risultato e con l'errore si è verificato:

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

Sembra essere che la gli elementi di input erano troppo piccoli che python li ha trasformati in zeri, ma in effetti la divisione ha il suo risultato.

Come risolvere questo tipo di problema?

Answer 1

non puoi risolverlo. Semplicemente answer1.sum()==0, e non è possibile eseguire una divisione per zero.

Questo succede perché answer1 è l'esponenziale di 2 numeri molto grandi, negativi, in modo che il risultato sia arrotondato a zero.

nan viene restituito in questo caso a causa della divisione per zero.

Ora per risolvere il problema si potrebbe:

• andare per una libreria per la matematica ad alta precisione, come mpmath. Ma è meno divertente.
• in alternativa a un'arma più grande, eseguire alcune manipolazioni matematiche, come descritto di seguito.
• andare per una funzione su misura scipy/numpy che fa esattamente quello che vuoi! Dai un'occhiata alla risposta di @Warren Weckesser.

Qui spiego come eseguire alcune manipolazioni matematiche che aiutano questo problema. Abbiamo che per il numeratore:

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))) 

cui sopra x=3* 1089 e y=3* 1093. Ora, l'argomento di questa esponenziale è

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

Per il denominatore si potrebbe procedere allo stesso modo, ma ottenere che log(1+exp(-z+k)) è già arrotondato a 0, in modo che l'argomento della funzione esponenziale al denominatore è semplicemente arrotondato a -z=-3000 . Devi quindi che il risultato è

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

che è già molto vicino al risultato che si otterrebbe se si dovesse mantenere solo i 2 termini principali (vale a dire il primo numero 1089 al numeratore e il primo numero 1000 al denominatore):

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

Per il gusto di farlo, vediamo quanto siamo vicini dalla soluzione di Wolfram alpha (link):

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

La differenza tra questo numero e l'esponente sopra è +1.7053025658242404e-13, quindi l'approssimazione che abbiamo fatto al denominatore andava bene.

Il risultato finale è

'exp(-266.99999385580668) = 1.1050349147204485e-116 

Da wolframio alfa è (link)

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

e nuovo, è sicuro da usare NumPy anche qui.

Answer 2

È possibile utilizzare np.logaddexp (che implementa l'idea di @ gg349 risposta):

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 

Oppure si può utilizzare scipy.misc.logsumexp:

In [52]: from scipy.misc import logsumexp 

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

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