2016-07-26 17 views
0

Ich erstelle eine Funktion, die Fläche unter der Kurve berechnet und wenn ich die 2 Teiltöne nehme und sie für den Zähler multipliziere, überschreite ich 2^31 und dann wird ein Wert wie verwendet die Berechnung.Rcpp gibt eine große negative Zahl zurück, wenn 2 große Positive multipliziert werden

Hier sind die CPP-Brocken

#include <Rcpp.h> 
using namespace Rcpp; 


// [[Rcpp::export]] 
NumericVector sort_rcpp(NumericVector x) { 
    std::vector<double> tmp = Rcpp::as< std::vector<double> > (x); 
    std::sort(tmp.begin(), tmp.end()); 
    return wrap(tmp); 
} 

// [[Rcpp::export]] 
IntegerVector rank(NumericVector x) { 
    return match(x, sort_rcpp(x)); 
} 

// [[Rcpp::export]] 
double auc_(NumericVector actual, NumericVector predicted) { 

    double n = actual.size(); 

    IntegerVector Ranks = rank(predicted); 
    int NPos = sum(actual == 1); 
    int NNeg = (actual.size() - NPos); 

    int sumranks = 0; 

    for(int i = 0; i < n; ++i) { 
    if (actual[i] == 1){ 
     sumranks = sumranks + Ranks[i]; 
    } 
    } 

    double p1 = (sumranks - NPos*(NPos + 1)/2); 
    long double p2 = NPos*NNeg; 

    double auc = p1/p2; 
    return auc ; 

} 

und dann das Testbeispiel, das das Problem hat

N = 100000 
Actual = as.numeric(runif(N) > .65) 
Predicted = as.numeric(runif(N)) 

actual = Actual 
predicted = Predicted 

auc_(Actual, Predicted) 

ich auch diese in einem R-Paket am Putting

devtools::install_github("JackStat/ModelMetrics") 

N = 100000 
Actual = as.numeric(runif(N) > .65) 
Predicted = as.numeric(runif(N)) 

actual = Actual 
predicted = Predicted 

ModelMetrics::auc(Actual, Predicted) 
+1

https://xkcd.com/571/ –

Antwort

5

Sie verwenden int intern in Ihrer Funktion, was zum Überlauf führt. Verwenden Sie ein double und Dinge sehen sonnigeren:

R> sourceCpp("/tmp/jackstat.cpp") 

R> N <- 100000 

R> Actual <- as.numeric(runif(N) > .65) 

R> Predicted <- as.numeric(runif(N)) 

R> auc1(Actual, Predicted) # your function 
[1] -0.558932 

R> auc2(Actual, Predicted) # my variant using double 
[1] 0.499922 
R> 

Die gesamte korrigierte Datei unter:

#include <Rcpp.h> 

using namespace Rcpp; 

// [[Rcpp::export]] 
NumericVector sort_rcpp(NumericVector x) { 
    std::vector<double> tmp = Rcpp::as< std::vector<double> > (x); 
    std::sort(tmp.begin(), tmp.end()); 
    return wrap(tmp); 
} 

// [[Rcpp::export]] 
IntegerVector rank(NumericVector x) { 
    return match(x, sort_rcpp(x)); 
} 

// [[Rcpp::export]] 
double auc1(NumericVector actual, NumericVector predicted) { 

    double n = actual.size(); 

    IntegerVector Ranks = rank(predicted); 
    int NPos = sum(actual == 1); 
    int NNeg = (actual.size() - NPos); 

    int sumranks = 0; 

    for(int i = 0; i < n; ++i) { 
    if (actual[i] == 1){ 
     sumranks = sumranks + Ranks[i]; 
    } 
    } 

    double p1 = (sumranks - NPos*(NPos + 1)/2); 
    long double p2 = NPos*NNeg; 

    double auc = p1/p2; 
    return auc ; 

} 

// [[Rcpp::export]] 
double auc2(NumericVector actual, NumericVector predicted) { 

    double n = actual.size(); 

    IntegerVector Ranks = rank(predicted); 
    double NPos = sum(actual == 1); 
    double NNeg = (actual.size() - NPos); 

    double sumranks = 0; 

    for(int i = 0; i < n; ++i) { 
    if (actual[i] == 1){ 
     sumranks = sumranks + Ranks[i]; 
    } 
    } 

    double p1 = (sumranks - NPos*(NPos + 1)/2); 
    long double p2 = NPos*NNeg; 

    double auc = p1/p2; 
    return auc ; 

} 

/*** R 
N <- 100000 
Actual <- as.numeric(runif(N) > .65) 
Predicted <- as.numeric(runif(N)) 

auc1(Actual, Predicted) 
auc2(Actual, Predicted) 
*/ 
+2

Dirk, du bist genial! – JackStat

+0

Immer ein Vergnügen. Das war leicht zu beheben, aber hin und wieder haben wir alle Scheuklappen. Ich war dort, habe das gemacht ... –

+0

Ich baue das Paket um zu lernen, wie man Rcpp und C++ im Allgemeinen benutzt, also ist es alles neu für mich :) – JackStat