查找多边形区域：垂直和水平几何约束下的积分

Question

我正在尝试计算零以上和曲线以下的面积。我的曲线具有离散的 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下图

我现在正在为两件事而苦恼：

我使用 which(diff(sign(diff(y))) == 2) + 1 确定了局部最小值的位置索引，但这并没有给我 C[= 的 x 上限值D 的 50=] 或更低的 x 值。如何获得曲线与零相交的那些点？

我想如果我能得到 1) 零以上的局部最小值的正确列表，2) 零处的那些交叉点，3) 零以上的局部最大值的正确列表，我知道所有A、B、C 和 D 的边界点, 因此可以计算它们的面积。但这在 R 中编写代码似乎并不简单。这真的是解决我的问题的最简单方法，还是有更好的方法？

Answer 1

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

分析方法

您需要的计算可以分两步完成：

每个线性段上的分段积分。如果你只想整合零以上的比例，没有困难（详见下文）；
为 A、B、C 和 C 区域适当聚合分段结果D（详见下文）。

第一步：零以上比例的分段积分

如果有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

查找多边形区域：垂直和水平几何约束下的积分

Find areas of polygons: integration under vertical and horizontal geometric constraints

r

polygon

area

integral

numerical-integration

分析方法

数值方法