2014-06-13 50 views
3

這是我對R社區的第一篇文章,所以如果它很愚蠢,請原諒我。我想使用ggplot2中的函數geom_density2d和stat_density2d來繪製內核密度估計值,但問題是他們無法處理加權數據。據我所知,這兩個函數從MASS包中調用函數kde2d來進行內核密度估計。而kde2d不會將數據權重作爲參數。ggplot2 - 修改geom_density2d以接受權重作爲參數?

現在,我發現kde2d http://www.inside-r.org/node/226757的這個改變版本,它以權重作爲參數,並且基於kde2d的源代碼。這個函數的代碼:

kde2d.weighted <- function (x, y, w, h, n = 25, lims = c(range(x), range(y))) { 
    nx <- length(x) 
    if (length(y) != nx) 
    stop("data vectors must be the same length") 
    if (length(w) != nx & length(w) != 1) 
    stop("weight vectors must be 1 or length of data") 
    gx <- seq(lims[1], lims[2], length = n) # gridpoints x 
    gy <- seq(lims[3], lims[4], length = n) # gridpoints y 
    if (missing(h)) 
    h <- c(bandwidth.nrd(x), bandwidth.nrd(y)); 
    if (missing(w)) 
    w <- numeric(nx)+1; 
    h <- h/4 
    ax <- outer(gx, x, "-")/h[1] # distance of each point to each grid point in x-direction 
    ay <- outer(gy, y, "-")/h[2] # distance of each point to each grid point in y-direction 
    z <- (matrix(rep(w,n), nrow=n, ncol=nx, byrow=TRUE)*matrix(dnorm(ax), n, nx)) %*% t(matrix(dnorm(ay), n, nx))/(sum(w) * h[1] * h[2]) # z is the density 
    return(list(x = gx, y = gy, z = z)) 
} 

我想提出的功能geom_density2d和stat_density2d通話kd2d.weighted代替kde2d,並通過讓他們接受加權數據。

我從來沒有改變現有的R包中的任何功能,所以我的問題是這樣做的最簡單的方法是什麼?

回答

4

你實際上可以將你自己的密度數據傳遞給geom_contour,這可能是最簡單的。讓我們從樣本數據集開始,爲間歇數據添加權重。

library("MASS") 
data(geyser, "MASS") 
geyserw <- transform(geyser, 
    weigh = sample(1:5, nrow(geyser), replace=T) 
) 

現在我們使用加權函數來計算密度,並把它變成一個data.frame

dens <- kde2d.weighted(geyserw$duration, geyserw$waiting, geyserw$weight) 
dfdens <- data.frame(expand.grid(x=dens$x, y=dens$y), z=as.vector(dens$z)) 

現在我們繪製數據

ggplot(geyserw, aes(x = duration, y = waiting)) + 
    geom_point() + xlim(0.5, 6) + ylim(40, 110) 
    geom_contour(aes(x=x, y=y, z=z), data= dfdens) 

而且應該這樣做

resulting weighted density

+0

謝謝!正是我在找的東西。雖然,是否有可能使用這個geom_contour做填充等值線圖? – AntonvSchantz

+0

stat_contour似乎做的伎倆:) – AntonvSchantz

相關問題