Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- f_x <- function(x){
- 1/((2*pi)*(1 + ((x-0)/2)^2))
- }
- n <- 10000
- x <- rep(0, n)
- q <- function(x){
- rnorm(1, x, 2)
- }
- # initial value
- x[1] <- 1
- for(i in 1:(n-1)){
- # pick new sample
- xnew <- q(x[i])
- # acceptance condition
- u <- runif(1)
- if (u < f_x(xnew)/f_x(x[i])){
- x[i+1] <- xnew
- }
- else{
- x[i+1] <- x[i]
- }
- }
- hist(x, prob = TRUE)
- y <- seq(-40, 40, 0.01)
- lines(y, f_x(y), col = "red")
- # Monte Carlo estimator of E(X)
- mean(x)
- p <- seq(0.025, 0.975, 0.01)
- Qhat <- quantile(x,p)
- Q <- qcauchy(Qhat, location = 0, scale = 2)
- round(rbind(Qhat, Q), 3)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement