1

I'm attempting to calculate the area between two plots using R's integrate function. I have two forecast curves that I'm able to place on the same plot and also shade in the area between the two curves, for visualisation: 

x1 = seq(1,200,1)
y1 = # 200 numbers from 0.02 - 1.000
y2 = # 200 numbers from 0.00 - 0.95

plot(x1,y1,type="l",bty="L",xlab="Forecast Days",ylab="Curve Actuals")
points(x1,y2,type="l",col="red")

polygon(c(x1,rev(x1)),c(y2,rev(y1)),col="skyblue")

Following the example here https://stat.ethz.ch/pipermail/r-help/2010-September/251756.html I've attempted to execute the same code for my data to calculate the are between the two curves. As stated in the example "the area between two curves is the same as the integral of the difference between these two curves (resp. its absolute values)":

f1 = approxfun(x1, y1-y2)     # piecewise linear function
f2 = function(x) abs(f1(x))   # take the positive value

integrate(f2, 1, 200)

However I get the following error:

Error in integrate(f2, 1, 200) : maximum number of subdivisions reached

Appreciate any clarity that can be shed on this. Thanks!

Qaribbean
  • 178
  • 2
  • 3
  • 17
  • 1
    You can increase the number of `subdivisions`, but then you are likely to get a roundoff error. Instead I'd suggest to calculate the integral without approximations. It's not that difficult using some geometry. – Roland Oct 17 '16 at 11:02

2 Answers2

1

Increasing the number of subdivisions, as suggested in an earlier comment by @Roland, executes correctly. It does produce an absolute error, but very minute.

f1 = approxfun(x1, y1-y2)     # piecewise linear function
f2 = function(x) abs(f1(x))   # take the positive value

area1 = integrate(f2, 1, 200, subdivisions = 300)
> area1
9.327692 with absolute error < 8.1e-05
Qaribbean
  • 178
  • 2
  • 3
  • 17
0

Is approxfun() returning a function? You need to have the difference of f1 and f2 saved as a function, not a set a point differences. For example

f1 <- function(x) { 2 * x - 1}
f2 <- function(x) { x^2 - 3 * x}
abs_dif <- function(x) { abs( f1(x) - f2(x) ) } 
integrate(abs_dif, -1, 1)

Try expressing the f1 and f2 as functions then try the second two lines of code.

BonStats
  • 358
  • 2
  • 9
  • Thanks for responding. In the example I included from the R Help forum and where I got the original code, it was executed exactly as I have written mine and the data itself is not too dissimilar to my own (the seq used on the x axis is smaller, the y-values similar to my own). When running the example used there, I get the exact same answer, so I'm not sure why I would need to make the changes to mine that you suggest. – Qaribbean Oct 17 '16 at 14:30