cmdstanR：对伯努利参数的推断

Question

我使用 cmdstanR 在 R 中使用伯努利分布构建了一个简单的模型。

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 次伯努利图，样本数为 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) 时，我得到以下结果

我有两个问题：

1. 当我从 c(0,1) 以概率 c(0.4, 0.6) 创建样本时。我预计 lambda 大约为 0.6，至少对于具有 10,000 个样本的数据集。然而后验模式是~0.4。

2. 当我从 c(0,1) 以概率 c(0.6, 0.4) 创建样本时。我预计 lambda 大约为 0.4，至少对于具有 10,000 个样本的数据集。后模接近于0。

Answer 1

那是因为您使用 logit-Bernoulli 分布。

那么，在第一种情况下，后验集中于：

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

在第二种情况下，有：

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

但您在 logit(p) 上的先前分布仅限于 (0,1) 范围。所以后验也限制在这个范围内，然后集中在0。

我不知道在 Stan 中是否存在由 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);
  }
}