用Jags估计dacay常量

我尝试估计我的数据的衰减。让我解释。想象一下，你观察事件X，并且每个事件都有一个影响y，随着时间exp（-t/Tau）衰减。你观察时间t和事件x以及预测它的影响y。让我给你看看我的JAGS代码。用Jags估计dacay常量

model{ 
for(j in 1:N){ 
    for(i in 1:p){ 
    td[j,i] <- exp(- t[j,i]/Tau[i]) 
    } 

    mu[j] <- inprod(X[j,]*td[j,] ,beta[]) 
    Y[j] ~ dnorm(mu[j], sigma) 
} 

    for(j in 1:p){ 
    bsigma[j] ~dgamma(0.001,0.001); 
    beta[j] ~ dnorm(0,bsigma[j]); 
    Tau[j] ~ dgamma(0.001,0.001); 
    } 
    sigma ~ dgamma(0.001,0.001) 
}

我产生R测试数据如下：

N = 1000; 

sigma = 0.1; 
beta = c(0.75,0.33) 
tau = c(5.7,1.3) 

X<-cbind(rnorm(N,1,1),rnorm(N,2,1)) 
t<-cbind(rnorm(N,1,1),rnorm(N,2,1)) 
t = abs(t); 
Y<- rnorm(N,(X*exp(- t/tau))%*%as.matrix(beta),sigma)

随着我的模型，我可以顺利找到测试值，但我无法估计头正确的价值观。
这里的完整代码：

N = 1000; 

sigma = 0.1; 
beta = c(0.75,0.33) 
tau = c(5.7,1.3) 

X<-cbind(rnorm(N,1,1),rnorm(N,2,1)) 
t<-cbind(rnorm(N,1,1),rnorm(N,2,1)) 
t = abs(t); 
Y<- rnorm(N,(X*exp(- t/tau))%*%as.matrix(beta),sigma) 

####JAGS 
################## 
library(mcmcplots) 
library(runjags) 
library(rjags) 

hmodel_jags<- function(X,Y,t){ 
    modelstring = " 
    model{ 
    for(j in 1:N){ 
     for(i in 1:p){ 
     td[j,i] <- exp(- t[j,i]/Tau[i]) 
     } 

    mu[j] <- inprod(X[j,]*td[j,] ,beta[]) 
     Y[j] ~ dnorm(mu[j], sigma) 
    } 

    for(j in 1:p){ 
     bsigma[j] ~dgamma(0.001,0.001); 
     beta[j] ~ dnorm(0,bsigma[j]); 
     Tau[j] ~ dgamma(0.001,0.001); 
    } 
    sigma ~ dgamma(0.001,0.001) 
    }" 
    writeLines(modelstring,con="dec.txt") 
    ######## 
    set.seed(123) 


    jags_data <- list(Y = Y, 
       t = t, 
       X = X, 
       p = ncol(X), 
       N=nrow(X) 
       ) 
    params <- c("Tau",'sigma','beta') 
    adapt <- 1000 
    burn <- 1000 
    iterations <- 1000 
    inits <- list() 

    sample <- run.jags(model="dec.txt", thin =2, monitor=params,data=jags_data, n.chains=2, inits=inits, adapt=adapt, burnin=burn,  sample=iterations, summarise=T, method="parallel") 

    sample 
} 
res_jags_het <- hmodel_jags(X,Y,t)

来源

2015-12-02 user4666604

我觉得很难理解你的代码。如果你从模拟开始解释什么样的数据需要处理，以及你想估计什么参数（如果你避免使用sigma作为JAGS代码的精度，那么也会有帮助， R脚本中的标准偏差）。一个直接的困惑：在你的R脚本'beta'中如果修复的话。但是在你的JAGS代码中，不是给出一个“beta”以前的无信息，而是假设它来自具有固定平均值和模拟精度的正态分布。为什么？ –

感谢您的评论，并指出我与西格玛的错误。观察到t，y，x。我想估计beta和tau。以前均值为0的正态伽马组合应该类似于岭回归或者像ARD那样考虑它。贝叶斯回归，基本上。 – user4666604

错就错在你的数据模拟。您有

Y<- rnorm(N,(X*exp(- t/tau))%*%as.matrix(beta),sigma)

专注于t/tau。你看，如果你做

t <- matrix(c(1,1,1,1,1,1,2,2,2,2,2,2), ncol=2) 
tau <- c(1,20) 
t/tau 
[,1] [,2] 
[1,] 1.00 2.0 
[2,] 0.05 0.1 
[3,] 1.00 2.0 
[4,] 0.05 0.1 
[5,] 1.00 2.0 
[6,] 0.05 0.1

有几种方法来解决这个问题，其中最直观的就是遍历T的行会发生什么。

tt <- matrix(data=NA, nrow=dim(t)[1], ncol=dim(t)[2]) 
for(i in 1:dim(t)[1]){ 
    tt[i,] <- t[i,]/tau 
} 
tt 
[,1] [,2] 
[1,] 1 0.1 
[2,] 1 0.1 
[3,] 1 0.1 
[4,] 1 0.1 
[5,] 1 0.1 
[6,] 1 0.1

虽然我没有时间重新运行JAGS模型以确认这是唯一的问题 - 我必须走出门！

来源

2015-12-03 15:11:42

是的，就是这样。非常感谢！愚蠢的我，我没有适当地调试它，我假设它做我认为它做的。你真的帮助我。谢啦！ – user4666604

用Jags估计dacay常量

回答

相关问题