## (x, y) data
y <- c(-1, 5, 2, 3.5, 1, 5.5, -2, 3, -1)
x <- 1:length(y)
Аналитический метод
Ваше желаемое вычисление может быть сделано в два шага:
- кусочный интеграл на каждом линейном сегменте. Нет проблем, если вы хотите интегрировать только нулевую пропорцию (см. Подробности ниже);
- агрегированные по частям результаты соответственно для областей A , B , C и D (подробности см. Ниже).
шаг 1: кусочный интеграл для пропорции выше нуля
Если есть данные n
(x, y), то будет (n - 1)
сегментов. Обозначим (xl, yl)
как левую точку сегмента, а (xr, yr)
как правую точку.
- если
(yl < 0) && (yr < 0)
, весь сегмент ниже нуля;
- если
(yl > 0) && (yr > 0)
, весь сегмент выше нуля;
- если
(yl < 0) && (yr > 0)
, сегмент увеличивается, пересекая ноль;
- , если
(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)
Поскольку вы хотите интегрировать пропорцию выше нуля каждого сегмента, мы имеем для каждого случая:
- интеграл равен 0;
- площадь трапеции и интеграл
(yl + yr) * (xr - xl) / 2
;
- область представляет собой треугольник между
(xm, 0)
и (xr, yr)
; интеграл составляет yr * (xr - xm) / 2
;
- область представляет собой треугольник между
(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