Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- set.seed(123)
- x <- rnorm(7)
- dx <- density(x)
- plot(dx)
- f <- approxfun(dx, yleft = 0, yright = 0)
- integrate(f, lower = -6, upper = 6)
- xx <- seq(-6, 6, by = 0.001)
- plot(xx, f(xx), type = "l")
- y <- sample(xx, 1e5, prob = f(xx), replace = TRUE)
- hist(y, freq = FALSE, 100)
- lines(xx, f(xx), col = "red")
- R <- 1e5
- y <- numeric(R)
- tmp <- 0
- for (i in 1:R) {
- prop <- tmp + rnorm(1)
- u <- runif(1)
- A <- dnorm(tmp - prop, log = TRUE) - dnorm(prop - tmp, log = TRUE)
- A <- A + log(f(prop)) - log(f(tmp))
- if (log(u) < A) {
- tmp <- prop
- }
- y[i] <- tmp
- }
- lines(xx, f(xx), col = "blue", lty = 2)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement