Estimating the area under a curve using random points in R

[This article was first published on The Beginner Programmer, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Even though this post is not the first one on this blog, the code I used to estimate the area under the curve x^2 and posted below is one of my first attempts to code something “complicated” with R. It has been more than 5 months since I wrote it, it does look a bit clunky now. Indeed. However I believe it is fun to keep an eye on our own written code and see how it has evolved. Even though the code is lengthy, the final result looks nice.

# Estimating the area under the curve x^2 using random points
# Given function
parabola <- function(x)
{
return(x^2)
}
plot(parabola)
# Definite integral over the interval [a,b]
areaintegrparab <- function(a,b)
{
return(b^3/3-a^3/3)
}
# Simulation approach to the calculation of the area
# a is the lower bound of [a,b]
# b is the upper bound of [a,b]
# c is the number of random points to be used
# e is the maximum allowed error (in percentage points)
areamt <- function(a,b,c,e)
{
randx<- runif(c,a,b)
randy<- runif(c,min(parabola(a),parabola(b)),max(parabola(a),parabola(b)))
c <- c(0)
h = 1
n = 0
m = 1
o = 1
vecy = c(0)
vecx = c(0)
for(number in randy)
{
if(number< parabola(randx[h]))
{
vecy[m]<-number
vecx[o]<-randx[h]
n = n + 1
h = h + 1
m = m + 1
o = o + 1
}else
{
h = h + 1
}
}
# Points ratio
pr=n/length(randx)
# Estimated area
mca = n/length(randx)*abs(b-a)*abs(max(parabola(a),parabola(b))-min(parabola(a),parabola(b)))
# Exact area
aint = areaintegrparab(a,b)
# Error
change = (mca-aint)/aint*100
# Vector containing the summary of the entire process
data = c(pr,mca,aint,change)
# Let's put some labels
if(change>=0)
{
namechange = "% Overestimated"
}else
{
namechange = "% Underestimated"
}
names(data)=c("Points ratio","Estimated area","Real area",namechange)
# Check if error is compatible with the maximum error allowed
# If TRUE, data is printed and plotted, else the process is reapeted with more points
if(abs(change) < e)
{
print(data)#Estimating the area under the curve x^2 using random points
# Given function
parabola <- function(x)
{
return(x^2)
}
plot(parabola)
# Definite integral over the interval [a,b]
areaintegrparab <- function(a,b)
{
return(b^3/3-a^3/3)
}
# Simulation approach to the calculation of the area
# a is the lower bound of [a,b]
# b is the upper bound of [a,b]
# c is the number of random points to be used
# e is the maximum allowed error (in percentage points)
areamt <- function(a,b,c,e)
{
randx<- runif(c,a,b)
randy<- runif(c,min(parabola(a),parabola(b)),max(parabola(a),parabola(b)))
c <- c(0)
h = 1
n = 0
m = 1
o = 1
vecy = c(0)
vecx = c(0)
for(number in randy)
{
if(number< parabola(randx[h]))
{
vecy[m]<-number
vecx[o]<-randx[h]
n = n + 1
h = h + 1
m = m + 1
o = o + 1
}else
{
h = h + 1
}
}
# Points ratio
pr=n/length(randx)
# Estimated area
mca = n/length(randx)*abs(b-a)*abs(max(parabola(a),parabola(b))-min(parabola(a),parabola(b)))
# Exact area
aint = areaintegrparab(a,b)
# Error
change = (mca-aint)/aint*100
# Vector containing the summary of the entire process
data = c(pr,mca,aint,change)
# Let's put some labels
if(change>=0)
{
namechange = "% Overestimated"
}else
{
namechange = "% Underestimated"
}
names(data)=c("Points ratio","Estimated area","Real area",namechange)
# Check if error is compatible with the maximum error allowed
# If TRUE, data is printed and plotted, else the process is reapeted with more points
if(abs(change) < e)
{
print(data)
plot(parabola,xlim=c(a,b),ylim=c(parabola(a),parabola(b)),col="red",main="Estimating the area under the curve")
points(vecx,vecy,col="green")
return(n/length(randx)*abs(b-a)*abs(max(parabola(a),parabola(b))-min(parabola(a),parabola(b))))
}else
{
return(areamt(a,b,c+1000,e))
}
}
# Calculate the area under the curve between 0 and 2 with 30000 points
# and a maximum allowed error of 2%
areamt(0,3,30000,2)

 

Here is the result you should get:



To leave a comment for the author, please follow the link and comment on their blog: The Beginner Programmer.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)