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) 时，我得到以下结果

我有两个问题：

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

Answer 1

那是因为您使用了 logit-伯努利分布。

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

> 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);
  }
}

cmdstanR：对伯努利参数的推断

cmdstanR: inference on a bernouli parameter

r

stan

rstan