【发布时间】：2017-06-15 15:44:00
【问题描述】：
去年我使用了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).
哪里出错了？
【问题讨论】：
-
我建议我编辑的代码使其更具可读性（例如，缩进，诚然是主观的）。我假设您使用的是
library(mvtnorm)，所以我也添加了它。如果我不正确，请更正。 （并且请确保您自己包含它。）如果您之前的 SO 问题/答案与上下文相关，则包含该链接可能会很有用，但最好这个问题是独立的。