Como determinar elegantemente a área de um ciclo de histerese (problema interno / externo)?

8

Eu medi dois parâmetros (carbono orgânico dissolvido DOC = y, e descarga = x). Quando essas duas variáveis ​​são plotadas uma contra a outra, obtemos um loop de histerese (veja o exemplo de código e a figura).

Agora, para uma análise mais aprofundada, eu gostaria de determinar a área desse loop histérico. Eu descobri que isso pode ser feito usando o método de dardo em Monte Carlo. Este método diz que a área de uma área desconhecida é proporcional à área de um retângulo conhecido vezes os acertos no campo interno (o loop).

Meu problema agora é: como resolver o problema interno / externo usando R. Como posso desenhar um retangular com uma área conhecida e como posso destacar os acertos aleatórios dentro e fora do loop histerético?

Observe que estou aberto a qualquer outro método ...

Pesquisei em vários sites estatísticos, mas não consegui encontrar uma resposta. Qualquer ajuda direta ou ligação a outros sites / postagens é muito apreciada.

Meu ciclo de histerese

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")
Strohmi
fonte

Respostas:

6

Uma maneira completamente diferente seria calcular diretamente a área do seu polígono:

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

Isso produz 0,606.

Stephan Kolassa
fonte
+1. Além disso, as áreas (diferentemente dos comprimentos, por exemplo) são notavelmente robustas a erros independentes nas coordenadas dos vértices. Observe que esta solução não pressupõe muito sobre o polígono; em particular, respeita sua falta (manifesta) de convexidade.
whuber
1

Uma possibilidade seria esta: parece-me que o ciclo de histerese deve ser convexo, certo? Assim, pode-se gerar pontos e testar para cada ponto se faz parte do casco convexo da união do seu conjunto de dados e o ponto aleatório - se sim, está fora do conjunto de dados original. Para acelerar as coisas, podemos trabalhar com um subconjunto do conjunto de dados original, ou seja, os pontos que compõem seu próprio casco convexo.

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

casco convexo

Obviamente, a maneira R de fazer as coisas não usaria meu forloop, mas isso é fácil de entender e não é muito lento. Recebo uma área estimada de 0,703.

EDIT: é claro, isso depende da suposta convexidade do relacionamento. Por exemplo, parece haver uma parte não-convexa no canto inferior direito. Em princípio, poderíamos estimar Monte-Carlo da área dessa área da mesma maneira e subtraí-la da estimativa original da área.

Stephan Kolassa
fonte
Seria muito mais rápido e preciso apenas rasterizar o polígono e estimar sua área contando as células que ele contém.
whuber