stanf_jfst {neodistr}

Stan function of Jones and Faddys Skew-t Distribution

Description

Stan code of JFST distribution for custom distribution in stan

Usage

stanf_jfst(vectorize = TRUE, rng = TRUE)

Arguments

vectorize

logical; if TRUE, Vectorize Stan code of Jones and faddy distribution are given The default value of this parameter is TRUE

rng

logical; if TRUE, Stan code of quantile and random number generation of Jones and faddy distribution are given The default value of this parameter is TRUE

Details

Jones-Faddy’s Skew-t distribution has density:

f(y |\mu,\sigma,\beta,\alpha)= \frac{c}{\sigma} {\left[{1+\frac{z}{\sqrt{\alpha+\beta+z^2}}}\right]}^{\alpha+\frac{1}{2}} {\left[{1-\frac{z}{\sqrt{\alpha+\beta+z^2}}}\right]}^{\beta+\frac{1}{2}}

where -\infty<y<\infty, -\infty<\mu<\infty, \sigma>0, \alpha>0, \beta>0, z =\frac{y-\mu}{\sigma} , c = {\left[2^{\left(\alpha+\beta-1\right)} {\left(\alpha+\beta\right)^{\frac{1}{2}}} B(a,b)\right]}^{-1} ,

This function gives stan code of log density, cumulative distribution, log of cumulatif distribution, log complementary cumulative distribution, quantile, random number of Jones-Faddy's Skew-t distribution

Value

jfst_lpdf gives stan's code of the log of density, jfst_cdf gives stan's code of the distribution function, jfst_lcdf gives stan's code of the log of distribution function and jfst_rng gives stan's code of generates random numbers.

Author(s)

Anisa' Faoziah and Achmad Syahrul Choir

References

Jones, M.C. and Faddy, M. J. (2003) A skew extension of the t distribution, with applications. Journal of the Royal Statistical Society, Series B, 65, pp 159-174

Rigby, R.A. and Stasinopoulos, M.D. and Heller, G.Z. and De Bastiani, F. (2019) Distributions for Modeling Location, Scale, and Shape: Using GAMLSS in R.CRC Press

Examples

library (neodistr)
library (rstan)

# inputting data
set.seed(400)
dt <- neodistr::rjfst(100,mu=0, sigma=1, alpha = 2, beta = 2) # random generating JFST data
dataf <- list(
 n = 100,
 y = dt
 )
 
 
#### not vector
## Calling the function of the neo-normal distribution that is available in the package.
func_code<-paste(c("functions{",neodistr::stanf_jfst(vectorize=FALSE),"}"),collapse="\n")

# Define Stan Model Code
model <-"
    data{
      int<lower=1> n;
      vector[n] y;
    }
    parameters{
      real mu;
      real <lower=0> sigma;
      real <lower=0> alpha;
      real <lower=0> beta;
    }
    model {
      for(i in 1 : n){
      y[i] ~ jfst(mu,sigma, alpha, beta);
      }
      mu ~ cauchy(0,1);
      sigma ~ cauchy(0, 2.5);
      alpha ~ lognormal(0,5);
      beta ~ lognormal(0,5);
      
    }
"

# Merge stan model code and selected neo-normal stan function
fit_code <- paste (c(func_code,model,"\n"), collapse = "\n")

# Create the model using Stan Function
fit1 <- stan(
    model_code = fit_code,  # Stan Program
    data = dataf,           # named list data
    chains = 2,             # number of markov chains
    warmup = 5000,          # total number of warmup iterarions per chain
    iter = 10000,           # total number of iterations iterarions per chain
    cores = 2,              # number of cores (could use one per chain)
    control = list(         # control sampel behavior
      adapt_delta = 0.99
    ),
    refresh = 1000          # progress has shown if refresh >=1, else no progress shown
)

# Showing the estimation result of the parameters that were executed using the Stan file
print(fit1, pars = c("mu", "sigma", "alpha", "beta", "lp__"), probs=c(.025,.5,.975))


#### Vector
## Calling the function of the neonormal distribution that is available in the package.
func_code_vector<-paste(c("functions{",neodistr::stanf_jfst(vectorize=TRUE),"}"),collapse="\n")

# Define Stan Model Code
model_vector <-"
    data{
      int<lower=1> n;
      vector[n] y;
    }
    parameters{
      real mu;
      real <lower=0> sigma;
      real <lower=0> alpha;
      real <lower=0>beta;
    }
    model {
      y ~ jfst(rep_vector(mu,n),sigma, alpha, beta);
      mu ~ cauchy (0,1);
      sigma ~ cauchy (0, 2.5);
      alpha ~ lognormal(0,5);
      beta ~ lognormal(0,5);
      
    }
 "
 
 # Merge stan model code and selected neo-normal stan function
fit_code_vector <- paste (c(func_code_vector,model_vector,"\n"), collapse = "\n")

# Create the model using Stan Function
fit2 <- stan(
    model_code = fit_code_vector,  # Stan Program
    data = dataf,                  # named list data
    chains = 2,                    # number of markov chains
    warmup = 5000,                 # total number of warmup iterarions per chain
    iter = 10000,                  # total number of iterations iterarions per chain
    cores = 2,                     # number of cores (could use one per chain)
    control = list(                # control sampel behavior
      adapt_delta = 0.99
    ),
    refresh = 1000                 # progress has shown if refresh >=1, else no progress shown
)

# Showing the estimation result of the parameters that were executed using the Stan file
print(fit2, pars = c("mu", "sigma", "alpha", "beta", "lp__"), probs=c(.025,.5,.975))