**(y), not y. To get the corresponding parameters for the original-scale data, you have to transform the estimated parameters. And the prior on the mean and precision have to be appropriate for log(y), as opposed to y. This post shows an example of a JAGS program (with rjags) for estimating the parameters of log-normal distributed data.**

*log***UPDATED CODE CAN BE FOUND AT THIS MORE RECENT POST.**

The data for the example are shown as the (pink) histogram in the top left panel. A smattering of credible log-normal distributions is superimposed (as blue curves) on the data:

The upper-middle panel shows that the mean of the log-normal distribution is estimated to be 5.0. Notice that a value of 5.0 is nowhere near the bulk of the data, because the mean is for the log data, not the original-scale data. The transformed parameters are shown in the lower panels. Thus, the mode of the log-normal distribution, on the original scale, is about 115, and the mean of the log-normal distribution, on the original scale, is about 169.

The complete program I used is included below. Notice, in particular, that the priors are set to be vague on the scale of the

*data.*

**log**graphics.off()

rm(list=ls(all=TRUE))

fileNameRoot="YmetricXsingleJagsDLNORM" # for constructing output filenames

source("openGraphSaveGraph.R")

source("plotPost.R")

require(rjags) # Kruschke, J. K. (2011). Doing Bayesian Data Analysis:

# A Tutorial with R and BUGS. Academic Press / Elsevier.

#------------------------------------------------------------------------------

# THE MODEL.

modelstring = "

model {

for( i in 1 : N ) {

y[i] ~ dlnorm( muOfLogY , 1/sigmaOfLogY^2 )

}

sigmaOfLogY ~ dunif( 0.001*sdOfLogY , 1000*sdOfLogY )

muOfLogY ~ dnorm( meanOfLogY , 0.001*1/sdOfLogY^2 )

muOfY <- exp(muOfLogY+sigmaOfLogY^2/2)

modeOfY <- exp(muOfLogY-sigmaOfLogY^2)

sigmaOfY <- sqrt(exp(2*muOfLogY+sigmaOfLogY^2)*(exp(sigmaOfLogY^2)-1))

}

" # close quote for modelstring

writeLines(modelstring,con="model.txt")

#------------------------------------------------------------------------------

# THE DATA.

# Generate random data from known parameter values:

set.seed(47405)

trueM = 5.0

trueSD = 0.5

y = rnorm( n=125 ) # R dnorm uses mean and SD

LogY = (y-mean(y))/sd(y)*trueSD + trueM

y = exp(LogY)

dataList = list(

y = y ,

N = length(LogY) ,

meanOfLogY = mean(LogY) ,

sdOfLogY = sd(LogY)

)

#------------------------------------------------------------------------------

# INTIALIZE THE CHAINS.

# Let JAGS do it

#------------------------------------------------------------------------------

# RUN THE CHAINS

parameters = c("muOfLogY" , "sigmaOfLogY" , "muOfY" , "modeOfY" , "sigmaOfY" )

adaptSteps = 1000 # Number of steps to "tune" the samplers.

burnInSteps = 1000 # Number of steps to "burn-in" the samplers.

nChains = 3 # Number of chains to run.

numSavedSteps=30000 # Total number of steps in chains to save.

thinSteps=1 # Number of steps to "thin" (1=keep every step).

nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.

# Create, initialize, and adapt the model:

jagsModel = jags.model( "model.txt" , data=dataList ,

n.chains=nChains , n.adapt=adaptSteps )

# Burn-in:

cat( "Burning in the MCMC chain...\n" )

update( jagsModel , n.iter=burnInSteps )

# The saved MCMC chain:

cat( "Sampling final MCMC chain...\n" )

codaSamples = coda.samples( jagsModel , variable.names=parameters ,

n.iter=nPerChain , thin=thinSteps )

# resulting codaSamples object has these indices:

# codaSamples[[ chainIdx ]][ stepIdx , paramIdx ]

#------------------------------------------------------------------------------

# EXAMINE THE RESULTS

checkConvergence = FALSE

if ( checkConvergence ) {

openGraph(width=7,height=7)

autocorr.plot( codaSamples[[1]] , ask=FALSE )

show( gelman.diag( codaSamples ) )

effectiveChainLength = effectiveSize( codaSamples )

show( effectiveChainLength )

}

# Convert coda-object codaSamples to matrix object for easier handling.

# But note that this concatenates the different chains into one long chain.

# Result is mcmcChain[ stepIdx , paramIdx ]

mcmcChain = as.matrix( codaSamples )

chainLength = NROW(mcmcChain)

openGraph(width=10,height=6)

layout(matrix(1:6,nrow=2,byrow=TRUE))

# posterior predictive

hist( dataList$y , xlab="y" , main="Data w. Post. Pred." , breaks=30 ,

col="pink" , border="white" , prob=TRUE , cex.lab=1.5)

pltIdx = floor(seq(1,chainLength,length=20))

xComb = seq( min(dataList$y) , max(dataList$y) , length=501 )

for ( chnIdx in pltIdx ) {

lines( xComb ,

dlnorm( xComb, mcmcChain[chnIdx,"muOfLogY"], mcmcChain[chnIdx,"sigmaOfLogY"] ),

col="skyblue" )

}

# param's of log(y)

postInfo = plotPost( mcmcChain[,"muOfLogY"] , xlab="mu of log(y)" )

postInfo = plotPost( mcmcChain[,"sigmaOfLogY"] , xlab="sigma of log(y)" )

# param's of y

postInfo = plotPost( mcmcChain[,"modeOfY"] , xlab="mode of y" )

postInfo = plotPost( mcmcChain[,"muOfY"] , xlab="mu of y" )

postInfo = plotPost( mcmcChain[,"sigmaOfY"] , xlab="sigma of y" , showMode=TRUE)

saveGraph(file=fileNameRoot,type="jpg")

#------------------------------------------------------------------------------