Comment déterminer avec élégance l'aire d'une boucle d'hystérésis (problème intérieur / extérieur)?

J'ai mesuré deux paramètres (carbone organique dissous DOC = y et décharge = x). Lorsque ces deux variables sont tracées l'une contre l'autre, nous obtenons une boucle d'hystérésis (voir l'exemple de code et l'image).

Maintenant, pour une analyse plus approfondie, je voudrais déterminer la zone de cette boucle hystérétique. J'ai compris que cela peut être fait en utilisant la méthode de fléchettes Monte Carlo. Cette méthode dit que l'aire d'une zone inconnue est proportionnelle à l'aire d'un rectangulaire connu multiplié par les coups dans le champ intérieur (la boucle).

Mon problème est maintenant, comment résoudre le problème intérieur / extérieur en utilisant R. Comment puis-je dessiner un rectangle avec une zone connue et comment puis-je exceller les coups aléatoires à l'intérieur et à l'extérieur de la boucle hystérétique?

Veuillez noter que je suis ouvert à toute autre méthode ...

J'ai recherché sur Google et fouillé divers sites statistiques, mais je n'ai pas trouvé de réponse. Toute aide directe ou tout lien vers d'autres sites Web / messages est grandement apprécié.

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

Une manière complètement différente serait de calculer directement l'aire de votre polygone:

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

Cela donne 0,606.

Une possibilité serait la suivante: il me semble que la boucle d'hystérésis devrait être convexe, non? Ainsi, on pourrait générer des points et tester pour chaque point s'il fait partie de la coque convexe de l'union de votre jeu de données et du point aléatoire - si oui, il se trouve en dehors du jeu de données d'origine. Pour accélérer les choses, nous pouvons travailler avec un sous-ensemble de l'ensemble de données d'origine, à savoir les points comprenant sa coque convexe elle-même.

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

Bien sûr, la manière R de faire les choses n'utiliserait pas ma forboucle, mais c'est facile à comprendre et pas trop lent. J'obtiens une superficie estimée à 0,703.

EDIT: bien sûr, cela repose sur la convexité présumée de la relation. Par exemple, il semble y avoir une partie non convexe en bas à droite. En principe, nous pourrions Monte-Carlo-estimer la surface d'une telle zone de la même manière et la soustraire de l'estimation de surface d'origine.