Estimating dacay constant with Jags

Question

I try to estimate decays for my data. let me explain. imagine, you observe events X, and each event has an effect y, which decays over time exp(-t/Tau). You observe the time t and the event x and what to predict its effect y. Let me show you my JAGS code.

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)
}

I generate test data in R as follows:

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)

With my model I can successfully find the values for beta, but I fail to estimate the correct values for Tau.
here the complete code:

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

The mistake is in your data simulation. You have

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

Focus on t/tau. Look what happens if you do

 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

There are several ways to fix this, the most intuitive of which is to loop over the rows of 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

I haven't had time to rerun the JAGS model though to confirm that this is the only problem--I gotta get out the door!