generate uncertainty for a prediction using JAGS via rjags in R Say I have a relationship between two variables, that I have successfully fit using JAGS via the rjags package for R. Below is the code to generate the data, and the fitted JAGS model.
data:
#generate data
age <- seq(1:200)
age.m1 <- 0.1    * exp(-age * .1) #mortality rate declines initially from juvenile stage
age.m2 <- 0.001 * exp(age* 0.015) #mortality rate creeps up as trees get older
x <- age.m1 + age.m2 #combine the two mortality rates
#generate some scatter in the data to plot
y <- age.m3 + rnorm(length(age.m3),0, 0.002)

JAGS model:
require(rjags)

jags.model = "
model {
    for (i in 1:N){
        y[i] ~ dnorm(y.hat[i], tau)
        y.hat[i] <- c*exp(-a * x[i]) + d*exp (b * x[i])
        }
    a ~ dnorm(0, .0001)I(0,10)
    b ~ dnorm(0, .0001)I(0,10)
    c ~ dnorm(0, .0001)I(0,10)
    d ~ dnorm(0, .0001)I(0,10)
    tau <- pow(sigma, -2)
    sigma ~ dunif(0, 100)
}
"

data = list(y=y, x=x, N = length(y))

j.model <- jags.model(file = textConnection(jags.model),
                          data=data,
                          n.chains=3)

#sample from the posterior
jags.out   <- coda.samples (model = j.model,
                            variable.names = c('a','b','c','d'),
                            n.iter = 5000)

I can then predict y as a function of x, based on known values of x, and parameter estimates from my JAGS model, like this:
#grab parameter values, estimate y over the range of x.
pars <- summary(jags.out)$statistics[,1]
m.out <- pars[3] * exp(-pars[1] * x) + pars[4] * exp(pars[2]* x)

#plot predicted y over the range of x
plot(m.out ~ x, pch=16, cex = 0, main = 'fitted JAGS model', ylim=c(0,0.1))
lines(smooth.spline(m.out ~ age), lwd = 2, col = 'purple')
par(new=T)
plot(scatter, pch=16, cex = 0.5, ylim=c(0,0.1))


However, this is only the estimate of the mean. Using R, how can I go about estimating +/- 1 standard deviation on this mean prediction estimate, accounting for uncertainty in all the parameter values used to generate it? If this was a fit done in lm/glm I would use the predict function to generate +/- 1 standard error of the fit.
 A: This can be done by estimating the tau parameter, the measurement of the precision of y, when calling coda while sampling the JAGS model output. Modify the coda line to include tau:
#sample from the posterior
jags.out   <- coda.samples (model = j.model,
                            variable.names = c('a','b','c','d','tau'),
                            n.iter = 5000)

Then estimate the upper and lower bounds of y by addding or substracting sqrt(1/tau), which is the standard deviation of the estimate. 
pars <- summary(jags.out)$statistics[,1]
m.out       <- pars[3] * exp(-pars[1] * age) + pars[4] * exp(pars[2]* age)
m.out.plus  <- pars[3] * exp(-pars[1] * age) + pars[4] * exp(pars[2]* age) + sqrt(1/pars[5])
m.out.minus <- pars[3] * exp(-pars[1] * age) + pars[4] * exp(pars[2]* age) - sqrt(1/pars[5])

You can then plot a shaded confidence interval that represents +/- 1 standard devation with the following code:
plot(m.out ~ age, pch=16, cex = 0, main = 'fitted JAGS model', ylim=c(0,0.1))
lines(smooth.spline(m.out ~ age), lwd = 2, col = 'purple')
polygon(c((1:length(x)), rev(1:length(x))),c(m.out.plus, rev(m.out.minus)), col=adjustcolor('purple', 0.3), lty=0)
par(new=T)
plot(scatter, pch=16, cex = 0.5, ylim=c(0,0.1))


