dyn.load("src/stable.so")
###########################################################################
#' Density of a stable distribution
#'
#' This function produces the stable density computed at y.
#' Note that it is obtained by integrating a function from 0 to Infinity.
#' This integral is approximated by the finite integral from 0 to UP
#' using the Simpson's method with npt points or Romberg's integration. Core of the code is based on Philipe Lambert's code https://github.com/cran/stable
#'
#'
#' @param input series of input
#' @param vec vector of class "stableFit" MLE estimates from mleFit
#' @return vector of densities of input for given 4 parameters (alpha-tail, beta-skew, gamma-scale, delta-location)
#' @examples denStable(rnorm(500), alpha=2,beta=0,gamma=1/sqrt(2), delta=0)
#'
#' @useDynLib stable stable
denStable <- function(input, vec=NULL, alpha=1.8, beta=0, gamma=1, delta=0,
npt=501, up=10, eps=1.0e-6){
if(!is.null(vec) && class(vec)=="stableFit" && !any(is.na(vec$estimate)) ){
alpha=c(vec$estimate[1])
beta=c(vec$estimate[2])
gamma=c(vec$estimate[3])
delta=c(vec$estimate[4])
}else{alpha=alpha;beta=beta;gamma=gamma;delta=delta}
#special cases
# if (alpha == 2 && beta==0) {
# return(stats::dnorm(y, mean = 0, sd = 1))
# } else{ if (alpha == 1 && beta == 0) {
# return(stats::dcauchy(y, log=log))}}
ly <- length(input)
z0 <- rep(0,ly)
skew <- beta+z0
tail <- alpha+z0
yy <- (input-delta)/gamma
z <- .C("stable",
as.integer(ly),
yy,
skew,
tail,
as.integer(npt),
up,
eps,
err=integer(1),
ffy=double(ly))
z$ffy/gamma
}
###########################################################################
#' Function computes
#' @useDynLib stable stable
probStable <- function(y, loc=0,disp=1/sqrt(2),skew=0,tail=2,eps=1.0e-6){
if (alpha == 2 && beta==0) {
return(pnorm(y, mean = 0, sd = 1, log=log))
} else if (alpha == 1 && beta == 0) {
return(pcauchy(y, log=log))}
yy <- (y-loc)/disp
ly <- length(yy)
z0 <- rep(0,ly)
skew <- skew+z0
tail <- tail+z0
eta <- skew*(1-abs(1-tail))*pi/2
if ((yy==0)&(eta==0)) return(0.5+z0)
z <- .C("pstable",
as.integer(ly),
yy,
skew,
tail,
eps,
err=integer(1),
ffy=double(ly))
z$ffy}
###########################################################################
# Quantile function of a stable random variable
#
quantStable <- function(p, alpha=1.8,beta=0,gamma=1,delta=0){
if(alpha==2){return(qnorm(p))}
if(alpha==1 && beta){return(qcauchy(p))}
q<-Vectorize(stabledist::qstable,vectorize.args = c("p","alpha","beta","gamma","delta"))
q(p,alpha,beta,gamma,delta,pm=1)
}
###########################################################################
# Generation of stable random deviates
#
randStable <- function(n=1,delta=0,gamma=1,beta=0,alpha=2){
return(quantStable(runif(n),delta=delta,gamma=gamma,beta=beta,alpha=alpha))}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.