В настоящее время я пробираюсь через превосходную книгу Ричарда Мак-Элирея Статистическое переосмысление .Я нахожусь в главе 8. У меня есть 2 источника проблем (объяснено ниже).Обратите внимание, что для этого вам потребуется установить пакет stan.
Сначала создайте данные для моделирования и просмотрите их.Эта простая регрессия имеет перехват 10 и наклон 2.
df <- data.frame(x = 1:100,
noise = rnorm(100, 0, 20))
df$y <- 10 + df$x*2 + df$noise
ggplot(df, aes(x, y)) + geom_point() + geom_smooth(method = "lm")
Теперь для модели.Функция map2stan является оболочкой для stan.Модель ниже
library(rethinking)
mod <- map2stan(
alist(
y ~ dnorm(mu, sigma),
mu <- a + b*x,
a ~ dnorm(0, 100),
b ~ dnorm(0, 10),
sigma ~ dunif(0,10)
), data = df, chains = 2, iter = 4000, warmup = 1000)
Убедитесь, что модель вернула правильные параметры
precis(mod)
# Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
# a 7.68 2.00 4.41 10.79 2612 1
# b 1.94 0.03 1.88 1.99 2543 1
# sigma 9.96 0.04 9.90 10.00 3980 1
Задача 1
Теперь в соответствии с книгой, если я запускаю
plot(mod)
Я должен получить прекрасный график трассировки для каждой из переменных, построив обе цепи и периода прогрева.Вместо этого я получаю ошибку
Error in as.double(y) :
cannot coerce type 'S4' to vector of type 'double'
Вопрос 1: Что я делаю не так?Почему график не работает?
Задача 2
Если мы напечатаем mod объект
mod
, он скажет нам
map2stan model fit
6000 samples from 2 chains
Мы можем обойти вышеуказанную проблему, извлекая и нанося образцы из модели map2stan (при этом извлекается только параметр перехвата a.
dfSamp1 <- data.frame(chain1 = extract.samples(mod, n=4e3)$a,
row = 1:4e3)
ggplot(dfSamp1, aes(x = row, y = chain1, group = 1)) + geom_line(colour = "skyblue")
Но не хватает нескольких цепочек. Отбор образцов из модели
dfSamp2 <- data.frame(chain1 = extract.samples(mod, n=4e3)$a,
chain2 = extract.samples(mod, n=4e3)$a,
row = 1:4e3)
Но выборки в двух цепях абсолютно одинаковы.
ggplot(dfSamp2) +
geom_line(aes(x = row, y = chain1, group = 1), colour = "skyblue") +
geom_line(aes(x = row, y = chain2, group = 1), colour = "black")
На графике трассировки должны быть видны два цвета, но только 1, потому что вторая черная линия идентична первой (синей).
Что мы можем проверить, если мы посмотрим на фрейм данных
dfSamp2[1:10,]
chain1 chain2 row
1 5.738354 5.738354 1
2 4.934791 4.934791 2
3 5.701355 5.701355 3
4 9.701134 9.701134 4
5 6.123426 6.123426 5
6 11.769974 11.769974 6
7 7.675946 7.675946 7
8 8.858765 8.858765 8
9 6.279705 6.279705 9
10 5.805226 5.805226 10
Точно то же самое: не то, что вы получите от двух отдельных случайных прогулок.
Вопрос 2: Итак, как мне извлечь 2 цепи из объекта map2stan или извлечь другие цепочки?оценки ples с использованием функции extract.samples()?