R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1 Patched (2005-08-04), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> invisible(options(echo = TRUE))
> 
>  ##### GOFMC toy problem.  X is bivariate normal with mean vector
> 
>  mu <- c(1, 1)
> 
>  ##### and covariance matrix
> 
>  V <- matrix(c(1, 0.5, 0.5, 1), 2, 2)
> 
>  ##### What is the probability that X lies in the first quadrant?
> 
>  ##### Note (as with many toy MC problems) can do by numerical integration
> 
>  detV <- prod(eigen(V, symmetric = TRUE, only.values = TRUE)$values)
>  const <- 1 / (2 * pi * sqrt(detV))
>  invV <- solve(V)
> 
>  dens <- function(x) {
+      stopifnot(is.numeric(x))
+      stopifnot(length(x) == 2)
+      delta <- x - mu
+      const * exp(- as.numeric(t(delta) %*% invV %*% delta) / 2)
+  }
> 
>  ##### check
> 
>  library(adapt)
> 
>  adapt(2, lower = c(-10, -10), upper = c(10, 10), functn = dens)
      value      relerr      minpts      lenwrk       ifail 
   1.100264 0.009757585         477          73           0 
> 
>  adapt(2, lower = c(-20, -20), upper = c(20, 20), functn = dens)
      value      relerr      minpts      lenwrk       ifail 
   1.024290 0.008149177        1201         283           0 
> 
>  ##### pretty sloppy! (by default)
> 
>  adapt(2, lower = c(-10, -10), upper = c(10, 10), functn = dens, eps = 1e-5)
       value       relerr       minpts       lenwrk        ifail 
    1.000096 7.308158e-06         1877          563            0 
> 
>  adapt(2, lower = c(-20, -20), upper = c(20, 20), functn = dens, eps = 1e-5)
       value       relerr       minpts       lenwrk        ifail 
    1.000098 8.066745e-06         2421          563            0 
> 
>  ##### much better! (know and use the optional arguments)
> 
>  true.answer <- adapt(2, lower = c(0, 0), upper = c(10, 10), functn = dens,
+      eps = 1e-5)$value
>  print(true.answer)
[1] 0.7452255
> 
>  ##### to about 5 sig. figs.
> 
>  ##### now for MCMC
> 
>  library(MASS)
> 
>  set.seed(42)
> 
>  nboot <- 1000
> 
>  x <- mvrnorm(nboot, mu = mu, Sigma = V)
> 
>  mc.answer <- mean(apply(x > 0, 1, prod))
>  print(mc.answer)
[1] 0.738
> 
>  ##### MC standard error (if you don't do this, you're just fooling around)
> 
>  mc.error <- sqrt(mc.answer * (1 - mc.answer) / nboot)
> 
>  (mc.answer - true.answer) / mc.error
[1] -0.5196212
> 
>  ##### closer than we expect to be on average, try more
> 
>  nfoo <- 30
>  z <- as.double(nfoo)
>  for (i in 1:nfoo) {
+      x <- mvrnorm(nboot, mu = mu, Sigma = V)
+      mc.answer <- mean(apply(x > 0, 1, prod))
+      mc.error <- sqrt(mc.answer * (1 - mc.answer) / nboot)
+      z[i] <- (mc.answer - true.answer) / mc.error
+  }
> 
>  stem(z)

  The decimal point is at the |

  -2 | 4
  -1 | 
  -1 | 3
  -0 | 9988775
  -0 | 443332211111
   0 | 3
   0 | 5679
   1 | 134
   1 | 
   2 | 2

> 
>  ##### Looks pretty standard normal to me.
> 
>  ##### COMMENTS #####
>  ##
>  ## 1. Always do (at least some) MC standard errors
>  ##
>  ## 2. The way we did here only works for probabilities where the sum
>  ##    of indicator random variables is binomial.  In general you have
>  ##    to calculate as in "z" procedures.
>  ##
>  ## 3. The calculation the dumb way (the nfoo for loop) is almost never
>  ##    necessary.  We're statisticians.  We know n p (1 - p) and the CLT
>  ##    and don't have to calculate them by simulation.  In general, we know
>  ##    we can plug in the sample standard deviation for the population
>  ##    standard deviation and use the CLT.
>  ##
>  ##    In short: we can use INTERNAL (one run) rather than EXTERNAL (multi-run,
>  ##    a. k. a., bootstrap) estimates of standard error.
>  ##
>  ##    This "in short" will be a constant theme throughout the course.
>  ##
>  ##### END COMMENTS #####
> 
> proc.time()
[1]  2.48  0.05 10.27  0.00  0.00
>