neodistr (version 0.1.2)

stanf_gmsnburr: Stan function of GMSNBurr Distribution

Description

Stan code of GMSNBurr distribution for custom distribution in stan

Usage

stanf_gmsnburr(vectorize = TRUE, rng = TRUE)

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 distribution function (1-gmsnburr_lcdf), and gmsnburr_rng the stans's code of generates random numbers.

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

Author

Achmad Syahrul Choir

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 cumulative distribution, log complementary cumulative distribution, quantile, random number of GMSNBurr distribution

References

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

Examples

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