neodistr (version 0.1.2)

stanf_jsep: Stan function of Jones Skew Exponential Power

Description

Stan code of jsep distribution for custom distribution in stan

Usage

stanf_jsep(vectorize = TRUE)

Value

jsep_lpdf gives stan's code of the log of density, jsep_cdf gives stan's code of the distribution function, and jsep_lcdf gives stan's code of the log of distribution function

Arguments

vectorize

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

Author

Meischa Zahra Nur Adhelia and Achmad Syahrul Choir

Details

The Jones Skew Exponential Power with parameters \(\mu\), \(\sigma\),\(\alpha\), and \(\beta\) has density: $$ f(y | \mu, \sigma, \alpha, \beta) = \left\{ \begin{array}{ll} \frac{c}{\sigma} \exp\left(-|z|^{\alpha}\right), & \text{if } y < \mu \\ \frac{c}{\sigma} \exp\left(-|z|^{\beta}\right), & \text{if } y \geq \mu \end{array} \right. $$ where: $$z = \frac{y - \mu}{\sigma},$$ $$c = \left[ \Gamma(1 + \beta^{-1}) + \Gamma(1 + \alpha^{-1}) \right]^{-1}.$$

References

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::rjsep(100,mu=0, sigma=1, alpha = 2, beta = 2) # random generating jsep 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_jsep(vectorize=TRUE),"}"),collapse="\n")

# Define Stan Model Code
model <-"
    data{
      int n;
      vector[n] y;
    }
    parameters{
      real mu;
      real  sigma;
      real  alpha;
      real  beta;
    }
    model {
      y ~ jsep(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 <- 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_jsep(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 ~ jsep(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))
 }