cmdstanR: вывод по параметру Бернулли

Question

Я построил простую модель, используя распределение Бернулли в R, используя cmdstanR.

Файл stan:

data {
  int<lower=0> N;
  int<lower=0, upper=1> obs_data[N];
}

parameters {
  real<lower=0, upper=1> lambda;
}

model {
  target += beta_lpdf(lambda | 1,1);
  for (n in 1:N) {
    target += bernoulli_logit_lpmf(obs_data[n] | lambda);
  }
}

Затем я создал 4 Бернулли др aws, с количеством образцов как 10, 100, 1000 и 10000. Я хотел заметить, что с увеличением количества точек данных неопределенность, связанная с параметром, уменьшается.

Код r выглядит следующим образом:

extract_lambda_draws <- function(mod, obs_data, iter = 1) {

  dl <- list(N = length(obs_data), obs_data = obs_data)
  print(paste("Model build iteration: ", iter))

  fit <- mod$sample(data = dl, num_chains = 4, num_cores = 4)

  print("Model build competed ...")
  draws <- fit$draws()[,,1] %>% as_tibble() 
  return(round(draws,3))
}

num_tosses <- c(10, 100, 1000, 10000)

results <- tibble()

m <- cmdstan_model("coin-flip.stan")

for (i in num_tosses) {
  coin_tosses <- sample(c(0,1), i, replace = T, prob = c(0.4, 0.6))
  d <- extract_lambda_draws(m, coin_tosses, i)
  d <- d %>% mutate(iter = i)
  results <- rbind(results, d)
}

results %>%
  pivot_longer(cols = c(ends_with("lambda")), names_to = "chains", values_to = "lambda" ) %>% 
  mutate(chains = gsub(".lambda", "", chains)) %>% 
  ggplot(aes(x = lambda)) + geom_density() + facet_wrap(iter~., nrow = 4, ncol = 5)

Я получаю следующее распределение плотности по параметру

Когда я изменяю вероятность для 0 и 1 на c (0,6, 0,4) , Я получаю следующее

У меня есть 2 вопроса:

1. Когда я создаю образцы из c (0,1) с вероятностью c (0,4, 0,6). Я ожидаю, что лямбда будет около 0,6, по крайней мере, для набора данных с 10000 образцов. Однако апостериорный режим составляет ~ 0,4.

2. Когда я создаю выборки из c (0,1) с вероятностью c (0,6, 0,4). Я ожидаю, что лямбда будет около 0,4, по крайней мере, для набора данных с 10000 образцов. Задний режим близок к 0.

Answer 1

Это потому, что вы используете lo git - распределение Бернулли.

Затем, в первой ситуации, задняя часть концентрируется примерно на:

> car::logit(0.6)
[1] 0.4054651

В вторая ситуация, одна имеет:

> car::logit(0.4)
[1] -0.4054651

Но ваше предыдущее распределение по lo git (p) ограничено диапазоном (0,1). Таким образом, апостериор также ограничен этим диапазоном, и затем он концентрируется на 0.

Я не знаю, существует ли функция для распределения Бернулли, параметризованная p в Стэне. Но вы могли бы сделать что-то подобное (я не уверен в синтаксисе):

parameters {
  real<lower=0, upper=1> p;
}
transformed_parameters {
  lambda = log(p/(1-p)) // not sure of the syntax here
}
model {
  target += beta_lpdf(p | 1,1);
  for (n in 1:N) {
    target += bernoulli_logit_lpmf(obs_data[n] | lambda);
  }
}