Улучшение алгоритмов сходимости с помощью числового итератора в R - PullRequest
0 голосов
/ 02 февраля 2019

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

Вот изображение y, построенное на фоне x Actual Problem Я создал рабочийНапример, определив асимптотическую функцию, которая адекватно отражает проблему

y <- (-1/x)+0.05

, которая при построении графика дает

x <- 1:100
y <- (-1/x)+0.05
DT <- data.frame(cbind(x=x,y=y))
ggplot(DT, aes(x, y)) + geom_point() + geom_hline(yintercept=0, color="red")

Working Example

Я разработал TWO итерационные алгоритмы для аппроксимации x-перехвата.

Решение 1 : x изначально очень маленькое, пошаговое 1...nраз.Размер шагов предопределен, начинайте с большого (10-кратное увеличение).После каждого шага вычисляется y.i.Если abs(y.i) < y[i-1], то этот большой шаг повторяется, если только y.i не изменит знак, который указывает, что этот шаг превысил x-intercept.Если алгоритм отклоняется, мы возвращаемся назад и делаем меньший шаг (увеличение в 2 раза).С каждым выбросом делаются все меньшие и меньшие шаги, начиная с 10 *, 2 *, 1,1 *, 1,05 *, 1,01 *, 1,005 *, 1,001 *.

x.i <- x <- runif(1,0.0001,0.001)
y.i <- y <- (-1/x.i)+0.05
i <- 2
while(abs(y.i)>0.0001){
  x.i <- x[i-1]*10
  y.i <- (-1/x.i)+0.05
  if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
    x <- c(x,x.i); y <- c(y,y.i)
  } else {
    x.i <- x[i-1]*2
    y.i <- (-1/x.i)+0.05
    if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
      x <- c(x,x.i); y <- c(y,y.i)
    } else {
      x.i <- x[i-1]*1.1
      y.i <- (-1/x.i)+0.05
      if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
        x <- c(x,x.i); y <- c(y,y.i)
      } else {
        x.i <- x[i-1]*1.05
        y.i <- (-1/x.i)+0.05
        if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
          x <- c(x,x.i); y <- c(y,y.i)
        } else {
          x.i <- x[i-1]*1.01
          y.i <- (-1/x.i)+0.05
          if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
            x <- c(x,x.i); y <- c(y,y.i)
          } else {
            x.i <- x[i-1]*1.005
            y.i <- (-1/x.i)+0.05
            if(abs(y.i)<abs(y[i-1]) & sign(y.i)==sign(y[i-1])){
              x <- c(x,x.i); y <- c(y,y.i)
            } else {
              x.i <- x[i-1]*1.001
              y.i <- (-1/x.i)+0.05
            }
          }
        }
      }
    }
  }
  i <- i+1
}

Решение 2 : Этот алгоритм основан на идеях метода Ньютона-Рафсона, где шаги основаны на скорости изменения y.Чем больше изменение, тем меньше предпринятых шагов.

x.i <- x <- runif(1,0.0001,0.001)
y.i <- y <- (-1/x.i)+0.05
i <- 2
d.i <- d <- NULL
while(abs(y.i)>0.0001){
  if(is.null(d.i)){
    x.i <- x[i-1]*10
    y.i <- (-1/x.i)+0.05
    d.i <- (y.i-y[i-1])/(x.i-x[i-1])
    x <- c(x,x.i); y <- c(y,y.i); d <- c(d,d.i)
  } else {
    x.i <- x.i-(y.i/d.i)
    y.i <- (-1/x.i)+0.05
    d.i <- (y.i-y[i-1])/(x.i-x[i-1])
    x <- c(x,x.i); y <- c(y,y.i); d <- c(d,d.i)
  }
  i <- i+1
}

Сравнение

  1. Решение 1 требует последовательно меньше итераций, чем Решение 2 (1/2, если не 1/3).
  2. Решение 2 более элегантно, не требует произвольного уменьшения размера шага.
  3. Я могу представить несколько сценариев, когда Решение 1 застревает (например, если даже на самом маленьком шаге цикл несходятся на достаточно малом значении для y.i)

Вопросы

  1. Есть ли лучший (меньшее количество итераций) способ аппроксимации x-intercept в таких сценариях?
  2. Может ли кто-нибудь указать мне какую-нибудь литературу, в которой рассматриваются такие проблемы (желательно написанные понятным образом для начинающего, такого как я)?
  3. Любые предложения по номенклатуре или ключевым словам, которые представляют этот класс проблемы / алгоритмадобро пожаловать.
  4. Можно ли улучшить представленные мной решения?
  5. Есть предложения о том, как сделать заголовок / вопрос более доступным?для более широкого сообщества или экспертов с потенциальными решениями приветствуются.

1 Ответ

0 голосов
/ 06 февраля 2019

Следуя рекомендациям и предложениям @Lyngbakr и @LutzL, алгоритм поиска корней, известный как Метод Брента , оказался эффективным и значительно более быстрым, чем моя реализация Ньютона-Рафсона (Решение 2).Этот алгоритм реализован uniroot в R.

f <- function(x){(-1/x)+0.05}
uniroot(f, interval=c(0,100), maxiter=100)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...