Advertisement
Guest User

Untitled

a guest
Aug 17th, 2017
64
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 0.59 KB | None | 0 0
  1. set.seed(123)
  2.  
  3. x <- rnorm(7)
  4.  
  5. dx <- density(x)
  6. plot(dx)
  7.  
  8. f <- approxfun(dx, yleft = 0, yright = 0)
  9.  
  10. integrate(f, lower = -6, upper = 6)
  11.  
  12. xx <- seq(-6, 6, by = 0.001)
  13. plot(xx, f(xx), type = "l")
  14. y <- sample(xx, 1e5, prob = f(xx), replace = TRUE)
  15.  
  16. hist(y, freq = FALSE, 100)
  17. lines(xx, f(xx), col = "red")
  18.  
  19.  
  20. R <- 1e5
  21. y <- numeric(R)
  22. tmp <- 0
  23.  
  24. for (i in 1:R) {
  25. prop <- tmp + rnorm(1)
  26. u <- runif(1)
  27. A <- dnorm(tmp - prop, log = TRUE) - dnorm(prop - tmp, log = TRUE)
  28. A <- A + log(f(prop)) - log(f(tmp))
  29. if (log(u) < A) {
  30. tmp <- prop
  31. }
  32. y[i] <- tmp
  33. }
  34.  
  35. lines(xx, f(xx), col = "blue", lty = 2)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement