18

次のような非常に小さな数値のリストを計算する必要があります

(0.1)^1000、0.2^(1200)、

そしてそれらを正規化して、合計が1になるようにします。

a1 = 0.1^1000、a2 = 0.2^1200

そして、a1' = a1/(a1+a2), a2'=a2(a1+a2) を計算したいと思います。

a1=0 になるため、アンダーフローの問題が発生しています。どうすればこれを回避できますか? 理論的にはログを処理することができ、log(a1) = 1000*log(0.l) はアンダーフローの問題なしに a1 を表す方法になります - しかし、正規化するには log(a1+a2) を取得する必要があります - a1 を直接表すことができないため、計算できません。

私はRでプログラミングしています-私が知る限り、C＃にはDecimalのような倍精度値よりも優れたデータ型はありません。

任意の提案をいただければ幸いです、ありがとう

4

4 に答える 4

16

数学的に言えば、これらの数字の 1 つが appx になります。ゼロ、そしてもう一つ。あなたの数字の違いは大きいので、これが理にかなっているのかどうかさえ疑問に思っています.

logspace_addしかし、一般的にそれを行うには、R のフードの下にある C 関数のアイデアを使用できます。 logxpy ( =log(x+y) )whenlx = log(x)およびly = log(y)asを定義できます。

logxpy <- function(lx,ly) max(lx,ly) + log1p(exp(-abs(lx-ly)))

つまり、次を使用できます。

> la1 <- 1000*log(0.1)
> la2 <- 1200*log(0.2)

> exp(la1 - logxpy(la1,la2))
[1] 5.807714e-162

> exp(la2 - logxpy(la1,la2))
[1] 1

より多くの数値がある場合は、この関数を再帰的に呼び出すこともできます。1 は 1 であり、 1 - ではありません5.807...e-162。本当にもっと精度が必要で、プラットフォームが long double 型をサポートしている場合は、すべてを C や C++ などでコーディングし、後で結果を返すことができます。しかし、私が正しければ、R は - 今のところ - 通常の double しか扱えないので、最終的には結果が表示されたときに再び精度を失うことになります。

編集 ：

あなたのために計算をする：

log(x+y) = log(exp(lx)+exp(ly))
         = log( exp(lx) * (1 + exp(ly-lx) )
         = lx + log ( 1 + exp(ly - lx)  )

ここで、最大のものを lx として、式 in に到達しlogxpy()ます。

EDIT 2：なぜ最大を取るのですか？exp(lx-ly) で負の数を使用することを確認するのは簡単です。lx-ly が大きくなりすぎると、exp(lx-ly) は Inf を返します。それは正しい結果ではありません。exp(ly-lx) は 0 を返すため、はるかに優れた結果が得られます。

lx=1 と ly=1000 とすると、次のようになります。

> 1+log1p(exp(1000-1))
[1] Inf
> 1000+log1p(exp(1-1000))
[1] 1000
于 2011-04-27T11:40:03.857 に答える
8

このBrobdingnagパッケージは、非常に大きな数または小さな数を扱い、基本的に Joris の回答を便利な形式にまとめています。

a1 <- as.brob(0.1)^1000
a2 <- as.brob(0.2)^1200
a1_dash <- a1 / (a1 + a2)
a2_dash <- a2 / (a1 + a2)
as.numeric(a1_dash)
as.numeric(a2_dash)
于 2011-04-27T12:59:14.443 に答える
2

任意精度のパッケージを試してください:

  • Rmpfr「R MPFR - 複数精度浮動小数点の信頼性」
  • Ryacas「「Yacas」コンピューター代数システムへのRインターフェイス」-任意の精度を実行できる場合もあります。
于 2011-04-27T12:29:40.103 に答える
1

a1 と a2 を分数として扱うことができるかもしれません。あなたの例では、

a1 = (a1num/a1denom)^1000  # 1/10
a2 = (a2num/a2denom)^1200  # 1/5

あなたが到着するだろう

a1' = (a1num^1000 * a2denom^1200)/(a1num^1000 * a2denom^1200 + a1denom^1000 * a2num^1200)
a2' = (a1denom^1000 * a2num^1200)/(a1num^1000 * a2denom^1200 + a1denom^1000 * a2num^1200)

gmp パッケージを使用して計算できます。

library(gmp)
a1 <- as.double(pow.bigz(5,1200) / (pow.bigz(5,1200)+ pow.bigz(10,1000)))
于 2011-04-27T11:49:57.053 に答える