Как я могу использовать RStan без петель?

Question

Есть ли способ более эффективно выполнять следующие вычисления в RStan?

Я предоставил только минимальное количество необходимого кода:

parameters {
  real beta_0;
  real beta_1;
}     
model {
  vector [n] p_i = exp(beta_0 + beta_1*x)/[1 + exp(beta_0 + beta_1*x)];
  y ~ bernoulli(p_i);
  /* Likelihood:
  for(i in 1:n){
    p_i[i] = exp(beta_0 + beta_1*x[i])/(1 + exp(beta_0 + beta_1*x[i]));
    y[i] ~ bernoulli(p_i[i]);
  */}
// Prior:
  beta_0 ~ normal(m_beta_0, s_beta_0);
  beta_1 ~ normal(m_beta_1, s_beta_1);
}

Я получаю следующую ошибкусообщение: «Элементы выражения матрицы должны иметь тип row_vector, а элементы выражения вектора строки должны быть целыми или действительными, но найдены элементы типа vector».Если я использую цикл for (который закомментирован), код работает нормально, но я бы хотел ограничить использование циклов for в моем коде.В приведенном выше коде x является вектором длины n.

Другой пример:

parameters {
  real gamma1;
  real gamma2;
  real gamma3;
  real gamma4;
}
model {
// Likelihood:
  real lambda;
  real beta;
  real phi;
  for(i in 1:n){
    lambda = exp(gamma1)*x[n_length[i]]^gamma2;
    beta = exp(gamma3)*x[n_length[i]]^gamma4;
    phi = lambda^(-1/beta);
    y[i] ~ weibull(beta, phi);
  }
  //y ~ weibull(exp(gamma1)*x^gamma2, exp(gamma3)*x^gamma4); //cannot raise a vector to a power
// Prior:
  gamma1 ~ normal(m_gamma1, s_gamma1);
  gamma2 ~ normal(m_gamma2, s_gamma2);
  gamma3 ~ normal(m_gamma3, s_gamma3);
  gamma4 ~ normal(m_gamma4, s_gamma4);
}

Приведенный выше код работает, но закомментированный расчет вероятности не работает, так как я "не могу поднять вектор до степени" (но вы можете в R),Я хотел бы еще раз, чтобы меня не заставляли использовать для циклов.В приведенном выше коде n_length является вектором длины n.

Последний пример.Если я хочу нарисовать 10000 выборок из нормального распределения в R, я могу просто указать

rnorm(10000, mu, sigma)

Но в RStan мне придется использовать цикл for, например

parameters {
      real mu;
      real sigma;
    }
generated quantities {
  vector[n] x;
  for(i in 1:n) {
    x[i] = normal_rng(mu, sigma);
  }
}

Что я могу сделать, чтобы ускорить мои примеры RStan?

Answer 1

Эта строка кода:

vector [n] p_i = exp(beta_0 + beta_1*x)/[1 + exp(beta_0 + beta_1*x)];

недопустимый синтаксис в языке Stan, поскольку квадратные скобки используются только для индексации. Вместо этого это может быть

vector [n] p_i = exp(beta_0 + beta_1*x) ./ (1 + exp(beta_0 + beta_1*x));

, который использует оператор поэлементного деления, или, что еще лучше,

vector [n] p_i = inv_logit(beta_0 + beta_1*x);

В этом случае y ~ bernoulli(p_i); будет работать как вероятность. А еще лучше, просто сделай

y ~ bernoulli_logit(beta_0 + beta_1 * x);

и он сделает преобразование для вас численно стабильным. Вы также можете использовать bernoulli_logit_glm, что немного быстрее, особенно с большими наборами данных.

В Stan 2.19.x, я думаю, вы можете извлечь N значений из распределения вероятностей в блоке сгенерированных количеств. Но вы слишком беспокоитесь о for петлях. Программа Stan переносится на C ++, где циклы бывают быстрыми, и почти все функции в языке Stan, которые принимают векторные входные данные и производят векторные выходные данные, на самом деле включают тот же цикл в C ++, как если бы вы сделали цикл самостоятельно.