R语言--蒙特卡洛计算定积分

　　参考上一篇蒙特卡洛计算圆周率

rm(list = ls())
x <- seq(0,1,0.001)
y <- x^2
d <- data.frame(x,y)
ggplot(d,aes(x,y))+geom_area(fill='brown1')

#求定积分从0到1 ，函数为x^2 
#求随机点在点（0，0）和点（1，1）这个正方形中，y< x^2 的随机点比例
#积分得：1/3*x^3当x在[0,1]的值==1/3
origin <- c(0,0)
distance <- function(a){
  a[2]-a[1]^2 
# 定义函数时考虑，用apply时，若按行，先取出了第一行，所以不能用a[,2]
}  #数据框只有一行时，不能选某一列
n <- 100000
A <- matrix(runif(2*n,0,1),ncol = 2,byrow = T)
b <- apply(A, 1, distance) 
mean(b<0)

par(bg='beige')
plot(A,col='azure3',xlab='',ylab='',main='MC',asp=1) #所有点
points(A[b<0,],col='aquamarine3')#积分下点
abline(h=0,col='goldenrod4',lty='dotdash',lwd=3) #画上下左右的直线
abline(h=1,col='goldenrod4',lty='dotdash',lwd=3)
abline(v=1,col='goldenrod4',lty='dotdash',lwd=3)
abline(v=0,col='goldenrod4',lty='dotdash',lwd=3)
lines(x,y,col='brown1',lwd=2)      #画曲线