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.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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.