top of page

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida

you are kindly asked to include the complete citation if you used this material in a publication

â€‹

â€‹

â€‹

Code 10.7 Lognormal model in R using JAGS to describe the initial mass function (IMF)

==================================================================

library(R2jags)

â€‹

N <-nrow(IMF)
x <- IMF\$Mass

# prepare data to JAGS
jags_data <- list(x = x,
N = N                                             # sample size
)

# Fit
LNORM <-" model{
# Uniform prior for standard deviation
tau <- pow(sigma, -2)                                          # precision
sigma ~ dunif(0, 100)                                          # standard deviation
mu ~ dnorm(0,1e-3)

â€‹

# Likelihood function
for (i in 1:N){
x[i] ~ dlnorm(mu,tau)
}
}"

â€‹

# Identify parameters
params <- c("mu", "sigma")

â€‹

# Run mcmc
LN <- jags(
data = jags_data,
parameters = params,
model = textConnection(LNORM),
n.chains = 3,
n.iter = 5000,
n.thin = 1,
n.burnin = 2500)

â€‹

# Output
print(LN, justify = "left", intervals=c(0.025,0.975), digits=2)

==================================================================

â€‹

â€‹

â€‹

â€‹

â€‹

â€‹

Code 10.8 Plotting routine, in R, for Figure 10.6

====================================

require(jagstools)
require(ggplot2)

â€‹

# Create new data for prediction
M = 750
xx = seq(from = 0.75*min(x),
to = 1.05*max(x),
length.out = M)

â€‹

# Extract results
mx <- jagsresults(x=LN, params=c('mu'))
sigmax <- jagsresults(x=LN, params=c('sigma'))

â€‹

# Estimate values for the Lognormal PDF
ymean <- dlnorm(xx,meanlog=mx[,"50%"],sdlog=sigmax[,"50%"])
ylwr1 <- dlnorm(xx,meanlog=mx[,"25%"],sdlog=sigmax[,"25%"])
ylwr2 <- dlnorm(xx,meanlog=mx[,"2.5%"],sdlog=sigmax[,"2.5%"])
yupr1 <- dlnorm(xx,meanlog=mx[,"75%"],sdlog=sigmax[,"75%"])
yupr2 <- dlnorm(xx,meanlog=mx[,"97.5%"],sdlog=sigmax[,"97.5%"])

â€‹

# Create a data.frame for ggplot2
gdata <- data.frame(x=xx, mean = ymean,lwr1=ylwr1 ,lwr2=ylwr2,upr1=yupr1,upr2=yupr2)

â€‹

ggplot(gdata,aes(x=xx))+
geom_histogram(data=IMF,aes(x=Mass,y = ..density..),
colour="red",fill="gray99",size=1,binwidth = 0.075,
linetype="dashed")+
geom_ribbon(aes(x=xx,ymin=lwr1, ymax=upr1,y=NULL),
alpha=0.45, fill=c("#00526D"),show.legend=FALSE) +
geom_ribbon(aes(x=xx,ymin=lwr2, ymax=upr2,y=NULL),
alpha=0.35, fill = c("#00A3DB"),show.legend=FALSE) +
geom_line(aes(x=xx,y=mean),colour="gray25",
linetype="dashed",size=0.75,
show.legend=FALSE)+
ylab("Density")+
xlab(expression(M/M['\u0298']))+
theme_bw() +
theme(legend.background = element_rect(fill = "white"),
legend.key = element_rect(fill = "white",color = "white"),
plot.background = element_rect(fill = "white"),
legend.position = "top",
axis.title.y = element_text(vjust = 0.1,margin = margin(0,10,0,0)),
axis.title.x = element_text(vjust = -0.25),
text = element_text(size = 25))

====================================

â€‹

â€‹

â€‹

Anchor 1

Output on screen:

â€‹

Inference for Bugs model at "4", fit using jags,

3 chains, each with 5000 iterations (first 2500 discarded)

n.sims = 7500 iterations saved

mu.vect      sd.vect        2.5%      97.5%        Rhat        n.eff

mu                 -1.26          0.07        -1.40       -1.12              1         1300

sigma             1.04           0.05         0.94        1.14               1         6300

deviance       81.22          2.01        79.27     86.63               1         3400

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)

pD = 2.0 and DIC = 83.2

DIC is an estimate of expected predictive error (lower deviance is better).

bottom of page