R中的STAN：线性回归中的尺寸误差

Question

您可以在下面找到一个带有 2 个预测变量的线性回归的简单示例。它运作良好。但是，当只有 1 个预测器时（参见第二个脚本），会出现以下错误消息：

例外：在上下文中声明和发现的数字维度不匹配；处理阶段=数据初始化；变量名=x;暗淡声明=（20,1）；发现暗淡=(20)

问题在于 1 行的矩阵会自动转换为向量，因此与声明的维度不匹配。一种解决方案是将 x 声明为向量，但问题是我正在运行具有不同数量的预测变量（可能是 1 或更多）的相同脚本。

STAN 脚本：

write("// Stan model for simple linear regression
data {
int<lower=0> N;  // number of data items
int<lower=0> K;// number of predictors
matrix[N, K] x;// predictor matrix
vector[N] y;// outcome vector
}
parameters {
real alpha;       // intercept
vector[K] beta;       // coefficients for predictors
real<lower=0> sigma;  // error scale
}
model {
y ~ normal(x * beta + alpha, sigma);  // likelihood
}", "ex_dimension.stan")

具有 2 个预测变量的 R 脚本（工作）：

N=20
K=2
x1=1:N+rnorm(N,0,0.5)
x2=rnorm(N,2,1)
x=cbind(x1,x2)
a=2
b=3
y=a*x1+b*x2+rnorm(N,0,1)


stan_data=list(N=N,
               K=K,
               x=x,
               y=y)
fit <- stan(file = "ex_dimension.stan",
            data = stan_data,
            warmup = 500,
            iter = 2000,
            chains = 4,
            cores = 4,
            thin = 1,
            control=list(adapt_delta=0.8))
fit

带有 1 个预测器的脚本（不工作）：

stan_data=list(N=N,
               K=1,
               x=x[,1],
               y=y)
fit <- stan(file = "ex_dimension.stan",
            data = stan_data,
            warmup = 500,
            iter = 2000,
            chains = 4,
            cores = 4,
            thin = 1,
            control=list(adapt_delta=0.8))
fit

Answer 1

使用drop = FALSE 对矩阵进行子集化，以避免将其折叠为向量，从而将错误的输入传递给 Stan 模型（另请参见 Advanced R - Subsetting Chapter）。

library(rstan)

stan_data <- list(N = N, K = 1, x = x[, 1, drop = FALSE], y = y)

fit <- stan(
    model_code = "// Stan model for simple linear regression
        data {
            int<lower=0> N;       // number of data items
            int<lower=0> K;       // number of predictors
            matrix[N, K] x;       // predictor matrix
            vector[N] y;          // outcome vector
        }
        parameters {
            real alpha;           // intercept
            vector[K] beta;       // coefficients for predictors
            real<lower=0> sigma;  // error scale
        }
        model {
            y ~ normal(x * beta + alpha, sigma);  // likelihood
        }",
    data = stan_data,
    chains = 1
)

fit
#> Inference for Stan model: 4f8ba0f0c644593f519910e9d2741995.
#> 1 chains, each with iter=2000; warmup=1000; thin=1; 
#> post-warmup draws per chain=1000, total post-warmup draws=1000.
#> 
#>           mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
#> alpha     6.26    0.06 1.20   3.93   5.49   6.25   7.04   8.68   470    1
#> beta[1]   2.00    0.00 0.10   1.81   1.94   2.00   2.06   2.19   453    1
#> sigma     2.70    0.02 0.50   1.87   2.35   2.62   2.97   3.88   458    1
#> lp__    -28.15    0.06 1.21 -31.12 -28.80 -27.84 -27.23 -26.74   366    1
#> 
#> Samples were drawn using NUTS(diag_e) at Thu Aug 15 12:41:19 2019.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).