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