Как построить трехмерный график с 3 входными переменными в R?

Question

Я пишу логарифмическую поверхность для функции:

LN [Pr (Y_A = 186, Y_B = 38, Y_ {АВ} = 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 при условии, что A + B + O = 1, а затем построить график логарифмической вероятности, обозначив A: ось x, B: ось y, срубы вероятность: ось г

.

I не может избавиться от "O", потому что команды команды указывают, что параметром функции должен быть трехмерный вектор: 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. То есть Я бы ожидал что-то вроде:

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 необходимо векторизовать.