<- 7^5
A <- 2^31-1
M
<- 100000
N <- rep (0, N)
rn 1] <- 10
rn[for (i in 2:length (rn))
{<- (A * rn[i-1] ) %% M
rn[i]
}
<- rn/(M-1)
nrn par(mfrow=c(1,3),mar=c(4,4,3,1))
plot (nrn[1:100])
acf (nrn)
hist (nrn)
<- 100000
n <- runif(n)
a par(mfrow=c(1,2),mar=c(4,4,3,1))
hist(a,xlab="Random Numbers",main="")
acf(a,main="")
# generate exponenail random numbers
#use method of inverse cdf to generate iid sample from exp(1)
<- function(n)
gen_exp
{#generate unif(0,) random numbers
<- runif(n)
u #transform the random numbers
-log(1-u)
}
<- gen_exp (10000)
b par(mfrow=c(1,3),mar=c(4,4,3,1))
hist (b)
acf (b)
## r built in generators
hist (rexp (10000))
<- function(n)
gen_normal
{#calculates size of random samples, which is greater than half of n
<- ceiling(n/2)
size_sample
<- sqrt(2*rexp(size_sample))
R <- runif(size_sample,0,2*pi)
theta
<- R*cos(theta)
X <- R*sin(theta)
Y
c(X,Y)[1:n]
}
<- gen_normal(10000)
normal_sample
par(mfrow=c(1,2),mar=c(4,4,1,1))
hist(normal_sample,main="")
qqnorm(normal_sample)
qqline(normal_sample)
<- rnorm (10000)
nsample2 hist (nsample2)
qqnorm (nsample2)
qqline(nsample2)
<- 200
n <- rgamma (n, shape = 2)
rn
par (mfrow = c(2,2))
plot (rn[1:100])
hist (rn)
<- replicate(10, cumsum(rgamma (n, shape = 2))/(1:n))
xbar.rep <- replicate(10, (cumsum(rgamma (n, shape = 2))/(1:n) - 2) / sqrt(2/(1:n)))
se.rep plot(xbar.rep[,1], main = "Demonstrating LLN", log = "x", type = "b", ylim = c(0, 6)); abline (h = 2)
for (i in 2:ncol(xbar.rep)) points(xbar.rep[,i], col = i, pch = i, type = "b")
plot(se.rep[,1], ylim = c(-3,3), type="b",log = "x", main="Demonstrating CLT")
for (i in 2:ncol(se.rep)) points(se.rep[,i], col = i, pch = i, type = "b")
abline (h = c(0, -2,2), lty = c(1,2,2))
#### an application of monte carlo method in estimating pi
# n is the number of samples drawn uniformly from the rectangle (-1,1) * (-1,1)
# an estimate of pi is returned
<- function(n)
pi_est_mc
{#X and Y are independent, each with marginal distribution unif(-1,1)
<- runif(n,-1,1)
X <- runif(n,-1,1)
Y
<- 4 * (X^2 + Y^2 <= 1)
Z <- mean (Z)
mu <- 1.96 * sd (Z) /sqrt (n)
error list (pi.est = mu, error.95perc = error, ci.95perc = mu + c(-error, error))
}
pi_est_mc (100)
## $pi.est
## [1] 3.12
##
## $error.95perc
## [1] 0.3264052
##
## $ci.95perc
## [1] 2.793595 3.446405
pi_est_mc (10000)
## $pi.est
## [1] 3.1344
##
## $error.95perc
## [1] 0.03228595
##
## $ci.95perc
## [1] 3.102114 3.166686
pi_est_mc (100000)
## $pi.est
## [1] 3.13892
##
## $error.95perc
## [1] 0.0101899
##
## $ci.95perc
## [1] 3.12873 3.14911
pi_est_mc (10000000)
## $pi.est
## [1] 3.142116
##
## $error.95perc
## [1] 0.00101761
##
## $ci.95perc
## [1] 3.141099 3.143134