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)      #畫曲線