r - stan 中的多项式模型 - 如何拟合狄利克雷分布参数？

Question

我想学习如何在 stan 中使用 Dirichlet 分布。

我有一个表格，其中包含因子变量的六个级别中每个级别的观察总数:

     counts    n     factor_var 
    -------- ------ ------------
      3710    4386    level 1   
      252     4386    level 2   
      332     4386    level 3   
       59     4386    level 4   
       29     4386    level 5   
       4      4386    level 6  

我只想估计 dirichlet 参数，即 a 值来自给定水平的概率。我知道它可以在 R 中轻松完成 (gtools::rdirichlet(number_of_samples, df$counts))，但我最终的目标是基于该估计的层次模型。

---

这是我目前拥有的一个 stan 模型。该模型无法编译，因为分类分布的结果是 <1, n_levels> 中的离散数字。

    data {
      int<lower=1> n_levels;
      int counts[n_levels];
    }

    parameters {
      vector<lower=0>[n_levels] priors;
      simplex[n_levels] level_p;
    }

    model {
      level_p ~ dirichlet(priors);
      counts ~ categorical_logit(level_p);
    }

这是 R 代码:

    df <-
      data.frame(
        counts=c(3710, 252, 332, 59, 29, 4),
        n = c(4386, 4386, 4386, 4386, 4386, 4386),
        factor_var = factor(1:6, labels=paste0('level ', 1:6)))


    data<-list(
      n_levels=6,
      counts<-df$counts
    )

    fit1<-stan(
      file='model.stan')

运行时，编译后出现如下错误:

    SAMPLING FOR MODEL 'model' NOW (CHAIN 1).
    Rejecting initial value:
      Error evaluating the log probability at the initial value.
    Exception: categorical_logit_lpmf: categorical outcome out of support is 3710, but must be in the interval [1, 6]  (in 'model6fffc8dd594_proba_stan' at line 13)

    ...

    Initialization between (-2, 2) failed after 100 attempts. 
     Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
    [1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
    error occurred during calling the sampler; sampling not done

Answer 1

您需要使用多项式分布。

categorical_logit 假定对数几率参数，而不是单纯形。

手册中有如何拟合 Dirichlet 的示例，包括一些广义先验。