neodistr (version 0.1.2)

stanf_msnburr: Stan function of MSNBurr Distribution

Description

Stan code of MSNBurr distribution for custom distribution in stan

Usage

stanf_msnburr(vectorize = TRUE, rng = TRUE)

Value

msnburr_lpdf gives the log of density, msnburr_cdf gives the distribution function, msnburr_lcdf gives the log of distribution function, msnburr_lccdf gives the complement of log distribution function (1-msnburr_lcdf), and msnburr_rng generates random deviates.

Arguments

vectorize

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

rng

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

Author

Achmad Syahrul Choir and Nur Iriawan

Details

MSNBurr Distribution has density: $$f(y |\mu,\sigma,\alpha)=\frac{\omega}{\sigma}\exp{\left(-\omega{\left(\frac{y-\mu}{\sigma}\right)}\right)}{{\left(1+\frac{\exp{\left(-\omega{(\frac{y-\mu}{\sigma})}\right)}}{\alpha}\right)}^{-(\alpha+1)}}$$ where \(-\infty < y < \infty, -\infty < \mu< \infty, \sigma>0, \alpha>0, \omega = \frac{1}{\sqrt{2\pi}} {\left(1+\frac{1}{\alpha}\right)^{\alpha+1}}\)

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

References

Iriawan, N. (2000). Computationally Intensive Approaches to Inference in Neo-Normal Linear Models. Curtin University of Technology. 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 <- neodistr::rmsnburr(100,0,1,0.5) # random generating MSNBurr 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_msnburr(vectorize=FALSE),"}"),collapse="\n")
#define stan model code
model<-"
 data {
 int n;
 vector[n] y;
 }
 parameters {
 real mu;
 real  sigma;
 real  alpha;
 
 }
 model {
 for(i in 1:n){
 y[i]~msnburr(mu,sigma,alpha);
 }
 mu~cauchy(0,1);
 sigma~cauchy(0,2.5);
 alpha~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)
)

# Showing the estimation results of the parameters that were executed using the Stan file
print(fit1, pars=c("mu", "sigma", "alpha", "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_msnburr(vectorize=TRUE),"}"),collapse="\n")
# define stan model as vector
model_vector<-"
data {
  int n;
  vector[n] y;
}
parameters {
  real mu;
  real  sigma;
  real  alpha;
}
model {
  y~msnburr(rep_vector(mu,n),sigma,alpha);
  mu~cauchy(0,1);
  sigma~cauchy(0,2.5);
  alpha~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)
)

# Showing the estimation results of the parameters that were executed using the Stan file
print(fit2, pars=c("mu", "sigma", "alpha",  "lp__"), probs=c(.025,.5,.975))

}