Advertisement
Guest User

Untitled

a guest
Feb 20th, 2020
97
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 0.58 KB | None | 0 0
  1. f_x <- function(x){
  2. 1/((2*pi)*(1 + ((x-0)/2)^2))
  3. }
  4.  
  5. n <- 10000
  6. x <- rep(0, n)
  7.  
  8.  
  9. q <- function(x){
  10. rnorm(1, x, 2)
  11. }
  12.  
  13. # initial value
  14. x[1] <- 1
  15.  
  16.  
  17. for(i in 1:(n-1)){
  18. # pick new sample
  19. xnew <- q(x[i])
  20. # acceptance condition
  21. u <- runif(1)
  22. if (u < f_x(xnew)/f_x(x[i])){
  23. x[i+1] <- xnew
  24. }
  25. else{
  26. x[i+1] <- x[i]
  27. }
  28. }
  29.  
  30. hist(x, prob = TRUE)
  31. y <- seq(-40, 40, 0.01)
  32. lines(y, f_x(y), col = "red")
  33.  
  34. # Monte Carlo estimator of E(X)
  35. mean(x)
  36.  
  37. p <- seq(0.025, 0.975, 0.01)
  38. Qhat <- quantile(x,p)
  39. Q <- qcauchy(Qhat, location = 0, scale = 2)
  40.  
  41. round(rbind(Qhat, Q), 3)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement