STAN 中的简单 Gamma GLM
simple Gamma GLM in STAN
我正在 STAN 和 R 中尝试一个简单的 Gamma GLM,但它立即崩溃
生成数据:
set.seed(1)
library(rstan)
N<-500 #sample size
dat<-data.frame(x1=runif(N,-1,1),x2=runif(N,-1,1))
#the model
X<-model.matrix(~.,dat)
K<-dim(X)[2] #number of regression params
#the regression slopes
betas<-runif(K,-1,1)
shape <- 10
#simulate gamma data
mus<-exp(X%*%betas)
y<-rgamma(500,shape=shape,rate=shape/mus)
这是我的 STAN 模型:
model_string <- "
data {
int<lower=0> N; //the number of observations
int<lower=0> K; //the number of columns in the model matrix
matrix[N,K] X; //the model matrix
vector[N] y; //the response
}
parameters {
vector[K] betas; //the regression parameters
real<lower=0, upper=1000> shape; //the shape parameter
}
model {
y ~ gamma(shape, (shape/exp(X * betas)));
}"
当我 运行 这个模型时,R 立即崩溃:
m <- stan(model_code = model_string, data = list(X=X, K=3, N=500, y=y), chains = 1, cores=1)
update :我认为问题出在向量化的某个地方,因为我可以获得一个 运行ning 模型,我将 X 的每一列作为向量传递。
update2:这也有效
for(i in 1:N)
y[i] ~ gamma(shape, (shape / exp(X[i,] * betas)));
原始代码的问题是目前在 Stan 中没有为标量除以向量定义运算符。在这种情况下,
shape / exp(X * betas)
你也许可以做
shape[1:N] ./ exp(X * betas)
或者做不到,
(shape * ones_vector) ./ exp(X * betas)
我正在 STAN 和 R 中尝试一个简单的 Gamma GLM,但它立即崩溃
生成数据:
set.seed(1)
library(rstan)
N<-500 #sample size
dat<-data.frame(x1=runif(N,-1,1),x2=runif(N,-1,1))
#the model
X<-model.matrix(~.,dat)
K<-dim(X)[2] #number of regression params
#the regression slopes
betas<-runif(K,-1,1)
shape <- 10
#simulate gamma data
mus<-exp(X%*%betas)
y<-rgamma(500,shape=shape,rate=shape/mus)
这是我的 STAN 模型:
model_string <- "
data {
int<lower=0> N; //the number of observations
int<lower=0> K; //the number of columns in the model matrix
matrix[N,K] X; //the model matrix
vector[N] y; //the response
}
parameters {
vector[K] betas; //the regression parameters
real<lower=0, upper=1000> shape; //the shape parameter
}
model {
y ~ gamma(shape, (shape/exp(X * betas)));
}"
当我 运行 这个模型时,R 立即崩溃:
m <- stan(model_code = model_string, data = list(X=X, K=3, N=500, y=y), chains = 1, cores=1)
update :我认为问题出在向量化的某个地方,因为我可以获得一个 运行ning 模型,我将 X 的每一列作为向量传递。
update2:这也有效
for(i in 1:N)
y[i] ~ gamma(shape, (shape / exp(X[i,] * betas)));
原始代码的问题是目前在 Stan 中没有为标量除以向量定义运算符。在这种情况下,
shape / exp(X * betas)
你也许可以做
shape[1:N] ./ exp(X * betas)
或者做不到,
(shape * ones_vector) ./ exp(X * betas)