Saturday, April 27, 2013

Bayesian estimation of log-normal parameters

Using the log-normal density can be confusing because it's parameterized in terms of the mean and precision of the log-scale data, not the original-scale data. Thus, if your data, y, are nicely described by a log-normal distribution, the estimated mean and precision are for log(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.

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 log data.


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")

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


10 comments:

  1. Professor Kruschke,

    Thank your for this example. At first glance, I thought that it would be incredibly helpful for my own work... I work with quite a bit of environmental data where a lot of it is lognormal. Unfortunately, as I began working through the example, I was stumped in your data section. Can you explain the reasoning and origin of your LogY = (y-mean(y))/sd(y)*trueSD + trueM?

    Is this just LogY = zscore*B1 + B0 where B1 = trueSD and B0 = trueM?

    Could this be done in a more straight forward way by just taking the log of an appropriately parameterized random draw?

    And, the subtle differences between y ~ rnorm( n=125 ) and y = exp(LogY)? The two y's are not the same.

    I am asking because I am trying to generate a lognormal prior using a known, but N < 5 previous samples. I am trying to illustrate the use of Bayesian techniques with small sample sets and when N approaches infinity due to Monte Carlo simulation runs.

    Thank you for any guidance you can offer. If the questions are too basic, I'd be happy to read the any appropriate background information you could offer.

    ReplyDelete
  2. Dear Anonymous:

    Sorry if the data-generation code was mysterious. I didn't bother explaining it because in real life the data come from actual measurements, not from a random data generator.

    It would have been easiest to generate log-normal data with the command rlnorm. For example:

    > y = rlnorm(100,mean=5,sd=0.5)

    Then, if we look at the generated y values, we see that the generated sample data have mean and sd close to the true parameter values, but not exactly:

    > mean(log(y))
    [1] 4.911132
    > sd(log(y))
    [1] 0.5223407

    I simply wanted the sample to have sample mean and sample sd that exactly match the true parameter values, so I first standardized normally generated data, then multiplied by the desired sd and added the mean, then exponentiated. Then mean(log(y)) is exactly the desired sample mean, and sd(log(y)) is exactly the desired sample sd.

    ReplyDelete
  3. Thank you! This helps clear things up and knock me out of my mental record skip... ;)

    ReplyDelete
  4. Hi,

    Thank you. I have been having trouble with estimating log normal parameter. However, how do I make a posterior predictive distribution from the log normal parameters and make a 95% CI plot. Here is what I did:
    logofnewY = Normal(muoflogY(sample_i),sigmaoflogY(sample_i))

    I tried the above method and back transform the logofnewY to newY. It fails.

    Something that I don't understand. If there any references that I can refer to ?

    ReplyDelete
  5. Trying to work through this example and it looks like the graph of the Data w. Post. Pred. command is not quite correct. There may be some commas missing. Here's my partial correction: 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 )

    There is still an error for printing the blue sample distributions that I haven't been able to fix (e.g. "showMode" is not a graphical parameter).

    ReplyDelete
  6. Sean:

    The R code was merely rendered improperly in your browser because of a change I made to the HTML template a few days ago. (That is, until a few days ago, the R code was shown correctly. The new template interpreted the "$" as a formatting command instead of as an R symbol.) I have corrected the HTML template and the R code should be rendered correctly now. Thanks for letting me know of the glitch, and sorry for any frustration it may have caused.

    ReplyDelete
  7. P.S.: This post uses old R code that is pre-second edition of the book! The plotPost() command has some new arguments, and in particular there is no longer an argument "showMode". Instead, use

    cenTend=c("mode","median","mean")[1]

    You can see the definition of the plotPost() in the script, DBDA2E-utilities.R

    ReplyDelete
  8. how Bayesian estimation of log-normal parameters with winbugs?
    prior normal, posterior log-normal.
    I still do not understand. thanks

    ReplyDelete
  9. Howdy,
    I am trying to find the mean of the data set below which is lognormal dist. The mode of the post. mean= 82. However, when I run as log(y)~dnorm(mu, 1/sigma^2) as per DBDA page 456 the mode of the mean=30. Obvioulsy I am wondering how to interpret the difference between these results. All help appreciated.

    timea<-c(19 13 8 1 44 121 87 88 131 33 12 48 144 24 329 39 20 46 42 26 25 42 10 4)

    ltimea<-log(timea)
    dataList = list(
    y = timea,
    N = length(timea),
    meanOfLogY = mean(ltimea),
    sdOfLogY = sd(ltimea) )
    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) #change to original data scale.
    modeOfY <- exp(muOfLogY-sigmaOfLogY^2) #change to original scale.
    sigmaOfY <- sqrt(exp(2*muOfLogY+sigmaOfLogY^2)*(exp(sigmaOfLogY^2)-1)) #change to original data scale.
    }
    "
    writeLines(modelstring,con="model.txt")

    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 = 1 # Number of chains to run.
    numSavedSteps=30000 # Total number of steps in chains to save.
    thinSteps=2 # Number of steps to "thin" (1=keep every step).
    nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
    jagsModel = jags.model( "model.txt" , data=dataList ,
    n.chains=nChains , n.adapt=adaptSteps )

    codaSamples = coda.samples( jagsModel , variable.names=parameters ,
    n.iter=nPerChain , thin=thinSteps )
    mcmcChain = as.matrix( codaSamples )
    mc<-as.mcmc(mcmcChain)
    head(mc)

    ReplyDelete
  10. Hi Graeme.

    I have updated the code at a new post:
    http://doingbayesiandataanalysis.blogspot.com/2016/04/bayesian-estimation-of-log-normal.html

    Notice in the code the highlighted line that indicates where to put in your own data. You will find that the results for your data are a bit weird because your data are not well described by a log-normal distribution. The data are more skewed than a log-normal accommodates.

    ReplyDelete