编写循环计算 Jensen-Shannon 散度的巧妙方法

Question

我尝试使用应用或聚合函数来计算一些分布之间的 Jensen-Shannon 散度 (JS.dist)。

我正在模拟四种不同模型下的一些数据，并为每个数据计算一系列统计数据。

所以想象我有以下data.frame：

dataframe1:
Model Factor1 Factor2 stats1 stats2
M1    0.0001  0.2     -1.0   0.9
M1    0.0001  0.2     -1.3   0.5
M1    0.0002  0.3     -1.9   0.2
M2    0.0001  0.2     -2.0   0.2
M2    0.0001  0.2     -2.0   0.2
M2    0.0002  0.3     -2.1   0.4
M3    0.0001  0.2      9.9   0.4
M3    0.0001  0.2      8.3   0.4
M3    0.0002  0.3      8.0   0.4
M4    0.0001  0.2      3.0   0.1
M4    0.0001  0.2      3.5   0.3
M4    0.0002  0.3      3.2   0.3

计算JS.dist的函数如下：

在日志中将 Inf 或 -Inf 更改为零的功能。它以数字的日志为参数

 test.logs <- function(log.num){

  log.num[log.num == -Inf | log.num == Inf] <- 0
  return (log.num)

}

# 计算 kl.dist（Kullback–Leibler 散度）的函数。它将两个分布（x.p 和 y.p）的概率向量（见下文）作为参数

kl.dist <- function(x.p, y.p) {
  # x.p, y.p: probability vectors for x and y distributions

  log.x <- test.logs(log(x.p))
  log.y <- test.logs(log(y.p))

  sum(x.p * (log.x - log.y))

}

# 计算 js.dist 的函数。它将 x、y 和 M 的概率向量作为参数。M 是中间分布

js.dist <- function(x.p, y.p, M.p){
  0.5 * kl.dist(x.p, M.p) + .5 * kl.dist(y.p, M.p)
}

要使用上述函数，我必须计算分布的密度曲线（按模型和因子计算统计 1 和统计 2）。

要计算这个，我必须设置一个最小值和最大值，以便为每个统计数据计算密度曲线。

例如：

x.d <- density(x, n=512, from=min, to=max)
y.d <- density(y, n=512, from=min, to=max)
M.d <- (x.d$y + y.d$y)/2

# width of the histogram
w <- x.d$x[2] - x.d$x[1]

# probability of x value in n-th bin
x.p <- x.d$y * w # (hist hight) * (bin width)
y.p <- y.d$y * w
M.p <- M.d * w

我尝试编写一个代码，其中有两个 for 循环（针对每个因素），并按模型对数据进行子集化，并计算每个统计数据的最小值和最大值。最后我计算密度和概率，只有在我可以计算 JS.dist 之后。

以R代码为例：

density_js.dist <- function(data.df){
# gets the unique values for mutation rate
factor1 <- unique(data.df$Factor1)
# gets the unique values for rate of new copies
factor2 <- unique(data.df$factor2)

# calculates the minimum and maximum value for each of the statistics
# showing only for stats1
stats1.min <- min(data.df$stats1)
stats1.max <- max(data.df$stats1)



# for loop to calculate the densities and probabilities and JS distance for each combination of factor1 and factor2

for (f1 in factor1){
  for (f2 in factor2){

  new.df <- subset(data.df, factor1 == f1 & factor2 == f2)

  # subsetting data. One data frame for each of the four models
  MM.df <- subset(new.df, Model == "M1")
  TM.df <- subset(new.df, Model == "M2")

  MI.df <- subset(new.df, Model == "M3")
  TI.df <- subset(new.df, Model == "M4")

  # densitiy and probability for each stats

  #1.stats1
  # calculating densities for model M1 and M2
  MM1.d <- density(MM.df$stats1, n=512, from=stats1.min, to=stats1.max)
  TM1.d <- density(TM.df$stats1, n=512, from=stats1.min, to=stats1.max)

  # Density for the middle distribution between models M1 and M2 
  Middle12.d <- (MM1.d$y + TM1.d$y)/2

  # width for models
  w12 <- MM1.d$x[2] - MM1.d$x[1]

  # calculating probabilities for each models
  MM1.p <- MM1.d$y * w12 # (hist hight) * (bin width)
  TM1.p <- TM1.d$y * w12
  Middle12.p <- Middle12.d * w12 

  # calculating densities for models M3 and M4
  MI1.d <- density(MI.df$stats1, n=512, from=stats1.min, to=stats1.max)
  TI1.d <- density(TI.df$stats1, n=512, from=stats1.min, to=stats1.max)
  Middle34.d <- (MI1.d$y + TI1.d$y)/2

  w34 <- MI1.d$x[2] - MI1.d$x[1]

  # calculating probabilities for M3 and M4 models
  MI1.p <- MM1.d$y * w34 
  TI1.p <- TM1.d$y * w34
  Middle34.p <- Middle34.d * w34 


 js.dist(MM1.p, TM1.p, Middle12.p)
 js.dist(MI1.p, TI1.p, Middle34.p)
  }
 }
}

我的问题是：

我尝试使用应用或聚合，但是我不知道如何将每个统计数据的最小值和最大值作为参数传递，以便能够创建密度曲线？ 请注意，此最小值和最大值是针对因素和模型的所有组合而不是针对每个子集计算的。例如，为了便于比较，我无法通过因子和模型计算子集的最小值和最大值。

我的数据实际上要复杂得多。我有 10 个不同的统计数据，我想按因子计算两个分布之间的 JS.dist。我的两个分布是 M1 和 M2，以及 M3 和 M4。 上面的代码有效，但它需要我写更多的 700 行，我真的认为这不是很聪明。

如果有人能帮我解决这个问题，我将不胜感激。

Answer 1

这是一种使用列表同时计算所有 10 个系列的 hack 方式。由于代码的长度和冗长，如果您想要一个功能解决方案，则需要完全重写。只能测试前两个系列的代码（甚至不能完全测试，因为不止一个 factor1:factor2 组合只有 1 个观察值，因此无法进行密度计算）。还删除了该功能，因为它对您完全没有任何作用。

library(dplyr)
L = list()

  # gets the unique values for mutation rate
  factor1 <- unique(data.df$Factor1)
  # gets the unique values for rate of new copies
  factor2 <- unique(data.df$Factor2)

  # calculates the minimum and maximum value for each of the statistics
  # Store all 10 min and max in a vector
  vector.min <- lapply(data.df %>% select(stats1:stats10), min)
  vector.max <- lapply(data.df %>% select(stats1:stats10), max)

  # for loop to calculate the densities and probabilities and JS distance for each combination of factor1 and factor2

  for (f1 in factor1){
    for (f2 in factor2){
      new.df <- subset(data.df, factor1 == f1 & factor2 == f2)    
      # subsetting data. One data frame for each of the four models
      MM.df <- subset(new.df, Model == "M1")
      TM.df <- subset(new.df, Model == "M2")
      MI.df <- subset(new.df, Model == "M3")
      TI.df <- subset(new.df, Model == "M4")

      # densitiy and probability for each stats

      # calculating densities for model M1 and M2
      MM.d = lapply(1:10, function(i) density(MM.df %>% select(i+3) %>% unlist, n = 512, from = vector.min[[i]], to = vector.min[[i]]))
      TM.d = lapply(1:10, function(i) density(TM.df %>% select(i+3) %>% unlist, n = 512, from = vector.min[[i]], to = vector.min[[i]]))

      # Density for the middle distribution between models M1 and M2 
      Middle12.d <- mapply(function(d1, d2) (d1$y+d2$y)/2, MM.d, TM.d, SIMPLIFY = F)

      # width for models
      w12 = lapply(MM.d, function(y) {y$x[2] - y$x[1]})

      # calculating probabilities for each models
      MM1.p = mapply(function(arg1, arg2) {arg1$y * arg2}, MM.d, w12)  # (hist hight) * (bin width)
      TM1.p = mapply(function(arg1, arg2) {arg1$y * arg2}, TM.d, w12)
      Middle12.p = mapply("*", Middle12.d, w12)

      # calculating densities for models M3 and M4
      MI.d = lapply(1:10, function(i) density(MI.df %>% select(i+3) %>% unlist, n = 512, from = vector.min[[i]], to = vector.min[[2]]))
      TI.d = lapply(1:10, function(i) density(TI.df %>% select(i+3) %>% unlist, n = 512, from = vector.min[[i]], to = vector.min[[2]]))
      Middle34.d <- mapply(function(d1, d2) (d1$y+d2$y)/2, MI.d, TI.d)

      w34 = lapply(MI.d, function(y) {y$x[2] - y$x[1]})      

      # calculating probabilities for M3 and M4 models
      MI1.p = mapply(function(arg1, arg2) {arg1$y * arg2}, MI.d, w34)  # (hist hight) * (bin width)
      TI1.p = mapply(function(arg1, arg2) {arg1$y * arg2}, TI.d, w34)
      Middle34.p = mapply("*", Middle34.d, w34)

      L = c(L, list(mapply(js.dist, MM1.p, TM1.p, Middle12.p), mapply(js.dist, MI1.p, TI1.p, Middle34.p)))
    }
  }