Rcpp函数：f（x） - f（x）不是0

Question

由于计算速度的原因，我使用Rcpp重写了一个纯R函数。

下面是函数：

set.seed(1) 
X <- cbind(1, rnorm(1000)) 
Y <- X %*% c(0, 1) + rnorm(1000) 
b0 <- rnorm(4) 
es_loss_c(b0, Y, X, 0.01, 100) 
[1] 64.22035 

到目前为止是这种情况，结果是一样的（最多1E-15）与：

// [[Rcpp::export]] 
double es_loss_c(NumericVector b, NumericVector Y, NumericMatrix X, double alpha, double delta) { 

    int n = X.nrow();       // Number of observations 
    int k = X.ncol();       // Number of parameters 

    NumericVector b_q = b[Range(0, k-1)];  // Quantile parameters 
    NumericVector b_e = b[Range(k, 2*k-1)];  // Expected shortfall parameters 

    // Initialize required variables 
    NumericVector Xi; 
    double out, Yi, Xb_q, Xb_e, u, exp_e, H; 

    for (int i = 0; i < n; i++) { 

    Yi = Y(i); 
    Xi = X(i,_); 

    // Pre-store some values 
    Xb_q = sum(Xi * b_q); 
    Xb_e = sum(Xi * b_e); 

    u = Yi - Xb_q; 
    exp_e = exp(Xb_e); 

    // Indicator function or its approximation 
    if (delta > 0) { 
     H = 1/(1 + exp(delta * u)); 
    } else { 
     H = u <= 0; 
    } 

    // Shortfall loss 
    out += (alpha - H) * u + exp_e/(1 + exp_e) * (-1/alpha * H * u + Xb_e - Xb_q) - log(1 + exp_e); 
    } 

    return out/n; 
} 

你可以用测试功能纯粹的R实现。

但随后，一些真正奇怪的事情发生：

es_loss_c(b0, Y, X, 0.01, 100) - es_loss_c(b0, Y, X, 0.01, 100) 
[1] -64.22035 
es_loss_c(b0, Y, X, 0.01, 100) + es_loss_c(b0, Y, X, 0.01, 100) 
[1] 192.6611 
es_loss_c(b0, Y, X, 0.01, 100)/es_loss_c(b0, Y, X, 0.01, 100) 
[1] 0.5 

你能向我解释这里发生了什么？

Answer 1

您在使用前忘记初始化out = 0;。只需在for循环前添加它，每次拨打es_loss_c将产生相同的结果。