Estimating the area under a curve using random points in R

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)

view raw area_under_curve_random_points.R hosted with ❤ by GitHub

 

Here is the result you should get: