Bin Fu的算法實現沒有給出正確的結果

Question

我試圖用真實語言實現Bin Fu'sapproximate sum algorithm 以更好地感受它是如何工作的。對於$ s（x）= \ sum_ {i = 1}的值，計算$ \ hat {s}（x）$，$（1+ε）$近似值的算法如下：^n x_i $ （例如，這意味着$ \ hat {s}（x）$滿足：$ \ hat {s}（x）/（1+ε）\ leq s（x）\ leq（1+ \小量）\帽子{S}（x）的$ [1]）。

但是，我一定在做錯事，因爲運行我的實現並沒有給出正確的結果，例如， $ \ hat {s}（x）$我離開它並不滿足[1]。

我懷疑在我的實現下，我現在太早了，但我不明白是什麼原因造成的。

ApproxRegion<-function(x,b,delta,n){ 
    if(n<=1 || x[n]<b)   return(NULL) 
    if(x[n-1]<b)    return(c(n,n)) 
    if(x[1]>=b)    return(c(1,n)) 
    m<-2 
    while(n-m**2>0 && x[n-m**2+1]>=b) m<-m**2 
    r<-m 
    while(m>=(1+delta)){ 
     m<-sqrt(m) 
     if(n-floor(m*r)>=0 && x[n-floor(m*r)+1]>=b) r=m*r 
    } 
    return(c(n-floor(m*r)+1,n)) 
}  
ApproxSum<-function(x,n,epsilon){ 
    if(x[n]==0)   return(0) 
    delta=3*epsilon/4 
    rp<-n 
    i<-0 
    s<-0 
    b<-x[n]/(1+delta) 
    while(b>=delta*x[n]/(3*n)){ 
     R<-ApproxRegion(x,b,delta,rp) 
     if(is.null(R))  break 
     rp<-R[1]-1; 
     b<-x[rp]/(1+delta) 
     si<-(R[2]-R[1]+1)*b 
     s<-s+si 
     i<-i+1 
    } 
    return(list(s=s,i=i)) 
} 

然而，當我運行它

n<-100; 
set.seed(123) 
x<-sort(rexp(n)); 
eps<-1/10 
y0<-ApproxSum(x=x,n=n,epsilon=eps); 
y0$s*(1+eps) 
sum(x) 

我得到y0$s*(1+eps)小於sum(x)

Answer 1

你好像失去跟蹤我的VS i + 1的兩個地方，第二次在ApproxRegion中循環，在ApproxSum中循環。這看起來像你的例子：

ApproxRegion<-function(x,b,delta,n){ 
    if(n<=1 || x[n]<b)   return(NULL) 
    if(x[n-1]<b)    return(c(n,n)) 
    if(x[1]>=b)    return(c(1,n)) 
    m<-2 
    while(n-m**2>0 && x[n-m**2+1]>=b) m<-m**2 
    r<-m 
    while(m>=(1+delta)){ 
     m<-sqrt(m) 
     if(n-floor(m*r)>=0 && x[n-floor(m*r)+1]>=b) r=m*r 
    } 
    return(c(n-floor(r)+1,n)) 
} 
ApproxSum<-function(x,n,epsilon){ 
    if(x[n]==0)   return(0) 
    delta=3*epsilon/4 
    rp<-n 
    i<-0 
    s<-0 
    b<-x[n]/(1+delta) 
    while(b>=delta*x[n]/(3*n)){ 
     R<-ApproxRegion(x,b,delta,rp) 
     if(is.null(R))  break 
     si<-(R[2]-R[1]+1)*b 
     s<-s+si 
     i<-i+1 
     rp<-R[1]-1; 
     b<-x[rp]/(1+delta) 
    } 
    return(list(s=s,i=i)) 
} 

n<-100; 
set.seed(123) 
x<-sort(rexp(n)); 
eps<-0.001 
y0<-ApproxSum(x=x,n=n,epsilon=eps); 


> y0$s*(1+eps) 
[1] 104.5955 

> sum(x) 
[1] 104.5719 

> y0$s/(1+eps) 
[1] 104.3866