Thursday, 28 May 2015

Estimation of PI using simulation

## Estimating the value of PI using simulation

library(psych)
mydata <- runif(200000,0,1)
dim(mydata) <- c(100000,2)
colnames(mydata) <- c("u", "v")
mydata <- transform(mydata, ss = (mydata[,1]^2 + mydata[,2]^2))
head(mydata)
##           u          v         ss
## 1 0.3124052 0.26868775 0.16979011
## 2 0.3438789 0.57032157 0.44351941
## 3 0.8312577 0.87588807 1.45816935
## 4 0.9655283 0.09276772 0.94085070
## 5 0.4588156 0.31979795 0.31278246
## 6 0.2589525 0.10603420 0.07829967
headTail(mydata,hlength = 4,tlength = 4)
##           u    v   ss
## 1      0.31 0.27 0.17
## 2      0.34 0.57 0.44
## 3      0.83 0.88 1.46
## 4      0.97 0.09 0.94
## ...     ...  ...  ...
## 99997  0.81  0.5 0.91
## 99998  0.17  0.5 0.28
## 99999  0.49 0.32 0.35
## 100000 0.88 0.72  1.3
p <- sum(mydata[,3] < 1) / nrow(mydata) #proportion of (x,y) falling inside or on 
                                        #the circle
(PI <- 4*p) #approximate value of PI
## [1] 3.14744

No comments:

Post a Comment