Оцениваем модель-игрушку

Остальные материалы mcmc-курса

Создаём мир!

n.obs <- 300
chain.len <- 10^4

x <- rnorm(n.obs, mean = 42, sd = 6)

Готовим нужные функции

ln_prior <- function(a) {
    # log of prior distribution at the point a=(mu,sigma^2)
    result <- -a[1]^2/2/10^6 - a[2]/2 + 499 * log(a[2])
    return(result)
}

ln_conditional <- function(a) {
    # log of conditional distribution at the point a=(mu,sigma^2)
    result <- (-1) * n.obs/2 * log(a[2]) - sum((x - a[1])^2)/2/a[2]
    return(result)
}


alpha <- function(a, b) {
    # probability density of transition from state a to state b
    res <- ln_prior(b) + ln_conditional(b) - ln_prior(a) - ln_conditional(a)
    res <- exp(res)
    return(min(c(1, res)))
}

Создаём исходную цепь

go_q <- function(a) {
    # y_{1,t} = y_{1,t-1}+e_{1,t} y_{2,t} = y_{2,t-1}*exp(e_{2,t})
    return(c(a[1] + rnorm(1), a[2] * exp(rnorm(1))))
}

Создаём MH-испорченную цепь

go_p <- function(a) {
    proposal <- go_q(a)
    p.accept <- alpha(a, proposal)

    res <- a
    if (runif(1) < p.accept) 
        res <- proposal

    return(res)
}

Много шагов по изменённой цепи

gogogo <- function(a, n.steps = 1) {
    res <- data.frame(mu = rep(a[1], times = n.steps), sigma2 = rep(a[2], times = n.steps))
    for (i in 2:n.steps) res[i, ] <- go_p(as.numeric(res[i - 1, ]))
    return(res)
}

Смотрим на результат

library(ggplot2)
library(gridExtra)

## Loading required package: grid

h <- gogogo(c(1, 1), chain.len)
h$obs <- 1:chain.len

g1 <- ggplot(h, aes(x = obs, y = mu)) + geom_line()
g2 <- ggplot(h, aes(x = mu)) + geom_histogram(binwidth = 0.05) + coord_flip()
g <- arrangeGrob(g1, g2, nrow = 1)

## Warning: position_stack requires constant width: output may be incorrect

g