Neo-Normal Distributions Family for MCMC Models in JAGS

if (requireNamespace("neojags", quietly = TRUE)){
      neojags::load.neojagsmodule()
} 
#> module neojags loaded
if (requireNamespace("neojags", quietly = TRUE)){
      library(rjags)
} 
#> Loading required package: coda
#> Linked to JAGS 4.3.2
#> Loaded modules: basemod,bugs,neojags

Generate data

Create model for JAGS

mod <- "
model {
  # Likelihood
  for (i in 1:100) {
    x[i] ~ djskew.ep(2,1,0.8,1)
  }
}
"

Compile the model

modelv <- jags.model(textConnection(mod), n.chains=1, inits = list(".RNG.name" = "base::Wichmann-Hill",".RNG.seed" = 314159))
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 0
#>    Unobserved stochastic nodes: 100
#>    Total graph size: 103
#> 
#> Initializing model

Generate samples

samplesv <- coda.samples(modelv, variable.names = c("x"), n.iter = 1)
gen_datav <- (as.data.frame(as.matrix(samplesv)))
x <- as.numeric(gen_datav[1,])

Parameter Estimation

Create model for JAGS

model_string <- "
model {
  # Likelihood
  for (i in 1:100) {
    x[i] ~ djskew.ep(mu, tau,nu1, nu2)
  }
  
  # Prior distributions
  mu ~ dnorm(2,10000)
  tau ~ dgamma(10,10)
  nu1 ~ dgamma(10,10)
  nu2 ~ dgamma(10,10)
}
"

Compile the model

model <- jags.model(textConnection(model_string), data = list(x=c(x)),n.chains=2)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 100
#>    Unobserved stochastic nodes: 4
#>    Total graph size: 107
#> 
#> Initializing model

Generate samples from the posterior distribution

samples<- coda.samples(model, variable.names = c("mu", "tau", "nu1", "nu2"), n.iter = 2000)

Summary Samples

summary(samples)
#> 
#> Iterations = 1001:3000
#> Thinning interval = 1 
#> Number of chains = 2 
#> Sample size per chain = 2000 
#> 
#> 1. Empirical mean and standard deviation for each variable,
#>    plus standard error of the mean:
#> 
#>       Mean     SD  Naive SE Time-series SE
#> mu  1.9986 0.0101 0.0001597      0.0002007
#> nu1 0.7446 0.0634 0.0010024      0.0021681
#> nu2 1.1688 0.1547 0.0024468      0.0049008
#> tau 0.9430 0.2540 0.0040161      0.0089956
#> 
#> 2. Quantiles for each variable:
#> 
#>       2.5%    25%    50%    75%  97.5%
#> mu  1.9786 1.9918 1.9985 2.0054 2.0183
#> nu1 0.6338 0.6993 0.7402 0.7847 0.8787
#> nu2 0.9033 1.0599 1.1553 1.2610 1.5043
#> tau 0.5380 0.7616 0.9160 1.0883 1.5210

Traceplot

traceplot(samples)

Probability Density Function (PDF), Cumulative Density Function (CDF), and Invers CDF (Quantile)

model_string1 <- "
model {
    d <- djskew.ep(0.5,2,2,2,2)
        p <- pjskew.ep(0.5,2,2,2,2)
        q <- qjskew.ep(0.5,2,2,2,2)
}
"

Compile the model

model1 <- jags.model(textConnection(model_string1),  n.chains=2)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 0
#>    Unobserved stochastic nodes: 0
#>    Total graph size: 5
#> 
#> Initializing model

Generate samples from the posterior distribution

samples1<- coda.samples(model1, variable.names = c("d","p","q"), n.iter = 2)

Summary samples

summary(samples1)
#> 
#> Iterations = 1:2
#> Thinning interval = 1 
#> Number of chains = 2 
#> Sample size per chain = 2 
#> 
#> 1. Empirical mean and standard deviation for each variable,
#>    plus standard error of the mean:
#> 
#>       Mean SD Naive SE Time-series SE
#> d 0.008864  0        0              0
#> p 0.001350  0        0              0
#> q 2.000000  0        0              0
#> 
#> 2. Quantiles for each variable:
#> 
#>       2.5%      25%      50%      75%    97.5%
#> d 0.008864 0.008864 0.008864 0.008864 0.008864
#> p 0.001350 0.001350 0.001350 0.001350 0.001350
#> q 2.000000 2.000000 2.000000 2.000000 2.000000