stanf_jfst: Stan function of Jones and Faddys Skew-t Distribution In neodistr: Neo-Normal Distribution

 stanf_jfst R Documentation

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

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
),
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