# Coding the normal PDF from scratch normal <- function(y, mu, sigma){ #y must be first! exp(-(y - mu) ^ 2 / (2 * sigma^2)) / (sigma * sqrt(2 * pi)) } # Checking that our function integrates to 1 integrate(normal, lower = -1000, upper = 1000, mu = 0, #put any extra arguments needed for your function here sigma = 3) # Simulating from our function values <- seq(-10, 10, .001) weights <- normal(values, mu = 0, sigma=1) draws <- sample(values, size = 10000, prob = weights, replace = T) # Plot the histogram of our simulated draws hist(draws, xlim = c(-4,4), prob = T, xlab = "y", main = "Histogram of simulations from our normal function", col = "red") curve(normal(x, mu = 0, sigma = 1), add = T, col = "blue", lwd = "2")