J'essaie de créer une fonction de mode ultra rapide pour R à utiliser pour agréger de grands ensembles de données catégoriques. La fonction doit prendre une entrée vectorielle de tous les types R pris en charge et retourner le mode. J'ai lu Ce message a >, Cette page d'aide et d'autres, mais je ne l'étais pas capable de faire en sorte que la fonction accepte tous les types de données R. Mon code fonctionne maintenant pour les vecteurs numériques, je m'appuie sur les fonctions de wrapper de sucre Rcpp:
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
int Mode(NumericVector x, bool narm = false)
{
if (narm) x = x[!is_na(x)];
NumericVector ux = unique(x);
int y = ux[which_max(table(match(x, ux)))];
return y;
}
De plus je me demandais si l'argument ' narm ' peut être renommé ' na.rm 'sans donner d'erreurs, et bien sûr s'il existe un moyen plus rapide de coder une fonction de mode en C ++, je serais reconnaissant de le savoir.
3 Réponses :
Dans votre fonction Mode , puisque vous appelez principalement des fonctions de wrapper de sucre, vous ne verrez pas beaucoup d'amélioration par rapport à la base R . En fait, en écrivant simplement une traduction de base R fidèle, nous avons:
## base R returned in 0.943s system.time(fastIntMode(s)) user system elapsed 0.217 0.006 0.224
Et l'analyse comparative, nous avons:
fastIntMode(s) ##[1] 85433 baseMode(s) ##[1] 85433 Mode(s) ##[1] 85433
Typiquement, quand nous entreprendre l'effort d'écrire notre propre code compilé, nous nous attendrions à des gains plus importants. Emballer simplement ces fonctions compilées déjà efficaces dans Rcpp ne vous rapportera pas comme par magie les gains que vous attendez. En fait, sur de plus grands exemples, la solution de base est plus rapide. Observez:
microbenchmark(Mode(s), baseMode(s), fastIntMode(s), times = 15, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
Mode(s) 6.428343 6.268131 6.622914 6.134388 6.881746 7.78522 15
baseMode(s) 9.757491 9.404101 9.454857 9.169315 9.018938 10.16640 15
fastIntMode(s) 1.000000 1.000000 1.000000 1.000000 1.000000 1.00000 15
Pour répondre à votre question d'écrire une fonction de mode plus rapide, nous pouvons utiliser std :: unordered_map , qui est très similaire à table sous le capot (c'est-à-dire qu'ils sont tous les deux des tables de hachage en leur cœur). De plus, puisque vous renvoyez un seul entier, nous pouvons supposer en toute sécurité que nous pouvons remplacer NumericVector par IntegerVector et aussi que vous n'êtes pas concerné par le renvoi tous les em> qui se produit le plus.
L'algorithme ci-dessous peut être modifié pour renvoyer le true mode , mais je laisserai cela comme un exercice (indice: vous aurez besoin de std :: vector avec une sorte d'action lorsque it-> seconde == myMax ). N.B. vous devrez également ajouter // [[Rcpp :: plugins (cpp11)]] en haut de votre fichier cpp pour std :: unordered_map et auto .
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::plugins(cpp11)]]
#include <unordered_map>
// [[Rcpp::export]]
int fastIntMode(IntegerVector x, bool narm = false) {
if (narm) x = x[!is_na(x)];
int myMax = 1;
int myMode = 0;
std::unordered_map<int, int> modeMap;
modeMap.reserve(x.size());
for (std::size_t i = 0, len = x.size(); i < len; ++i) {
auto it = modeMap.find(x[i]);
if (it != modeMap.end()) {
++(it->second);
if (it->second > myMax) {
myMax = it->second;
myMode = x[i];
}
} else {
modeMap.insert({x[i], 1});
}
}
return myMode;
}
Et les benchmarks:
set.seed(1234) sBig <- sample(1e6, replace = TRUE) system.time(Mode(sBig)) user system elapsed 1.410 0.036 1.450 system.time(baseMode(sBig)) user system elapsed 0.915 0.025 0.943
Maintenant, nous parlons ... environ 6 fois plus rapide que l'original et 9x plus rapide que la base. Ils renvoient tous la même valeur:
set.seed(1234)
s <- sample(1e5, replace = TRUE)
library(microbenchmark)
microbenchmark(Mode(s), baseMode(s), times = 10, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
Mode(s) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 10
baseMode(s) 1.490765 1.645367 1.571132 1.616061 1.637181 1.448306 10
Et pour notre exemple plus large:
baseMode <- function(x, narm = FALSE) {
if (narm) x <- x[!is.na(x)]
ux <- unique(x)
ux[which.max(table(match(x, ux)))]
}
En plus de l'implicite Licence CC-BY-SA Par la présente, j'accepte le code de cette réponse sous la GPL> = 2 . p >
Il est également recommandé d'utiliser boost :: unordered_map au lieu de std :: unordered_map pour plus d'amélioration de la vitesse.
Merci @ F.Privé, je vais examiner cela
Merci beaucoup @JosephWood pour cette réponse étendue! Je pensais déjà que le sucre Rcpp ne mènerait pas loin. Je suppose qu'avec le mode vrai vous voulez dire retourner plusieurs modes le cas échéant ... J'ai juste besoin de l'un d'entre eux. Dans le cas où vous savez également comment faire la méthode dispatch, c'est-à-dire travailler avec des facteurs, des caractères et des vecteurs logiques, ce serait super, car c'est vraiment ce que je passe beaucoup de temps à essayer de faire sans réussir.
Merci aussi @ F.Privé pour cette suggestion, j'ai essayé de la coder: #include Mais ce n'était pas plus rapide: environ 5% plus lent que le code ci-dessus.
Merci pour les commentaires.
Afin de faire fonctionner la fonction pour n'importe quelle entrée vectorielle, vous pouvez implémenter l'algorithme de @ JosephWood pour tout type de données que vous souhaitez prendre en charge et l'appeler à partir d'un commutateur (TYPEOF (x)) . Mais ce serait beaucoup de duplication de code. Au lieu de cela, il est préférable de créer une fonction générique qui peut fonctionner sur n'importe quel argument Vector. Si nous suivons le paradigme de R selon lequel tout est un vecteur et laissons la fonction renvoyer également un Vector, alors nous pouvons utiliser RCPP_RETURN_VECTOR . Notez que nous avons besoin de C ++ 11 pour pouvoir passer des arguments supplémentaires à la fonction appelée par RCPP_RETURN_VECTOR . Une chose délicate est que vous avez besoin du type de stockage pour Vector afin de créer un std :: unordered_map approprié. Ici, Rcpp :: traits :: storage_type vient à la rescousse. Cependant, std :: unordered_map ne sait pas comment gérer les nombres complexes de R. Pour simplifier, je désactive ce cas particulier.
Pour tout rassembler:
> set.seed(1234) > s <- sample(1e5, replace = TRUE) > fastMode(s) [1] 85433 > fastMode(s + 0.1) [1] 85433.1 > l <- sample(c(TRUE, FALSE), 11, replace = TRUE) > fastMode(l) [1] TRUE > c <- sample(letters, 1e5, replace = TRUE) > fastMode(c) [1] "z" > f <- as.factor(c) > fastMode(f) [1] z Levels: a b c d e f g h i j k l m n o p q r s t u v w x y z
Résultat:
#include <Rcpp.h>
using namespace Rcpp ;
// [[Rcpp::plugins(cpp11)]]
#include <unordered_map>
template <int RTYPE>
Vector<RTYPE> fastModeImpl(Vector<RTYPE> x, bool narm){
if (narm) x = x[!is_na(x)];
int myMax = 1;
Vector<RTYPE> myMode(1);
// special case for factors == INTSXP with "class" and "levels" attribute
if (x.hasAttribute("levels")){
myMode.attr("class") = x.attr("class");
myMode.attr("levels") = x.attr("levels");
}
std::unordered_map<typename Rcpp::traits::storage_type<RTYPE>::type, int> modeMap;
modeMap.reserve(x.size());
for (std::size_t i = 0, len = x.size(); i < len; ++i) {
auto it = modeMap.find(x[i]);
if (it != modeMap.end()) {
++(it->second);
if (it->second > myMax) {
myMax = it->second;
myMode[0] = x[i];
}
} else {
modeMap.insert({x[i], 1});
}
}
return myMode;
}
template <>
Vector<CPLXSXP> fastModeImpl(Vector<CPLXSXP> x, bool narm) {
stop("Not supported SEXP type!");
}
// [[Rcpp::export]]
SEXP fastMode( SEXP x, bool narm = false ){
RCPP_RETURN_VECTOR(fastModeImpl, x, narm);
}
/*** R
set.seed(1234)
s <- sample(1e5, replace = TRUE)
fastMode(s)
fastMode(s + 0.1)
l <- sample(c(TRUE, FALSE), 11, replace = TRUE)
fastMode(l)
c <- sample(letters, 1e5, replace = TRUE)
fastMode(c)
f <- as.factor(c)
fastMode(f)
*/
Comme indiqué ci-dessus, l'algorithme utilisé provient de Réponse de Joseph Wood , qui a fait l'objet d'une double licence explicite sous CC-BY-SA et GPL> = 2. Je suis Joseph et j'autorise par la présente le code de cette réponse sous GPL (version 2 ou ultérieure) en plus de la licence implicite CC-BY-SA.
Merci beaucoup @RalfStubner !! Je n'aurais jamais réussi cela moi-même. Votre fonction est environ 25 fois plus rapide sur de grandes données catégorielles que la solution de base, c'est parfait! La seule chose que j'ajouterais au code est par défaut la première valeur si le mode ne peut pas être trouvé, c'est-à-dire initialiser `Vector
@Sebastian Vous êtes les bienvenus. Le mérite de la vitesse revient à Joseph Wood, cependant.
@RalfStubner, excellent travail comme toujours! J'ai en fait pensé à vous pendant que j'écrivais ma réponse car je sais que vous travaillez sur l'implémentation d'une fonction de hachage pour une fonction d'échantillonnage dans dqrng . J'ai hâte de voir le résultat final!
C'est une solution vraiment rapide. Comment devrions-nous adopter le code pour renvoyer la fréquence du mode comme attribut "freq"?
@AndriSignorell Vous pouvez ajouter myMode.attr ("freq") = myMax; avant de renvoyer myMode .
Marche parfaitement! Merci, j'utiliserai votre code pour la fonction DescTools :: Mode et je ferai référence à vous dans cet article. J'espère que ça ira.
@AndriSignorell C'est plus compliqué qu'il n'y paraît :-( L'algorithme initial a été fourni par Joseph Wood (malheureusement je ne peux pas @ -mention deux personnes) et (implicitement) sous licence CC-BY-SA 3.0. J'ai adapté l'algorithme pour utiliser arbitraire entrée, toujours avec la même licence (implicite). Récemment, SO a changé la licence pour tout le contenu en CC-BY-SA 4.0. Qu'ils aient ou non le droit de le faire est une autre question. Désormais, DescTools est sous licence GPL> = 2. Il serait probablement préférable que Joseph et moi obtenions une nouvelle licence pour le code.
@JosephWood Andri souhaite utiliser notre code dans un package sous licence GPL. Ce serait plus facile si nous devions ré-licencier le code sous GPL> = 2. Je serais d'accord avec ça, mais vous devriez commencer car mon travail est basé sur le vôtre ;-)
@RalfStubner, ça sonne bien! Je ne sais pas vraiment comment procéder. Je sais que je peux obtenir une licence pour mon code dans un dépôt git. Pour ce que cela vaut, je suivrai volontiers votre exemple.
@JosephWood Il devrait suffire de lui indiquer la licence supplémentaire sur SO. De préférence dans une modification de votre réponse, car les commentaires peuvent être supprimés / modifiés. Quelque chose comme ce qui suit devrait faire (IANAL!): En plus de la licence implicite CC-BY-SA, j'accepte par la présente le code de cette réponse sous GPL (version 2 ou ultérieure).
@RalfStubner, comment est-ce? :)
Mieux encore, merci d'avoir clarifié la situation juridique complexe. Je vais donc vous citer tous les deux en tant qu'auteurs avec un lien vers l'article.
Pour faire suite à une auto-promotion éhontée, j'ai maintenant publié un package collapse sur CRAN qui comprend un ensemble complet de Fonctions statistiques rapides , parmi lesquelles le générique fonction fmode . L'implémentation est basée sur le hachage d'index et encore plus rapide que la solution ci-dessus. fmode peut être utilisé pour effectuer des calculs en mode simple, groupé et / ou pondéré sur des vecteurs, des matrices, des data.frames et des tibbles groupés dplyr. Syntaxe:
collap(data, ~ id1 + id2, FUN = fmean, catFUN = fmode)
où x est un vecteur, une matrice, data.frame ou grouped_df, g est un vecteur de regroupement ou une liste de regroupement vecteurs, et w est un vecteur de poids. Une solution compacte aux problèmes d'agrégation catégorielle et mixte est en outre fournie par la fonction collap . Le code
fmode(x, g = NULL, w = NULL, ...)
agrège les données de type mixte data.frame data en appliquant fmean au numérique et fmode code> aux colonnes catégorielles. Des appels plus personnalisés sont également possibles. Associé aux Fonctions statistiques rapides , collap est tout aussi rapide que data.table sur des données numériques volumineuses, et les agrégations catégorielles et pondérées sont nettement plus rapides que tout ce qui peut actuellement être fait avec data.table.
@JosephWood a répondu à votre question sur les performances. Qu'avez-vous essayé concernant la prise en charge des types de vecteurs de tout type?
J'ai essayé d'utiliser un modèle avec une macro de type basée sur ce message . Avec le code de Joseph, je commencerais comme ce modèle
Vector Mode_temp (Vector x, bool narm = false) {if (narm) x = x [! Is_na (x)]; int myMax = 1; monMode = x [0]; std :: unordered_map <, > modeMap; // [[Rcpp :: export]] SEXP Mode (SEXP x) {RCPP_RETURN_VECTOR (Mode_temp, x); }, mais ne comprend pas la syntaxe, et si le monde fonctionne avec unordered_map.Une meilleure option me semble la syntaxe 'switch' décrite sur cette page mais ici j'ai aussi eu des problèmes avec l'entrée vectorielle, la sortie scalaire et la gestion du code en interne (c'est-à-dire en créant des vecteurs et des salars d'un certain type dans chaque cas de commutateur.)