2016-02-02 4 views
2

У меня есть набор переменных временного ряда различной длины и Я пытаюсь оценить иерархическую модель авторегрессии с использованием языка СТ (используя 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")) 

ответ

1

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

+0

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

+1

В руководстве также есть глава о том, как использовать длинные структуры данных вместо широких, если вы хотите, чтобы решение было более жестким в памяти. Независимо от того, что вы делаете, вы, вероятно, не должны быть параметрами дополнений. Существует также одна разреженная функция векторного умножения, которую мы используем больше всего, что может помочь в некоторых ситуациях. –

+0

@BobCarpenter, спасибо. Я полагаю, вы указали главу 13 «Редкие и рваные структуры данных» руководства Стэна. Я еще не добрался туда, но это определенно обсуждает, как я могу работать! –