stanf_gmsnburr {neodistr}

Stan function of GMSNBurr Distribution

Description

Stan code of GMSNBurr distribution for custom distribution in stan

Usage

stanf_gmsnburr(vectorize = TRUE, rng = TRUE)

Arguments

vectorize

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

rng

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

Details

GMSNBurr Distribution has density:

f(y |\mu,\sigma,\alpha,\beta) = {\frac{\omega}{{B(\alpha,\beta)}\sigma}}{{\left(\frac{\beta}{\alpha}\right)}^\beta} {{\exp{\left(-\beta \omega {\left(\frac{y-\mu}{\sigma}\right)}\right)} {{\left(1+{\frac{\beta}{\alpha}} {\exp{\left(-\omega {\left(\frac{y-\mu}{\sigma}\right)}\right)}}\right)}^{-(\alpha+\beta)}}}}

where -\infty<y<\infty, -\infty<\mu<\infty, \sigma>0, \alpha>0, \beta>0 and \omega = {\frac{B(\alpha,\beta)}{\sqrt{2\pi}}}{{\left(1+{\frac{\beta}{\alpha}}\right)}^{\alpha+\beta}}{\left(\frac{\beta}{\alpha}\right)}^{-\beta}

This function gives stan code of log density, cumulative distribution, log of cumulatif distribution, log complementary cumulative distribution, quantile, random number of GMSNBurr distribution

Value

msnburr_lpdf gives the stans's code of log of density, msnburr_cdf gives the stans's code of distribution function, gmsnburr_lcdf gives the stans's code of log of distribution function, gmsnburr_lccdf gives the stans's code of complement of log ditribution function (1-gmsnburr_lcdf), and gmsnburr_rng the stans's code of generates random numbers.

Author(s)

Achmad Syahrul Choir

References

Choir, A. S. (2020). The New Neo-Normal DDistributions and their Properties. Disertation. Institut Teknologi Sepuluh Nopember.

Examples

library(neodistr)
library(rstan)
#inputting data
set.seed(136)
dt <- rgmsnburr(100,0,1,0.5,0.5) # random generating MSNBurr-IIA 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_gmsnburr(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]~gmsnburr(mu,sigma,alpha,beta);
  }
  mu~cauchy(0,1);
  sigma~cauchy(0,2.5);
  alpha~cauchy(0,1);
  beta~cauchy(0,1);
  }
   "
#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 of data
  chains = 2,             # number of Markov chains
  #warmup = 5000,          # number of warmup iterations per chain
  iter = 10000,           # total number of iterations per chain
  cores = 2,              # number of cores (could use one per chain)
  control = list(         #control samplers behavior
    adapt_delta=0.9
  )
)

# Showing the estimation results 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 neo-normal distribution that is available in the package.
func_code_vector<-paste(c("functions{",neodistr::stanf_gmsnburr(vectorize=TRUE),"}"),collapse="\n")
# define stan model as vector
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~gmsnburr(rep_vector(mu,n),sigma,alpha,beta);
   mu~cauchy(0,1);
   sigma~cauchy(0,2.5);
   alpha~cauchy(0,1);
   beta~cauchy(0,1);
 }
 "
#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 of data
  chains = 2,             # number of Markov chains
  #warmup = 5000,          # number of warmup iterations per chain
  iter = 10000,           # total number of iterations per chain
  cores = 2,              # number of cores (could use one per chain)
  control = list(         #control samplers behavior
    adapt_delta=0.9
  )
)

# Showing the estimation results of the parameters 
print(fit2, pars=c("mu", "sigma", "alpha","beta",  "lp__"), probs=c(.025,.5,.975))