# MAS115 Class 4 code/solutions # section 1 While/repeat loops # A while loop set.seed(12) nobs <- 10 x <- runif(nobs) while((mean(x) < 0.49) || (mean(x) > 0.51)) { nobs <- 2 * nobs cat("Mean is", mean(x), "\n") cat("Need to try number = ", nobs, "\n") x <- runif(nobs) } nobs mean(x) # A loop which never ends while(TRUE) { cat("Trapped inside a never ending loop", "\n") } # A repeat loop i <- 0 repeat { i <- i+1 x <- rnorm(1) cat("Try", i, "realisation is", x, "\n") if(x > 1.2) break } i ### section 2 Lab class tasks (and homework) ##### ######### Cumulative sums ############ # While loop to count sum of exponential random variables sum <- 0 count <- 0 while(sum < 100) { sum <- sum + rexp(1, rate = 0.2) count <- count + 1 } count # While loop to count sum of exponential random variables sum <- 0 count <- 0 while(sum <= 100) { sum <- sum + rexp(1, rate = 0.2) count <- count + 1 } count # Putting it inside a for loop n <- rep(NA, 1000) for(i in 1:1000) { sum <- 0 count <- 0 while(sum < 100) { sum <- sum + rexp(1, rate = 0.2) count <- count + 1 } n[i] <- count } pdf("SumRV.pdf") hist(n, breaks = 30, xlab = "Number of RVs needed", main = "Histogram of number of Exponential RVs needed to sum to 100") dev.off() ######### Death... ############ ## Using a for loop ## y <- 100 z <- 40 N <- 10000 # Number of realisation to create C <- rep(NA, N) # Vector to store created values Q <- 100 h <- 50 for(i in 1:N) { # Sample values for u, sy and sz u<-rlnorm(1,2,0.1^0.5) sy<-rlnorm(1,10,0.2^0.5) sz<-rlnorm(1,5,0.05^0.5) # Store in C[i] value of conc. at site of interest C[i] <-Q/(2*pi*u*(sz*sy)^0.5)*exp(-0.5*(y^2/sy+(z-h)^2/sz)) } head(C) # Plot histogram hist(C) # Find 95th quantile quantile(C, 0.95) ## Using vectorisation N <- 10000 # Create values for u, sy and sz u <- rlnorm(N,2,0.1^0.5) sy <- rlnorm(N,10,0.2^0.5) sz <- rlnorm(N,5,0.05^0.5) # Create a C for each C<-100/(2*pi*u*(sz*sy)^0.5)*exp(-0.5*(100^2/sy+(40-50)^2/sz)) # Plot histogram hist(C) # Find 95th quantile quantile(C,0.95) library(MASS) pdf("Conc.pdf") truehist(C, xlab = "Concentration") dev.off() ######### Triangular numbers ############ target <- 7 k <- 0 T_k <- 0 # kth triangular number countSquares <- 0 # count of cases for T_k = perfect square sols <- matrix(NA,nrow=target,ncol=3, dimnames=list(NULL,c("k","T_k","sqroot"))) while(countSquares