# Univariate distributed data with multiple parameters ----------------------
# Parameters
set.seed(1231)
n <- 1e3
pars <- c(mean = 2.6, sd = 3)
# Generating data and writing the log likelihood function
D <- rnorm(n, pars[1], pars[2])
fun <- function(x) {
x <- log(dnorm(D, x[1], x[2]))
sum(x)
}
# Calling MCMC, but first, loading the coda R package for
# diagnostics
library(coda)
ans <- MCMC(
fun, initial = c(mu=1, sigma=1), nsteps = 2e3,
kernel = kernel_normal_reflective(scale = .1, ub = 10, lb = 0)
)
# Ploting the output
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1,2))
boxplot(as.matrix(ans),
main = expression("Posterior distribution of"~mu~and~sigma),
names = expression(mu, sigma), horizontal = TRUE,
col = blues9[c(4,9)],
sub = bquote(mu == .(pars[1])~", and"~sigma == .(pars[2]))
)
abline(v = pars, col = blues9[c(4,9)], lwd = 2, lty = 2)
plot(apply(as.matrix(ans), 1, fun), type = "l",
main = "LogLikelihood",
ylab = expression(L("{"~mu,sigma~"}"~"|"~D))
)
par(oldpar)
# In this example we estimate the parameter for a dataset with ----------------
# With 5,000 draws from a MVN() with parameters M and S.
# \donttest{
# Loading the required packages
library(mvtnorm)
library(coda)
# Parameters and data simulation
S <- cbind(c(.8, .2), c(.2, 1))
M <- c(0, 1)
set.seed(123)
D <- rmvnorm(5e3, mean = M, sigma = S)
# Function to pass to MCMC
fun <- function(pars) {
# Putting the parameters in a sensible way
m <- pars[1:2]
s <- cbind( c(pars[3], pars[4]), c(pars[4], pars[5]) )
# Computing the unnormalized log likelihood
sum(log(dmvnorm(D, m, s)))
}
# Calling MCMC
ans <- MCMC(
initial = c(mu0=5, mu1=5, s0=5, s01=0, s2=5),
fun,
kernel = kernel_normal_reflective(
lb = c(-10, -10, .01, -5, .01),
ub = 5,
scale = 0.01
),
nsteps = 1e3,
thin = 10,
burnin = 5e2
)
# Checking out the outcomes
plot(ans)
summary(ans)
# Multiple chains -----------------------------------------------------------
# As we want to run -fun- in multiple cores, we have to
# pass -D- explicitly (unless using Fork Clusters)
# just like specifying that we are calling a function from the
# -mvtnorm- package.
fun <- function(pars, D) {
# Putting the parameters in a sensible way
m <- pars[1:2]
s <- cbind( c(pars[3], pars[4]), c(pars[4], pars[5]) )
# Computing the unnormalized log likelihood
sum(log(mvtnorm::dmvnorm(D, m, s)))
}
# Two chains
ans <- MCMC(
initial = c(mu0=5, mu1=5, s0=5, s01=0, s2=5),
fun,
nchains = 2,
kernel = kernel_normal_reflective(
lb = c(-10, -10, .01, -5, .01),
ub = 5,
scale = 0.01
),
nsteps = 1e3,
thin = 10,
burnin = 5e2,
D = D
)
summary(ans)
# }
Run the code above in your browser using DataLab