ゼロとの浮動小数点の比較

Question

私はGEV分布の分位数を計算する関数を書いています。この質問に関連する態様は、パラメータ（形状パラメータまたはカッパ）の一方がゼロ

これは一般にアドレス指定されるプログラムである場合、機能の異なる形態が必要とされることです次のように（：qgevとlmomcoが似ている:: quagevこれはEVDからの抜粋です）：

（編集：lmomcoのバージョン2.2.2、この問題の識別問題に対処しています）

if (shape == 0) 
     return(loc - scale * log(-log(p))) 
    else return(loc + scale * ((-log(p))^(-shape) - 1)/shape) 

これは、shape/kappaがゼロに正確に等しいがゼロ付近で奇妙な振る舞いがある場合にうまくいきます。同じ答えがゼロに近く、ゼロで返されるので、これは正常に見える

Qgev_zero <- function(shape){ 
    # p is an exceedance probability 
    p= 0.01 
    location=0 
    scale=1 

    if(shape == 0) return(location - scale*(log(-log(1-p)))) 

    location + (scale/shape)*((-log(1-p))^-shape - 1) 

} 

Qgev_zero(0) 
#[1] 4.600149 

Qgev_zero(1e-8) 
#[1] 4.600149 

：

は例を見てみましょう。しかし、ゼロに近いものを見てください。

k.seq <- seq(from = -4e-16, to = 4e-16, length.out = 1000) 
plot(k.seq, sapply(k.seq, Qgev_zero), type = 'l') 

関数によって返される値は、多くの場合、間違っています振動します。

ゼロとの直接比較をall.equalに置き換えると、これらの問題はなくなります。 all.equalのヘルプを見ると

if(isTRUE(all.equal(shape, 0))) return(location - scale*(log(-log(1-p)))) 

は、デフォルト値については、1.5E-8よりも小さいものをゼロとして扱われることを示唆しています。

もちろん、ゼロ付近の奇妙な振る舞いは一般的には問題ではありませんが、私の場合は既知の分位数​​からパラメータを決定するために最適化/根検索を使用しています。

質問には：all.equal(target, 0)この問題に対処する適切な方法を使用していますか？このアプローチが日常的に使われていないのはなぜですか？

Answer 1

一部の関数は、浮動小数点表現で明白な方法で実装されたとき、特定のポイントで不正な動作をします。これは特に、関数を1つのポイントで手動で定義しなければならない場合に当てはまります。あるポイントで事態が絶対に定義されていない場合、近い将来親愛なる人生にぶつかっている可能性が高いです。

この場合、それはκの負の指数と戦うκの分母からのものです。どちらが勝利するかは、ビットごとに決定され、それぞれが時には「より大きな規模への丸め」コンテストに勝つこともあります。

これらの問題を解決するには、さまざまなアプローチがありますが、すべてがケースバイケースで設計されています。しばしば欠陥があり実装が容易なアプローチの1つは、問題のあるポイントの近くで、より適切な表現（例えば、カッパに関するテイラー展開）に切り替えることです。境界線に不連続性が導入されます。必要に応じて、2つの間で補間を試すことができます。

Answer 2

Sneftelの提案に従うと、k = -1e-7とk = 1e-7で分位数を計算し、kの引き数がこれらの限度の間に入るかどうかを補間します。これはうまくいくようです。

私はlmomcoからGEV分位関数のパラメータ化を使用しています。このコードでは

:: quagev

（編集：lmomcoのバージョン2.2.2はこの質問で特定された問題に対処してい）

関数Qgevは問題のあるバージョン（プロットの黒線）ですが、Qgev_interpはゼロに近い値を補間します（プロットの緑線）。

Qgev <- function(K, f, XI, A){ 

# K = shape 
# f = probability 
# XI = location 
# A = scale 

    Y <- -log(-log(f)) 
    Y <- (1-exp(-K*Y))/K 
    x <- XI + A*Y 
    return(x) 
} 

Qgev_interp <- function(K, f, XI, A){ 

    .F <- function(K, f, XI, A){ 
    Y <- -log(-log(f)) 
    Y <- (1-exp(-K*Y))/K 
    x <- XI + A*Y 
    return(x) 
    } 


    k1 <- -1e-7 
    k2 <- 1e-7 
    y1 <- .F(k1, f, XI, A) 
    y2 <- .F(k2, f, XI, A) 

    F_nearZero <- approxfun(c(k1, k2), c(y1, y2)) 

    if(K > k1 & K < k2) { 
    return(F_nearZero(K)) 
    } else { 
    return(.F(K, f, XI, A)) 
    } 
} 



k.seq <- seq(from = -1.1e-7, to = 1.1e-7, length.out = 1000) 
plot(k.seq, sapply(k.seq, Qgev, f = 0.01, XI = 0, A = 1), col=1, lwd = 1, type = 'l') 
lines(k.seq, sapply(k.seq, Qgev_interp, f = 0.01, XI = 0, A = 1), col=3, lwd = 2)