2017-01-11 9 views
1

residuesという2つのベクトルと、scoresという2つのベクトルがあることを考慮すると、各残差に1つの正の数が31点あります。説明するために、以下のように2つのベクターを得た。R - 複数のシナリオで曲線下の面積を最大にする

例示のためにランダム配列を検討している。私がやりたいことは以下の通りです

enter image description here

:私はスキップし、15の残基(以後、15量体と呼ばれる)の特定の組み合わせを検討します、私はプロットするとresidues X scoresは、私は次のグラフィックを持つことになります 1つの残基(すなわち、1:15、2:16、3:17までに17:31まで)、これらの17通りの組み合わせすべてについて曲線下面積(AUC)を計算したいと考えています。私の最終的な目標は、最も高いAUCを有する15量体を選択することです。

this questionに示すように、動物園パッケージのロールメン機能を使用してAUCを計算することができます。しかし、この例では17通りの組み合わせがありますが、私はプロセスを自動化するスクリプトを見つけようとしています。 ありがとうございます。ここで

答えて

2
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 
:あなたのデータに応じて、あなたは私たちが結果を次のようにループ

for (i in 1:length(N)) { 
    output[i, ] <- dat[N[i]:(N[i] + n -1)] 
} 

でそれらのすべてを返却することを決定した最大を得る複数のシーケンスを持っているかもしれません

またはカスタム機能付き

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) 
+0

非常に便利で、私が簡単に望んでいたものとまったく同じでした。ありがとう – BCArg

0

は、私はあなたが述べたように動物園パッケージから

rollmean(dat, n) 

裏返しそれも実行

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) 

ことができますを使用していた以下の機能は、私たちに私たち

のローリング平均を与えています
max(rollmean(dat, n)) 

転がり軸受の最大値を見つける

max(rollmean(dat, n)) == rollmean(dat, n) 

最大

N <- which(max(rollmean(dat, n)) == rollmean(dat, n)) 

等しい圧延手段のうちTRUE/FALSEベクトルが最大値のインデックスを返します。

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 
関連する問題