###################################################################### ### Simulating from the exponential distribution with lambda = 0.5 ### ### using the U(0,1)-distribution ### ###################################################################### xseq <- seq(-5, 10, len = 200) plot(xseq, pexp(xseq, 1/2), type = "l", xlab = "Exponential outcomes", ylab = "Uniform outcomes", cex.lab = 2) set.seed(123) us <- runif(5) arrows(rep(-5,5), us, qexp(us, 1/2), us, lty = 3, length = 0.2) arrows(qexp(us, 1/2), us, qexp(us, 1/2), 0, lty = 3, length = 0.2) ################################## ### Simulation: Area of plates ### ################################## set.seed(345) k = 10000 # Number of simulations X = rnorm(k, 2, 0.01) Y = rnorm(k, 3, 0.02) A = X*Y mean(A) var(A) mean(abs(A - 6) > 0.1) ######################################################################## ### Simulations: Mean of 10 exponential distributed random variables ### ######################################################################## set.seed(9876) # Number of simulations k <- 100000 # Simulate 10 exponentials with the 'right' mean k times sim_samples <- replicate(k, rexp(10, 1/26.08)) # Compute the mean of the 10 simulated observations k times sim_means <- apply(sim_samples, 2, mean) # Find relevant quantiles of the k simulated means quantile(sim_means, c(0.025, 0.975)) # Make histogram of simulated means hist(sim_means, col = "blue", nclass = 30, main = "", prob = TRUE, xlab = "Simulated means") ########################################################################## ### Simulations: Median of 10 exponential distributed random variables ### ########################################################################## set.seed(9876) # Number of simulations k <- 100000 # Simulate 10 exponentials with the 'right' mean k times sim_samples <- replicate(k, rexp(10, 1/26.08)) # Compute the median of the 10 simulated observations k times sim_medians <- apply(sim_samples, 2, median) # Find relevant quantiles of the k simulated medians quantile(sim_medians, c(0.025, 0.975)) # Make histogram of simulated medians hist(sim_medians, col = "blue", nclass = 30, main = "", prob = TRUE, xlab = "Simulated medians") ###################################################### ### Simulation: CI for Q3 of a normal distribution ### ###################################################### set.seed(9876) # Heights data x <- c(168, 161, 167, 179, 184, 166, 198, 187, 191, 179) n <- length(x) # Define a Q3-function Q3 <- function(x){ quantile(x, 0.75)} # Set number of simulations k <- 100000 # Simulate k samples of n = 10 normals with the 'right' mean and variance sim_samples <- replicate(k, rnorm(n, mean(x), sd(x))) # Compute the Q3 of the n = 10 simulated observations k times simQ3s <- apply(sim_samples, 2, Q3) # Find the two relevant quantiles of the k simulated Q3s quantile(simQ3s, c(0.005, 0.995)) ################################################################################# ### Simulations: Confidence interval for difference between exponential means ### ################################################################################# set.seed(9876) # Day 1 data x <- c(32.6, 1.6, 42.1, 29.2, 53.4, 79.3, 2.3 , 4.7, 13.6, 2.0) n1 <- length(x) # Day 2 data y <- c(9.6, 22.2, 52.5, 12.6, 33.0, 15.2, 76.6, 36.3, 110.2, 18.0, 62.4, 10.3) n2 <- length(y) # Set number of simulations: k <- 100000 # Simulate k samples of each n1 = 10 and n2 = 12 exponentials # with the 'right' means simX_samples <- replicate(k, rexp(n1, 1/mean(x))) simY_samples <- replicate(k, rexp(n2, 1/mean(y))) # Compute the difference between the simulated means k times sim_dif_means <- apply(simX_samples, 2, mean) - apply(simY_samples, 2, mean) # Find the relevant quantiles of the k simulated differences of means: quantile(sim_dif_means, c(0.025, 0.975)) ############################################## ### Example: Womens' cigarette consumption ### ############################################## # Data x1 <- c(8, 24, 7, 20, 6, 20, 13, 15, 11, 22, 15) x2 <- c(5, 11, 0, 15, 0, 20, 15, 19, 12, 0, 6) # Compute differences dif <- x1-x2 dif # Compute average difference mean(dif) # 95% CI for mean by non-parametric bootstrap k = 100000 sim_samples = replicate(k, sample(dif, replace = TRUE)) sim_means = apply(sim_samples, 2, mean) quantile(sim_means, c(0.025,0.975)) # 95% CI for median by non-parametric bootstrap k = 100000 sim_samples = replicate(k, sample(dif, replace = TRUE)) sim_medians = apply(sim_samples, 2, median) quantile(sim_medians, c(0.025,0.975)) ############################# ### Example: Tooth health ### ############################# # Reading in data x <- c(9, 10, 12, 6, 10, 8, 6, 20, 12) y <- c(14,15,19,12,13,13,16,14,9,12) # 95% CI for mean difference by non-parametric bootstrap k <- 100000 simx_samples <- replicate(k, sample(x, replace = TRUE)) simy_samples <- replicate(k, sample(y, replace = TRUE)) sim_mean_difs <- apply(simx_samples, 2, mean)- apply(simy_samples, 2, mean) quantile(sim_mean_difs, c(0.025,0.975)) # 99% CI for median difference by non-parametric bootstrap k <- 100000 simx_samples <- replicate(k, sample(x, replace = TRUE)) simy_samples <- replicate(k, sample(y, replace = TRUE)) sim_median_difs <- apply(simx_samples, 2, median)- apply(simy_samples, 2, median) quantile(sim_median_difs, c(0.005,0.995))