#buffon's needle buf<-function(n) { ttt <- NULL ttt[1] <- 0 x <- runif(n) th <- runif(n, 0, pi) st <- sin(th) for(i in 1:n) { if(st[i] > x[i]) ttt[i + 1] <- ttt[i] + 1 else ttt[i + 1] <- ttt[i] } ttt if(ttt[n + 1] > 0) plot((0:n)[ttt > 0], 2*(0:n)[ttt > 0]/ttt[ttt > 0], type = "l", xlab = "number of simulations", ylab = "proportion of hits") else print("no successes") abline(pi, 0) } #Simulation by inversion #------------------------- "dsim"<- function(n, inv.df) { u <- runif(n) inv.df(u) } --------------------- "inv.df.exp" <- function(x, lam = 1) { -(log(1-x))/lam} --------------------------- #Rejection sampling --------------------- "rej.sim"<- function(n) { r <- rep(-99,n) for(i in 1:n) { t <- -1 while(t < 0) { x <- rexp(1, 1) y <- runif(1, 0, exp( - x)) if(x > 1) t <- - y else t <- x^2 * exp( - x) - y } r[i] <- x } r } ------------------ "rej.sim"<-function(n) { for(i in 1:n) { t <- -1 while(t < 0) { x <- rexp(1, 1) y <- runif(1, 0, exp( - x)) if(x > 1) t <- - y else t <- x^2 * exp( - x) - y } r[i] <- x } r } ------------------ #simulation of a standard Cauchy ----------------- ru.sim<-function(n) { r<-NULL for(i in 1:n) { t<-1 while(t>0) { u<-runif(1,0,1) v<-runif(1,-1,1) t<-u^2 + v^2 -1 } r[i] <-u/v } r } #Importance sampling example #----------------------------- "i.s"<- function(n) { x <- 2/runif(n) psi <- x^2/(2 * pi * (1 + x^2)) mean(psi) } ---------------------