effizientere Art und Weise ein spezifische Bandmatrix

Question

Ich verwende den folgenden Code zu erstellen:

var <- c(rep(4,4),rep(9,5)) 
cov <- diag(var) 
n <- length(var) 
rho <- 0.2 
for(i in 1:(n-1)){ 
    for(j in (i+1):n){ 
    if (j <= i + 4) cov[i,j] <- rho/(j-i)* sqrt(var[i] * var[j]) 
    } 
} 

So erstellen Sie diesen gewünschten Matrixausgang:

 [,1] [,2] [,3]  [,4] [,5] [,6] [,7] [,8] [,9] 
[1,] 4 0.8 0.4 0.2666667 0.3 0.0 0.0 0.0 0.00 
[2,] 0 4.0 0.8 0.4000000 0.4 0.3 0.0 0.0 0.00 
[3,] 0 0.0 4.0 0.8000000 0.6 0.4 0.3 0.0 0.00 
[4,] 0 0.0 0.0 4.0000000 1.2 0.6 0.4 0.3 0.00 
[5,] 0 0.0 0.0 0.0000000 9.0 1.8 0.9 0.6 0.45 
[6,] 0 0.0 0.0 0.0000000 0.0 9.0 1.8 0.9 0.60 
[7,] 0 0.0 0.0 0.0000000 0.0 0.0 9.0 1.8 0.90 
[8,] 0 0.0 0.0 0.0000000 0.0 0.0 0.0 9.0 1.80 
[9,] 0 0.0 0.0 0.0000000 0.0 0.0 0.0 0.0 9.00 

Allerdings ist dieser Code zu langsam, um den Fall zu berechnen von großem n. Haben Sie effiziente Lösungen?

Answer 1

Ich gehe davon aus Ihrem n recht groß ist, so dass Sie eine spärliche Bandmatrix mit einer Bandbreite von 5

Zunächst wird eine Hilfsfunktion wie diff, die eine beliebige Funktion, anstatt nur Subtraktion erlaubt (-).

fdiff <- function(x,lag,f) { 
    i1 <- -seq_len(lag) 
    f(x[i1],x[-length(x):-(length(x)-lag+1L)]) 
} 

In diesem Fall ist die Funktion, die wir wollen, ist

gm <- function(x,y) sqrt(x*y) 

Damit der erste superdiagonal ist gegeben durch

x <- c(rep(4,4),rep(9,5)) 
0.2*fdiff(x,1,gm)/1 
# [1] 0.8 0.8 0.8 1.2 1.8 1.8 1.8 1.8 

Um eine Ersatzbandmatrix zu füllen, verwenden wir 'bandSparse' von der Matrix Bibliothek

library(Matrix) 
x <- c(rep(4,4),rep(9,5)) 
bandSparse(n,k=0:4,diagonals= 
    c(list(x),lapply(1:4,function(lag) 0.2*fdiff(x,lag,gm)/lag))) 

Output

:

9 x 9 sparse Matrix of class "dgCMatrix" 

[1,] 4 0.8 0.4 0.2666667 0.3 . . . . 
[2,] . 4.0 0.8 0.4000000 0.4 0.3 . . . 
[3,] . . 4.0 0.8000000 0.6 0.4 0.3 . . 
[4,] . . . 4.0000000 1.2 0.6 0.4 0.3 . 
[5,] . . . .   9.0 1.8 0.9 0.6 0.45 
[6,] . . . .   . 9.0 1.8 0.9 0.60 
[7,] . . . .   . . 9.0 1.8 0.90 
[8,] . . . .   . . . 9.0 1.80 
[9,] . . . .   . . . . 9.00 

Answer 2

Was auch immer Sie verwenden, müssen Sie die obere/untere Dreieck der Matrix füllen. In diesem Fall ist das Muster sehr einfach, d. H. mat1 <- rho/(col(mat)-row(mat));diag(mat1)=mat;#Cov mat. Um die Korrelationsmatrix zu berechnen, beachten Sie einfach R Speicherelemente in der Spalte und tun Sie v=sqrt(diag(mat1));mat1=mat1/v;mat1=t(t(mat1)/v);. Dies könnte in einer Zeile kombiniert werden und mat1=mat1/v Art von Kopie in den Kommentaren anderer zu vermeiden.

Aber wenn das Muster in jedem Eintrag viel komplizierter ist, könnten Sie die Verwendung von inline::cxxfunction in Betracht ziehen, die wenig Modifikation des ursprünglichen R-Codes und viel viel schneller benötigt. Um eine Speicherkopie zu vermeiden, könnten Sie auch den Index der Matrix (in der Tat einen Doppelvektor) verwenden.

library(inline) 
include=" 
#include <math.h> 
#include <vector> 
" 

body=" 
NumericMatrix x(X); 
int nrow = x.nrow(); 
int ncol = x.ncol(); 
std::vector<double> diag(nrow); 
for (int i=0;i<nrow;i++){ 
    diag[i] = sqrt(x(i,i)); 
} 
double rho = .2; 
for(int j=1;j<ncol;j++){ 
     for(int i=0; i<(j-1);i++){ 
       if (j < i + nrow-4){// Change to your version 
        x(i,j) = rho/double(j-i)*diag[i]*diag[j]; 
        x(j,i) = x(i,j); 
} 
    } 
} 
return(x); 
" 
f1 <- cxxfunction(signature(X='matrix'),body,plugin='Rcpp',include=include) 

Verbrauch:

> dim(C) 
[1] 3000 3000 
> system.time(f1(C)) 
    user system elapsed 
    0.086 0.000 0.087 

daher die Kompilierung Zeit ignorieren, es ist ziemlich schnell, auch auf mein Chrom Buch.