r - 如何在 R 中绘制具有 3 个输入变量的 3-D 图？

Question

我正在为函数编写对数似然曲面:

ln[Pr(Y_A=186,Y_B=38,Y_{AB}=13,Y_O=284)]= ln(G+186*ln(A^2+2*A*O)+38*ln(B^2+2*B*O)+13*ln(2*A*B)+284*ln( O^2))

由A+B+O=1约束

A = seq(0.0001, .9999,length=50)
B = A
O = A
G = 1.129675e-06   
f = function(A,B,O){F = ifelse(A+B+O==1,
                    G+186*log(A*A+2*A*O)+38*log(B*B+2*B*O)+13*log(2*A*B)+284*log(O*O), O)}
Z <- outer(A, B, O, f)
png()
persp(A,B,Z, theta=60, phi=30 )
dev.off()

错误告诉我没有对象“O”。

Error in get(as.character(FUN), mode = "function", envir = envir)

我的意思是在A+B+O=1的约束条件下输入A，B，O，然后画出对数似然面让A:x轴，B:y轴，对数似然:z 轴。

我不能去掉“O”，因为指令命令函数的参数应该是一个 3 维向量:A、B、O。

那么我应该怎么做才能改进我当前的代码呢？如果我需要更改功能，有人可以建议使用功能吗？(我想也许我可以使用重心坐标，但我认为这是我最不想做的事情。)

Answer 1

outer功能无法按照您尝试使用的方式工作。 outer接受两个数字参数 X和 Y和一个函数参数 FUN应用了前两个参数。参见 ?outer .所以不是没有对象O根本。而是

中的错误消息

Z <- outer(A, B, O, f)
#Error in get(as.character(FUN), mode = "function", envir = envir) : 
#   object 'O' of mode 'function' was not found

表示函数 O没找到。确实没有这个功能。

你的f也有一些问题定义。首先，它不返回任何东西。它将结果保存为 F但永远不会归还它。其次，即使它返回F ，输出并不总是满足您的约束。当你的约束不满足时，它只是输出 O 的值。 .最后是对比A+B+O==1是一个糟糕的测试，因为它可能不会评估为 TRUE即使你期望它是由于浮点精度(尝试运行 3 - 2.9 == 0.1 )。基于网格的评估使事情变得更糟。也许，你应该测试 abs(A+B+O-1) < epsilon相反，如果你坚持为 f 提供三个参数. IE。我会期待这样的事情:

f <- function(A, B, O){
  G <- 1.129675e-06
  epsilon <- 1e-3
  ifelse(abs(A+B+O-1) < epsilon,
         G+186*log(A*A+2*A*O)+38*log(B*B+2*B*O)+13*log(2*A*B)+284*log(O*O),
         NA)
}

然后你可以这样做:

dat <- expand.grid(A = A, B = B, O = O) # All combinations of A, B, O
dat$Z <- f(dat$A, dat$B, dat$O) # Apply function
head(dat)
#           A     B     O  Z
#1 0.00010000 1e-04 1e-04 NA
#2 0.02050408 1e-04 1e-04 NA
#3 0.04090816 1e-04 1e-04 NA
#4 0.06131224 1e-04 1e-04 NA
#5 0.08171633 1e-04 1e-04 NA
#6 0.10212041 1e-04 1e-04 NA

但我看不出您如何轻松地将 Z 绘制为 A 和 B 的函数。您需要子集来删除 NA，这在计算上似乎非常浪费。另请注意 any(dat$A + dat$B + dat$O == 1)返回 FALSE , 所以你原来的约束测试确实总是在这个网格上失败。

话虽如此，你为什么不确定O给出A和 B在函数中使用你的约束？

A <- seq(0.0001, .9999,length=50)
B <- A

f <- function(A, B){
  G <- 1.129675e-06
  O <- 1 - A - B
  out <-  G+186*log(A*A+2*A*O)+38*log(B*B+2*B*O)+13*log(2*A*B)+284*log(O*O)
  return(out)
}

Z <- outer(A, B, f)
#Warning messages:
#1: In log(A * A + 2 * A * O) : NaNs produced
#2: In log(B * B + 2 * B * O) : NaNs produced

Z[is.infinite(Z)] <- NA

persp(A, B, Z, theta=60, phi=30, zlim = range(Z, na.rm = TRUE))

这样看起来对吗？就是这样persp和 outer旨在至少使用。

当然可以修改f因此避免了警告消息。请记住 f需要矢量化。