Cómo elegantemente determinar el área de un ciclo de histéresis (dentro/fuera problema)?

Question

He medido dos parámetros (carbono orgánico Disuelto DOC = y, y de la descarga = x). Cuando estas dos variables se representan el uno contra el otro, tenemos un bucle de histéresis (véase el ejemplo de código y la imagen).

Ahora, para su posterior análisis, me gustaría para determinar el área de este bucle de histéresis. Me di cuenta de que esto se puede hacer utilizando el Monte Carlo dardos método. Este método dice que el área de una zona desconocida es proporcional al área de un conocido rectangular veces los golpes en el interior del campo (el lazo).

Mi problema ahora es, ¿cómo resolver el dentro / fuera problema mediante R. ¿Cómo puedo dibujar un rectángulo con un área conocida y cómo puedo excel los golpes al azar dentro y fuera del bucle de histéresis?

Por favor, tenga en cuenta, que estoy abierto a cualquier otro método...

Busqué en google, y busqué estadísticas varias webs, pero no podía encontrar una respuesta. Cualquier ayuda directa o enlaces a otros sitios web/blogs es muy apreciado.

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

Answer 1

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

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

Esto produce 0.606.

Answer 2

Una posibilidad sería este: a mí me parece que el bucle de histéresis debe ser convexo, ¿verdad? Por lo que se podría generar puntos de prueba y para cada punto de si es parte de el casco convexo de la unión de un conjunto de datos y el punto al azar - 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 componen su casco convexo de sí mismo.

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, el R manera de hacer las cosas que no uso mi for bucle, pero esto es fácil de entender y no demasiado lento. Puedo obtener una superficie estimada de 0.703.

EDIT: por supuesto, esto se basa en la presunción de la convexidad de la relación. Por ejemplo, no parece ser un nonconvex parte en la parte inferior derecha. En principio, se podría Monte-Carlo-calcular el área de un área de este tipo en la misma forma y restar este de la zona original de la estimación.