R - 最大化多个场景的曲线下面积
R - maximize area under curve for multiple scenarios
考虑到我有两个向量,一个叫residues
,另一个叫scores
,它们有31个分数,每个残基一个,都是正数。为了说明,得到两个向量如下图:
residues <- 1:31
scores <- runif(n = 31, min = 0.35, max = 3.54)
我正在考虑一个随机序列来举例说明。
如果我绘制 residues
x scores
我将得到以下图形:
我想做的是:我将考虑 15 个残基的特定组合(以下称为 15mer),跳过一个残基(即 1:15、2:16、3:17 一直到 17:31),我想计算所有这 17 种组合的曲线下面积 (AUC)。我的最终目标是 select 具有最高 AUC 的 15mer。
可以使用 zoo 包中的 rollmean 函数计算 AUC,如 this question 所示。然而,在这个例子中,我有 17 种可能的组合,我试图找到一个脚本来自动化这个过程。
提前致谢。
下面是我用过的函数
scores <- runif(n = 31, min = 0.35, max = 3.54)
fun <- function(dat, n) {
require(zoo)
N <- which(max(rollmean(dat, n)) == rollmean(dat, n))
output <- matrix(0, length(N), n)
for (i in 1:length(N)) {
output[i, ] <- dat[N[i]:(N[i] + n - 1)]
}
output
}
fun(scores, 15)
让运行彻底翻过来
rollmean(dat, n)
您提到的 zoo 包为我们提供了滚动平均值,我们
max(rollmean(dat, n))
找到滚动平均值的最大值
max(rollmean(dat, n)) == rollmean(dat, n)
returns 一个TRUE/FALSE向量,其滚动均值等于最大值
N <- which(max(rollmean(dat, n)) == rollmean(dat, n))
returns 最大值的索引。根据您的数据,您可能有多个序列获得最大值,我们决定 return 所有这些都使用以下循环
for (i in 1:length(N)) {
output[i, ] <- dat[N[i]:(N[i] + n -1)]
}
结果:
set.seed(12345)
scores <- runif(n = 31, min = 0.35, max = 3.54)
fun(scores, 15)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.588179 1.633928 0.9208938 3.385791 1.797393 1.39234
[,7] [,8] [,9] [,10] [,11] [,12]
[1,] 3.429675 2.606867 2.406091 1.593553 2.578354 2.085545
[,13] [,14] [,15]
[1,] 1.07243 1.895739 2.879693
fun(rpois(1000, 1), 10)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 4 2 1 1 3 3 2 2
[2,] 1 4 2 1 1 3 3 2 2 1
[3,] 4 2 1 1 3 3 2 2 1 1
library(zoo)
set.seed(555)
residues <- 1:31
scores <- runif(n = 31, min = 0.35, max = 3.54)
which.max(sapply(1:17, function(x){sum(diff(residues[x:(x+14)])*rollmean(scores[x:(x+14)],2))}))
# result 7 i.e. 7:21
或
sapply(1:17, function(x){sum(diff(residues[x:(x+14)])*rollmean(scores[x:(x+14)],2))}) # gives you the AUCs
# result [1] 28.52530 29.10203 28.52847 27.65325 27.19925 28.77782 29.29373 28.13133 28.23705 27.68724 25.75294 25.27226 25.44963 25.81201 25.49907 23.48632
#[17] 22.45763
或使用自定义函数
f_AUC <- function(x, y, lngth){
sapply(1:(length(x)-lngth+1), function(z) sum(diff(x[z:(z+lngth-1)])*rollmean(y[z:(z+lngth-1)],2)))
}
f_AUC(x=residues, y=scores, lngth=15)
考虑到我有两个向量,一个叫residues
,另一个叫scores
,它们有31个分数,每个残基一个,都是正数。为了说明,得到两个向量如下图:
residues <- 1:31
scores <- runif(n = 31, min = 0.35, max = 3.54)
我正在考虑一个随机序列来举例说明。
如果我绘制 residues
x scores
我将得到以下图形:
我想做的是:我将考虑 15 个残基的特定组合(以下称为 15mer),跳过一个残基(即 1:15、2:16、3:17 一直到 17:31),我想计算所有这 17 种组合的曲线下面积 (AUC)。我的最终目标是 select 具有最高 AUC 的 15mer。
可以使用 zoo 包中的 rollmean 函数计算 AUC,如 this question 所示。然而,在这个例子中,我有 17 种可能的组合,我试图找到一个脚本来自动化这个过程。 提前致谢。
下面是我用过的函数
scores <- runif(n = 31, min = 0.35, max = 3.54)
fun <- function(dat, n) {
require(zoo)
N <- which(max(rollmean(dat, n)) == rollmean(dat, n))
output <- matrix(0, length(N), n)
for (i in 1:length(N)) {
output[i, ] <- dat[N[i]:(N[i] + n - 1)]
}
output
}
fun(scores, 15)
让运行彻底翻过来
rollmean(dat, n)
您提到的 zoo 包为我们提供了滚动平均值,我们
max(rollmean(dat, n))
找到滚动平均值的最大值
max(rollmean(dat, n)) == rollmean(dat, n)
returns 一个TRUE/FALSE向量,其滚动均值等于最大值
N <- which(max(rollmean(dat, n)) == rollmean(dat, n))
returns 最大值的索引。根据您的数据,您可能有多个序列获得最大值,我们决定 return 所有这些都使用以下循环
for (i in 1:length(N)) {
output[i, ] <- dat[N[i]:(N[i] + n -1)]
}
结果:
set.seed(12345)
scores <- runif(n = 31, min = 0.35, max = 3.54)
fun(scores, 15)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.588179 1.633928 0.9208938 3.385791 1.797393 1.39234
[,7] [,8] [,9] [,10] [,11] [,12]
[1,] 3.429675 2.606867 2.406091 1.593553 2.578354 2.085545
[,13] [,14] [,15]
[1,] 1.07243 1.895739 2.879693
fun(rpois(1000, 1), 10)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 4 2 1 1 3 3 2 2
[2,] 1 4 2 1 1 3 3 2 2 1
[3,] 4 2 1 1 3 3 2 2 1 1
library(zoo)
set.seed(555)
residues <- 1:31
scores <- runif(n = 31, min = 0.35, max = 3.54)
which.max(sapply(1:17, function(x){sum(diff(residues[x:(x+14)])*rollmean(scores[x:(x+14)],2))}))
# result 7 i.e. 7:21
或
sapply(1:17, function(x){sum(diff(residues[x:(x+14)])*rollmean(scores[x:(x+14)],2))}) # gives you the AUCs
# result [1] 28.52530 29.10203 28.52847 27.65325 27.19925 28.77782 29.29373 28.13133 28.23705 27.68724 25.75294 25.27226 25.44963 25.81201 25.49907 23.48632
#[17] 22.45763
或使用自定义函数
f_AUC <- function(x, y, lngth){
sapply(1:(length(x)-lngth+1), function(z) sum(diff(x[z:(z+lngth-1)])*rollmean(y[z:(z+lngth-1)],2)))
}
f_AUC(x=residues, y=scores, lngth=15)