R中的二元蒙特卡罗积分器无法运行

Question

我有这个代码来计算二元正态分布X~N(3,16)， Y~N(1,1)的概率，相关性-(2/3)落在矩形0<4和0><1，在r中有100,000个样本>

for(i in 1:100000){
X[i] = runif(100000)
Y[i] = runif(100000)
Z[i] = runif(100000)
if((3/(8*(pi)*sqrt(5)))*exp((-9/10)*(((X[i])^2/16)-((17/24)*(X[i]))+(X[i]*Y[i])/3)+((Y[i])^2)-(3*(Y[i]))+41/16) < Z[i])
count = count + 1}

count/100000表示X, Y落在上述矩形内的概率。

当我运行这个时，R返回一个错误，说

Error: object 'Z' not found
In addition: Warning messages:
1: In X[i] <- runif(1e+05) :
number of items to replace is not a multiple of replacement length
2: In Y[i] <- runif(1e+05) :
number of items to replace is not a multiple of replacement length

我该如何解决这个问题?

我尝试运行类似结构的蒙特卡罗积分器对x = 1的积分，但这返回了关于替换长度的类似错误。

Answer 1

X <- rep(NA, 100000)
Y <- rep(NA, 100000)
Z <- rep(NA, 100000)

X[1:100000] <- runif(100000)
Y[1:100000] <- runif(100000)
Z[1:100000] <- runif(100000)

定义X, Y和Z，然后分配runif。runif是矢量化的，所以不需要循环。这将消除你的错误和警告。

count = 0
for (i in 1:100000){
if((3/(8*(pi)*sqrt(5)))*exp((-9/10)*(((X[i])^2/16)-((17/24)*(X[i]))+(X[i]*Y[i])/3)+((Y[i])^2)-(3*(Y[i]))+41/16) < Z[i])
count = count + 1}

这是假设您出于其他原因想要保留X, Y和Z。你可以使用X, Y和Z作为长度为1的向量，并在每个循环中用新值替换它们。

你的计数计算也可以矢量化:

counter[1:100000] <- rep(NA, 100000)
counter <- ((3/(8*(pi)*sqrt(5)))*exp((-9/10)*((X^2/16)-((17/24)*X)+(X*Y)/3)+(Y^2)-(3*Y)+41/16) < Z)
sum(counter)/length(counter)