Матрица с симплексными столбцами в stan - PullRequest
3 голосов
/ 01 октября 2019

Есть ли способ построить матрицу с симплексными столбцами в Stan? Модель, которую я хочу построить, похожа на следующую, где моя модель считается dirichlet-multinomial:

data {
  int g;
  int c;
  int<lower=0> counts[g, c];
}

parameters {
  simplex [g] p;
}

model {
  for (j in 1:c) {
    p ~ dirichlet(rep_vector(1.0, g));
    counts[, j] ~ multinomial(p);
  }
}

Однако я хотел бы использовать скрытую матрицу [g, c] для других слоев иерархической модели, аналогичнойк следующему:

parameters {
  // simplex_matrix would have columns which are each a simplex.
  simplex_matrix[g, c] p;
}
model {
  for (j in 1:c) {
    p[, j] ~ dirichlet(rep_vector(1.0, g));
    counts[, j] ~ multinomial(p[, j]);
  }
}

Если есть другой способ создать эту скрытую переменную, это, конечно, также было бы здорово! Я не очень хорошо знаком с тем, что Стэн реализовал только несколько иерархических моделей.

1 Ответ

4 голосов
/ 01 октября 2019

Чтобы ответить на вопросы, которые вы задали, вы можете объявить массив симплексов в блоке параметров программы Stan и использовать их для заполнения матрицы. Например,

parameters {
  simplex[g] p[c];
}
model {
  matrix[g, c] col_stochastic_matrix;
  for (i in 1:c) col_stochastic_matrix[,c] = p[c];
}

Однако на самом деле вам не нужно формировать стохастическую матрицу столбцов в приведенном вами примере, поскольку вы можете создать модель с многочленом-Дирихле, индексируя массив симплексов, таких как

data {
  int g;
  int c;
  int<lower=0> counts[g, c];
}
parameters {
  simplex [g] p[c];
}
model {
  for (j in 1:c) {
    p[j] ~ dirichlet(rep_vector(1.0, g));
    counts[, j] ~ multinomial(p[j]);
  }
}

Наконец, вам вообще не нужно объявлять массив симплексов, поскольку они могут быть интегрированы вне апостериорного распределения и восстановлены в блоке сгенерированных количеств программы Stan. См. wikipedia для деталей, но суть этого дается этой функцией Stan

functions {
  real DM_lpmf(int [] n, vector alpha) {
    int N = sum(n);
    real A = sum(alpha);
    return lgamma(A) - lgamma(N + A) 
           + sum(lgamma(to_vector(n) + alpha)
           - sum(lgamma(alpha));
  }
}
...