Appliquer une fonction sur chaque élément d'une matrice

Question

L'année dernière, j'ai utilisé un code qui produit les différentes valeurs de probabilité cumulative d'une distribution normale trivariée, lorsque les paramètres changent de valeur simultanément selon une fonction. Je le code suivant:

library(mvtnorm) 
Y <- mapply(function(x,y,z) 
    pmvnorm(mean = c(18, 12.72, (18*(x+y) +12.72*z)), 
      sigma = { 
      s1 <- matrix(c(5.7, 0, 5.7*(x+y), 
          0, 30.38, 30.38*z, 
          5.7*(x+y), 30.38*z, 5.7*(x+y)^2+30.38*(z)^2), 
         3) 
      replace(s1,s1==0, 1e-20) 
      }, 
      lower = c(15, -Inf, p+15*y), 
      upper = c(Inf, 15, Inf)), 
    m1, m2, m3) 

où m1, m2, m3 (associée à x, y et z) sont les vecteurs déjà définis. Comme ce sont des vecteurs 1x10, le code produit un vecteur 1x10.

Maintenant, si m1, m2 et m3 sont des matrices Je produis 10x36 trois matrices avec ce code:

a <- 5 
b <- 5 
vals <- rep(0:a, (b+1)) 
x <- rep(1:(a+b), ((a+1)*(b+1))) 
c <- matrix(pmin(x, vals[rep(1:((a+1)*(b+1)), each = (a+b))]), nrow = (a+b)) 
m1 <- c 
d1 <- matrix(rep(c(1:(a+b)), ((a+1)*(b+1))), ncol=((a+1)*(b+1)), nrow=(a+b),  byrow=F) 
d2 <- matrix(rep(rowSums(expand.grid(0:a, 0:b)), (a+b)), ncol=((a+1)*(b+1)), nrow=(a+b), byrow=T) 
d3 <- pmax((d1-d2),0) 
d4 <- matrix(rep(0:a, b+1), ncol=((a+1)*(b+1)), nrow=(a+b), byrow=T) 
d5 <- a-d4 
d6 <- pmin(d3,d5) 
d7 <- matrix(nrow = (a+b), ncol = ((a+1)*(b+1))) 
for (i in 1:(a+b)) { 
    for (j in 1:((a+1)*(b+1))) { 
    d7[i,j] <- a 
    } 
} 
d8 <- pmin(d7,d1) 
d <- d6-d8 
m2 <- d 
e1 <- pmin(d1,d2) 
e2 <- e1-d4 
e3 <- matrix(nrow = (a+b), ncol = ((a+1)*(b+1))) 
for (i in 1:(a+b)) { 
    for (j in 1:((a+1)*(b+1))) { 
    e3[i,j] <- 0 
    } 
} 
e <- pmax(e2, e3) 
m3 <- e 
r <- 0.04 
k <- rep(rowSums(expand.grid(0:a, 0:b))) 
p <- k*(15*(1-r)^(k-1)) 
p <- t(matrix(rep(p,(a+b)), ncol=(a+b), nrow=((a+1)*(b+1)))) 

Maintenant, je veux produire une matrice 10x36 où chaque élément est la valeur de la probabilité cumulative d'une distribution normale trivariée, où la matrice de moyenne et de covariance dépend des éléments des matrices m1 m2 et m3. J'ai essayé à nouveau avec le code:

Y <- mapply(function(x,y,z) 
    pmvnorm(mean = c(18, 12.72, (18*(x+y) +12.72*z)), 
      sigma = { 
      s1 <- matrix(c(5.7, 0, 5.7*(x+y), 
          0, 30.38, 30.38*z, 
          5.7*(x+y), 30.38*z, 5.7*(x+y)^2+30.38*(z)^2), 
         3) 
      replace(s1, s1==0, 1e-20) 
      }, 
      lower = c(15, -Inf, p+15*y), 
      upper = c(Inf, 15, Inf)), 
    m1, m2, m3) 

Mais j'obtiens l'erreur suivante:

Error in checkmvArgs(lower = lower, upper = upper, mean = mean, corr = corr,  : 
    at least one element of ‘lower’ is larger than ‘upper’ 
Warning message: 
In cbind(lower, upper, mean) : 
number of rows of result is not a multiple of vector length (arg 2). 

Où est l'erreur?

Answer 1

Depuis m1, m2, m3, etp sont toutes dimensionnées 10x36, j'inférant que vous avez l'intention d'utiliser chaque valeur individuelle de p au lieu de l'ensemble de la matrice à chaque fois. Par exemple, si vous entrez dans l'appel à pmvnorm, vous verrez:

Y <- mapply(function(x,y,z) { 
    browser() 
    pmvnorm(mean = c(18, 12.72, (18*(x+y) +12.72*z)), 
      sigma = { 
      s1 <- matrix(c(5.7, 0, 5.7*(x+y), 
          0, 30.38, 30.38*z, 
          5.7*(x+y), 30.38*z, 5.7*(x+y)^2+30.38*(z)^2), 
         3) 
      replace(s1, s1==0, 1e-20) 
      }, 
      lower = c(15, -Inf, p+15*y), 
      upper = c(Inf, 15, Inf)) 
}, m1, m2, m3) 
# debug at c:/Users/r2/AppData/Local/Temp/foo.R!12268ZNM#3: pmvnorm(mean = c(18, 12.72, (18 * (x + y) + 12.72 * z)), sigma = { 
#  s1 <- matrix(c(5.7, 0, 5.7 * (x + y), 0, 30.38, 30.38 * z, 
#   5.7 * (x + y), 30.38 * z, 5.7 * (x + y)^2 + 30.38 * (z)^2), 
#   3) 
#  replace(s1, s1 == 0, 1e-20) 
# }, lower = c(15, -Inf, p + 15 * y), upper = c(Inf, 15, Inf)) 
# Browse[2]> 
x 
# [1] 0 

Cela devrait.

# Browse[2]> 
dim(p) 
# [1] 10 36 

Je pense que cela devrait plutôt être:

# Browse[2]> 
p 
# [1] 0 

Dans ce cas, vous pouvez commencer à résoudre le problème en ajoutant p à la liste des arguments:

pm <- p # 'pm' is 'p-matrix' (?) 
rm(p) 

(I a fait cela pour supprimer l'ambiguïté d'avoir p à la fois à l'intérieur et à l'extérieur de l'appel à mapply. Cela fonctionne de toute façon, mais le dépannage d'une fermeture où vous n'êtes pas certain qui instance de la variable est utilisée peut être frustrant et certainement pas intuitive.)

set.seed(2) 
Y <- mapply(function(x,y,z,p) 
    pmvnorm(mean = c(18, 12.72, (18*(x+y) +12.72*z)), 
      sigma = { 
      s1 <- matrix(c(5.7, 0, 5.7*(x+y), 
          0, 30.38, 30.38*z, 
          5.7*(x+y), 30.38*z, 5.7*(x+y)^2+30.38*(z)^2), 
         3) 
      replace(s1, s1==0, 1e-20) 
      }, 
      lower = c(15, -Inf, p+15*y), 
      upper = c(Inf, 15, Inf)), 
    m1, m2, m3, pm) 

Un problème que vous verrez est que mapply a retourné un vecteur. Ceci peut être corrigé en affectant la même dimensionnalité que la variable d'entrée. (Depuis que vous avez appelé matrix avec la valeur par défaut byrow=FALSE, cela fonctionne en l'état. Si vous aviez fait byrow=TRUE, vous auriez besoin d'ajuster ici.)

dim(Y) <- dim(m1) 
dim(Y) 
# [1] 10 36 
Y[1:4,1:5] 
#   [,1]  [,2]  [,3]  [,4]  [,5] 
# [1,] 0.2957254 0.2957254 0.0000000 0.0000000 0.0000000 
# [2,] 0.2957254 0.2957254 0.5914508 0.0000000 0.0000000 
# [3,] 0.2957254 0.2957254 0.5914508 0.5914508 0.0000000 
# [4,] 0.2957254 0.2957254 0.5914508 0.5914508 0.5914508