【问题标题】:How to measure area between 2 distribution curves in R / ggplot2如何测量R / ggplot2中2条分布曲线之间的面积
【发布时间】:2014-09-04 17:34:21
【问题描述】:

具体的例子是假设x是0到10之间的一些连续变量,红线是“商品”的分布,蓝线是“坏”的分布,我想看看合并这个是否有价值变量来检查“善良”,但我想首先量化蓝色>红色区域中的东西数量

因为这是一个分布图,比例看起来是一样的,但实际上我的样本中有 98 倍的好样本,这使事情变得复杂,因为它实际上不仅仅是测量曲线下的面积,而是测量坏样本它的分布是沿着大于红色的线。

我一直在努力学习 R,但我什至不知道如何接近这个,任何帮助表示赞赏。

编辑 样本数据: http://pastebin.com/7L3Xc2KU

图表是用

创建的
graph <- qplot(sample_x, bad_is_1, data=sample_data, geom="density", color=bid_is_1)

【问题讨论】:

  • 这里有一些sample data 会很有帮助。
  • 感谢@MrFlick,添加了一个链接来显示我正在使用的数据类型
  • base 图形中,polygon 是绘制这种差异的方式。要计算它,找到两条曲线的交点,或者只找到“红色”
  • @CarlWitthoft 我最终计算了曲线之间的梯形面积。当你有一个点列表时,有没有更好的方法来用 R 进行数值积分?

标签: r ggplot2


【解决方案1】:

这是一种对两个密度图之间的区域进行着色并计算该区域大小的方法。

# Create some fake data
set.seed(10)
dat = data.frame(x=c(rnorm(1000, 0, 5), rnorm(2000, 0, 1)), 
                 group=c(rep("Bad", 1000), rep("Good", 2000)))

# Plot densities
# Use y=..count.. to get counts on the vertical axis
p1 = ggplot(dat) +
       geom_density(aes(x=x, y=..count.., colour=group), lwd=1)

一些额外的计算来遮蔽两个密度图之间的区域 (改编自this SO question):

pp1 = ggplot_build(p1)

# Create a new data frame with densities for the two groups ("Bad" and "Good")
dat2 = data.frame(x = pp1$data[[1]]$x[pp1$data[[1]]$group==1],
                 ymin=pp1$data[[1]]$y[pp1$data[[1]]$group==1],
                 ymax=pp1$data[[1]]$y[pp1$data[[1]]$group==2])

# We want ymax and ymin to differ only when the density of "Good" 
# is greater than the density of "Bad"
dat2$ymax[dat2$ymax < dat2$ymin] = dat2$ymin[dat2$ymax < dat2$ymin]

# Shade the area between "Good" and "Bad"
p1a = p1 +  
    geom_ribbon(data=dat2, aes(x=x, ymin=ymin, ymax=ymax), fill='yellow', alpha=0.5)

这是两个图:

要获取GoodBad 特定范围内的面积(值的数量),请对每个组使用density 函数(或者您可以继续使用从ggplot 提取的数据,如上,但这样您可以更直接地控制密度分布的生成方式):

## Calculate densities for Bad and Good. 
# Use same number of points and same x-range for each group, so that the density 
# values will line up. Use a higher value for n to get a finer x-grid for the density
# values. Use a power of 2 for n, because the density function rounds up to the nearest 
# power of 2 anyway.
bad = density(dat$x[dat$group=="Bad"], 
             n=1024, from=min(dat$x), to=max(dat$x))
good = density(dat$x[dat$group=="Good"], 
             n=1024, from=min(dat$x), to=max(dat$x))

## Normalize so that densities sum to number of rows in each group

# Number of rows in each group
counts = tapply(dat$x, dat$group, length)

bad$y = counts[1]/sum(bad$y) * bad$y
good$y = counts[2]/sum(good$y) * good$y

## Results

# Number of "Good" in region where "Good" exceeds "Bad"
sum(good$y[good$y > bad$y])
[1] 1931.495  # Out of 2000 total in the data frame

# Number of "Bad" in region where "Good" exceeds "Bad"
sum(bad$y[good$y > bad$y])
[1] 317.7315  # Out of 1000 total in the data frame

【讨论】:

    【解决方案2】:

    我能想到的唯一方法是使用简单的梯形计算曲线之间的面积。首先我们手动计算密度

    d0 <- density(sample$sample_x[sample$bad_is_1==0])
    d1 <- density(sample$sample_x[sample$bad_is_1==1])
    

    现在我们创建将在观察到的密度点之间进行插值的函数

    f0 <- approxfun(d0$x, d0$y)
    f1 <- approxfun(d1$x, d1$y)
    

    接下来我们找到密度重叠的x范围

    ovrng <- c(max(min(d0$x), min(d1$x)), min(max(d0$x), max(d1$x)))
    

    并将其分成 500 个部分

    i <- seq(min(ovrng), max(ovrng), length.out=500)
    

    现在我们计算密度曲线之间的距离

    h <- f0(i)-f1(i)
    

    并使用梯形面积公式,我们将 d1>d0 区域的面积相加

    area<-sum( (h[-1]+h[-length(h)]) /2 *diff(i) *(h[-1]>=0+0))
    # [1] 0.1957627
    

    我们可以绘制区域使用

    plot(d0, main="d0=black, d1=green")
    lines(d1, col="green")
    jj<-which(h>0 & seq_along(h) %% 5==0); j<-i[jj]; 
    segments(j, f1(j), j, f1(j)+h[jj])
    

    【讨论】:

    • 好吧,你可以总是使用抛物线拟合,或切比雪夫多项式等:-)。哎呀-你问这个:integrate 使用自适应正交;辛普森法则有多种实现方式。可能其他我还没有找到。
    • @CarlWitthoft 当然,有很多方法可以计算这样的面积。我只是不确定是否有任何我刚刚丢失的内置内容。
    猜你喜欢
    • 2022-11-24
    • 2019-08-15
    • 2013-02-27
    • 1970-01-01
    • 2018-02-11
    • 2017-10-12
    相关资源
    最近更新 更多