根据 R 中的 2D 密度图计算值的概率答案

【问题标题】：Calculate probability of value based on 2D density plot in R根据 R 中的 2D 密度图计算值的概率
【发布时间】：2020-04-17 20:24:35
【问题描述】：

我正在寻找一个函数来计算 B 和 R 的某个组合的可能性。当前的数据说明如下所示：

ggplot(df, aes(R,B)) +
geom_bin2d(binwidth = c(1,1))

有没有办法根据这两个正偏态的离散相关变量计算每个组合的概率（例如 R = 23，B = 30）？

是否可以使用 stat_density_2d 来解决或者是否有更好的方法？

谢谢。

【问题讨论】：

我在赞成票和反对票之间摇摆不定，因为 1) 并不完全清楚您想要什么 - 只是计算值或实际绘图？ 2）您没有提供样本数据。最后我决定投票，因为这让@JonSpring 给出了一个有趣的答案。

标签： r ggplot2 probability-density

【解决方案1】：

stat_density_2d 在后台使用MASS::kde2d。我想有更巧妙的方法可以做到这一点，但我们可以将数据输入该函数并将其转换为整洁的数据，以获得该类型估计的平滑版本。

首先，像你这样的一些数据：

library(tidyverse)
set.seed(42)
df <- tibble(
  R = rlnorm(1E4, 0, 0.2) * 100,
  B = R * rnorm(1E4, 1, 0.2)
)

ggplot(df, aes(R,B)) +
  geom_bin2d(binwidth = c(1,1))

这里运行密度并转换为与数据具有相同坐标的小标题。（有更好的方法吗？）

n = 201 # arbitrary grid size, chosen to be 1 more than the range below 
        #   so the breaks are at integers
smooth <- MASS::kde2d(df$R, df$B, lims = c(0, 200, 0, 200),
                      # h = c(20,20),  # could tweak bandwidth here 
                      n = n) 
df_smoothed <- smooth$z %>% 
  as_tibble() %>%
  pivot_longer(cols = everything(), names_to = "col", values_to = "val") %>% 
  mutate(R = rep(smooth$x, each = n), # EDIT: fixed, these were swapped
         B = rep(smooth$y, n))

df_smoothed 现在包含 R 和 B 维度中从 0:200 开始的所有坐标，val 列中的每个组合的概率。这些加起来几乎是 1（在这种情况下为 99.6%）。我认为剩下的 smidgen 是坐标超出指定范围的概率。

sum(df_smoothed$val)
#[1] 0.9960702

任何特定组合的机会只是那个点的密度值。所以 R = 70 和 B = 100 的几率是 0.013%。

df_smoothed %>%
  filter(R == 70, B == 100)
## A tibble: 1 x 4
#  col        val     R     B
#  <chr>    <dbl> <int> <int>
#1 V101   0.0000345    70   100

R 在 50-100 之间和 B 在 50-100 之间的机会是 36.9%：

df_smoothed %>%
  filter(R %>% between(50, 100),
         B %>% between(50, 100)) %>%
  summarize(total_val = sum(val))
## A tibble: 1 x 1
#total_val
#<dbl>
#  1     0.369

以下是平滑数据和原始数据的外观：

ggplot() +
  geom_tile(data = df_smoothed, aes(R, B, alpha = val), fill = "red") +
  geom_point(data = df %>% sample_n(500), aes(R, B), size = 0.2, alpha = 1/5)

【讨论】：

这很有趣。我认为“错误”必须以某种方式存在于矩阵到数据框的转换中？解释 OP 的一种方法是，他们实际上不想绘制密度估计，而是要检索每个坐标的预测值 - 那么您能否详细说明如何将预测的密度矩阵转换为坐标？
感谢@JonSpring 的回复，非常感谢！我做了一些检查，当我扩大sum(df_smoothed$val) 的范围时，它出现了 1。我还发现由于某种原因斜率也不同。感谢您提供从那里计算的方法，我将尝试找出梯度为何如此，但感谢您让我走上正轨！
好的，如果你翻转kde2d中的df$R和df$B，它看起来就像我们想要的形状！
感谢您的提示 - 看起来我搞砸了并交换了 mutate(R = rep(smooth$x, each = n), B = rep(smooth$y, n)) 部分，并且混淆了 rep 替代品。这修复了后面的绘图。

【解决方案2】：

如果只是关于绘图，可以简单地关闭轮廓并使用geom = raster，就像建议的in the ggplot2 reference 一样。

感谢@JonSpring 提供示例数据！

library(tidyverse)

df <- tibble(
  R = rlnorm(1E4, 0, 0.2) * 100,
  B = R * rnorm(1E4, 1, 0.2)
)

ggplot(df, aes(R,B)) +
  stat_density2d(geom = 'raster', aes(fill = stat(density)), contour = FALSE)

^{由reprex package (v0.3.0) 于 2019 年 12 月 28 日创建}

【讨论】：