Поиск областей полигонов: интегрирование при вертикальных и горизонтальных геометрических ограничениях

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

y <- c(-1, 5, 2, 3.5, 1, 5.5, -2, 3, -1)
plot(1:length(y), y, type = "l")
abline(h = 0)

Я пытаюсь вычислить области, подверженные вертикальным и горизонтальным геометрическим ограничениям:

  • (вертикально) под кривой, но выше нуля;
  • (горизонтальный) между двумя соседними локальными минимумами.

То есть мне нужны площади полигонов A, B, C и D на изображении ниже.

пытается вычислить площади A, B, C и D

Я сейчас борюсь с двумя вещами:

  1. Я определил индекс положения локальных минимумов, используя: which(diff(sign(diff(y))) == 2) + 1, но это не дало мне верхнее значение x для C или нижнее значение x для D. Как получить те точки, где кривая пересекает нули?

  2. Я думаю, что если я смогу получить 1) правильный список локальных минимумов выше нуля, 2) эти точки пересечения в нуле, 3) правильный список локальных максимумов выше нуля, я знаю все граничные точки A, B, C и D, чтобы можно было вычислить их площади. Но это не кажется простым кодированием в R. Действительно ли это самый простой способ решить мою проблему или есть методы получше?


r area polygon integral numerical-integration
person Jennifer    schedule 03.09.2018    source источник

Ответы (1)


arrow_upward
2
arrow_downward 

## (x, y) data
y <- c(-1, 5, 2, 3.5, 1, 5.5, -2, 3, -1)
x <- 1:length(y)

Аналитический метод

Желаемый расчет можно выполнить в два этапа:

  1. кусочный интеграл на каждом линейном отрезке. Нет никаких трудностей, если вы хотите интегрировать только пропорцию выше нуля (подробности см. ниже);
  2. агрегировать по частям результаты соответствующим образом для областей A, B, C и D (подробности см. ниже).

шаг 1: кусочный интеграл для положительной доли

Если есть n (x, y) данных, будет (n - 1) сегментов. Обозначим (xl, yl) как левую точку сегмента, а (xr, yr) - как правую точку.

  1. если (yl < 0) && (yr < 0), весь сегмент ниже нуля;
  2. если (yl > 0) && (yr > 0), весь сегмент выше нуля;
  3. если (yl < 0) && (yr > 0), сегмент увеличивается, пересекая ноль;
  4. если (yl > 0) && (yr < 0), сегмент уменьшается, пересекая ноль.

В случаях 3 и 4 обозначьте (xm, 0) как точку пересечения. xm легко определить. Уравнение для отрезка прямой

y = yl + (yr - yl) * (x - xl) / (xr - xl)

Установив y в 0, вы получите

xm = xl - yl * (xr - xl) / (yr - yl)

Поскольку вы хотите интегрировать положительную долю каждого сегмента, у нас есть для каждого случая:

  1. интеграл равен 0;
  2. площадь представляет собой трапецию, а интеграл равен (yl + yr) * (xr - xl) / 2;
  3. область представляет собой треугольник между (xm, 0) и (xr, yr); интеграл равен yr * (xr - xm) / 2;
  4. область представляет собой треугольник между (xl, yl) и (xm, 0); интеграл равен yl * (xm - xl) / 2.

Поскольку в конечном итоге вы захотите применить вычисления к длинным векторам, я бы представил вычисления в функции Rcpp.

library(Rcpp)

cppFunction('NumericVector foo_cpp (NumericVector x, NumericVector y) {
  int n_segments = x.size() - 1;
  NumericVector integral(n_segments);
  double xl, xr, yl, yr, xm; int i;
  for (i = 0; i < n_segments; i++) {
    xl = x[i]; xr = x[i + 1];
    yl = y[i]; yr = y[i + 1];
    if (yl < 0 && yr < 0) integral[i] = 0.0;
    if (yl > 0 && yr > 0) integral[i] = 0.5 * (yl + yr) * (xr - xl);
    if (yl < 0 && yr > 0) {
      xm = xl - yl * (xr - xl) / (yr - yl);
      integral[i] = 0.5 * yr * (xr - xm);
      }
    if (yl > 0 && yr < 0) {
      xm = xl - yl * (xr - xl) / (yr - yl);
      integral[i] = 0.5 * yl * (xm - xl);
      }
    }
  return integral;
  }')

z <- foo_cpp(x, y)
#[1] 2.083333 3.500000 2.750000 2.250000 3.250000 2.016667 0.900000 1.125000

Мне лень заниматься дальнейшей оптимизацией кода. Его скорость достаточно хороша для вашего практического использования.

Шаг 2: объединение

Вы фактически разрезаете сегменты на куски по локальным минимумам и стремитесь вычислить интеграл для каждого куска.

Индекс позиции для локальных минимумов (как вы определили в своем вопросе):

which(diff(sign(diff(y))) == 2) + 1
#[1] 3 5 7

Это означает, что сегменты должны быть разделены точками останова:

b <- which(diff(sign(diff(y))) == 2)
#[1] 2 4 6

То есть,

## number of segments per chunk
n_chunks <- length(x) - 1
n_segments_per_chunk <- diff(c(0, b, n_chunks))
#[1] 2 2 2 2

## grouping index for each chunk
grp <- rep.int(seq_along(n_segments_per_chunk), n_segments_per_chunk)
#[1] 1 1 2 2 3 3 4 4

Таким образом, области A, B, C и D:

sapply(split(z, grp), sum)
#       1        2        3        4 
#5.583333 5.000000 5.266667 2.025000

Численный метод

## original linear interpolation function
f <- approxfun(x, y)

## a function zeroing out below-zero part of `f`
g <- function (x) {
  fx <- f(x)
  ifelse(fx > 0, fx, 0)
  }

## local minima
x_minima <- x[which(diff(sign(diff(y))) == 2) + 1]

## break points for numerical integration
xx <- c(x[1], x_minima, x[length(x)])

## integration will happen on:
# cbind(xx[-length(xx)], xx[-1])
#     [,1] [,2]
#[1,]    1    3  ## A
#[2,]    3    5  ## B
#[3,]    5    7  ## C
#[4,]    7    9  ## D

## use `mapply`
mapply(function (lwr, upr) integrate(g, lower = lwr, upper = upr)$value,
       xx[-length(xx)], xx[-1])
#[1] 5.583333 5.000000 5.266679 2.025000
person Zheyuan Li    schedule 03.09.2018