Mnożenie macierzy wielowymiarowych w R

Question

Chciałbym wykonać pewne złożone wielowymiarowe mnożenie macierzy gdzie I pomnożyć przez określone marginesy z tablic.

Rozważmy następujący przykład, gdzie mam występowania cechą grup (A i B) przez niektórych marginesie populacji:

# setup data 

random=runif(4) 

group.prevalence <- aperm (array(c(random,1-random), 
        dim=c(2,2,2), 
        dimnames=list(age=c("young","old"), 
           gender=c("male","female"), 
           group=c("A","B"))) , c(3,1,2)) 

group.prevalence 
# A + B = 1 

Przypuśćmy teraz, że mam populację zainteresowanie ...

population <- round(array(runif(4, min=100,max=200) %o% c(1,1*(1+random[1]),1*(1+random[1])^2), 
          dim=c(2,2,3), dimnames=list(age=c("young","old"), 
                 gender=c("male","female"), 
                 year=c("year1","year2","year3")))) 

population 

... dla których chciałbym obliczyć występowanie "A" i "B".

Zła rozwiązaniem byłoby, aby wypełnić to wszystko w pętli:

# bad solution 
grouped.population <- array(NA, dim=c(2,2,2,3), 
          dimnames=list(group=c("A","B"), 
              age=c("young","old"), 
              gender=c("male","female"), 
              year=c("year1","year2","year3"))) 

for (group in c("A","B")) 
    for(gender in c("male","female")) 
    for (age in c("young","old")) 
     grouped.population[group,age,gender,] <- group.prevalence[group,age,gender] * population[age,gender,] 

Ale przypuszczam, że jakieś zastosowanie mogłoby się przydać, ewentualnie aaply plyr, ponieważ wymiary wynikiem powinno być zachowane. Próbowałem:

library(plyr) 
aaply(population, c(1,2), function(x) x * group.prevalence) 
# too many dimensions 

Z zadowoleniem przyjmuję wszelkie sugestie.

Answer 1

dla konkretnego przypadku, można obliczyć:

out <- rep(group.prevalence, times=last(dim(population))) * 
     rep(population, each=first(dim(group.prevalence))) 

a następnie można ustawić wymiary tego array:

array(out, dim=c(2,2,2,3), 
     dimnames=list(group=c("A","B"), 
        age=c("young","old"), 
        gender=c("male","female"), 
        year=c("year1","year2","year3"))) 

Kluczem jest align wymiary dwóch tablice za pomocą transpozycji wymiarów i rozszerzenia/replikacji w celu uzupełnienia brakujących wymiarów, które są w drugiej tablicy. Ogólnie rzecz biorąc, procedura:

1. Określenie przecinające wymiary. Tutaj jest (age,gender).
2. Dla argumentu po lewej stronie pomnożenia, group.prevalence, zmień wymiary (używając aperm) tak, aby wszystkie niepodzielące wymiary (tj. group) były pierwsze. Następnie replikuj tę tablicę razy (używając), gdziejest rozmiarem nieprzecinających się wymiarów (tj. year) z argumentu po prawej stronie, population.
3. W przypadku argumentu po prawej stronie mnożenia, population, należy przerobić wymiary tak, aby wszystkie nieprzecinające się wymiary (tj. year) były ostatnie. Następnie replikuj każdy element z tablicy M razy (przy użyciu each), gdzie M jest rozmiarem nieprzecinających się wymiarów (tj. group) z lewej strony argumentu, group.prevalence.
4. Następnie wystarczy (tablica) pomnożyć, co jest wektoryzowane i szybkie.
5. Wymiary wyniku to po prostu nie przecinające się wymiary argumentu po lewej stronie, po których następują przecinające się wymiary, a następnie nieprzecinające się wymiary prawej strony (tj. (group, age, gender, year)).Następnie możesz przesłać te wymiary w razie potrzeby do wyjścia, aby uzyskać to, co chcesz.

Jako sprawdzenia:

# bad solution 
grouped.population <- array(NA, dim=c(2,2,2,3), 
          dimnames=list(group=c("A","B"), 
              age=c("young","old"), 
              gender=c("male","female"), 
              year=c("year1","year2","year3"))) 

for (group in c("A","B")) 
    for(gender in c("male","female")) 
    for (age in c("young","old")) 
     grouped.population[group,age,gender,] <- group.prevalence[group,age,gender] * population[age,gender,] 

# another approach 
grouped.population2 <- array(rep(group.prevalence, times=last(dim(population))) * 
          rep(population, each=first(dim(group.prevalence))), 
          dim=c(2,2,2,3), 
          dimnames=list(group=c("A","B"), 
              age=c("young","old"), 
              gender=c("male","female"), 
              year=c("year1","year2","year3"))) 

# check 
all.equal(grouped.population,grouped.population2) 
##[1] TRUE 

---

zaktualizowany o wartość odniesienia:

library(microbenchmark) 

f1 <- function(group.prevalence, population) { 
    grouped.population <- array(NA, dim=c(2,2,2,3), 
           dimnames=list(group=c("A","B"), 
              age=c("young","old"), 
              gender=c("male","female"), 
              year=c("year1","year2","year3"))) 
    for (group in c("A","B")) { 
    for(gender in c("male","female")) { 
     for (age in c("young","old")) { 
     grouped.population[group,age,gender,] <- group.prevalence[group,age,gender] * population[age,gender,]}}} 
} 

f2 <- function(group.prevalence, population) { 
    grouped.population2 <- array(rep(group.prevalence, times=last(dim(population))) * 
           rep(population, each=first(dim(group.prevalence))), 
           dim=c(2,2,2,3), 
           dimnames=list(group=c("A","B"), 
              age=c("young","old"), 
              gender=c("male","female"), 
              year=c("year1","year2","year3"))) 
} 

print(microbenchmark(f1(group.prevalence, population))) 
##Unit: microseconds 
##        expr  min  lq  mean median  uq  max neval 
## f1(group.prevalence, population) 101.473 103.998 149.2562 106.8865 115.372 1185.32 100 
print(microbenchmark(f2(group.prevalence, population))) 
##Unit: microseconds 
##        expr min  lq  mean median  uq  max neval 
## f2(group.prevalence, population) 66.392 67.672 70.19873 68.454 69.4205 173.284 100 

Uważam, że wydajność będzie rozchodzą się jeszcze bardziej jak liczba wymiarów i wielkości w każdym wymiarze wzrasta.