library(HDInterval) #For calculation of HDI use hdi() from HDInterval
library(rjags) #Run JAGS
library(rstan) #Run STAN

1 Typical response distribution

Scale y Dist. \(y \sim f(\mu, ...)\) Inverse link \(\mu = g^{-1}(\eta)\)
Metric \(y \sim dnorm(\mu,\sigma)\) \(\mu = \eta\)
Binary \(y \sim dbern(\mu)\) \(\mu = logistic(\eta)\)
Nominal \(y \sim dmultinom(\mu_1, ...,\mu_p)\) \(\mu_k = \frac{exp(\eta_k)}{\sum_j exp(\eta_j)}\)
Ordinal \(y \sim dmultinom(\mu_1, ...,\mu_p)\) \(\mu = \Phi(\frac{\theta_k - \eta}{\sigma}) - \Phi(\frac{\theta_{k-1} - \eta}{\sigma})\)
Count \(y \sim dpois(\mu)\) \(\mu = exp(\eta)\)

2 Fitting normal model for 1 group with no predictors

We adapted an example in (Kruschke 2015), Chapter 16.

In this example we just need to estimate both parameters of normal distribution for random variable \(y\).

The diagram of the model structure (adapted from Textbook)

# preparing data
set.seed(03172021)
y <- rnorm(100, mean = 6, sd = 3)
Ntotal = length(y)

Create a data list, include sample mean and standard deviation in it.

dataList = list(
  y = y ,
  Ntotal = Ntotal ,
  mean_mu = mean(y) ,
  sd_mu = sd(y)
)

2.1 Running in JAGS

Specify the model. Recall that in JAGS (Web Page 2017) normal distribution is specified by precision \(\frac{1}{\sigma^2}\) instead of standard deviation or variance. Select an uninformative prior for σ and normal (conjugate) prior for μ.

What do we suggest as parameters of the normal distribution based on the sample? \[ y_i \sim N(\mu,\tau), \ \text{where} \ \tau = \frac{1}{\sigma^2} \\ \mu \sim N(M,S) \\ \sigma \sim Uniform[Lw,Mx] \ \text{OR} \ Gamma(\alpha, \beta) \]

modelString = "
  model {
    for (i in 1:Ntotal) {
      y[i] ~ dnorm(mu , 1/sigma^2) #JAGS uses precision
    }
    #mu ~ dnorm(mean_mu , 1/(100*sd_mu)^2)
    mu ~ dnorm(mean_mu , Ntotal/sd_mu^2) #JAGS uses precision
    sigma ~ dunif(sd_mu/1000, sd_mu*1000)
  }
  " # close quote for modelString
# Write out modelString to a text file
writeLines(modelString, con="TEMPmodel.txt")

Initialize chains.

initsList <-function() {
  upDown<-sample(c(1,-1),1)
  m <- mean(y)*(1+upDown*.05)
  s <- sd(y)*(1-upDown*.1) 
  list(mu = m, sigma = s)
}

Run the chains.

  1. Set parameters:
# Create, initialize, and adapt the model:
parameters = c( "mu" , "sigma")     # The parameters to be monitored
adaptSteps = 500               # Number of steps to "tune" the samplers
burnInSteps = 1000
numSavedSteps=50000
nChains = 4 
thinSteps = 1
nIter = ceiling((numSavedSteps*thinSteps)/nChains )
  1. Send model to JAGS:
jagsModel = jags.model("TEMPmodel.txt" , data=dataList , inits=initsList , 
                          n.chains=nChains , n.adapt=adaptSteps )
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 100
##    Unobserved stochastic nodes: 2
##    Total graph size: 114
## 
## Initializing model
  1. Burn in and run:
# Burn-in:
update(jagsModel, n.iter=burnInSteps )

# Run it
# The saved MCMC chain:
codaSamples = coda.samples(jagsModel, variable.names=parameters , 
                          n.iter=nIter, thin=thinSteps)

Check how chains converged.

summary(codaSamples)
## 
## Iterations = 1501:14000
## Thinning interval = 1 
## Number of chains = 4 
## Sample size per chain = 12500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean     SD Naive SE Time-series SE
## mu    5.922 0.2259 0.001010       0.001013
## sigma 3.206 0.2321 0.001038       0.001372
## 
## 2. Quantiles for each variable:
## 
##        2.5%   25%   50%   75% 97.5%
## mu    5.479 5.770 5.923 6.075 6.366
## sigma 2.790 3.044 3.194 3.354 3.704

The parameters are estimated close to what we simulated and very similar to what point estimation would give.

mean(y)
## [1] 5.922201
sd(y)
## [1] 3.173288

The plot of the samples and the densities of the parameters.

plot(codaSamples)

Plot autocorrelations.

autocorr.plot(codaSamples,ask=F)

Autocorrelation function plot shows that standard deviation effective size must be pretty small.
In fact:

effectiveSize(codaSamples)
##       mu    sigma 
## 49728.17 28610.27

Shrink factor shows that even with long memory for standard deviation distributions converged:

gelman.diag(codaSamples)
## Potential scale reduction factors:
## 
##       Point est. Upper C.I.
## mu             1          1
## sigma          1          1
## 
## Multivariate psrf
## 
## 1
gelman.plot(codaSamples)

Observed HDIs of the chains:

lapply(codaSamples,function(z) hdi(as.matrix(z)))
## [[1]]
##             mu    sigma
## lower 5.492515 2.774088
## upper 6.379380 3.670242
## attr(,"credMass")
## [1] 0.95
## 
## [[2]]
##             mu    sigma
## lower 5.474003 2.743650
## upper 6.359025 3.652929
## attr(,"credMass")
## [1] 0.95
## 
## [[3]]
##             mu    sigma
## lower 5.474012 2.755222
## upper 6.354005 3.664931
## attr(,"credMass")
## [1] 0.95
## 
## [[4]]
##             mu    sigma
## lower 5.478756 2.769952
## upper 6.360833 3.673049
## attr(,"credMass")
## [1] 0.95

2.2 Running in Stan

3 Robust estimation using t-distribution

Create a sample with heavy tails in order to check robust estimation of parameters of normal distribution. Refer this simulation of Leptokurtic Distributions.

nSample<-1000
sd.Values<-c(2,3.4,.8,2.6)
sd.process<-rep(c(rep(sd.Values[1],50),
                  rep(sd.Values[2],75),
                  rep(sd.Values[3],75),
                  rep(sd.Values[4],50)), 4)
            
plot(sd.process,type="l")

3.1 Generalized Student distribution

3.2 JAGS

3.3 Stan

References

Kruschke, John K. 2015. Doing Bayesian Data Analysis : A Tutorial with r, JAGS, and Stan. Book. 2E [edition]. Amsterdam: Academic Press is an imprint of Elsevier.
Web Page. 2017. https://sourceforge.net/projects/mcmc-jags/files/Manuals/.