如何在R中使用3个输入变量绘制三维图形？

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轴。

由于函数的参数应为三维向量A、B、O，所以我不能去掉“O”。

那么，我该如何改进当前的代码呢？如果需要更改函数，有人能推荐一个要使用的函数吗？（我认为可能可以使用重心坐标，但我考虑将其作为最后一件事来完成。）

- Dianafreedom

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，它也可能不会评估为TRUE（尝试运行3 - 2.9 == 0.1）。基于网格的评估使事情更加糟糕。如果你坚持让f有三个参数，那么你应该测试abs(A+B+O-1) < epsilon。我希望看到类似于以下的内容：

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的函数。您需要对子集进行操作以删除NAs，这在计算上似乎非常浪费。此外，请注意any(dat$A + dat$B + dat$O == 1)返回FALSE，因此您原始的约束测试确实在此网格上始终失败。

话虽如此，为什么您不使用函数内的约束条件确定给定A和B时的O呢？

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需要矢量化。

- Anders Ellern Bilgrau