Come faccio a calcolare il logaritmo di 1 meno l'esponente di un dato numero piccolo in python

Question

Sto eseguendo un calcolo di probabilità. Ho molti numeri molto piccoli, tutti dei quali voglio sottrarre da 1, e farlo in modo accurato. Posso calcolare con precisione il logaritmo di questi piccoli numeri. La mia strategia è stata finora in questo modo (usando NumPy):

dato un array del registro dei piccoli numeri x, calcolare:

y = numpy.logaddexp.reduce(x) 

Ora voglio calcolare qualcosa come 1-exp(y) o meglio ancora log(1-exp(y)) , ma non sono sicuro di come farlo senza perdere la precisione.

Infatti, anche la funzione logaddexp è in esecuzione in problemi di precisione. I valori nel vettore x possono variare da -2 a -800 o anche più negativi. Il vettore y di cui sopra avrebbe fondamentalmente un'intera sezione di numeri intorno a 1e-16, che è il eps del tipo di dati. Così, ad esempio, i dati con precisione calcolati potrebbero apparire così:

In [358]: x 
Out[358]: 
[-5.2194676211172837, 
-3.9050377656308362, 
-3.1619783292449615, 
-2.71289594096134, 
-2.4488395891021639, 
-2.3129210706827568, 
-2.2709987626652346, 
-2.3007776073511259, 
-2.3868404149802434, 
-2.5180718876609163, 
-2.68619816583087, 
-2.8849022632856958, 
-3.1092603032627686, 
-3.3553673369747834, 
-3.6200806272462351, 
-3.9008385919463073, 
-4.1955300857178379, 
-4.5023981074719899, 
-4.8199676154248081, 
-5.1469905756384904, 
-5.4824035553480428, 
-5.8252945959126876, 
-6.174877049340779, 
-6.5304687083067563, 
-6.8914750074202473, 
-7.25737538919104, 
-7.6277121540338797, 
-8.0020812775389558, 
-8.3801247986220773, 
-8.7615244716292437, 
-9.1459964426584435, 
-9.5332867613176404, 
-9.9231675781398394, 
-10.315433907978701, 
-10.709900863130784, 
-11.106401278287066, 
-11.50478366390567, 
-11.904910436107656, 
-12.30665638039909, 
-12.709907313918777, 
-13.114558916892051, 
-13.52051570882999, 
-13.927690148982549, 
-14.336001843810081, 
-14.745376846921289, 
-15.155747039147968, 
-15.567049578271309, 
-15.979226409456359, 
-16.39222382873956, 
-16.805992092998878, 
-17.22048507074976, 
-17.63565992888303, 
-18.051476851117201, 
-18.467898784496384, 
-18.884891210740903, 
-19.302421939667397, 
-19.720460922243518, 
-20.138980081145718, 
-20.557953156947775, 
-20.977355568292495, 
-21.397164284594595, 
-21.817357709992422, 
-22.237915577412224, 
-22.658818851739369, 
-23.080049641202237, 
-23.501591116172762, 
-23.923427434676114, 
-24.345543673975158, 
-24.767925767665417, 
-25.190560447772668, 
-25.61343519140047, 
-26.036538171518259, 
-26.459858211524278, 
-26.883384743252066, 
-27.307107768123842, 
-27.731017821180984, 
-28.155105937748402, 
-28.579363622513654, 
-29.003782820820732, 
-29.428355891997484, 
-29.853075584553352, 
-30.27793501309668, 
-30.702927636836705, 
-31.128047239545907, 
-31.553287910869187, 
-31.978644028878307, 
-32.404110243774596, 
-32.82968146265631, 
-33.255352835270173, 
-33.681119740674262, 
-34.106977774747804, 
-34.532922738484046, 
-34.958950627012712, 
-35.385057619298891, 
-35.811240068471022, 
-36.237494492735493, 
-36.663817566835519, 
-37.090206114019054, 
-37.516657098479527, 
-37.943167618239784, 
-38.369734898447348, 
-38.796356285056333, 
-39.223029238868548, 
-39.64975132991276, 
-40.076520232137909, 
-40.5033337184027, 
-40.930189655741344, 
-41.357086000888444, 
-41.784020796047173, 
-42.210992164885965, 
-42.637998308748706, 
-43.065037503066776, 
-43.492108093959985, 
-43.919208495015312, 
-44.346337184233221, 
-44.773492701130749, 
-45.200673643993753, 
-45.627878667267964, 
-46.055106479082156, 
-46.482355838895614, 
-46.909625555262096, 
-47.336914483704675, 
-47.764221524695017, 
-48.191545621730768, 
-48.618885759506213, 
-49.04624096217151, 
-49.473610291673936, 
-49.900992846179292, 
-50.328387758566748, 
-50.755794194994508, 
-51.183211353532613, 
-51.610638462858901, 
-52.0380747810147, 
-52.46551959421754, 
-52.892972215728378, 
-53.320431984769073, 
-53.747898265489198, 
-54.175370445978274, 
-54.602847937323247, 
-55.030330172705362, 
-55.457816606538813, 
-55.885306713645889, 
-56.312799988467418, 
-56.740295944308855, 
-57.167794112617116, 
-57.59529404228897, 
-58.02279529900909, 
-58.450297464615232, 
-58.877800136490578, 
-59.305302926981085, 
-59.732805462838542, 
-60.160307384683506, 
-60.587808346493375, 
-61.015308015110463, 
-61.442806069768608, 
-61.87030220164138, 
-62.297796113406662, 
-62.725287518829532, 
-63.15277614236129, 
-63.580261718755196, 
-64.007743992695964, 
-64.435222718445743, 
-64.862697659501919, 
-65.290168588270035, 
-65.717635285748088, 
-66.14509754122389, 
-66.572555151982783, 
-67.000007923029216, 
-67.427455666815376, 
-67.854898202982099, 
-68.282335358110231, 
-68.709766965479957, 
-69.137192864839108, 
-69.564612902180784, 
-69.992026929530198, 
-70.419434804735829, 
-70.8468363912732, 
-71.274231558051156, 
-71.701620179229167, 
-72.129002134037705, 
-72.556377306608397, 
-72.983745585807242, 
-73.411106865077045, 
-73.838461042282461, 
-74.265808019561746, 
-74.693147703185559, 
-75.120480003416901, 
-75.547804834380145, 
-75.97512211393132, 
-76.402431763534764, 
-76.829733708143749, 
-77.257027876085431, 
-77.684314198948414, 
-78.111592611476681, 
-78.538863051464546, 
-78.966125459656723, 
-79.393379779652037, 
-79.820625957809625, 
-80.24786394315754, 
-80.675093687306912, 
-81.102315144366912] 

poi cerco di calcolare la somma registro degli esponenti:

In [359]: np.logaddexp.accumulate(x) 
Out[359]: 
array([ -5.21946762e+00, -3.66710221e+00, -2.68983273e+00, 
     -2.00815067e+00, -1.51126604e+00, -1.14067818e+00, 
     -8.60829425e-01, -6.48188808e-01, -4.86276416e-01, 
     -3.63085873e-01, -2.69624488e-01, -1.99028599e-01, 
     -1.45996863e-01, -1.06408884e-01, -7.70565672e-02, 
     -5.54467248e-02, -3.96506186e-02, -2.81859503e-02, 
     -1.99225261e-02, -1.40061296e-02, -9.79701394e-03, 
     -6.82045164e-03, -4.72733966e-03, -3.26317960e-03, 
     -2.24396350e-03, -1.53767347e-03, -1.05026994e-03, 
     -7.15209142e-04, -4.85690052e-04, -3.28980607e-04, 
     -2.22305294e-04, -1.49890553e-04, -1.00858788e-04, 
     -6.77380054e-05, -4.54139175e-05, -3.03974537e-05, 
     -2.03154477e-05, -1.35581905e-05, -9.03659252e-06, 
     -6.01552344e-06, -3.99984336e-06, -2.65671945e-06, 
     -1.76283376e-06, -1.16860435e-06, -7.73997496e-07, 
     -5.12213574e-07, -3.38706792e-07, -2.23809375e-07, 
     -1.47785898e-07, -9.75226648e-08, -6.43149957e-08, 
     -4.23904687e-08, -2.79246430e-08, -1.83858489e-08, 
     -1.20995365e-08, -7.95892319e-09, -5.23300609e-09, 
     -3.43929670e-09, -2.25953475e-09, -1.48391255e-09, 
     -9.74194956e-10, -6.39351406e-10, -4.19466218e-10, 
     -2.75121795e-10, -1.80397409e-10, -1.18254918e-10, 
     -7.74993004e-11, -5.07775611e-11, -3.32619009e-11, 
     -2.17835737e-11, -1.42634249e-11, -9.33764336e-12, 
     -6.11190167e-12, -3.99989955e-12, -2.61737204e-12, 
     -1.71253165e-12, -1.12043465e-12, -7.33052079e-13, 
     -4.79645919e-13, -3.13905885e-13, -2.05519681e-13, 
     -1.34650094e-13, -8.83173582e-14, -5.80300378e-14, 
     -3.82338678e-14, -2.52963381e-14, -1.68421145e-14, 
     -1.13181549e-14, -7.70918073e-15, -5.35155125e-15, 
     -3.81152630e-15, -2.80565548e-15, -2.14872312e-15, 
     -1.71971577e-15, -1.43957518e-15, -1.25665732e-15, 
     -1.13722927e-15, -1.05925916e-15, -1.00835857e-15, 
     -9.75131524e-16, -9.53442707e-16, -9.39286186e-16, 
     -9.30046550e-16, -9.24016349e-16, -9.20080954e-16, 
     -9.17512772e-16, -9.15836886e-16, -9.14743318e-16, 
     -9.14029759e-16, -9.13564174e-16, -9.13260398e-16, 
     -9.13062204e-16, -9.12932898e-16, -9.12848539e-16, 
     -9.12793505e-16, -9.12757603e-16, -9.12734183e-16, 
     -9.12718905e-16, -9.12708939e-16, -9.12702438e-16, 
     -9.12698198e-16, -9.12695432e-16, -9.12693627e-16, 
     -9.12692451e-16, -9.12691683e-16, -9.12691183e-16, 
     -9.12690856e-16, -9.12690643e-16, -9.12690504e-16, 
     -9.12690414e-16, -9.12690355e-16, -9.12690316e-16, 
     -9.12690291e-16, -9.12690275e-16, -9.12690264e-16, 
     -9.12690257e-16, -9.12690252e-16, -9.12690249e-16, 
     -9.12690248e-16, -9.12690246e-16, -9.12690245e-16, 
     -9.12690245e-16, -9.12690245e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16, 
     -9.12690244e-16, -9.12690244e-16, -9.12690244e-16]) 

che alla fine porta a:

In [360]: np.logaddexp.reduce(x) 
Out[360]: -9.1269024387687033e-16 

quindi la mia precisione è già stata cancellata. Qualche idea su come aggirare questo?

Answer 1

In Python 2.7, abbiamo aggiunto math.expm1() per questo caso d'uso:

>>> from math import exp, expm1 
>>> exp(1e-5) - 1 # gives result accurate to 11 places 
1.0000050000069649e-05 
>>> expm1(1e-5) # result accurate to full precision 
1.0000050000166668e-05 

Inoltre, v'è math.fsum() per la fase di somma, senza perdita di precisione:

>>> sum([.1, .1, .1, .1, .1, .1, .1, .1, .1, .1]) 
0.9999999999999999 
>>> fsum([.1, .1, .1, .1, .1, .1, .1, .1, .1, .1]) 
1.0 

Infine, se non altro aiuta , è possibile utilizzare il decimal module che supporta l'aritmetica di altissima precisione:

>>> from decimal import * 
>>> getcontext().prec = 200 
>>> (1 - 1/Decimal(7000000)).ln() 
Decimal('-1.4285715306122546161332083855139723669559469615692284955124609122046580004888309867906750714454869716398919778588515625689415322789136206397998627088895481989036005482451668027002380442299229191323673E-7') 

Answer 2

Esattamente come Raymond Hettinger detto, ma si sarebbe sicuramente bisogno di moltiplicare per -1 dopo, perché si vuole 1 - exp invece di exp - 1

Answer 3

suggerisco di sostituire exp() e log() con il loro Taylor series nelle vicinanze di 0 e 1, corrispondentemente. In questo modo, non perderai la precisione usando numeri grandi (mio, ho appena chiamato 1 un numero elevato: ^). Utilizza la formula Lagrange remainder o solo l'espressione del membro con qualche riserva per determinare da quando la discrepanza andrà oltre la tua precisione.

Aggiornamento:

Python 2.7 di math.expm1 (exp(x)-1) e math.log1p (log(1+x)) fare questo per voi se la precisione della libreria C della piattaforma (tipicamente doppia) è sufficiente. (in caso contrario, dovrai ricorrere a uno speciale software matematico (la FPU di x86 può calcolare con precisione estesa)).

Answer 4

io non sono sicuro se questo è ciò che si vuole

numpy.expm1(x[, out]) 
Calculate exp(x) - 1 for all elements in the array. 



>>> import numpy as np 
>>> np.expm1(x).sum() 
-200.0 
>>> (-np.expm1(x)).sum() 
200.0 
>>> from scipy import special 
>>> (-special.expm1(x)).sum() 
200.0 
>>> np.log((-special.expm1(x)).sum()) 
5.2983173665480363 

Edit:

Siamo spiacenti, non mi rendevo conto che questa è solo la versione NumPy della risposta di Raymond Hettinger.

(non una risposta al problema numerico)

io non sono ancora sicuro di cosa esattamente la domanda è, però, invece di buttare decimale o mpmath a lui, forse una riformulazione del problema sarà di aiuto. Se sommate le probabilità in Poisson, per esempio, alla fine arriverete sempre "vicini" a 1. Ma per alcune domande possiamo evitare alcuni problemi di lavoro con la funzione di sopravvivenza invece del cdf.

Answer 5

Ho utilizzato mpmath per problemi simili. È una libreria Python 100% molto buona e ben documentata. Se la velocità è un problema per te; considera l'utilizzo di mpmath con gmpy. Vedi questo link per farlo.

Answer 6

Non so molto di Python, faccio gran parte del mio lavoro in Java ma mi sembra che sarebbe meglio eseguire prima il trucco log-sum-exp su tutti i valori, piuttosto che farlo due -by-two utilizzando numpy.logaddexp.accumulate.

Una rapida ricerca in google mostra un candidato tra le librerie python: scipy.misc.logsumexp.

In ogni caso che non sarebbe difficile da programmare te stesso:

logsumexp(probs) := max(probs) + log(sum[i](exp(probes[i]-max(probs)))) 

così qualcosa di simile:

maxValue = -Inf; 
    for x in probs 
    if x > maxValue then maxValue = x 
    expSum = 0 
    for x in probs 
    expSum += exp(x - maxValue) 
    return log(expSum) 

Il valore singolo restituito, dicono p, è solo il registro del somma di tutte le probabilità nei problemi . Si noti che se vi è una grande scala diversa tra il valore più grande e quello minore nelle probabilità di input, i piccoli saranno ignorati, ma solo se il loro contributo è molto piccolo rispetto ai grandi numeri che nella maggior parte delle applicazioni dovrebbero andare bene.

È possibile utilizzare strategie più complesse per rendere conto di quei piccoli valori nel caso in cui vi sia un numero elevato di numeri piccoli, ad esempio probe = 0,5 + 1E-7 + 1E-7 ... un milione di volte, quindi aggiungere fino a 0,1 . Quello che puoi fare è dividere le singole somme in diverse in cui i valori approssimativamente della stessa scala vengono sommati prima di essere combinati. oppure puoi utilizzare un valore di pivot intermedio anziché il massimo, ma devi assicurarti che il valore più grande non sia troppo grande perché exp (probs [i] - pivot) causerebbe un overflow in quel caso.

Una volta fatto questo è ancora necessario fare il calcolo log (1-exp (p))

Per questo ho trovato questo documento che descrive un modo per evitare la perdita di precisione, per quanto possibile utilizzando funzioni logistiche che puoi trovare nella maggior parte delle librerie comuni di matematica della lingua.

Maechler M, Accurately Computing log(1 − exp(− |a|)) Assessed by the Rmpfr package, 2012

La chiave è quella di utilizzare uno dei due possibili approcci a seconda del valore di ingresso un.

Definizioni:

log1mexp(a) := log(1-exp(a)) ### function that we seek to implement. 
log1p(a) := log(1+a) # function provided by common math libraries. 
exp1m(a) := exp(a) - 1 # function provided by common math libraries. 

C'è un modo ovvio per implementare log1mexp utilizzando log1p:

log1mexp(a) := log1p(-exp(a)) 

Con exp1m si può fare:

log1mexp(a) := log(-expm1(a)) 

Si dovrebbe quindi utilizzare il log1p approccio quando a < log (.5) e l'expm1 quando a> = log (.5).

log1mexp(a) := (a < log(.5)) ? log1p(-exp(a)) : log(-expm1(a)). 

Fare riferimento al collegamento esterno per ulteriori informazioni.