# Difference of two distributions in JAGS: How to construct likelihood?

I'd like to get the parameters of these "unobserved variables":

X1 <- rgamma(100,shape=5, rate=0.5)
X2 <- rexp(100,rate=0.1)


from this "observed variable" that is their difference:

Z <- X1 - X2


using JAGS, mostly because then I could use the posteriors of the parameters as estimates.

It works if I do their sum as observed (i.e., Z <- X1 + X2) using

library(jagsUI)
library(rjags)
library(R2jags)
library(R2WinBUGS)

sink("gamma.txt")

cat("
model{
shap1 ~ dgamma(0.001,0.001) ## flat
rate1 ~ dgamma(0.001,0.001)
rate2 ~ dgamma(0.001,0.001)
for(i in 1:n){
G[i] ~ dgamma(shap1,rate1)
D[i] ~ dexp(rate2)
Z[i] ~ sum(G[i],D[i])
}
}",
fill=T)
sink()  #turn off connection

#data
n <- length(Z)
more.data <- c("Z","n")

#inits function
inits <- function(){list(shap1=1, rate1=1, rate2=1)}

#parameters estimated
params <- c("shap1",
"rate1",
"rate2")

#MCMC settings
nc=1
ni=1000
nb=100
nt=2

#gibbs sampler
out.put <- jags(data=more.data,
model="gamma.txt",
parameters=params,
inits=inits,
n.chains=nc, n.iter=ni, n.burnin=nb,n.thin=nt,
working.directory=getwd())


The toy above is shortened to run quickly and gives poor parameter estimates. If I bulk up the gamma and exponential samples for X1 and X2 it does quite well.

I can also solve this problem using method of moments, but I do not know how to get standard error estimates using that technique.

• Could you please decide whether you are asking about the sum or the difference? – whuber Jul 2 '18 at 3:46
• Sure, the difference. The code works for the sum. That part is easy. The difference is eluding me. – Dial Jul 2 '18 at 3:47
• And for simplicity, let's assume the difference (as in the title) between the gamma and the exponential (as in the code snippet generating the simulation toy) is between two distributions that are independent of each other. – Dial Jul 2 '18 at 3:52
• And again, as I seem to be somewhat confusing in my request, I can solve this difference problem using the first three moments to get the three parameters (two for gamma and one for exponential), but I do not know how to get their standard errors that way. I also have provided the code that gives me posteriors (a kind of substitute for standard errors) for the three parameters but only for the sum, which I really don't care about but started with because, as I have said, it was straightforward. – Dial Jul 2 '18 at 3:55
• @whuber see comments above and many thanks for responding. Sorry to be opaque. – Dial Jul 2 '18 at 4:06

In principle you can fake a negative exponential distribution by using a truncated double exponential:

    negD[i] ~ ddexp(0, rate2) T(,0)


This has the same distribution as -D[i]

Then the sum observable function should work

    Z[i] ~ sum(G[i],negD[I])


This is of course a work-around and not a generic solution to this kind of problem. Roman and I have corresponded privately on this question and here is my more generic answer:

Your example highlights one of the issues in modelling with the BUGS language. There is no general mechanism for handling posterior deterministic constraints on the unobserved nodes. Handling such cases requires me to define an "observable function" in the BUGS language (e.g. dsum, sum, dinterval) and create a custom sampling mechanism to make sure the constraints are respected at each iteration. It is the custom sampling mechanism that takes most of the work.

As you have seen, it is relatively easy to express a constraint in the BUGS language that the custom sampler cannot handle.