Блог

У меня есть уравнение ln(1 - P) * ln(1 - X2/d) = ln(1 - Q) * ln(1 - X1/d) . Я хочу определить d влауэ. Я написал следующий r код.

 Q = 0.75
P = 0.5
X2 = 5.57665
X1 = 1.473618
fun = function(d) (log(1-P)*log(1 - X2/d) - log(1 - Q)*log(1 - X1/d))^2
optimize(fun, c(0, 1000), maximum = TRUE, tol = 1e-10)
 

Верен ли приведенный выше код? Я установил интервал c(0, 1000) . Как выбрать подходящий вариант? Кроме того, есть ли какой-нибудь другой способ найти d ? Спасибо.

Ну, для начала, проверьте, есть ли вообще решение этого уравнения, которое, к счастью, находится в этом диапазоне x:

 library(ggplot2)

f1 <- function(d) log(1-P)*log(1-X2/d)
f2 <- function(d) log(1-Q)*log(1-X1/d)

(p <- ggplot()  
  xlim(-12, 12)   
  ylim(-10, 10)  
  geom_function(fun = f1, n = 1E4, aes(color = "f1"))  
  geom_function(fun = f2, n = 1E4, aes(color = "f2")))
 

Обратите внимание, что существует вертикальная асимптота при x = 0 и горизонтальная асимптота при y = 0. Это будет полезно для настройки интервала поиска.

Далее, вам нужно, чтобы ваше уравнение было правильным. Я не уверен, откуда ^2 это взялось. Перемещая LHS уравнения в RHS, вы получаете эту функцию:

 f <- function(d) log(1-P)*log(1-X2/d) - (log(1-Q)*log(1-X1/d))
 

Это мы можем решить с uniroot помощью функции из базы R :

 (eq <- uniroot(f, lower = -10, upper = -0.01))
$root
[1] -0.8258847

$f.root
[1] -5.59396e-06

$iter
[1] 10

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05
 

И подтвердить сюжетом:

 p   
  geom_point(data = data.frame(x = eq$root, y = f1(eq$root)), aes(x = x, y = y))
 

Примечание y могло быть найдено с любым f1 или f2 , поскольку на данный момент они равны. Также это значение было найдено с помощью оптимизации, поэтому есть некоторая ошибка, как показано на eq$estim.prec рисунке .