在矩陣的每個元素上應用函數

Question

去年，我使用了一個code，該函數根據函數同時更改值，從而產生三元正態分佈的累積概率的不同值。我使用此代碼：

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) 

其中m1，m2，m3（與x，y和z相關）是已定義的向量。由於它們是1x10矢量，代碼會生成一個矢量1x10。

現在，如果M1，M2和M3的矩陣我產生三個矩陣10x36與此代碼：

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

現在我想產生一個矩陣10x36每個元素的累積概率的價值的三元正態分佈，其中平均值和協方差矩陣依賴於矩陣m1 m2和m3的元素。 我與代碼再次嘗試：

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) 

，但我得到了以下錯誤：

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). 

哪裏錯誤？

Answer 1

由於m1，m2，m3，和p全部尺寸10x36，我推斷你打算每次使用的p每個單獨的值，而不是整個矩陣的。例如，如果你步入調用pmvnorm，你會看到：

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 

這是預期。

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

我認爲這應該改爲：

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

在這種情況下，你就可以開始加入p到的參數列表來解決這個問題：

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

（我這樣做是爲了消除在mapply的呼叫之內和之外都有p的不明確性。它無論如何都可以工作，但是如果你不確定該變量的情況下被使用，可令人沮喪，當然也不直觀。）

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) 

一個問題，你會看到的是，mapply返回的向量。

head(Y) 
# [1] 0.2957254 0.2957254 0.2957254 0.2957254 0.2957254 0.2957254 
length(Y) 
# [1] 360 
dim(Y) 
# NULL 

這可以通過分配與輸入變量相同的維度來解決。 （既然你叫matrix與byrow=FALSE默認的，這個工作原樣。如果你做了byrow=TRUE，你需要在此進行調整。）

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