如何优雅地确定磁滞回线的面积（内部/外部问题）？

机器算法验证 r 数据可视化 几何学
2022-03-10 19:52:32

我测量了两个参数（溶解有机碳 DOC = y，排放量 = x）。当这两个变量相互绘制时，我们得到一个滞后回线（参见代码示例和图片）。

现在，为了进一步分析，我想确定这个滞后回路的面积。我发现这可以使用蒙特卡洛飞镖方法来完成。该方法表示未知区域的面积与已知矩形的面积乘以内场（循环）中的命中数成正比。

我现在的问题是，如何使用 R 解决内部/外部问题。如何绘制一个具有已知面积的矩形，以及如何在滞后循环内部和外部超越随机命中？

请注意，我对任何其他方法持开放态度......

我用谷歌搜索，搜索了各种统计网站，但找不到答案。非常感谢任何直接帮助或与其他网站/帖子的链接。

我的磁滞回线

Data <- read.table("http://dl.dropbox.com/u/2108381/DOC_Q_hystersis.txt", sep = ";",
               header = T)

head(Data)
plot(Data$Q, Data$DOC, type = "o", xlab = "Discharge (m3 s-1)", ylab = "DOC (mg C l-1)",
 main = "Hystersis loop of the C/Q relationship")
2个回答

一种完全不同的方法是直接计算多边形的面积：

library(geometry)
polyarea(x=Data$Q, y=Data$DOC)

这产生 0.606。

一种可能性是：在我看来，磁滞回线应该是凸的，对吧？因此，可以为每个点生成点并测试它是否是数据集和随机点联合的凸包的一部分 - 如果是，它位于原始数据集之外。为了加快速度，我们可以使用原始数据集的一个子集，即构成其凸包本身的点。

Data.subset <- Data[chull(Data$Q, Data$DOC),c("Q","DOC")]

x.min <- min(Data.subset$Q)
x.max <- max(Data.subset$Q)
y.min <- min(Data.subset$DOC)
y.max <- max(Data.subset$DOC)

n.sims <- 1000
random.points <- data.frame(Q=runif(n=n.sims,x.min,x.max),
  DOC=runif(n=n.sims,y.min,y.max))
hit <- rep(NA,n.sims)
for ( ii in 1:n.sims ) {
  hit[ii] <- !((nrow(Data.subset)+1) %in%
    chull(c(Data.subset$Q,random.points$Q[ii]),
      c(Data.subset$DOC,random.points$DOC[ii])))
}

points(random.points$Q[hit],random.points$DOC[hit],pch=21,bg="black",cex=0.6)
points(random.points$Q[!hit],random.points$DOC[!hit],pch=21,bg="red",col="red",cex=0.6)

estimated.area <- (y.max-y.min)*(x.max-x.min)*sum(hit)/n.sims

凸包

当然，R 的做事方式不会使用我的for循环，但这很容易理解，也不会太慢。我得到的估计面积为 0.703。

编辑：当然，这取决于假定的关系凸性。例如，右下角似乎有一个非凸部分。原则上，我们可以用相同的方式蒙特卡罗估计这样一个区域的面积，然后从原始面积估计中减去它。

其它你可能感兴趣的问题
上一篇如何转换时间序列数据，以便使用更简单的技术进行故障预测？ 下一篇使用 Spearman 的 rho 报告确定系数