Алгоритм классификации разрыва головы / хвоста в R

Question

Я пытаюсь реализовать алгоритм классификации разрыва головы / хвоста в R ( см. Здесь ). Этот относительно новый алгоритм является менее дорогой вычислительной альтернативой другим методам классификации, используемым в картографии для сильно искаженных данных.

До сих пор я искал в качестве шаблона код в Python ( см. Здесь ) с относительно успешным результатом. Вот моя реализация в R:

# fake data to classify
pareto_data <- c()
for (i in 1:100){
  pareto_data[i] <- (1.0/i)^1.16
}

# head/tail breaks algorithm
ht <- function(data){
  ln <- length(data)
  mn <- mean(data)
  res <- append(c(),mn) # this is where I was hopping to store my output
  head <- subset(data,data>=mn)
  while (length(head)>=1 & length(head)/ln <= 0.40){
    print(res)
    return(ht(head))
  }
  #return(res)
}

ht(pareto_data)

В результате выполнения кода выше, я смог напечатать следующее:

[1] 0.03849691
[1] 0.1779904
[1] 0.4818454 

Этот вывод, скорее всего, совпадает с запустить оригинальный код Python, который я использовал в качестве шаблона. Однако мне не удалось сохранить его ни в векторе, ни в списке.

Я был бы очень благодарен, если бы вы могли дать подсказки, чтобы преодолеть эту проблему, а также улучшить мой код (что не совсем такой же, как и в Python, особенно в условиях оператора while.

Answer 1

Возможная рекурсивная версия алгоритма может быть следующей.

ht_breaks <- function(x){
  ht_inner <- function(x, mu){
    n <- length(x)
    mu <- c(mu, mean(x))
    h <- x[x > mean(x)]
    if(length(h) > 1 && length(h)/n <= 0.4){
      ht_inner(h, mu)
    } else mu
  }
  ht_inner(x, NULL)
}

pareto_data <- (1.0/(1:100))^1.16
ht_breaks(pareto_data)
#[1] 0.03849691 0.17799039 0.48184535