¿Cómo determinar con elegancia el área de un ciclo de histéresis (problema interno / externo)?

He medido dos parámetros (carbono orgánico disuelto DOC = y, y descarga = x). Cuando estas dos variables se trazan una contra la otra, obtenemos un ciclo de histéresis (vea el ejemplo de código y la imagen).

Ahora, para un análisis más detallado, me gustaría determinar el área de este circuito histérico. Descubrí que esto se puede hacer usando el método de lanzamiento de Monte Carlo. Este método dice que el área de un área desconocida es proporcional al área de un rectángulo conocido multiplicado por los golpes en el campo interior (el bucle).

Mi problema ahora es cómo resolver el problema interno / externo con R. ¿Cómo puedo dibujar un rectángulo con un área conocida y cómo puedo sobresalir los golpes aleatorios dentro y fuera del circuito histérico?

Tenga en cuenta que estoy abierto a cualquier otro método ...

Busqué en Google y busqué en varios sitios estadísticos, pero no pude encontrar una respuesta. Cualquier ayuda directa o enlace a otros sitios web / publicaciones es muy apreciada.

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

Una forma completamente diferente sería calcular directamente el área de su polígono:

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

Esto produce 0,606.

Una posibilidad sería esta: me parece que el ciclo de histéresis debería ser convexo, ¿verdad? Por lo tanto, se podrían generar puntos y probar para cada punto si es parte del casco convexo de la unión de su conjunto de datos y el punto aleatorio; en caso afirmativo, se encuentra fuera del conjunto de datos original. Para acelerar las cosas, podemos trabajar con un subconjunto del conjunto de datos original, es decir, los puntos que comprenden su propio 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

Por supuesto, la forma R de hacer las cosas no usaría mi forciclo, pero esto es fácil de entender y no demasiado lento. Obtengo un área estimada de 0.703.

EDITAR: por supuesto, esto se basa en la supuesta convexidad de la relación. Por ejemplo, parece haber una parte no convexa en la esquina inferior derecha. En principio, podríamos estimar Montecarlo el área de tal área de la misma manera y restar esto de la estimación del área original.