У меня есть набор переменных временного ряда различной длины и Я пытаюсь оценить иерархическую модель авторегрессии с использованием языка СТ (используя rstan
). Я изучил основы Стэна, но не уверен, как я должен выразить набор векторов с различной длиной.Как определить набор векторов с различной длиной в модели Стэна
Я закончил тем, что увеличил более короткие векторы, заполнив произвольные числа, чтобы все векторы имели одинаковую длину и игнорировали эти заполненные значения в разделе модели.
Хотя это, кажется, работает, я хотел бы знать, есть ли лучшие способы определения этого типа данных. Любое предложение приветствуется!
Я нахожусь на Ubuntu 14.04
с R version 3.2.3
и rstan_2.9.0
.
Данные
Мои данные выглядят следующим образом. Следует отметить, что существуют три вектора в y
, каждый из которых имеет длину 3, 5 и 10. Я заполнил короткие векторы с 100
так, что все векторы имеют размер 10.
$N
[1] 3 5 10
$K
[1] 3
$maxN
[1] 10
$y
$y[[1]]
[1] -0.626453811 -0.069571811 -0.004434404 100.000000000 100.000000000 100.000000000 100.000000000
[8] 100.000000000 100.000000000 100.000000000
$y[[2]]
[1] 1.5952808 0.9305912 0.4832488 0.3903673 0.3690161 100.0000000 100.0000000 100.0000000 100.0000000
[10] 100.0000000
$y[[3]]
[1] 0.5757814 0.5876644 0.7800761 0.8410528 0.7948234 0.5938711 0.7469771 0.7677860 0.7893884 0.9048332
Стан код
Мои Стан код выглядит следующим образом:
data {
int<lower=0> K; // number of series
int<lower=0> N[K]; // length of each series
int<lower=0> maxN; // maximum length of series
vector[maxN] y[K]; // time series (define as the same length)
}
parameters {
// For simplicity, assume only beta varies across series
real alpha;
real beta[K];
real<lower=0> sigma;
real beta0;
real<lower=0> beta_sig;
}
model {
for (k in 1:K) {
beta[k] ~ normal(beta0, beta_sig);
y[k][2:N[k]] ~ normal(alpha + beta[k]*y[k][1:(N[k]-1)], sigma);
}
}
R сценарий
код ниже воспроизводит мои данные, Стан код и фитинга в R.
library(rstan)
dat <- structure(list(N = c(3, 5, 10), K = 3, maxN = 10, y = list(c(-0.626453810742332,
-0.0695718108004915, -0.00443440448115216, 100, 100, 100, 100,
100, 100, 100), c(1.59528080213779, 0.930591178250432, 0.483248750713414,
0.390367280599556, 0.3690161108127, 100, 100, 100, 100, 100),
c(0.575781351653492, 0.587664377772507, 0.780076056840342,
0.841052774797451, 0.794823439263525, 0.593871106619423,
0.746977087771791, 0.767786018093089, 0.789388389973885,
0.904833172045027))), .Names = c("N", "K", "maxN", "y"))
stancode <- structure("data {\n int<lower=0> K; // number of series\n int<lower=0> N[K]; // length of each series\n\n int<lower=0> maxN; // maximum length of series\n vector[maxN] y[K]; // time series\n}\n\nparameters {\n // For simplicity, assume only beta varies across series\n real alpha;\n real beta[K];\n real<lower=0> sigma;\n\n real beta0;\n real<lower=0> beta_sig;\n}\n\nmodel {\n for (k in 1:K) {\n beta[k] ~ normal(beta0, beta_sig);\n y[k][2:N[k]] ~ normal(alpha + beta[k]*y[k][1:(N[k]-1)], sigma);\n }\n}", model_name2 = "HierarchicalAR")
fit <- stan(model_code = stancode, data = dat,
pars = c("alpha", "beta", "sigma", "beta0", "beta_sig"))
Спасибо за ваш ответ! Приятно знать, что моя настройка не такая уж плохая идея. Кажется, это ответ, но я сохраню этот поток на пару дней на всякий случай, если у кого-то будут разные предложения. –
В руководстве также есть глава о том, как использовать длинные структуры данных вместо широких, если вы хотите, чтобы решение было более жестким в памяти. Независимо от того, что вы делаете, вы, вероятно, не должны быть параметрами дополнений. Существует также одна разреженная функция векторного умножения, которую мы используем больше всего, что может помочь в некоторых ситуациях. –
@BobCarpenter, спасибо. Я полагаю, вы указали главу 13 «Редкие и рваные структуры данных» руководства Стэна. Я еще не добрался туда, но это определенно обсуждает, как я могу работать! –