Probleme mit RCPP Präzision

Question

Betrachten Sie die folgende C++ Funktion in R RCPP mit:

cppFunction('long double statZn_cpp(NumericVector dat, double kn) { 
    double n = dat.size(); 
    // Get total sum and sum of squares; this will be the "upper sum" 
    // (i.e. the sum above k) 
    long double s_upper, s_square_upper; 
    // The "lower sums" (i.e. those below k) 
    long double s_lower, s_square_lower; 
    // Get lower sums 
    // Go to kn - 1 to prevent double-counting in main 
    // loop 
    for (int i = 0; i < kn - 1; ++i) { 
    s_lower += dat[i]; 
    s_square_lower += dat[i] * dat[i]; 
    } 
    // Get upper sum 
    for (int i = kn - 1; i < n; ++i) { 
    s_upper += dat[i]; 
    s_square_upper += dat[i] * dat[i]; 
    } 
    // The maximum, which will be returned 
    long double M = 0; 
    // A candidate for the new maximum, used in a loop 
    long double M_candidate; 

    // Compute the test statistic 
    for (int k = kn; k <= (n - kn); ++k) { 
    // Update s and s_square for both lower and upper 
    s_lower += dat[k-1]; 
    s_square_lower += dat[k-1] * dat[k-1]; 
    s_upper -= dat[k-1]; 
    s_square_upper -= dat[k-1] * dat[k-1]; 

    // Get estimate of sd for this k 
    long double sdk = sqrt((s_square_lower - pow(s_lower, 2.0)/k + 
         s_square_upper - 
         pow(s_upper, 2.0)/(n - k))/n); 
    M_candidate = abs(s_lower/k - s_upper/(n - k))/sdk; 
    // Choose new maximum 
    if (M_candidate > M) { 
     M = M_candidate; 
    } 
    } 

    return M * sqrt(kn); 
}') 

Sie den Befehl statZn_cpp(1:20,4), und Sie werden 6.963106 bekommen, was die richtige Antwort ist. Skalierung sollte keine Rolle spielen; statZn_cpp(1:20*10,4) ergibt auch die richtige Antwort von 6.963106. Aber statZn_cpp(1:20/10,4) ergibt die falsche Antwort von 6.575959, und statZn_cpp(1:20/100,4) gibt Ihnen wieder die offensichtlich falsche Antwort von 0. Mehr auf den Punkt (und relevant für meine Forschung, die Simulationsstudien beinhaltet), wenn ich statZn_cpp(rnorm(20),4) versuche, ist die Antwort fast immer 0, die falsch ist.

Offensichtlich hat das Problem mit Rundungsfehlern zu tun, aber ich weiß nicht, wo sie sind oder wie man sie behebt (ich bin brandneu in C++). Ich habe versucht, die Präzision so weit wie möglich zu erweitern. Gibt es eine Möglichkeit, das Rundungsproblem zu beheben? (Eine R-Wrapper-Funktion ist zulässig, wenn ich versuchen sollte, was einem Vorverarbeitungsschritt entspricht, aber robust sein muss, um für allgemeine Genauigkeitsgrade zu arbeiten.) EDIT: Hier ist ein "äquivalenter" R-Code:

statZn <- function(dat, kn = function(n) {floor(sqrt(n))}) { 
    n = length(dat) 
    return(sqrt(kn(n))*max(sapply(
     floor(kn(n)):(n - floor(kn(n))), function(k) 
     abs(1/k*sum(dat[1:k]) - 
       1/(n-k)*sum(dat[(k+1):n]))/sqrt((sum((dat[1:k] - 
       mean(dat[1:k]))^2)+sum((dat[(k+1):n] - 
       mean(dat[(k+1):n]))^2))/n)))) 
} 

auch die R nachfolgende Code grundsätzlich repliziert die Methode, die durch den C++ Code verwendet werden soll. Es ist in der Lage, die richtige Antwort zu erreichen.

n = length(dat) 
    s_lower = 0 
    s_square_lower = 0 
    s_upper = 0 
    s_square_upper = 0 
    for (i in 1:(kn-1)) { 
    s_lower = s_lower + dat[i] 
    s_square_lower = s_square_lower + dat[i] * dat[i] 
    } 
    for (i in kn:n) { 
    s_upper = s_upper + dat[i] 
    s_square_upper = s_square_upper + dat[i] * dat[i] 
    } 
    M = 0 

    for (k in kn:(n-kn)) { 
    s_lower = s_lower + dat[k] 
    s_square_lower = s_square_lower + dat[k] * dat[k] 
    s_upper = s_upper - dat[k] 
    s_square_upper = s_square_upper - dat[k] * dat[k] 

    sdk = sqrt((s_square_lower - (s_lower)^2/k + 
         s_square_upper - 
         (s_upper)^2/(n-k))/n) 
    M_candidate = sqrt(kn) * abs(s_lower/k - s_upper/(n - k))/sdk 

    cat('k', k, '\n', 
     "s_lower", s_lower, '\n', 
     's_square_lower', s_square_lower, '\n', 
     's_upper', s_upper, '\n', 
     's_square_upper', s_square_upper, '\n', 
     'sdk', sdk, '\n', 
     'M_candidate', M_candidate, '\n\n') 

    if (M_candidate > M) { 
     M = M_candidate 
    } 
    } 

Answer 1

1: Sie sollten nicht long double verwenden, da R alle numerischen Werte in der double Typ darstellt. Die Verwendung eines präziseren Typs für Zwischenberechnungen ist äußerst unwahrscheinlich und führt eher zu seltsamen Inkonsistenzen zwischen den Plattformen.

2: Sie sind nicht s_upper Initialisierung s_square_upper, s_lower und s_square_lower. (Sie sind Initialisierung sie tatsächlich in der R Implementierung, sondern Sie vergessen haben, in der C++ Implementierung.)

3: Minor Punkt, aber ich würde die pow(x,2.0) Anrufe mit x*x ersetzen. Obwohl das nicht wirklich wichtig ist.

4: Das ist, was es für mich behoben: Sie müssen Aufrufe an C++ - Standardbibliotheksfunktionen mit ihrem enthaltenden Namespace qualifizieren. IOW, std::sqrt() anstelle von nur sqrt(), std::abs() anstelle von nur abs(), und std::pow() anstelle von nur pow(), wenn Sie es weiterhin verwenden.

---

cppFunction('double statZn_cpp(NumericVector dat, double kn) { 
    int n = dat.size(); 
    double s_upper = 0, s_square_upper = 0; // Get total sum and sum of squares; this will be the "upper sum" (i.e. the sum above k) 
    double s_lower = 0, s_square_lower = 0; // The "lower sums" (i.e. those below k) 
    for (int i = 0; i < kn - 1; ++i) { s_lower += dat[i]; s_square_lower += dat[i] * dat[i]; } // Get lower sums; Go to kn - 1 to prevent double-counting in main 
    for (int i = kn - 1; i < n; ++i) { s_upper += dat[i]; s_square_upper += dat[i] * dat[i]; } // Get upper sum 
    double M = 0; // The maximum, which will be returned 
    double M_candidate; // A candidate for the new maximum, used in a loop 
    // Compute the test statistic 
    for (int k = kn; k <= (n - kn); ++k) { 
    // Update s and s_square for both lower and upper 
    s_lower += dat[k-1]; 
    s_square_lower += dat[k-1] * dat[k-1]; 
    s_upper -= dat[k-1]; 
    s_square_upper -= dat[k-1] * dat[k-1]; 
    // Get estimate of sd for this k 
    double sdk = std::sqrt((s_square_lower - s_lower*s_lower/k + s_square_upper - s_upper*s_upper/(n - k))/n); 
    M_candidate = std::abs(s_lower/k - s_upper/(n - k))/sdk; 
    if (M_candidate > M) M = M_candidate; // Choose new maximum 
    } 
    return std::sqrt(kn) * M; 
}'); 

statZn_cpp(1:20,4); ## you will get 6.963106, which is the correct answer 
## [1] 6.963106 
statZn_cpp(1:20*10,4); ## Scaling should not matter; will also yield the correct answer of 6.963106 
## [1] 6.963106 
statZn_cpp(1:20/10,4); ## yields the wrong answer of 6.575959 
## [1] 6.963106 
statZn_cpp(1:20/100,4); ## again gives you the obviously wrong answer of 0. 
## [1] 6.963106 
set.seed(1L); statZn_cpp(rnorm(20),4); ## More to the point (and relevant to my research, which involves simulation studies), the answer is almost always 0, which is wrong. 
## [1] 1.270117