Ndtyjky

I am writing a log-likelihood surface for the function:

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

constraint by 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()

The error told me that there isn't object "O".

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

What I mean to do is to input A, B and O under the constraint that A+B+O=1, and then to plot the log-likelihood surface letting A:x-axis, B:y-axis, log-likelihood:z-axis.

I cannot get rid of "O" cause the instruction commands that the parameter of the function should be a 3-dimensional vector: A,B,O.

So what should I do to improve my current code?
If I need to change a function, can anyone suggest a function to use?
(I think maybe I can use barycentric coordinates but I consider it as the last thing I want to do.)

The outer function does not work the way you try to use it. outer takes two numeric arguments X and Y and a function argument FUN to which the first two arguments are applied. See ?outer. So it is not that there is no object O at all. Rather the error message in

Z <- outer(A, B, O, f)

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

#   object 'O' of mode 'function' was not found

means that the function O was not found. Indeed there is no such function.

There are also some problems with your f definition. First, it does not return anything. It saves a result as F but never returns it. Secondly, even if it returned F, the output would not always satisfy your constraint. When your constraint is not met, it simply outputs the value of O. Lastly, the comparison A+B+O==1 is a bad test as it may not evaluate to TRUE even if you'd expect it to due to floating point precision (try to run 3 - 2.9 == 0.1). The grid based evaluation makes things worse. Probably, you should be testing abs(A+B+O-1) < epsilon instead if you insist have three arguments for f. I.e. I would have expect something like:

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)

}

Then you can do something like:

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

But I do not see how you readily plot Z as a function of A and B from this. You'd need to subset to remove NAs and it seems very wasteful computationally. Also note any(dat$A + dat$B + dat$O == 1) returns FALSE, so your original constraint test indeed always fails on this grid.

With that said, why do you not determine O given A and B using your constraint within the function?

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

Does that look right? That is how persp and outer is intended to be used at least.

Of course, you can modify f so the warning messages are avoided. Just keep in mind that f needs to be vectorized.