2017-06-03 41 views
1

我試圖計算指數的自然對數加上數組的中值,但浮點數的精度必須在2000或答案將始終爲0.來自BigFloats數組的Python Numpy指數值

這是我到目前爲止有:

import bigfloat  
x = np.array([-15349.79, -15266.66, -15242.86]) 
answer = np.median(x) - bigfloat.log(np.mean(bigfloat.exp(x, precision(2000)) + np.median(x))) 

此代碼返回下面的錯誤,因爲BigFloat.exp不能在列表類型的工作。

__main__:1: RuntimeWarning: invalid value encountered in log 
    array([ nan, nan, nan]) 
    >>> np.median - bigfloat.log(np.mean(bigfloat.exp(x, precision(2000)) + np.median(x))) 
    Traceback (most recent call last): 
     File "<stdin>", line 1, in <module> 
     File "/usr/local/lib/python2.7/dist-packages/bigfloat/core.py", line 1446, in exp 
     (BigFloat._implicit_convert(x),), 
     File "/usr/local/lib/python2.7/dist-packages/bigfloat/core.py", line 800, in _implicit_convert 
     "to BigFloat" % (arg, type(arg))) 
    TypeError: Unable to convert argument [-15349.79 -15266.66 -15242.86] of type <type 'numpy.ndarray'> to BigFloat 

然後,我試着列表理解建立精度= 2000指數,但輸出失去了精度。

>>>[bigfloat.exp(num, precision(2000)) + BigFloat(np.median(x), context=precision(2000)) for num in x] 

[BigFloat.exact('-15266.660000000000', precision=53), BigFloat.exact('-15266.660000000000', precision=53), BigFloat.exact('-15266.660000000000', precision=53)] 

我可以從x使用列表理解 這樣得到的指數值精度:

>>> y = np.array([bigfloat.exp(num, precision(2000)) for num in x]) 
array([ BigFloat.exact('4.687104391719151845495683439738733373338442839115295629901427982078731632624437676931490944129698360820181432182432247030174990996132679553420878801399467284014932435857939862423137788809352433385343338286781879066456873472082636971652587852379587642184738259777920214145750081457830886095059600271300801524086458081526778407096875577469007859748418787409261916565428551066536598870494431653463691369520685259389041691328858076222360659135651318185915247243253437425313718584237418827954902711646850362440592578566737383074556346759332830833101196990845567845361685702120838400980658468192376607029699380e-6667', precision=2000), 
     BigFloat.exact('5.940252515613784074556309936460240719247289162721506036367173706446251980105141516508485353816420203092949808160441872102510858074752258842471232457791645630942356791030979819957875784926363292788234347364203226596687259399429746917920741954706671428927952815041593801046941161417184730521922391569252658359604561007344640221209247257690437938341582528542446104212627600044421405229891732857592357725820529732658234033631318425476205945557122313924830225909708184718714804450362854280687393958953216361511561704836206669636791590655443171684811683402330053964681285971765503804198242627954802184910024902e-6631', precision=2000), 
     BigFloat.exact('1.288289823457680769007768910906908351089465066737359292851136781233384628626343352992636825492031571595103435334370758196286825224816434202324210124540257467624499933926757723584914581286627071375803686797647455160335628799775858142489656885909172399293492295645703203222218182084249732700966821340431452575815746282651823925677898549539269454644804048695031225101971491985645951419388454738446335860070391719831039721150295437452237572831818257568221336263814212095143032293694959570063677473370798577031762607527221992020661837810811705007892502559603984057933079724597470162879821292957122400421187875e-6620', precision=2000)], dtype=object) 

所以,我怎麼指數陣列中添加np.median(x)到每一個項目,我如何獲得最終數組中元素的日誌?有沒有更簡單的方法來計算上面第一個代碼片段中變量answer給出的方程?

基本上,我想轉換該R代碼到Python:

llMed <- stats::median(x) 
metric <- as.double(
    llMed - log(Rmpfr::mean(exp(-Rmpfr::mpfr(x, prec=2000L) + llMed))) 
) 

回答

0

我想我是能夠轉以上將R代碼,使用gmpy2代替BigFloat但可能也可以使用兩種:

import gmpy2 as gp 
import numpy as np  
with gp.local_context(gp.context(), precision=2000) as ctx: 
    x = np.array([gp.mpfr(-15349.79), gp.mpfr(-15266.66), gp.mpfr(-15242.86)]) 
    answer = np.median(x) - gp.log(np.mean([gp.exp(gp.add(np.median(x), -gp.mpfr(item))) for item in x])) 

輸出的回報:

mpfr('-15348.69138771133276342351845677418587273364155071266454924792376077085597654860712492594599688058199377288639137666410888137048513058096745776113277020531535761588878042016412651412489608285582998650481415959482636259292554205272401992800167437149224493900214813760731711328834885525224938584036572786202764947113186177391441168267495103157033223800214492991771147327262071249020907408433551462034872024117419365924381891832161483692689444158313503155805562703358104122358438183824459422779020196496507161021965435781332319270106503099589290788685588200038739438059696970511568363022088057740627040541967',2000)