R言語での極小値の指数形式での出力

質問

R言語での極小値の指数形式での出力に関して質問です。
Windows vista, Rのバージョンは x64 2.12.1　です
統計に関する計算処理で exp関数を使うスクリプトを作っています。

入力する値を　a とすると
計算する式は
1 - exp(-exp( a )) となります(R script上)

a <- -10 とすると計算結果は
4.53988992012988e-5
となりますが

aが小さい数になってくると
a <- -40
計算結果は 0 となってしまいます。

本来指数表示の極小数が　1.0e-50のように表示されるべきなのですが切り捨てられて0として扱われてしまっているようです。
オプションで変更できるのかと思いoptions の　scipen　をいじってみましたがうまくいきません。
このような計算の結果の極小数を指数表示（1.0e-50のように）するにはどうすればいいのでしょうか？

------------スクリプト--------------
a <- -4
b <- -38
result_a <- 1- exp(-exp( a ))
result_b <- 1- exp(-exp( b ))
cat(result_a , "\n", result_b )

-----------出力---------------
0.01814893
0

回答No.4

Tacosan · ベストアンサー

expm1 なんてのもあるか?

お礼

expm1初めて知りました。
説明を読んで、ただexp(a)-1というだけかと思いましたが、小さい数字を扱う場合（今回のよう場合）に正確な値を算出してくれるようですね。

(a <- -4 の場合)
(exp使用)
> 1 - exp(- exp(-4))
[1] 0.01814893
(expm1使用)
> - expm1(- exp(-4))
[1] 0.01814893

(a <- -38 の場合)
(exp使用)
> 1 - exp(- exp(-38))
[1] 0
(expm1使用)
> - expm1(- exp(-38))
[1] 3.139133e-17

のようにすることで特にオプションを設定することなく求める値を得ることが出来ました。（元の式を変えているのでちょっと不思議な感じですが。。。）

回答No.3

ki073

なぜa=-40のときに答えは0.0になるか分かりましたか？
Rの計算は倍精度の浮動小数点演算で計算されますので有効桁数が15桁程度です。
-exp(-exp( a ))を計算しますと1.00000と0が続き、18桁目に0でない数字が現れます。有効桁数15桁だと1.0になるわけです。

もっと有効桁数をとれれば計算は可能です。Rではできるのかどうか分かりません。

Rubyではかなり大きな桁数をとることができ、手間はかかりますがこんな式で計算できます。有効桁数を100桁にした例です。
require "bigdecimal"
require "bigdecimal/math"
include BigMath
prec=100 # 有効桁数
a=-40
b=BigDecimal.new("1.0", prec)-exp(-exp(BigDecimal.new(a.to_s,prec), prec),prec)
puts b.to_f
答えは4.24835425529159e-18
となります。
Rも統計計算やグラフ作成に使っていますが、データ構造が複雑でとてもついていけません。
Rubyで前処理をしてそこからRを呼び出して統計計算やグラフ作成をしています。

ANo.2の理由ですが、ヒントはテイラー展開です。

お礼

回答ありがとうございます。

内部演算桁数が十分でなく、計算途中で0に丸められているのかな。。。というのが現在の理解です。

>もっと有効桁数をとれれば計算は可能です。Rではできるのかどうか分かりません。
出力用に有効桁数を変えることはできるようですが、Rの内部演算桁数に関しては今も調べている状態です。

>Rubyではかなり大きな桁数をとることができ、手間はかかりますがこんな式で計算できます。
普段はC++も使っているので、C++上で計算させようかなとも思っていました。なるべくR上で簡潔に済ませたいのですが。

>ANo.2の理由ですが、ヒントはテイラー展開です。
早速調べ始めましたが、なかなか難しそうですね。。。

回答No.2

f272

それくらいaが小さいときなら
1 - exp(-exp( a )) = exp( a )
で十分と思う。

お礼

実際にR上で計算させると
1 - exp(-exp( a ))
exp( a )
は非常に近い値になりますね。
exp(-40)も4.24835425529159e-18と値が返されます。

ただ、なぜ　1 - exp(-exp( a )) = exp( a )　と言えるのか理由が思いつきません。もしよければ簡単に説明してもらえないでしょうか？

回答No.1

Tacosan

計算精度は?

補足

resultの値はe-100くらいになることを想定しています。
ですのでその桁を切り捨て、丸め込みを行わないよう計算させたいのですが。。。

これで答えになっているでしょうか？