2

我正在尝试计算指数的自然对数加上数组的中位数，但浮点数的精度必须为 2000，否则答案将始终为 0。

这是我到目前为止所拥有的：

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

1 回答 1

0

我想我能够使用 gmpy2 而不是 BigFloat 来转置上面的 R 代码，但可能也可以使用：

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)
于 2017-06-03T21:13:34.097 回答