2011-03-03 128 views
1

我需要寫兩個函數toBand()replaceBand()r功能來創建帶矩陣

test3 <- toBand(test1,3) 
test4 <- replaceBand(test1, toBand(test2,3)) 

得到下面的輸出。

第一個返回對角線和對角線元素的矩陣,使用X的下三角形。第二個從所提供的頻帶矩陣中的數據中替換X中的對應元素。該函數然後返回修改的X.

是否有任何可用於執行此操作的軟件包?有關如何做到這一點的任何建議?

謝謝

> test1 
    [,1] [,2] [,3] [,4] [,5] [,6] 
[1,] "a11" "a21" "a31" "a41" "a51" "a61" 
[2,] "a21" "a22" "a32" "a42" "a52" "a62" 
[3,] "a31" "a32" "a33" "a43" "a53" "a63" 
[4,] "a41" "a42" "a43" "a44" "a54" "a64" 
[5,] "a51" "a52" "a53" "a54" "a55" "a65" 
[6,] "a61" "a62" "a63" "a64" "a65" "a66" 
> test2 
    [,1] [,2] [,3] [,4] [,5] [,6] 
[1,] "*a11*" "*a21*" "*a31*" "*a41*" "*a51*" "*a61*" 
[2,] "*a21*" "*a22*" "*a32*" "*a42*" "*a52*" "*a62*" 
[3,] "*a31*" "*a32*" "*a33*" "*a43*" "*a53*" "*a63*" 
[4,] "*a41*" "*a42*" "*a43*" "*a44*" "*a54*" "*a64*" 
[5,] "*a51*" "*a52*" "*a53*" "*a54*" "*a55*" "*a65*" 
[6,] "*a61*" "*a62*" "*a63*" "*a64*" "*a65*" "*a66*" 
> test3 
    [,1] [,2] [,3] [,4] [,5] [,6] 
[1,] "a11" "a22" "a33" "a44" "a55" "a66" 
[2,] "a21" "a32" "a43" "a54" "a65" NA 
[3,] "a31" "a42" "a53" "a64" NA NA 
[4,] "a41" "a52" "a63" NA NA NA 
> test4 
    [,1] [,2] [,3] [,4] [,5] [,6] 
[1,] "*a11*" "*a21*" "*a31*" "*a41*" "a51" "a61" 
[2,] "*a21*" "*a22*" "*a32*" "*a42*" "*a52*" "a62" 
[3,] "*a31*" "*a32*" "*a33*" "*a43*" "*a53*" "*a63*" 
[4,] "*a41*" "*a42*" "*a43*" "*a44*" "*a54*" "*a64*" 
[5,] "a51" "*a52*" "*a53*" "*a54*" "*a55*" "*a65*" 
[6,] "a61" "a62" "*a63*" "*a64*" "*a65*" "*a66*" 

回答

2

這很簡單,如果你不需要使用來自toBand輸出:

replaceBand <- function(a, b, k) { 
    swap <- abs(row(a) - col(a)) <= k 
    a[swap] <- b[swap] 
    a 
} 

使得矩陣證明:

test1 <- matrix(ncol=6, nrow=6) 
test1 <- matrix(paste("a", row(test1), col(test1), sep=""), nrow=6) 
test1b <- matrix(paste("a", col(test1), row(test1), sep=""), nrow=6) 
test1[upper.tri(test1)] <- test1b[upper.tri(test1b)] 
test2 <- matrix(paste("*", test1, "*", sep=""), nrow=6) 

的輸出完全符合要求:

> replaceBand(test1, test2, 3) 

    [,1] [,2] [,3] [,4] [,5] [,6] 
[1,] "*a11*" "*a21*" "*a31*" "*a41*" "a51" "a61" 
[2,] "*a21*" "*a22*" "*a32*" "*a42*" "*a52*" "a62" 
[3,] "*a31*" "*a32*" "*a33*" "*a43*" "*a53*" "*a63*" 
[4,] "*a41*" "*a42*" "*a43*" "*a44*" "*a54*" "*a64*" 
[5,] "a51" "*a52*" "*a53*" "*a54*" "*a55*" "*a65*" 
[6,] "a61" "a62" "*a63*" "*a64*" "*a65*" "*a66*" 

以下是版本toBandreplaceBand,其工作方式如上所述。我認爲用算術來確定如何填充矩陣會更清晰,但這是一種不需要非常努力就可以完成的方法。也許別人會這樣回答。

toBand <- function(x,k) { 
    n <- nrow(x) 
    out <- matrix(nrow=n, ncol=n) 
    out[row(out) + col(out) - 1 <= n] <- x[lower.tri(x, diag=TRUE)] 
    out[1:(k+1),] 
} 

replaceBand <- function(a, b) { 
    b[row(b)+col(b)-1 <= ncol(b)] 
    swap <- abs(row(a) - col(a)) <= nrow(b) - 1 
    a[swap & lower.tri(a, diag=TRUE)] <- b[row(b)+col(b)-1 <= ncol(b)] 
    a[upper.tri(a)] <- t(a)[upper.tri(a)] 
    a 
} 
+0

何亞倫,謝謝您的回覆。我還有一個天真的問題。如何爲需要矩陣和k值的S4「Band」類編寫initialize()方法?它應該被定義爲具有形式參數(x,k)的函數,並且應該將X的下三角形元素放入對象中。 – user635034 2011-03-04 03:40:36

+1

@ user635034:我只使用S3方法。如果通常的S4文檔不清楚,可以試試Hadley的wiki:https://github.com/hadley/devtools/wiki/S4 – Aaron 2011-03-04 05:12:16