8
$\begingroup$

I have measured two parameters (Dissolved organic carbon DOC = y, and discharge = x). When these two variables are plotted against each other, we get a hysteresis loop (see code example and picture).

Now, for further analysis, I would like to determine the area of this hysteretic loop. I figured out that this can be done using the Monte Carlo darting method. This method says that the area of an unknown area is proportional to the area of a known rectangular times the hits in the inside field (the loop).

My problem now is, how to solve the inside / outside problem using R. How can I draw a rectangular with a known area and how can I excel the random hits inside and outside the hysteretic loop?

Please note, that I am open to any other method...

I googled, and searched various statistical sites but could not find an answer. Any direct help or linkage to other websites/posts is greatly appreciated.

My hysteresis loop

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")
$\endgroup$
6
$\begingroup$

A completely different way would be to directly calculate the area of your polygon:

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

This yields 0.606.

$\endgroup$
  • $\begingroup$ +1. Moreover, areas (unlike lengths, for instance) are remarkably robust to independent errors in the vertex coordinates. Note that this solution does not assume much about the polygon; in particular, it respects its (manifest) lack of convexity. $\endgroup$ – whuber Nov 27 '12 at 17:05
1
$\begingroup$

One possibility would be this: it looks to me like the hysteresis loop should be convex, right? So one could generate points and test for each point whether it is part of the convex hull of the union of your dataset and the random point - if yes, it lies outside the original dataset. To speed up things, we can work with a subset of the original dataset, namely the points comprising its convex hull itself.

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

convex hull

Of course, the R way of doing things would not use my for loop, but this is easy to understand and not too slow. I get an estimated area of 0.703.

EDIT: of course, this relies on the presumed convexity of the relationship. For instance, there appears to be a nonconvex part in the lower right. In principle, we could Monte-Carlo-estimate the area of such an area in the same way and subtract this from the original area estimate.

$\endgroup$
  • $\begingroup$ It would be much faster and more accurate just to rasterize the polygon and estimate its area by counting the cells it contains. $\endgroup$ – whuber Nov 27 '12 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.