2015-08-08 13 views
7

मुझे लगता है कि ggplot2 में एक उत्तल पतवार एक scatterplot के समूह द्वारा के रूप मेंआर: एक 2d या 3 डी scatterplot

library(ggplot2) 
library(plyr) 
data(iris) 
df<-iris 
find_hull <- function(df) df[chull(df$Sepal.Length, df$Sepal.Width), ] 
hulls <- ddply(df, "Species", find_hull) 
plot <- ggplot(data = df, aes(x = Sepal.Length, y = Sepal.Width, colour=Species, fill = Species)) + 
    geom_point() + 
    geom_polygon(data = hulls, alpha = 0.5) + 
    labs(x = "Sepal.Length", y = "Sepal.Width") 
plot 

enter image description here

में जोड़ सकते हैं पता करने के लिए अल्फा बैग जोड़ने मैं हालांकि कैसे एक सोच रहा था इसके बजाय अल्फा बैग की गणना और जोड़ सकते हैं, यानी सबसे बड़ा उत्तल झुकाव जिसमें कम से कम अनुपात के सभी अनुपात 1-अल्फा शामिल हैं? या तो 2 डी में (ggplot2 के साथ प्रदर्शित करने के लिए) या 3 डी (आरजीएल के साथ प्रदर्शित करने के लिए)।

संपादित करें: मेरा प्रारंभिक विचार कम से कम दिए गए अंकों के मानदंड के साथ-साथ कन्वेक्स हुल को "छीलने" पर रखा जाना चाहिए, हालांकि कागज here में ऐसा लगता है कि वे एक अलग उपयोग करते हैं एल्गोरिदम (आईएसओडीपीएचटी, जो कि आर पैकेज पैकेज गहराई में कार्यान्वित प्रतीत होता है, फ़ंक्शन isodepth और aplpack::plothulls में जो भी मैं चाहता हूं उसके करीब भी लगता है (हालांकि यह केवल समोच्च के विपरीत एक पूर्ण साजिश उत्पन्न करता है), इसलिए मुझे लगता है कि इन्हें हल किया जा सकता है हालांकि ये फ़ंक्शन केवल 2 डी में काम करता है, और मुझे 3 डी एक्सटेंशन (आरजीएल में प्लॉट करने के लिए) में दिलचस्पी होगी। अगर किसी के पास कोई पॉइंटर्स है तो मुझे बताएं!

EDIT2: फ़ंक्शन depth::isodepth मुझे 2 डी समाधान मिला (se ई पोस्ट नीचे), हालांकि मैं अभी भी एक 3 डी समाधान की तलाश में हूं - अगर किसी को यह पता चल जाएगा कि ऐसा कैसे करें, तो कृपया मुझे बताएं!

उत्तर

3

हम के लिए एक पैरामीटर स्वीकार करने के लिए aplpack::plothulls समारोह संशोधित कर सकते हैं संलग्न करने के लिए अंक का अनुपात (एप्लैक में यह 50% पर सेट है)। फिर हम ggplot के लिए एक geom कस्टम बनाने के लिए इस संशोधित फ़ंक्शन का उपयोग कर सकते हैं।

यहां कस्टम geom है:

library(ggplot2) 
StatBag <- ggproto("Statbag", Stat, 
        compute_group = function(data, scales, prop = 0.5) { 

        ################################# 
        ################################# 
        # originally from aplpack package, plotting functions removed 
        plothulls_ <- function(x, y, fraction, n.hull = 1, 
              col.hull, lty.hull, lwd.hull, density=0, ...){ 
         # function for data peeling: 
         # x,y : data 
         # fraction.in.inner.hull : max percentage of points within the hull to be drawn 
         # n.hull : number of hulls to be plotted (if there is no fractiion argument) 
         # col.hull, lty.hull, lwd.hull : style of hull line 
         # plotting bits have been removed, BM 160321 
         # pw 130524 
         if(ncol(x) == 2){ y <- x[,2]; x <- x[,1] } 
         n <- length(x) 
         if(!missing(fraction)) { # find special hull 
         n.hull <- 1 
         if(missing(col.hull)) col.hull <- 1 
         if(missing(lty.hull)) lty.hull <- 1 
         if(missing(lwd.hull)) lwd.hull <- 1 
         x.old <- x; y.old <- y 
         idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx] 
         for(i in 1:(length(x)/3)){ 
          x <- x[-idx]; y <- y[-idx] 
          if((length(x)/n) < fraction){ 
          return(cbind(x.hull,y.hull)) 
          } 
          idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx]; 
         } 
         } 
         if(missing(col.hull)) col.hull <- 1:n.hull 
         if(length(col.hull)) col.hull <- rep(col.hull,n.hull) 
         if(missing(lty.hull)) lty.hull <- 1:n.hull 
         if(length(lty.hull)) lty.hull <- rep(lty.hull,n.hull) 
         if(missing(lwd.hull)) lwd.hull <- 1 
         if(length(lwd.hull)) lwd.hull <- rep(lwd.hull,n.hull) 
         result <- NULL 
         for(i in 1:n.hull){ 
         idx <- chull(x,y); x.hull <- x[idx]; y.hull <- y[idx] 
         result <- c(result, list(cbind(x.hull,y.hull))) 
         x <- x[-idx]; y <- y[-idx] 
         if(0 == length(x)) return(result) 
         } 
         result 
        } # end of definition of plothulls 
        ################################# 


        # prepare data to go into function below 
        the_matrix <- matrix(data = c(data$x, data$y), ncol = 2) 

        # get data out of function as df with names 
        setNames(data.frame(plothulls_(the_matrix, fraction = prop)), nm = c("x", "y")) 
        # how can we get the hull and loop vertices passed on also? 
        }, 

        required_aes = c("x", "y") 
) 

#' @inheritParams ggplot2::stat_identity 
#' @param prop Proportion of all the points to be included in the bag (default is 0.5) 
stat_bag <- function(mapping = NULL, data = NULL, geom = "polygon", 
        position = "identity", na.rm = FALSE, show.legend = NA, 
        inherit.aes = TRUE, prop = 0.5, alpha = 0.3, ...) { 
    layer(
    stat = StatBag, data = data, mapping = mapping, geom = geom, 
    position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
    params = list(na.rm = na.rm, prop = prop, alpha = alpha, ...) 
) 
} 


geom_bag <- function(mapping = NULL, data = NULL, 
        stat = "identity", position = "identity", 
        prop = 0.5, 
        alpha = 0.3, 
        ..., 
        na.rm = FALSE, 
        show.legend = NA, 
        inherit.aes = TRUE) { 
    layer(
    data = data, 
    mapping = mapping, 
    stat = StatBag, 
    geom = GeomBag, 
    position = position, 
    show.legend = show.legend, 
    inherit.aes = inherit.aes, 
    params = list(
     na.rm = na.rm, 
     alpha = alpha, 
     prop = prop, 
     ... 
    ) 
) 
} 

#' @rdname ggplot2-ggproto 
#' @format NULL 
#' @usage NULL 
#' @export 
GeomBag <- ggproto("GeomBag", Geom, 
        draw_group = function(data, panel_scales, coord) { 
        n <- nrow(data) 
        if (n == 1) return(zeroGrob()) 

        munched <- coord_munch(coord, data, panel_scales) 
        # Sort by group to make sure that colors, fill, etc. come in same order 
        munched <- munched[order(munched$group), ] 

        # For gpar(), there is one entry per polygon (not one entry per point). 
        # We'll pull the first value from each group, and assume all these values 
        # are the same within each group. 
        first_idx <- !duplicated(munched$group) 
        first_rows <- munched[first_idx, ] 

        ggplot2:::ggname("geom_bag", 
             grid:::polygonGrob(munched$x, munched$y, default.units = "native", 
                 id = munched$group, 
                 gp = grid::gpar(
                  col = first_rows$colour, 
                  fill = alpha(first_rows$fill, first_rows$alpha), 
                  lwd = first_rows$size * .pt, 
                  lty = first_rows$linetype 
                 ) 
            ) 
        ) 


        }, 

        default_aes = aes(colour = "NA", fill = "grey20", size = 0.5, linetype = 1, 
            alpha = NA, prop = 0.5), 

        handle_na = function(data, params) { 
        data 
        }, 

        required_aes = c("x", "y"), 

        draw_key = draw_key_polygon 
) 

और यहाँ है कि यह कैसे किया जा सकता है की एक उदाहरण है:

ggplot(iris, aes(Sepal.Length, Petal.Length, colour = Species, fill = Species)) + 
    geom_point() + 
    stat_bag(prop = 0.95) + # enclose 95% of points 
    stat_bag(prop = 0.5, alpha = 0.5) + # enclose 50% of points 
    stat_bag(prop = 0.05, alpha = 0.9) # enclose 5% of points 

enter image description here

+0

उस के लिए बहुत-बहुत धन्यवाद! यह वास्तव में उपयोगी होगा! –

6

समारोह depth::isodepth की मदद से हा मैं निम्नलिखित समाधान के साथ आया था - यहाँ मैं अल्फा बैग है कि सभी बिंदुओं का कम से कम 1-अल्फा के अनुपात होता लगता है:

library(mgcv) 
library(depth) 
library(plyr) 
library(ggplot2) 
data(iris) 
df=iris[,c(1,2,5)] 
alph=0.05 
find_bag = function(x,alpha=alph) { 
    n=nrow(x) 
    target=1-alpha 
    propinside=1 
    d=1 
    while (propinside>target) { 
     p=isodepth(x[,1:2],dpth=d,output=T, mustdith=T)[[1]] 
     ninside=sum(in.out(p,as.matrix(x[,1:2],ncol=2))*1) 
     nonedge=sum(sapply(1:nrow(p),function (row) 
      nrow(merge(round(setNames(as.data.frame(p[row,,drop=F]),names(x)[1:2]),5),as.data.frame(x[,1:2])))>0)*1)-3 
     propinside=(ninside+nonedge)/n 
     d=d+1 
    } 
    p=isodepth(x[,1:2],dpth=d-1,output=T, mustdith=T)[[1]] 
    p } 
bags <- ddply(df, "Species", find_bag,alpha=alph) 
names(bags) <- c("Species",names(df)[1:2]) 
plot <- ggplot(data = df, aes(x = Sepal.Length, y = Sepal.Width, colour=Species, fill = Species)) + 
    geom_point() + 
    geom_polygon(data = bags, alpha = 0.5) + 
    labs(x = "Sepal.Length", y = "Sepal.Width") 
plot 

enter image description here

EDIT2: उत्तल छील छीलने के अपने मूल विचार का उपयोग करके मैं निम्नलिखित समाधान के साथ आया जो अब 2 डी & 3 डी में काम करता है; परिणाम काफी एक ही isodepth एल्गोरिथ्म के साथ है नहीं है, लेकिन यह बहुत करीब है:

# in 2d 
library(plyr) 
library(ggplot2) 
data(iris) 
df=iris[,c(1,2,5)] 
alph=0.05 
find_bag = function(x,alpha=alph) { 
    n=nrow(x) 
    propinside=1 
    target=1-alpha 
    x2=x 
    while (propinside>target) { 
    propinside=nrow(x2)/n 
    hull=chull(x2) 
    x2old=x2 
    x2=x2[-hull,] 
    } 
    x2old[chull(x2old),] } 
bags <- ddply(df, "Species", find_bag, alpha=alph) 
plot <- ggplot(data = df, aes(x = Sepal.Length, y = Sepal.Width, colour=Species, fill = Species)) + 
    geom_point() + 
    geom_polygon(data = bags, alpha = 0.5) + 
    labs(x = "Sepal.Length", y = "Sepal.Width") 
plot 

enter image description here

# in 3d 
library(plyr) 
library(ggplot2) 
data(iris) 
df=iris[,c(1,2,3,5)] 
levels=unique(df[,"Species"]) 
nlevels=length(levels) 
zoom=0.8 
cex=1 
aspectr=c(1,1,0.7) 
pointsalpha=1 
userMatrix=matrix(c(0.80,-0.60,0.022,0,0.23,0.34,0.91,0,-0.55,-0.72,0.41,0,0,0,0,1),ncol=4,byrow=T) 
windowRect=c(0,29,1920,1032) 
cols=c("red","forestgreen","blue") 
alph=0.05 
plotbag = function(x,alpha=alph,grp=1,cols=c("red","forestgreen","blue"),transp=0.2) { 
    propinside=1 
    target=1-alpha 
    x2=x 
    levels=unique(x2[,ncol(x2)]) 
    x2=x2[x2[,ncol(x2)]==levels[[grp]],] 
    n=nrow(x2) 
    while (propinside>target) { 
    propinside=nrow(x2)/n 
    hull=unique(as.vector(convhulln(as.matrix(x2[,1:3]), options = "Tv"))) 
    x2old=x2 
    x2=x2[-hull,] 
    } 
    ids=t(convhulln(as.matrix(x2old[,1:3]), options = "Tv")) 
    rgl.triangles(x2old[ids,1],x2old[ids,2],x2old[ids,3],col=cols[[grp]],alpha=transp,shininess=50) 
} 

open3d(zoom=zoom,userMatrix=userMatrix,windowRect=windowRect,antialias=8) 
for (i in 1:nlevels) { 
    plot3d(x=df[df[,ncol(df)]==levels[[i]],][,1], 
     y=df[df[,ncol(df)]==levels[[i]],][,2], 
     z=df[df[,ncol(df)]==levels[[i]],][,3], 
     type="s", 
     col=cols[[i]], 
     size=cex, 
     lit=TRUE, 
     alpha=pointsalpha,point_antialias=TRUE, 
     line_antialias=TRUE,shininess=50, add=TRUE) 
    plotbag(df,alpha=alph, grp=i, cols=c("red","forestgreen","blue"), transp=0.3) } 
axes3d(color="black",drawfront=T,box=T,alpha=1) 
title3d(color="black",xlab=names(df)[[1]],ylab=names(df)[[2]],zlab=names(df)[[3]],alpha=1) 
aspect3d(aspectr) 

enter image description here