1

我正在尝试在 R 中拟合贝叶斯模型。我已经尝试了所有可用的软件包（rstanarm 和 brms），但它们需要几个小时。这就是我尝试编写 stan 代码的原因。为了让您了解模型，我正在尝试拟合具有随机效应（弹性）的对数模型。因此，在使用 log1p 转换我的数据后，我运行以下代码：

 data {
int<lower=0> N; // num individuals
int<lower=1> K; // num ind predictors
int<lower=1> J; // num groups
// int<lower=1> L; // num group predictors
int<lower=1,upper=J> jj[N]; // group for individual
matrix[N, K] x; // individual predictors
 // matrix[17,1] u; // group predictors
vector[N] y; // outcomes
}
parameters {
matrix[K, J] z;
cholesky_factor_corr[K] L_Omega;
vector<lower=0,upper=pi()/2>[K] tau_unif;
// matrix[1, K] gamma; // group coeffs
real<lower=0> sigma; // prediction error scale
}
transformed parameters {
matrix [J,K] beta;
vector<lower=0>[K] tau; // prior scale
for (k in 1:K) tau[k] = 2.5 * tan(tau_unif[k]);
// beta = u * gamma + (diag_pre_multiply(tau,L_Omega) * z)';
beta = (diag_pre_multiply(tau,L_Omega) * z)';
}
model {
to_vector(z) ~ normal(0, 1);
L_Omega ~ lkj_corr_cholesky(2);
// to_vector(gamma) ~ normal(0, 5);
y ~ normal(rows_dot_product(beta[jj] , x), sigma);
}

即使我只包含一个截距和一个回归量，拟合模型也需要 1 个小时以上。有任何想法吗 ？

X          <- mydata[,c("Price", "Distribution", "Distribution on Display","Distribution on Leaflet","Distribution Leaflet + Display"]
X          <- cbind(intercept=1,X)
stanDat    <- list()
stanDat$N  <- nrow(mydata)                           # num of observations   
stanDat$K  <- ncol(X)                                 # num of predictors
stanDat$J  <- length(unique(mydata$Brand))                # number of groups
stanDat$jj <- as.numeric(factor(mydata$Brand))         # group for individual
stanDat$x  <- X;                                      # individual predictors
stanDat$y  <- log1p(mydata$`number of units in the package`);                # outcomes

以及我如何称呼斯坦：

fit <- stan(data=stanDat, file= "blah blah\mymodel2.stan", cores=3, chains = 3,warmup = 500,iter = 1200,control = list(max_treedepth=15,adapt_delta = .99))
4

0 回答 0