#
# Calculating areas of complicated regions by Monte Carlo method
#

region <- function(x,y) {
    x^2 + y^2 < 1 && 9*(x-1/5)^2 + y^2/4 - 5*x*y < 1;
}

monte.carlo.area <- function(fun = region, npoints = 100, xmin = -1, ymin = -1, xmax = 1, ymax = 1, nruns = 5) {
  scores <- NULL;
  plot.new()                            #Erase graphics
  ## Setup ranges in the plot                                        
  plot.window(xlim = c(xmin, xmax), ylim = c(ymin, ymax))
  axis(1);   axis(2);   axis(3);   axis(4);
  for( i in 1:nruns) {
    x <- runif(npoints, min = xmin, max = xmax);
    y <- runif(npoints, min = ymin, max = ymax);
    s <- mapply(fun, x, y);
    points(x[s], y[s]);                  #Add the plot
    area = (length(s[s]) / npoints) * (xmax-xmin) / (ymax-ymin)
    scores = c(scores, area)
  };
  result <- c(mean(scores), sd(scores));
  names(result) <- c("mean", "std");
  result;
}