neodistr (version 0.1.2)

stanf_jfst: Stan function of Jones and Faddy's Skew-t Distribution

Description

Stan code of JFST distribution for custom distribution in stan

Usage

stanf_jfst(vectorize = TRUE, rng = TRUE)

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.

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

Author

Anisa' Faoziah and Achmad Syahrul Choir

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 cumulative distribution, log complementary cumulative distribution, quantile, random number of Jones-Faddy's Skew-t distribution

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

if (FALSE) {
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 n;
      vector[n] y;
    }
    parameters{
      real mu;
      real  sigma;
      real  alpha;
      real  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 n;
      vector[n] y;
    }
    parameters{
      real mu;
      real  sigma;
      real  alpha;
      real 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))
 }