用 Jags 估计 dacay 常数

Question

我尝试估计我的数据的衰减。让我解释。想象一下，你观察到事件 X，每个事件都有一个影响 y，它随时间衰减 exp(-t/Tau)。您观察时间 t 和事件 x 以及预测其影响 y 的内容。让我给你看看我的 JAGS 代码。

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)

使用我的模型，我可以成功找到 beta 的值，但我无法估计 Tau 的正确值。
这里是完整的代码：

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) 

Answer 1

错误在于您的数据模拟。你有

 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 模型以确认这是唯一的问题——我得出门了！