Найти root с помощью uniroot

Question

Я пытаюсь найти корень следующей функции (на основе гаммы ( gamma()) function) с помощью функции:

      cv = 0.056924/1.024987^2

fx2 = function(theta, eta){
  p1 = 1 - 2/(theta*(1-eta))
  p2 = 1 - 1/(theta*(1-eta))
  return(( gamma(p1)/(gamma(p2))^2 ) - (cv+1) )
}

Эта функция дает мне следующий график:

      v = seq(0, 1, 0.01)
plot(v, fx2(3.0, v), type='l' )

Мне кажется, что корень этой функции близок к 0,33, но uniroot() функция не находит корень, возвращая следующий результат:

      uniroot(fx2, interval = c(0,0.3), theta=3 )

Ошибка в uniroot (fx2, interval = c (0, 0,3), theta = 3): значения f () в конечных точках не противоположного знака

Как мне найти корень этой функции? Есть ли другие пакеты с более точным алгоритмом?

Answer 1

Сначала я переписал вашу функцию, чтобы (необязательно) выразить с точки зрения вычислений, которые сначала выполняются в логарифмической шкале (через ), а затем возведен в степень. Это более стабильно численно, и последствия станут ясны ниже ... (Возможно, я облажался с вычислением в логарифмическом масштабе - вам следует перепроверить его.)

      cv = 0.056924/1.024987^2
fx3 <- function(theta, eta, lgamma = FALSE) {
  p1 <- 1 - 2/(theta*(1-eta))
  p2 <- 1 - 1/(theta*(1-eta))
  if (lgamma) {
    val <- exp(lgamma(p1) - 2*lgamma(p2)) - (cv+1)
  } else {
    val <- ( gamma(p1)/(gamma(p2))^2 ) - (cv+1)
  }
}

Вычислите функцию с масштабированием журнала и без него:

      x <- seq(0, 1, length.out = 20001)
v <- sapply(x, fx3, theta = 3.0, lgamma =  TRUE)
v2 <- sapply(x, fx3, theta = 3.0, lgamma =  FALSE)

Найдите корень (более подробное объяснение ниже):

      uu <- uniroot(function(eta) fx3(3.0, eta, lgamma = TRUE),
        c(0.4, 0.5))

Постройте это:

      par(las=1, bty="l")
plot(x, abs(v), col = as.numeric(v<0) + 1, type="p", log="y",
     pch=".", cex=3)
abline(v = uu$root, lty=2)
cvec <- sapply(c("blue","magenta"), adjustcolor, alpha.f = 0.2)
points(x, abs(v2), col=cvec[as.numeric(v2<0) + 1], pch=".", cex=3)

Здесь я изображаю абсолютное значение в логарифмической шкале со знаком, обозначенным цветом (черный / синий>0, красный> пурпурный <0). Черный / красный - расчет в логарифмической шкале, синий / пурпурный - исходный расчет. Я также построил график функции в очень высоком разрешении, чтобы попытаться избежать пропуска или искажения характеристик.

Здесь происходит много странных вещей.

• обе версии функции делают что-то интересное около x=1/3; исходная версия выглядит как полюс (значение расходится до +∞, «возвращается» из -∞), в то время как вычисление в логарифмической шкале увеличивается до +∞ и возвращается без изменения знака.
• вычисление в логарифмическом масштабе имеет корень около x = 0,45 (абсолютное значение становится маленьким, когда знак переворачивается), но исходное вычисление нет - предположительно из-за какой-то катастрофической потери точности? Если мы дадим границы, не включающие полюс, он может найти этот корень.
• есть другие полюса и / или корни при больших значениях x, которые я не исследовал.

Все это в основном говорит о том, что довольно опасно возиться с этой функцией, не зная, каковы ее математические свойства. Я обнаружил кое-что путем численного исследования, но было бы лучше проанализировать функцию, чтобы вы действительно знали, что происходит; любое числовое исследование можно обмануть, если функция ведет себя достаточно странно.