5

I'm writing a JAGS script (hierarchical bayesian model) where the times of events are modelled as a race between two processes.

Observations: time is the measured times of events.

Model: two processes with gaussian rates - whichever process finishes first triggers the event.

Goal: Estimate the rates of the two processes. 

model{
  # Priors
  mu1 ~ dnorm( 0,1 )                      # rate of one process
  mu2 ~ dnorm( 0,1 )                      # rate of other process
  sigma1 <- 1                             # variability in rate
  sigma2 <- 0.1                           # variability in rate

  # Observations
  for (i in 1:N)
    rate1[i] ~ dnorm( mu1, sigma1 )        # Sample the two
    rate2[i] ~ dnorm( mu2, sigma2 )        # racing processes.

    rmax[i] <- max( rate1[i], rate2[i] )  # which was faster?

    time[i] ~ 1/rmax[i]      #### This is wrong!
  }
}

Question: How can I indicate that the times are sampled from the larger of two rates, each of which is sampled from a different distribution?

Example histogram of simulated time data, using mu1=3, mu2=3 with different standard deviations for the two processes (fixed at 1 and 0.1)

enter image description here

Sanjay Manohar
  • 6,920
  • 3
  • 35
  • 58
  • Just to clarify: do you know which of the processes finished first when the event is triggered? Or do you just know the time at which the first process finished without knowing which process won? I guess the latter but it does strongly affect the answer! – Matt Denwood Oct 09 '18 at 11:42
  • In this situation we don't know, hence the need to get the 'max' distribution. I'll update with an example event-time histogram. – Sanjay Manohar Oct 09 '18 at 11:45
  • I am wondering if there's a way to do this with censoring / `dinterval`? – Sanjay Manohar Oct 09 '18 at 12:16
  • You could certainly use dinterval if you knew which process finished first. Otherwise if this is information is latent then you need something like dnormmix but where the PDF represents the max/min of normals (with different parameterisations) rather than a mixture. Then you could use the zeros/ones trick to fit the distribution to your data. Otherwise you would need an intermediate layer (so that the row-wise expected failure time for each process is modelled explicitly) but this would necessitate another error distribution, which would not be identifiable with your data structure. – Matt Denwood Oct 09 '18 at 12:27
  • Dear @MattDenwood I really appreciate your advice on this - the bounty expires tomorrow but if you had an attempt at coding and answer (if you think this problem has a solution) that would be amazing - or even an answer with a solid reason why this sort of model is not possible in JAGS, that would also be acceptable. – Sanjay Manohar Oct 15 '18 at 14:27

0 Answers0