### R commands for section 6.1 & 6.2
### G. Givens & J. Hoeting, Computational Statistics, 2nd edition, Wiley, 2013.

## generating from standard parameter families
# available R functions
# rbeta -- The Beta Distribution 
# rbinom -- The Binomial Distribution 
# rcauchy -- The Cauchy Distribution 
# rchisq -- The (non-central) Chi-Squared Distribution 
# rexp -- The Exponential Distribution 
# rf -- The F Distribution 
# rgamma -- The Gamma Distribution 
# rgeom -- The Geometric Distribution 
# rhyper -- The Hypergeometric Distribution 
# rlnorm -- The Log Normal Distribution 
# rlogis -- The Logistic Distribution 
# rmultinom -- The Multinomial Distribution 
# rnbinom -- The Negative Binomial Distribution 
# rnorm -- The Normal Distribution 
# rpois -- The Poisson Distribution 
# rsignrank -- Distribution of the Wilcoxon Signed Rank Statistic 
# rt -- The Student t Distribution 
# runif -- The Uniform Distribution 
# rweibull -- The Weibull Distribution 
# rwilcox -- Distribution of the Wilcoxon Rank Sum Statistic 

#########################################################
## R tips: how to use R helps?
# want to know how to use function "runif"
?runif
# want to know how to generate random numbers from normal distribution
help.search("normal distribution")
?Normal
#########################################################

## simulation examples
# uniform distribution: runif(n, min=0, max=1)
x <- runif(1000, min=-1, max=1)
hist(x, freq=F)
curve(dunif(x, min=-1, max=1), lwd=2, lty=2, add=T)
# normal distribution: rnorm(n, mean=0, sd=1)
x <- rnorm(1000, mean=1, sd=2)
hist(x, freq=F)
curve(dnorm(x, mean=1, sd=2), lwd=2, lty=2, add=T)

#########################################################
## R tips: how to print or save R graphs?
# if use menu operations: [1] click the R graph window; 
#                         [2] (menu)File --> print or File --> Save as
# if use R command lines: 
x <- rnorm(1000, mean=1, sd=2)
postscript("normal.ps",paper="special",width=8,height=7,horizontal=F)
par(mfrow=c(1,1))
hist(x, freq=F)
curve(dnorm(x, mean=1, sd=2), lwd=2, lty=2, add=T)
dev.off()
#########################################################

#########################################################
## R tips: how to count the computation time (in seconds)
tempt <- proc.time()
temp <- runif(10000000)
proc.time()-tempt
#########################################################

## examples
## [1] int_2^4 ln(x) dx
4*log(4)-2*log(2)-2  #truth: 2.158883
# use loops
n <- 100000
ans <- 0
for(i in 1:n) {
  ans <- ans + 2*log(runif(1, min=2, max=4));
}
ans/n
# use vector 
x <- runif(n, min=2, max=4)
2*mean(log(x))

## [2] int_{-inf}^{inf} e^{-x^2} dx
sqrt(pi)  # truth: 1.772454
f <- function(x) {
  sqrt(2*pi)*exp(-x^2/2);
}
n <- 100000
x <- rnorm(n)
mean(f(x))

## [3] int_0^1 int_0^1 exp(x+y) dx dy
(exp(1)-1)^2      # truth: 2.952492
n <- 100000
x <- runif(n)
y <- runif(n)
mean(exp(x+y))

## [4] E[X(X-1)], Var[X(X-1)], X ~ N(0,1)
n <- 100000
x <- rnorm(n)
mu <- mean(x*(x-1))              # E[X(X-1)]
s2 <- mean(x^2*(x-1)^2) - mu^2   # Var[X(X-1)]
mu
s2