Hi, John. Long-time listener, first-time caller... I have a model that says X is a function of three (independent) random variables:

X ~ normal(mu, sigma) / uniform(a,a+b) - beta(v,w)

and I also have N random samples of X. Can I use JAGS to estimate mu, sigma, a, b, v and w? Or is this problem outside its scope? Not sure if it's even possible.

My reply:

In other words, you want

x <- r / u - b

where

r ~ normal( rMean , rSD )

u ~ uniform( uLow , uHigh )

b ~ beta( bA , bB )

but the problem is that JAGS won't accept that sort of model specification.

One solution is to change to an equivalent hierarchical formalization:

x ~ normal( rMean/u-b , rSD/u )

u ~ uniform( uLow , uHigh )

b ~ beta( bA , bB )

This hierarchical form is directly implementable in JAGS. But be careful that in JAGS the normal is parameterized by precision, not SD, so you'll need to write x ~ dnorm( rMean/u-b , 1/(rSD/u)^2 ).

To demonstrate that the two forms are equivalent, here is a Monte Carlo sample generated in R from each formalization:

# Specify parameter values and sample size:

rMean=7 ; rSD=2

uLow=3 ; uHigh=5

bA=4 ; bB=3

N = 50000

# Sample x values your original way:

b = rbeta(N, bA,bB)

u = runif(N, uLow,uHigh)

r = rnorm(N, rMean,rSD)

x1 = r/u-b

# Sample x values a new way:

x2 = rep(0,N)

for ( i in 1:N ) {

b = rbeta(1, bA,bB)

u = runif(1, uLow,uHigh)

x2[i] = rnorm(1, rMean/u-b , rSD/u )

}

# Plot results:

layout(matrix(1:2,nrow=2))

xLim=range(c(x1,x2))

hist(x1,xlim=xLim,breaks=21)

text( mean(x1) , 0 , "+" , cex=2 )

hist(x2,xlim=xLim,breaks=21)

text( mean(x2) , 0 , "+" , cex=2 )

The result is shown below, where you can see that the original and new distributions are identical:

ReplyDeleteGreat job done keep posting more. I will share your link with my friends.

Statistical Data analysis