From 327a184c76b28bfea2cd6fafbf63392666ee10de Mon Sep 17 00:00:00 2001 From: Alexios Galanos Date: Thu, 22 Aug 2024 21:51:51 +0300 Subject: [PATCH] translated the c code to c++ --- .Rbuildignore | 0 .github/.gitignore | 0 .github/workflows/R-CMD-check.yaml | 0 .gitignore | 5 + DESCRIPTION | 8 +- NAMESPACE | 7 +- NEWS.md | 6 + R/RcppExports.R | 127 +++ R/distribution.R | 0 R/estimation.R | 0 R/methods.R | 0 R/moments.R | 0 R/pdqr.R | 317 +++--- R/profile.R | 0 R/sandwich.R | 0 R/spd.R | 2 +- R/specification.R | 0 R/tsdistributions-package.R | 3 +- R/utils.R | 0 README.Rmd | 0 README.md | 4 +- cran-comments.md | 4 +- inst/REFERENCES.bib | 0 man/AIC.tsdistribution.estimate.Rd | 0 man/BIC.tsdistribution.estimate.Rd | 0 man/authorized_domain.Rd | 0 man/bread.tsdistribution.estimate.Rd | 0 man/coef.tsdistribution.estimate.Rd | 0 man/ddist.Rd | 0 man/distribution_bounds.Rd | 0 man/distribution_modelspec.Rd | 0 man/dskewness.Rd | 0 man/estfun.tsdistribution.estimate.Rd | 0 man/estimate.tsdistribution.spdspec.Rd | 0 man/estimate.tsdistribution.spec.Rd | 0 man/ged.Rd | 0 man/gh.Rd | 2 +- man/ghst.Rd | 0 man/ghyp.Rd | 67 ++ man/ghyptransform.Rd | 0 man/jsu.Rd | 0 man/logLik.tsdistribution.estimate.Rd | 0 man/nig.Rd | 0 man/print.summary.tsdistribution.Rd | 0 man/print.tsdistribution.profile.Rd | 0 man/sged.Rd | 0 man/snorm.Rd | 0 man/spd.Rd | 0 man/spd_modelspec.Rd | 2 +- man/sstd.Rd | 0 man/std.Rd | 0 man/summary.tsdistribution.estimate.Rd | 0 man/summary.tsdistribution.profile.Rd | 0 man/summary.tsdistribution.spdestimate.Rd | 0 man/tsdistributions-package.Rd | 0 man/tsmoments.tsdistribution.estimate.Rd | 0 man/tsprofile.tsdistribution.spec.Rd | 0 man/vcov.tsdistribution.estimate.Rd | 0 src/Makevars | 0 src/Makevars.win | 0 src/RcppExports.cpp | 520 ++++++++++ src/TMB/compile.R | 0 src/TMB/distfun.h | 0 src/TMB/distmodel.hpp | 0 src/TMB/tsdistributions_TMBExports.cpp | 2 - src/distributions.c | 848 ---------------- src/distributions.cpp | 1009 ++++++++++++++++++++ src/distributions.h | 207 ++-- src/tsdistributions_init.c | 79 -- tests/testthat.R | 0 tests/testthat/test-liktest.R | 0 tests/testthat/test-quantile.R | 0 vignettes/custom.css | 0 vignettes/estimation_demo.Rmd | 0 vignettes/location_scale_distributions.Rmd | 0 vignettes/profile_demo.Rmd | 0 vignettes/references.bib | 0 vignettes/spd_demo.Rmd | 0 78 files changed, 1980 insertions(+), 1239 deletions(-) mode change 100644 => 100755 .Rbuildignore mode change 100644 => 100755 .github/.gitignore mode change 100644 => 100755 .github/workflows/R-CMD-check.yaml mode change 100644 => 100755 .gitignore mode change 100644 => 100755 DESCRIPTION mode change 100644 => 100755 NAMESPACE mode change 100644 => 100755 NEWS.md create mode 100644 R/RcppExports.R mode change 100644 => 100755 R/distribution.R mode change 100644 => 100755 R/estimation.R mode change 100644 => 100755 R/methods.R mode change 100644 => 100755 R/moments.R mode change 100644 => 100755 R/pdqr.R mode change 100644 => 100755 R/profile.R mode change 100644 => 100755 R/sandwich.R mode change 100644 => 100755 R/spd.R mode change 100644 => 100755 R/specification.R mode change 100644 => 100755 R/tsdistributions-package.R mode change 100644 => 100755 R/utils.R mode change 100644 => 100755 README.Rmd mode change 100644 => 100755 README.md mode change 100644 => 100755 cran-comments.md mode change 100644 => 100755 inst/REFERENCES.bib mode change 100644 => 100755 man/AIC.tsdistribution.estimate.Rd mode change 100644 => 100755 man/BIC.tsdistribution.estimate.Rd mode change 100644 => 100755 man/authorized_domain.Rd mode change 100644 => 100755 man/bread.tsdistribution.estimate.Rd mode change 100644 => 100755 man/coef.tsdistribution.estimate.Rd mode change 100644 => 100755 man/ddist.Rd mode change 100644 => 100755 man/distribution_bounds.Rd mode change 100644 => 100755 man/distribution_modelspec.Rd mode change 100644 => 100755 man/dskewness.Rd mode change 100644 => 100755 man/estfun.tsdistribution.estimate.Rd mode change 100644 => 100755 man/estimate.tsdistribution.spdspec.Rd mode change 100644 => 100755 man/estimate.tsdistribution.spec.Rd mode change 100644 => 100755 man/ged.Rd mode change 100644 => 100755 man/gh.Rd mode change 100644 => 100755 man/ghst.Rd create mode 100644 man/ghyp.Rd mode change 100644 => 100755 man/ghyptransform.Rd mode change 100644 => 100755 man/jsu.Rd mode change 100644 => 100755 man/logLik.tsdistribution.estimate.Rd mode change 100644 => 100755 man/nig.Rd mode change 100644 => 100755 man/print.summary.tsdistribution.Rd mode change 100644 => 100755 man/print.tsdistribution.profile.Rd mode change 100644 => 100755 man/sged.Rd mode change 100644 => 100755 man/snorm.Rd mode change 100644 => 100755 man/spd.Rd mode change 100644 => 100755 man/spd_modelspec.Rd mode change 100644 => 100755 man/sstd.Rd mode change 100644 => 100755 man/std.Rd mode change 100644 => 100755 man/summary.tsdistribution.estimate.Rd mode change 100644 => 100755 man/summary.tsdistribution.profile.Rd mode change 100644 => 100755 man/summary.tsdistribution.spdestimate.Rd mode change 100644 => 100755 man/tsdistributions-package.Rd mode change 100644 => 100755 man/tsmoments.tsdistribution.estimate.Rd mode change 100644 => 100755 man/tsprofile.tsdistribution.spec.Rd mode change 100644 => 100755 man/vcov.tsdistribution.estimate.Rd mode change 100644 => 100755 src/Makevars mode change 100644 => 100755 src/Makevars.win create mode 100644 src/RcppExports.cpp mode change 100644 => 100755 src/TMB/compile.R mode change 100644 => 100755 src/TMB/distfun.h mode change 100644 => 100755 src/TMB/distmodel.hpp mode change 100644 => 100755 src/TMB/tsdistributions_TMBExports.cpp delete mode 100644 src/distributions.c create mode 100644 src/distributions.cpp delete mode 100644 src/tsdistributions_init.c mode change 100644 => 100755 tests/testthat.R mode change 100644 => 100755 tests/testthat/test-liktest.R mode change 100644 => 100755 tests/testthat/test-quantile.R mode change 100644 => 100755 vignettes/custom.css mode change 100644 => 100755 vignettes/estimation_demo.Rmd mode change 100644 => 100755 vignettes/location_scale_distributions.Rmd mode change 100644 => 100755 vignettes/profile_demo.Rmd mode change 100644 => 100755 vignettes/references.bib mode change 100644 => 100755 vignettes/spd_demo.Rmd diff --git a/.Rbuildignore b/.Rbuildignore old mode 100644 new mode 100755 diff --git a/.github/.gitignore b/.github/.gitignore old mode 100644 new mode 100755 diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml old mode 100644 new mode 100755 diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 index 0005d16..a4b0800 --- a/.gitignore +++ b/.gitignore @@ -40,3 +40,8 @@ vignettes/*.pdf *.o *.so *.Rproj +tests/.DS_Store +src/TMB/.DS_Store +src/.DS_Store +.github/.DS_Store +.DS_Store diff --git a/DESCRIPTION b/DESCRIPTION old mode 100644 new mode 100755 index e5f443e..b1e64ee --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,17 +1,17 @@ Package: tsdistributions Type: Package Title: Location Scale Standardized Distributions -Version: 1.0.1 +Version: 1.0.2 Authors@R: c(person("Alexios", "Galanos", role = c("aut", "cre","cph"), email = "alexios@4dscape.com")) Maintainer: Alexios Galanos Depends: R (>= 3.5.0), methods, tsmethods -LinkingTo: TMB, RcppEigen -Imports: TMB (>= 1.7.20), Rdpack, GeneralizedHyperbolic, KernSmooth, SkewHyperbolic, mev, stats, utils, data.table, zoo, Rsolnp, sandwich, future.apply, future, progressr +LinkingTo: Rcpp, TMB, RcppEigen +Imports: Rcpp, TMB (>= 1.7.20), Rdpack, GeneralizedHyperbolic, KernSmooth, SkewHyperbolic, mev, stats, utils, data.table, zoo, Rsolnp, sandwich, future.apply, future, progressr Description: Location-Scale based distributions parameterized in terms of mean, standard deviation, skew and shape parameters and estimation using automatic differentiation. Distributions include the Normal, Student and GED as well as their skewed variants ('Fernandez and Steel'), the 'Johnson SU', and the Generalized Hyperbolic. Also included is the semi-parametric piece wise distribution ('spd') with Pareto tails and kernel interior. License: GPL-2 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 VignetteBuilder: knitr RdMacros: Rdpack URL: https://www.nopredict.com/packages/tsdistributions, https://github.com/tsmodels/tsdistributions diff --git a/NAMESPACE b/NAMESPACE old mode 100644 new mode 100755 index 3e04aa5..70e0a1d --- a/NAMESPACE +++ b/NAMESPACE @@ -28,6 +28,7 @@ export(ddist) export(dged) export(dgh) export(dghst) +export(dghyp) export(distribution_bounds) export(distribution_modelspec) export(djsu) @@ -45,6 +46,7 @@ export(pdist) export(pged) export(pgh) export(pghst) +export(pghyp) export(pjsu) export(pnig) export(psged) @@ -56,6 +58,7 @@ export(qdist) export(qged) export(qgh) export(qghst) +export(qghyp) export(qjsu) export(qnig) export(qsged) @@ -67,6 +70,7 @@ export(rdist) export(rged) export(rgh) export(rghst) +export(rghyp) export(rjsu) export(rnig) export(rsged) @@ -79,8 +83,9 @@ import(data.table) import(methods) import(tsmethods) importFrom(GeneralizedHyperbolic,ghypMom) -importFrom(GeneralizedHyperbolic,rghyp) importFrom(KernSmooth,bkde) +importFrom(Rcpp,evalCpp) +importFrom(Rcpp,sourceCpp) importFrom(Rdpack,reprompt) importFrom(Rsolnp,solnp) importFrom(SkewHyperbolic,pskewhyp) diff --git a/NEWS.md b/NEWS.md old mode 100644 new mode 100755 index 2055dd2..c62b503 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# tsdistributions 1.0.2 + +* Replaced PI with M_PI to pass strict header checks +* Translated the distribution.c c functions to Rcpp. + + # tsdistributions 1.0.1 * Added the semi-parametric piece-wise distribution (`spd`). As this is a special diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..20fa631 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,127 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +c_dghst <- function(x, mu, sigma, skew, shape, logr) { + .Call(`_tsdistributions_c_dghst`, x, mu, sigma, skew, shape, logr) +} + +c_rghst <- function(n, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_rghst`, n, mu, sigma, skew, shape) +} + +c_dghyp <- function(x, alpha, beta, delta, mu, lambda, logr) { + .Call(`_tsdistributions_c_dghyp`, x, alpha, beta, delta, mu, lambda, logr) +} + +c_dgh <- function(x, mu, sigma, skew, shape, lambda, logr) { + .Call(`_tsdistributions_c_dgh`, x, mu, sigma, skew, shape, lambda, logr) +} + +c_dnig <- function(x, mu, sigma, skew, shape, logr) { + .Call(`_tsdistributions_c_dnig`, x, mu, sigma, skew, shape, logr) +} + +c_rstd <- function(n, mu, sigma, shape) { + .Call(`_tsdistributions_c_rstd`, n, mu, sigma, shape) +} + +c_dstd <- function(x, mu, sigma, shape, logr) { + .Call(`_tsdistributions_c_dstd`, x, mu, sigma, shape, logr) +} + +c_pstd <- function(q, mu, sigma, shape) { + .Call(`_tsdistributions_c_pstd`, q, mu, sigma, shape) +} + +c_qstd <- function(p, mu, sigma, shape) { + .Call(`_tsdistributions_c_qstd`, p, mu, sigma, shape) +} + +c_rsstd <- function(n, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_rsstd`, n, mu, sigma, skew, shape) +} + +c_dsstd <- function(x, mu, sigma, skew, shape, logr) { + .Call(`_tsdistributions_c_dsstd`, x, mu, sigma, skew, shape, logr) +} + +c_psstd <- function(q, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_psstd`, q, mu, sigma, skew, shape) +} + +c_qsstd <- function(p, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_qsstd`, p, mu, sigma, skew, shape) +} + +c_djsu <- function(x, mu, sigma, skew, shape, logr) { + .Call(`_tsdistributions_c_djsu`, x, mu, sigma, skew, shape, logr) +} + +c_qjsu <- function(p, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_qjsu`, p, mu, sigma, skew, shape) +} + +c_pjsu <- function(q, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_pjsu`, q, mu, sigma, skew, shape) +} + +c_rjsu <- function(n, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_rjsu`, n, mu, sigma, skew, shape) +} + +c_rsnorm <- function(n, mu, sigma, skew) { + .Call(`_tsdistributions_c_rsnorm`, n, mu, sigma, skew) +} + +c_dsnorm <- function(x, mu, sigma, skew, logr) { + .Call(`_tsdistributions_c_dsnorm`, x, mu, sigma, skew, logr) +} + +c_psnorm <- function(q, mu, sigma, skew) { + .Call(`_tsdistributions_c_psnorm`, q, mu, sigma, skew) +} + +c_qsnorm <- function(p, mu, sigma, skew) { + .Call(`_tsdistributions_c_qsnorm`, p, mu, sigma, skew) +} + +c_rged <- function(n, mu, sigma, shape) { + .Call(`_tsdistributions_c_rged`, n, mu, sigma, shape) +} + +c_dged <- function(x, mu, sigma, shape, logr) { + .Call(`_tsdistributions_c_dged`, x, mu, sigma, shape, logr) +} + +c_pged <- function(q, mu, sigma, shape) { + .Call(`_tsdistributions_c_pged`, q, mu, sigma, shape) +} + +c_qged <- function(p, mu, sigma, shape) { + .Call(`_tsdistributions_c_qged`, p, mu, sigma, shape) +} + +c_rsged <- function(n, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_rsged`, n, mu, sigma, skew, shape) +} + +c_dsged <- function(x, mu, sigma, skew, shape, logr) { + .Call(`_tsdistributions_c_dsged`, x, mu, sigma, skew, shape, logr) +} + +c_psged <- function(q, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_psged`, q, mu, sigma, skew, shape) +} + +c_qsged <- function(p, mu, sigma, skew, shape) { + .Call(`_tsdistributions_c_qsged`, p, mu, sigma, skew, shape) +} + +c_dhyp <- function(x, mu, sigma, skew, shape, logr) { + .Call(`_tsdistributions_c_dhyp`, x, mu, sigma, skew, shape, logr) +} + +c_rghyp <- function(n, mu = 0, delta = 1, alpha = 1, beta = 0, lambda = 1) { + .Call(`_tsdistributions_c_rghyp`, n, mu, delta, alpha, beta, lambda) +} + diff --git a/R/distribution.R b/R/distribution.R old mode 100644 new mode 100755 diff --git a/R/estimation.R b/R/estimation.R old mode 100644 new mode 100755 diff --git a/R/methods.R b/R/methods.R old mode 100644 new mode 100755 diff --git a/R/moments.R b/R/moments.R old mode 100644 new mode 100755 diff --git a/R/pdqr.R b/R/pdqr.R old mode 100644 new mode 100755 index ace1dc2..d376559 --- a/R/pdqr.R +++ b/R/pdqr.R @@ -48,15 +48,13 @@ dghst <- function(x, mu = 0, sigma = 1, skew = 1, shape = 8, log = FALSE) if (val_length[3] != max_n) sigma <- rep(sigma[1], max_n) if (val_length[4] != max_n) skew <- rep(skew[1], max_n) if (val_length[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_dghst", x = as.double(x), mu = as.double(mu), + sol <- try(c_dghst(x = as.double(x), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - logr = as.integer(log), PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape), logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -72,15 +70,13 @@ rghst <- function(n, mu = 0, sigma = 1, skew = 1, shape = 8) if (val_length[2] != n) sigma <- rep(sigma[1], n) if (val_length[3] != n) skew <- rep(skew[1], n) if (val_length[4] != n) shape <- rep(shape[1], n) - ans <- double(n) - sol <- try(.C("c_rghst", n = as.integer(n), mu = as.double(mu), + sol <- try(c_rghst(n = as.integer(n), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, - PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -185,15 +181,13 @@ dsnorm <- function(x, mu = 0, sigma = 1, skew = 1.5, log = FALSE) if (val_length[2] != max_n) mu <- rep(mu[1], max_n) if (val_length[3] != max_n) sigma <- rep(sigma[1], max_n) if (val_length[4] != max_n) skew <- rep(skew[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_dsnorm", x = as.double(x), mu = as.double(mu), - sigma = as.double(sigma), skew = as.double(skew), ans = ans, - n = as.integer(max_n), logr = as.integer(log), - PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_dsnorm(x = as.double(x), mu = as.double(mu), + sigma = as.double(sigma), skew = as.double(skew), + logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -207,16 +201,14 @@ psnorm <- function(q, mu = 0, sigma = 1, skew = 1.5, lower_tail = TRUE, log = FA if (val_length[2] != max_n) mu <- rep(mu[1], max_n) if (val_length[3] != max_n) sigma <- rep(sigma[1], max_n) if (val_length[4] != max_n) skew <- rep(skew[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_psnorm", q = as.double(q), mu = as.double(mu), - sigma = as.double(sigma), skew = as.double(skew), ans = ans, - n = as.integer(max_n), PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_psnorm(q = as.double(q), mu = as.double(mu), + sigma = as.double(sigma), skew = as.double(skew)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (!lower_tail) sol$ans <- 1 - sol$ans - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (!lower_tail) sol <- 1 - sol + if (log) sol <- log(sol) + return(sol) } } @@ -233,15 +225,13 @@ qsnorm <- function(p, mu = 0, sigma = 1, skew = 1.5, lower_tail = TRUE, log = FA if (val_length[2] != max_n) mu <- rep(mu[1], max_n) if (val_length[3] != max_n) sigma <- rep(sigma[1], max_n) if (val_length[4] != max_n) skew <- rep(skew[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_qsnorm", p = as.double(p), mu = as.double(mu), - sigma = as.double(sigma), skew = as.double(skew), ans = ans, - n = as.integer(max_n), PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_qsnorm(p = as.double(p), mu = as.double(mu), + sigma = as.double(sigma), skew = as.double(skew)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (log) sol <- log(sol) + return(sol) } } @@ -253,14 +243,12 @@ rsnorm <- function(n, mu = 0, sigma = 1, skew = 1.5) if (val_length[1] != n) mu <- rep(mu[1], n) if (val_length[2] != n) sigma <- rep(sigma[1], n) if (val_length[3] != n) skew <- rep(skew[1], n) - ans <- double(n) - sol <- try(.C("c_rsnorm", n = as.integer(n), mu = as.double(mu), - sigma = as.double(sigma), skew = as.double(skew), ans = ans, - PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_rsnorm(n = as.integer(n), mu = as.double(mu), + sigma = as.double(sigma), skew = as.double(skew)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -293,15 +281,13 @@ dged <- function(x, mu = 0, sigma = 1, shape = 2, log = FALSE) if (value_len[2] != max_n) mu <- rep(mu[1], max_n) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_dged", x = as.double(x), mu = as.double(mu), + sol <- try(c_dged(x = as.double(x), mu = as.double(mu), sigma = as.double(sigma), shape = as.double(shape), - ans = ans, n = as.integer(max_n), logr = as.integer(log), - PACKAGE = "tsdistributions"), silent = TRUE) + logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -315,17 +301,14 @@ pged <- function(q, mu = 0, sigma = 1, shape = 2, lower_tail = TRUE, log = FALSE if (value_len[2] != max_n) mu <- rep(mu[1], max_n) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_pged", q = as.double(q), mu = as.double(mu), - sigma = as.double(sigma), shape = as.double(shape), - ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_pged(q = as.double(q), mu = as.double(mu), + sigma = as.double(sigma), shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (!lower_tail) sol$ans <- 1 - sol$ans - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (!lower_tail) sol <- 1 - sol + if (log) sol <- log(sol) + return(sol) } } @@ -340,16 +323,13 @@ qged <- function(p, mu = 0, sigma = 1, shape = 2, lower_tail = TRUE, log = FALSE if (value_len[2] != max_n) mu <- rep(mu[1], max_n) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_qged", p = as.double(p), mu = as.double(mu), - sigma = as.double(sigma), shape = as.double(shape), - ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_qged(p = as.double(p), mu = as.double(mu), + sigma = as.double(sigma), shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (log) sol <- log(sol) + return(sol) } } @@ -361,14 +341,12 @@ rged <- function(n, mu = 0, sigma = 1, shape = 2) if (value_len[1] != n) mu <- rep(mu[1], n) if (value_len[2] != n) sigma <- rep(sigma[1], n) if (value_len[3] != n) shape <- rep(shape[1], n) - ans <- double(n) - sol <- try(.C("c_rged", n = as.integer(n), mu = as.double(mu), - sigma = as.double(sigma), shape = as.double(shape), - ans = ans, PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_rged(n = as.integer(n), mu = as.double(mu), + sigma = as.double(sigma), shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -403,15 +381,13 @@ dsged <- function(x, mu = 0, sigma = 1, skew = 1.5, shape = 2, log = FALSE) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_dsged", x = as.double(x), mu = as.double(mu), + sol <- try(c_dsged(x = as.double(x), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - logr = as.integer(log), PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape), logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -426,17 +402,15 @@ psged <- function(q, mu = 0, sigma = 1, skew = 1.5, shape = 2, lower_tail = TRUE if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_psged", q = as.double(q), mu = as.double(mu), + sol <- try(c_psged(q = as.double(q), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, - n = as.integer(max_n), PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (!lower_tail) sol$ans <- 1 - sol$ans - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (!lower_tail) sol <- 1 - sol + if (log) sol <- log(sol) + return(sol) } } @@ -452,16 +426,14 @@ qsged <- function(p, mu = 0, sigma = 1, skew = 1.5, shape = 2, lower_tail = TRUE if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol = try(.C("c_qsged", p = as.double(p), mu = as.double(mu), + sol = try(c_qsged(p = as.double(p), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (log) sol <- log(sol) + return(sol) } } @@ -474,15 +446,13 @@ rsged <- function(n, mu = 0, sigma = 1, skew = 1.5, shape = 2) if (value_len[2] != n) sigma <- rep(sigma[1], n) if (value_len[3] != n) skew <- rep(skew[1], n) if (value_len[4] != n) shape <- rep(shape[1], n) - ans <- double(n) - sol <- try(.C("c_rsged", n = as.integer(n), mu = as.double(mu), + sol <- try(c_rsged(n = as.integer(n), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), - ans = ans, PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -521,15 +491,13 @@ dstd <- function(x, mu = 0, sigma = 1, shape = 5, log = FALSE) if (value_len[2] != max_n) mu <- rep(mu[1], max_n) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_dstd", x = as.double(x), mu = as.double(mu), + sol <- try(c_dstd(x = as.double(x), mu = as.double(mu), sigma = as.double(sigma), shape = as.double(shape), - ans = ans, n = as.integer(max_n), logr = as.integer(log), - PACKAGE = "tsdistributions"), silent = TRUE) + logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -543,17 +511,14 @@ pstd <- function(q, mu = 0, sigma = 1, shape = 5, lower_tail = TRUE, log = FALSE if (value_len[2] != max_n) mu <- rep(mu[1], max_n) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_pstd", q = as.double(q), mu = as.double(mu), - sigma = as.double(sigma), shape = as.double(shape), - ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_pstd(q = as.double(q), mu = as.double(mu), + sigma = as.double(sigma), shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (!lower_tail) sol$ans <- 1 - sol$ans - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (!lower_tail) sol <- 1 - sol + if (log) sol <- log(sol) + return(sol) } } @@ -568,16 +533,13 @@ qstd <- function(p, mu = 0, sigma = 1, shape = 5, lower_tail = TRUE, log = FALSE if (value_len[2] != max_n) mu <- rep(mu[1], max_n) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_qstd", p = as.double(p), mu = as.double(mu), - sigma = as.double(sigma), shape = as.double(shape), - ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_qstd(p = as.double(p), mu = as.double(mu), + sigma = as.double(sigma), shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (log) sol <- log(sol) + return(sol) } } @@ -589,14 +551,12 @@ rstd <- function(n, mu = 0, sigma = 1, shape = 5) if (value_len[1] != n) mu <- rep(mu[1], n) if (value_len[2] != n) sigma <- rep(sigma[1], n) if (value_len[3] != n) shape <- rep(shape[1], n) - ans <- double(n) - sol <- try(.C("c_rstd", n = as.integer(n), mu = as.double(mu), - sigma = as.double(sigma), shape = as.double(shape), - ans = ans, PACKAGE = "tsdistributions"), silent = TRUE) + sol <- try(c_rstd(n = as.integer(n), mu = as.double(mu), + sigma = as.double(sigma), shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -631,15 +591,13 @@ dsstd <- function(x, mu = 0, sigma = 1, skew = 1.5, shape = 5, log = FALSE) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_dsstd", x = as.double(x), mu = as.double(mu), + sol <- try(c_dsstd(x = as.double(x), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - logr = as.integer(log), PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape), logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -654,17 +612,15 @@ psstd <- function(q, mu = 0, sigma = 1, skew = 1.5, shape = 5, lower_tail = TRUE if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol = try(.C("c_psstd", q = as.double(q), mu = as.double(mu), + sol = try(c_psstd(q = as.double(q), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (!lower_tail) sol$ans <- 1 - sol$ans - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (!lower_tail) sol <- 1 - sol + if (log) sol <- log(sol) + return(sol) } } @@ -680,16 +636,14 @@ qsstd <- function(p, mu = 0, sigma = 1, skew = 1.5, shape = 5, lower_tail = TRUE if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_qsstd", p = as.double(p), mu = as.double(mu), + sol <- try(c_qsstd(p = as.double(p), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (log) sol <- log(sol) + return(sol) } } @@ -702,15 +656,13 @@ rsstd <- function(n, mu = 0, sigma = 1, skew = 1.5, shape = 5) if (value_len[2] != n) sigma <- rep(sigma[1], n) if (value_len[3] != n) skew <- rep(skew[1], n) if (value_len[4] != n) shape <- rep(shape[1], n) - ans <- double(n) - sol <- try(.C("c_rsstd", n = as.integer(n), mu = as.double(mu), + sol <- try(c_rsstd(n = as.integer(n), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), - ans = ans, PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } # ------------------------------------------------------------------------------ @@ -762,6 +714,31 @@ rsstd <- function(n, mu = 0, sigma = 1, skew = 1.5, shape = 5) return(c(alpha = alpha, beta = beta, delta = delta, mu = mu)) } + +#' Generalized Hyperbolic Distribution (alpha-beta-delta-mu parameterization) +#' +#' @description Density, distribution, quantile function and random number +#' generation for the generalized hyperbolic distribution +#' using the alpha-beta-delta-mu-lambda parameterization. +#' @param x,q vector of quantiles. +#' @param p vector of probabilities. +#' @param n number of observations. +#' @param alpha tail parameter. +#' @param beta skewness parameter. +#' @param delta scale parameter. +#' @param mu location parameter. +#' @param lambda additional shape parameter determining subfamilies of this +#' distributions. +#' @param log (logical) if TRUE, probabilities p are given as log(p). +#' @param lower_tail if TRUE (default), probabilities are \eqn{P[X \le x]} otherwise, \eqn{P[X > x]}. +#' @return d gives the density, p gives the distribution function, q gives the quantile function +#' and r generates random deviates. Output depends on x or q length, or n for the random number +#' generator. +#' @rdname ghyp +#' @export +#' +#' +#' dghyp <- function(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1, log = FALSE) { value_len <- c(length(x), length(alpha), length(beta), length(delta), length(mu), length(lambda)) @@ -775,19 +752,19 @@ dghyp <- function(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1, log = F if (any(alpha <= 0)) stop("alpha must be greater than zero") if (any(delta <= 0)) stop("delta must be greater than zero") if (any(abs(beta) >= alpha)) stop("abs value of beta must be less than alpha") - ans <- double(max_n) - sol <- try(.C("c_dghyp", x = as.double(x), alpha = as.double(alpha), + sol <- try(c_dghyp(x = as.double(x), alpha = as.double(alpha), beta = as.double(beta), delta = as.double(delta), mu = as.double(mu), lambda = as.double(lambda), - ans = ans, n = as.integer(max_n), logr = as.integer(log), - PACKAGE = "tsdistributions"), silent = TRUE) + logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } +#' @rdname ghyp +#' @export pghyp <- function(q, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1, lower_tail = TRUE, log = FALSE) { value_len <- c(length(q), length(alpha), length(beta), length(delta), length(mu), length(lambda)) @@ -813,6 +790,8 @@ pghyp <- function(q, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1, lower_t return(ans) } +#' @rdname ghyp +#' @export qghyp <- function(p, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1, lower_tail = TRUE, log = FALSE) { if (!lower_tail) p <- 1 - p @@ -853,8 +832,17 @@ qghyp <- function(p, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1, lower_t return(ans) } +#' @rdname ghyp +#' @export +rghyp <- function(n, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1) +{ + if (alpha <= 0) stop("alpha must be greater than zero") + if (delta <= 0) stop("delta must be greater than zero") + if (abs(beta) >= alpha) stop("abs value of beta must be less than alpha") + return(c_rghyp(n, mu = mu[1], delta = delta[1], alpha = alpha[1], beta = beta[1], lambda = lambda[1])) +} -#' Generalized Hyperbolic Distribution +#' Generalized Hyperbolic Distribution (rho-zeta parameterization) #' #' @description Density, distribution, quantile function and random number #' generation for the generalized hyperbolic distribution parameterized in @@ -889,16 +877,14 @@ dgh <- function(x, mu = 0, sigma = 1, skew = 0, shape = 1, lambda = 1, log = FAL if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) if (value_len[6] != max_n) lambda <- rep(lambda[1], max_n) - ans <- double(max_n) - sol = try(.C("c_dgh", x = as.double(x), mu = as.double(mu), + sol = try(c_dgh(x = as.double(x), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), lambda = as.double(lambda), - ans = ans, n = as.integer(max_n), logr = as.integer(log), - PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape), lambda = as.double(lambda), + logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -942,7 +928,7 @@ qgh <- function(p, mu = 0, sigma = 1, skew = 0, shape = 1, lambda = 1, lower_tai rgh <- function(n, mu = 0, sigma = 1, skew = 0, shape = 1, lambda = 1) { params <- .paramGH(rho = skew, zeta = shape, lambda = lambda) - out <- rghyp(n = n, mu = params[4], delta = params[3], alpha = params[1], beta = params[2], lambda = lambda) + out <- rghyp(n = n, mu = params[4], delta = params[3], alpha = params[1], beta = params[2], lambda = lambda[1]) out <- mu + sigma * out return(out) } @@ -982,15 +968,13 @@ dnig <- function(x, mu = 0, sigma = 1, skew = 0, shape = 1, log = FALSE) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_dnig", x = as.double(x), mu = as.double(mu), + sol <- try(c_dnig(x = as.double(x), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - logr = as.integer(log), PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape), logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -1047,15 +1031,13 @@ djsu <- function(x, mu = 0, sigma = 1, skew = 1, shape = 0.5, log = FALSE) if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_djsu", x = as.double(x), mu = as.double(mu), + sol <- try(c_djsu(x = as.double(x), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - logr = as.integer(log), PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape), logr = as.integer(log)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } @@ -1071,16 +1053,15 @@ pjsu <- function(q, mu = 0, sigma = 1, skew = 1, shape = 0.5, lower_tail = TRUE, if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) ans <- double(max_n) - sol <- try(.C("c_pjsu", q = as.double(q), mu = as.double(mu), + sol <- try(c_pjsu(q = as.double(q), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (!lower_tail) sol$ans <- 1 - sol$ans - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (!lower_tail) sol <- 1 - sol + if (log) sol <- log(sol) + return(sol) } } @@ -1096,16 +1077,14 @@ qjsu <- function(p, mu = 0, sigma = 1, skew = 1, shape = 0.5, lower_tail = TRUE, if (value_len[3] != max_n) sigma <- rep(sigma[1], max_n) if (value_len[4] != max_n) skew <- rep(skew[1], max_n) if (value_len[5] != max_n) shape <- rep(shape[1], max_n) - ans <- double(max_n) - sol <- try(.C("c_qjsu", p = as.double(p), mu = as.double(mu), + sol <- try(c_qjsu(p = as.double(p), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), ans = ans, n = as.integer(max_n), - PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - if (log) sol$ans <- log(sol$ans) - return(sol$ans) + if (log) sol <- log(sol) + return(sol) } } @@ -1118,15 +1097,13 @@ rjsu <- function(n, mu = 0, sigma = 1, skew = 1, shape = 0.5) if (value_len[2] != n) sigma <- rep(sigma[1], n) if (value_len[3] != n) skew <- rep(skew[1], n) if (value_len[4] != n) shape <- rep(shape[1], n) - ans <- double(n) - sol <- try(.C("c_rjsu", n = as.integer(n), mu = as.double(mu), + sol <- try(c_rjsu(n = as.integer(n), mu = as.double(mu), sigma = as.double(sigma), skew = as.double(skew), - shape = as.double(shape), - ans = ans, PACKAGE = "tsdistributions"), silent = TRUE) + shape = as.double(shape)), silent = TRUE) if (inherits(sol, 'try-error')) { return(sol) } else { - return(sol$ans) + return(sol) } } diff --git a/R/profile.R b/R/profile.R old mode 100644 new mode 100755 diff --git a/R/sandwich.R b/R/sandwich.R old mode 100644 new mode 100755 diff --git a/R/spd.R b/R/spd.R old mode 100644 new mode 100755 index 63e58cc..0cec34a --- a/R/spd.R +++ b/R/spd.R @@ -3,7 +3,7 @@ #' @param y a numeric vector #' @param lower the probability for the lower GPD tail. #' @param upper the probability for the upper GPD tail. -#' @param kernel_type the choice of the kernel to use from the \code{\link{bkde}} +#' @param kernel_type the choice of the kernel to use from the \code{\link[KernSmooth]{bkde}} #' function. #' @param ... not currently used #' @returns An object of class \dQuote{tsdistribution.spd_spec}. diff --git a/R/specification.R b/R/specification.R old mode 100644 new mode 100755 diff --git a/R/tsdistributions-package.R b/R/tsdistributions-package.R old mode 100644 new mode 100755 index f3afb53..776c53d --- a/R/tsdistributions-package.R +++ b/R/tsdistributions-package.R @@ -5,13 +5,14 @@ #' @import methods #' @importFrom zoo coredata #' @importFrom SkewHyperbolic pskewhyp qskewhyp -#' @importFrom GeneralizedHyperbolic ghypMom rghyp +#' @importFrom GeneralizedHyperbolic ghypMom #' @importFrom KernSmooth bkde #' @importFrom Rsolnp solnp #' @importFrom sandwich estfun bwNeweyWest vcovHAC vcovOPG bread #' @importFrom stats dnorm integrate optim pnorm qnorm rnorm runif sd uniroot spline nlminb na.omit coef logLik printCoefmat ar lm logLik vcov AIC BIC ppoints approx rexp var #' @importFrom future.apply future_lapply #' @importFrom future %<-% +#' @importFrom Rcpp evalCpp sourceCpp #' @importFrom progressr handlers progressor #' @importFrom mev fit.gpd dgp pgp qgp rgp #' @importFrom utils tail diff --git a/R/utils.R b/R/utils.R old mode 100644 new mode 100755 diff --git a/README.Rmd b/README.Rmd old mode 100644 new mode 100755 diff --git a/README.md b/README.md old mode 100644 new mode 100755 index 6a897d4..1762116 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ [![R-CMD-check](https://github.com/tsmodels/tsdistributions/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/tsmodels/tsdistributions/actions/workflows/R-CMD-check.yaml) -[![Last-changedate](https://img.shields.io/badge/last%20change-2024--05--10-yellowgreen.svg)](/commits/master) -[![packageversion](https://img.shields.io/badge/Package%20version-1.0.1-orange.svg?style=flat-square)](commits/master) +[![Last-changedate](https://img.shields.io/badge/last%20change-2024--08--22-yellowgreen.svg)](/commits/master) +[![packageversion](https://img.shields.io/badge/Package%20version-1.0.2-orange.svg?style=flat-square)](commits/master) [![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/tsdistributions)](https://cran.r-project.org/package=tsdistributions) # tsdistributions diff --git a/cran-comments.md b/cran-comments.md old mode 100644 new mode 100755 index b3b6783..3da4835 --- a/cran-comments.md +++ b/cran-comments.md @@ -2,5 +2,5 @@ 0 errors | 0 warnings | 1 note -* Added a new distribution (`spd`). -* Fixed 'noRemap' in Additional Issues. +* Replaced PI to M_PI (strict headers warning) +* Translated c code to c++ (Rcpp) diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib old mode 100644 new mode 100755 diff --git a/man/AIC.tsdistribution.estimate.Rd b/man/AIC.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/man/BIC.tsdistribution.estimate.Rd b/man/BIC.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/man/authorized_domain.Rd b/man/authorized_domain.Rd old mode 100644 new mode 100755 diff --git a/man/bread.tsdistribution.estimate.Rd b/man/bread.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/man/coef.tsdistribution.estimate.Rd b/man/coef.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/man/ddist.Rd b/man/ddist.Rd old mode 100644 new mode 100755 diff --git a/man/distribution_bounds.Rd b/man/distribution_bounds.Rd old mode 100644 new mode 100755 diff --git a/man/distribution_modelspec.Rd b/man/distribution_modelspec.Rd old mode 100644 new mode 100755 diff --git a/man/dskewness.Rd b/man/dskewness.Rd old mode 100644 new mode 100755 diff --git a/man/estfun.tsdistribution.estimate.Rd b/man/estfun.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/man/estimate.tsdistribution.spdspec.Rd b/man/estimate.tsdistribution.spdspec.Rd old mode 100644 new mode 100755 diff --git a/man/estimate.tsdistribution.spec.Rd b/man/estimate.tsdistribution.spec.Rd old mode 100644 new mode 100755 diff --git a/man/ged.Rd b/man/ged.Rd old mode 100644 new mode 100755 diff --git a/man/gh.Rd b/man/gh.Rd old mode 100644 new mode 100755 index cf08ec7..b8cf0b6 --- a/man/gh.Rd +++ b/man/gh.Rd @@ -5,7 +5,7 @@ \alias{pgh} \alias{qgh} \alias{rgh} -\title{Generalized Hyperbolic Distribution} +\title{Generalized Hyperbolic Distribution (rho-zeta parameterization)} \usage{ dgh(x, mu = 0, sigma = 1, skew = 0, shape = 1, lambda = 1, log = FALSE) diff --git a/man/ghst.Rd b/man/ghst.Rd old mode 100644 new mode 100755 diff --git a/man/ghyp.Rd b/man/ghyp.Rd new file mode 100644 index 0000000..2c5df07 --- /dev/null +++ b/man/ghyp.Rd @@ -0,0 +1,67 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pdqr.R +\name{dghyp} +\alias{dghyp} +\alias{pghyp} +\alias{qghyp} +\alias{rghyp} +\title{Generalized Hyperbolic Distribution (alpha-beta-delta-mu parameterization)} +\usage{ +dghyp(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1, log = FALSE) + +pghyp( + q, + alpha = 1, + beta = 0, + delta = 1, + mu = 0, + lambda = 1, + lower_tail = TRUE, + log = FALSE +) + +qghyp( + p, + alpha = 1, + beta = 0, + delta = 1, + mu = 0, + lambda = 1, + lower_tail = TRUE, + log = FALSE +) + +rghyp(n, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1) +} +\arguments{ +\item{x, q}{vector of quantiles.} + +\item{alpha}{tail parameter.} + +\item{beta}{skewness parameter.} + +\item{delta}{scale parameter.} + +\item{mu}{location parameter.} + +\item{lambda}{additional shape parameter determining subfamilies of this +distributions.} + +\item{log}{(logical) if TRUE, probabilities p are given as log(p).} + +\item{lower_tail}{if TRUE (default), probabilities are \eqn{P[X \le x]} otherwise, \eqn{P[X > x]}.} + +\item{p}{vector of probabilities.} + +\item{n}{number of observations.} +} +\value{ +d gives the density, p gives the distribution function, q gives the quantile function +and r generates random deviates. Output depends on x or q length, or n for the random number +generator. +} +\description{ +Density, distribution, quantile function and random number +generation for the generalized hyperbolic distribution +using the alpha-beta-delta-mu-lambda parameterization. +} diff --git a/man/ghyptransform.Rd b/man/ghyptransform.Rd old mode 100644 new mode 100755 diff --git a/man/jsu.Rd b/man/jsu.Rd old mode 100644 new mode 100755 diff --git a/man/logLik.tsdistribution.estimate.Rd b/man/logLik.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/man/nig.Rd b/man/nig.Rd old mode 100644 new mode 100755 diff --git a/man/print.summary.tsdistribution.Rd b/man/print.summary.tsdistribution.Rd old mode 100644 new mode 100755 diff --git a/man/print.tsdistribution.profile.Rd b/man/print.tsdistribution.profile.Rd old mode 100644 new mode 100755 diff --git a/man/sged.Rd b/man/sged.Rd old mode 100644 new mode 100755 diff --git a/man/snorm.Rd b/man/snorm.Rd old mode 100644 new mode 100755 diff --git a/man/spd.Rd b/man/spd.Rd old mode 100644 new mode 100755 diff --git a/man/spd_modelspec.Rd b/man/spd_modelspec.Rd old mode 100644 new mode 100755 index 5636318..939f154 --- a/man/spd_modelspec.Rd +++ b/man/spd_modelspec.Rd @@ -19,7 +19,7 @@ spd_modelspec( \item{upper}{the probability for the upper GPD tail.} -\item{kernel_type}{the choice of the kernel to use from the \code{\link{bkde}} +\item{kernel_type}{the choice of the kernel to use from the \code{\link[KernSmooth]{bkde}} function.} \item{...}{not currently used} diff --git a/man/sstd.Rd b/man/sstd.Rd old mode 100644 new mode 100755 diff --git a/man/std.Rd b/man/std.Rd old mode 100644 new mode 100755 diff --git a/man/summary.tsdistribution.estimate.Rd b/man/summary.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/man/summary.tsdistribution.profile.Rd b/man/summary.tsdistribution.profile.Rd old mode 100644 new mode 100755 diff --git a/man/summary.tsdistribution.spdestimate.Rd b/man/summary.tsdistribution.spdestimate.Rd old mode 100644 new mode 100755 diff --git a/man/tsdistributions-package.Rd b/man/tsdistributions-package.Rd old mode 100644 new mode 100755 diff --git a/man/tsmoments.tsdistribution.estimate.Rd b/man/tsmoments.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/man/tsprofile.tsdistribution.spec.Rd b/man/tsprofile.tsdistribution.spec.Rd old mode 100644 new mode 100755 diff --git a/man/vcov.tsdistribution.estimate.Rd b/man/vcov.tsdistribution.estimate.Rd old mode 100644 new mode 100755 diff --git a/src/Makevars b/src/Makevars old mode 100644 new mode 100755 diff --git a/src/Makevars.win b/src/Makevars.win old mode 100644 new mode 100755 diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..a10dd09 --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,520 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include +#include + +using namespace Rcpp; + +#ifdef RCPP_USE_GLOBAL_ROSTREAM +Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +#endif + +// c_dghst +NumericVector c_dghst(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr); +RcppExport SEXP _tsdistributions_c_dghst(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dghst(x, mu, sigma, skew, shape, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_rghst +NumericVector c_rghst(int n, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_rghst(SEXP nSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_rghst(n, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_dghyp +NumericVector c_dghyp(NumericVector x, NumericVector alpha, NumericVector beta, NumericVector delta, NumericVector mu, NumericVector lambda, int logr); +RcppExport SEXP _tsdistributions_c_dghyp(SEXP xSEXP, SEXP alphaSEXP, SEXP betaSEXP, SEXP deltaSEXP, SEXP muSEXP, SEXP lambdaSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type alpha(alphaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type beta(betaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type lambda(lambdaSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dghyp(x, alpha, beta, delta, mu, lambda, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_dgh +NumericVector c_dgh(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, NumericVector lambda, int logr); +RcppExport SEXP _tsdistributions_c_dgh(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP, SEXP lambdaSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< NumericVector >::type lambda(lambdaSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dgh(x, mu, sigma, skew, shape, lambda, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_dnig +NumericVector c_dnig(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr); +RcppExport SEXP _tsdistributions_c_dnig(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dnig(x, mu, sigma, skew, shape, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_rstd +NumericVector c_rstd(int n, NumericVector mu, NumericVector sigma, NumericVector shape); +RcppExport SEXP _tsdistributions_c_rstd(SEXP nSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_rstd(n, mu, sigma, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_dstd +NumericVector c_dstd(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector shape, int logr); +RcppExport SEXP _tsdistributions_c_dstd(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP shapeSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dstd(x, mu, sigma, shape, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_pstd +NumericVector c_pstd(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector shape); +RcppExport SEXP _tsdistributions_c_pstd(SEXP qSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type q(qSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_pstd(q, mu, sigma, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_qstd +NumericVector c_qstd(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector shape); +RcppExport SEXP _tsdistributions_c_qstd(SEXP pSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type p(pSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_qstd(p, mu, sigma, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_rsstd +NumericVector c_rsstd(int n, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_rsstd(SEXP nSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_rsstd(n, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_dsstd +NumericVector c_dsstd(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr); +RcppExport SEXP _tsdistributions_c_dsstd(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dsstd(x, mu, sigma, skew, shape, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_psstd +NumericVector c_psstd(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_psstd(SEXP qSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type q(qSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_psstd(q, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_qsstd +NumericVector c_qsstd(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_qsstd(SEXP pSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type p(pSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_qsstd(p, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_djsu +NumericVector c_djsu(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr); +RcppExport SEXP _tsdistributions_c_djsu(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_djsu(x, mu, sigma, skew, shape, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_qjsu +NumericVector c_qjsu(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_qjsu(SEXP pSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type p(pSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_qjsu(p, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_pjsu +NumericVector c_pjsu(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_pjsu(SEXP qSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type q(qSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_pjsu(q, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_rjsu +NumericVector c_rjsu(int n, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_rjsu(SEXP nSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_rjsu(n, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_rsnorm +NumericVector c_rsnorm(int n, NumericVector mu, NumericVector sigma, NumericVector skew); +RcppExport SEXP _tsdistributions_c_rsnorm(SEXP nSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + rcpp_result_gen = Rcpp::wrap(c_rsnorm(n, mu, sigma, skew)); + return rcpp_result_gen; +END_RCPP +} +// c_dsnorm +NumericVector c_dsnorm(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, int logr); +RcppExport SEXP _tsdistributions_c_dsnorm(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dsnorm(x, mu, sigma, skew, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_psnorm +NumericVector c_psnorm(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector skew); +RcppExport SEXP _tsdistributions_c_psnorm(SEXP qSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type q(qSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + rcpp_result_gen = Rcpp::wrap(c_psnorm(q, mu, sigma, skew)); + return rcpp_result_gen; +END_RCPP +} +// c_qsnorm +NumericVector c_qsnorm(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector skew); +RcppExport SEXP _tsdistributions_c_qsnorm(SEXP pSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type p(pSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + rcpp_result_gen = Rcpp::wrap(c_qsnorm(p, mu, sigma, skew)); + return rcpp_result_gen; +END_RCPP +} +// c_rged +NumericVector c_rged(int n, NumericVector mu, NumericVector sigma, NumericVector shape); +RcppExport SEXP _tsdistributions_c_rged(SEXP nSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_rged(n, mu, sigma, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_dged +NumericVector c_dged(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector shape, int logr); +RcppExport SEXP _tsdistributions_c_dged(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP shapeSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dged(x, mu, sigma, shape, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_pged +NumericVector c_pged(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector shape); +RcppExport SEXP _tsdistributions_c_pged(SEXP qSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type q(qSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_pged(q, mu, sigma, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_qged +NumericVector c_qged(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector shape); +RcppExport SEXP _tsdistributions_c_qged(SEXP pSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type p(pSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_qged(p, mu, sigma, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_rsged +NumericVector c_rsged(int n, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_rsged(SEXP nSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_rsged(n, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_dsged +NumericVector c_dsged(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr); +RcppExport SEXP _tsdistributions_c_dsged(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dsged(x, mu, sigma, skew, shape, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_psged +NumericVector c_psged(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_psged(SEXP qSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type q(qSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_psged(q, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_qsged +NumericVector c_qsged(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape); +RcppExport SEXP _tsdistributions_c_qsged(SEXP pSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type p(pSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + rcpp_result_gen = Rcpp::wrap(c_qsged(p, mu, sigma, skew, shape)); + return rcpp_result_gen; +END_RCPP +} +// c_dhyp +NumericVector c_dhyp(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr); +RcppExport SEXP _tsdistributions_c_dhyp(SEXP xSEXP, SEXP muSEXP, SEXP sigmaSEXP, SEXP skewSEXP, SEXP shapeSEXP, SEXP logrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericVector >::type x(xSEXP); + Rcpp::traits::input_parameter< NumericVector >::type mu(muSEXP); + Rcpp::traits::input_parameter< NumericVector >::type sigma(sigmaSEXP); + Rcpp::traits::input_parameter< NumericVector >::type skew(skewSEXP); + Rcpp::traits::input_parameter< NumericVector >::type shape(shapeSEXP); + Rcpp::traits::input_parameter< int >::type logr(logrSEXP); + rcpp_result_gen = Rcpp::wrap(c_dhyp(x, mu, sigma, skew, shape, logr)); + return rcpp_result_gen; +END_RCPP +} +// c_rghyp +NumericVector c_rghyp(int n, double mu, double delta, double alpha, double beta, double lambda); +RcppExport SEXP _tsdistributions_c_rghyp(SEXP nSEXP, SEXP muSEXP, SEXP deltaSEXP, SEXP alphaSEXP, SEXP betaSEXP, SEXP lambdaSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type n(nSEXP); + Rcpp::traits::input_parameter< double >::type mu(muSEXP); + Rcpp::traits::input_parameter< double >::type delta(deltaSEXP); + Rcpp::traits::input_parameter< double >::type alpha(alphaSEXP); + Rcpp::traits::input_parameter< double >::type beta(betaSEXP); + Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP); + rcpp_result_gen = Rcpp::wrap(c_rghyp(n, mu, delta, alpha, beta, lambda)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_tsdistributions_c_dghst", (DL_FUNC) &_tsdistributions_c_dghst, 6}, + {"_tsdistributions_c_rghst", (DL_FUNC) &_tsdistributions_c_rghst, 5}, + {"_tsdistributions_c_dghyp", (DL_FUNC) &_tsdistributions_c_dghyp, 7}, + {"_tsdistributions_c_dgh", (DL_FUNC) &_tsdistributions_c_dgh, 7}, + {"_tsdistributions_c_dnig", (DL_FUNC) &_tsdistributions_c_dnig, 6}, + {"_tsdistributions_c_rstd", (DL_FUNC) &_tsdistributions_c_rstd, 4}, + {"_tsdistributions_c_dstd", (DL_FUNC) &_tsdistributions_c_dstd, 5}, + {"_tsdistributions_c_pstd", (DL_FUNC) &_tsdistributions_c_pstd, 4}, + {"_tsdistributions_c_qstd", (DL_FUNC) &_tsdistributions_c_qstd, 4}, + {"_tsdistributions_c_rsstd", (DL_FUNC) &_tsdistributions_c_rsstd, 5}, + {"_tsdistributions_c_dsstd", (DL_FUNC) &_tsdistributions_c_dsstd, 6}, + {"_tsdistributions_c_psstd", (DL_FUNC) &_tsdistributions_c_psstd, 5}, + {"_tsdistributions_c_qsstd", (DL_FUNC) &_tsdistributions_c_qsstd, 5}, + {"_tsdistributions_c_djsu", (DL_FUNC) &_tsdistributions_c_djsu, 6}, + {"_tsdistributions_c_qjsu", (DL_FUNC) &_tsdistributions_c_qjsu, 5}, + {"_tsdistributions_c_pjsu", (DL_FUNC) &_tsdistributions_c_pjsu, 5}, + {"_tsdistributions_c_rjsu", (DL_FUNC) &_tsdistributions_c_rjsu, 5}, + {"_tsdistributions_c_rsnorm", (DL_FUNC) &_tsdistributions_c_rsnorm, 4}, + {"_tsdistributions_c_dsnorm", (DL_FUNC) &_tsdistributions_c_dsnorm, 5}, + {"_tsdistributions_c_psnorm", (DL_FUNC) &_tsdistributions_c_psnorm, 4}, + {"_tsdistributions_c_qsnorm", (DL_FUNC) &_tsdistributions_c_qsnorm, 4}, + {"_tsdistributions_c_rged", (DL_FUNC) &_tsdistributions_c_rged, 4}, + {"_tsdistributions_c_dged", (DL_FUNC) &_tsdistributions_c_dged, 5}, + {"_tsdistributions_c_pged", (DL_FUNC) &_tsdistributions_c_pged, 4}, + {"_tsdistributions_c_qged", (DL_FUNC) &_tsdistributions_c_qged, 4}, + {"_tsdistributions_c_rsged", (DL_FUNC) &_tsdistributions_c_rsged, 5}, + {"_tsdistributions_c_dsged", (DL_FUNC) &_tsdistributions_c_dsged, 6}, + {"_tsdistributions_c_psged", (DL_FUNC) &_tsdistributions_c_psged, 5}, + {"_tsdistributions_c_qsged", (DL_FUNC) &_tsdistributions_c_qsged, 5}, + {"_tsdistributions_c_dhyp", (DL_FUNC) &_tsdistributions_c_dhyp, 6}, + {"_tsdistributions_c_rghyp", (DL_FUNC) &_tsdistributions_c_rghyp, 6}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_tsdistributions(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/TMB/compile.R b/src/TMB/compile.R old mode 100644 new mode 100755 diff --git a/src/TMB/distfun.h b/src/TMB/distfun.h old mode 100644 new mode 100755 diff --git a/src/TMB/distmodel.hpp b/src/TMB/distmodel.hpp old mode 100644 new mode 100755 diff --git a/src/TMB/tsdistributions_TMBExports.cpp b/src/TMB/tsdistributions_TMBExports.cpp old mode 100644 new mode 100755 index a07d44f..5dce1fd --- a/src/TMB/tsdistributions_TMBExports.cpp +++ b/src/TMB/tsdistributions_TMBExports.cpp @@ -1,5 +1,3 @@ -// Generated by TMBtools: do not edit by hand - #define TMB_LIB_INIT R_init_tsdistributions_TMBExports #include #include "distfun.h" diff --git a/src/distributions.c b/src/distributions.c deleted file mode 100644 index cdf648f..0000000 --- a/src/distributions.c +++ /dev/null @@ -1,848 +0,0 @@ -# include "distributions.h" -# include -# include -# include -# include - -/* - * Key Functions - */ -double dnorm_std(const double x) -{ - double pdf; - pdf = exp (-0.5 * x * x) / sqrt (2.0 * PI); - if(pdf == 0.0) pdf = 0.0 + 2.22507e-24; - return pdf; -} - -double signum(const double x) -{ - double res = -(x < 0.0) + (x > 0.0); - return res; -} - -double dhuge(void) -{ - return HUGE_VAL; -} - -double heaviside(const double x, const double a){ - double res = (signum(x - a) + 1.0) / 2.0; - return res; -} - -double depsilon(void) -{ - double r; - r = 1.0; - while( 1.0 < (double)(1.0 + r)) - { - r = r / 2.0; - } - double res = 2.0 * r; - return res; -} - -double kappagh(const double x, const double lambda) -{ - double kappa = 0.0; - if(lambda == -0.5) { - kappa = 1.0 / x; - } else { - kappa = (bessel_k(x, lambda + 1.0, 2.0) / bessel_k(x, lambda, 2.0)) / x; - } - return kappa; -} - -double deltakappagh(const double x, const double lambda) -{ - double deltakappa = 0.0; - deltakappa = kappagh(x, lambda + 1.0) - kappagh(x, lambda); - return deltakappa; -} - -double* paramgh(const double rho, const double zeta, const double lambda) -{ - double *param = malloc(4 * sizeof(double)); - double rho2 = 1.0 - pow(rho, 2.0); - double alpha = pow(zeta, 2.0) * kappagh(zeta, lambda) / rho2; - alpha = alpha * (1.0 + pow(rho, 2.0) * pow(zeta, 2.0) * deltakappagh(zeta, lambda) / rho2); - alpha = sqrt(alpha); - double beta = alpha * rho; - double delta = zeta / (alpha * sqrt(rho2)); - double mu = -beta * pow(delta, 2.0) * kappagh(zeta, lambda); - param[0] = (double) alpha; - param[1] = (double) beta; - param[2] = (double) delta; - param[3] = (double) mu; - return param; -} - -double* paramghst(const double betabar, const double nu) -{ - double *param = malloc(4 * sizeof(double)); - double delta = sqrt(1.0 / (((2.0 * betabar * betabar) / ((nu - 2.0) * (nu - 2.0) * (nu - 4.0))) + (1.0/(nu - 2.0)))); - double beta = betabar / delta; - double mu = -1.0 * ((beta * (delta * delta)) / (nu - 2.0)); - param[0] = (double) nu; - param[1] = (double) beta; - param[2] = (double) delta; - param[3] = (double) mu; - return param; -} - -/* - * GH Skew Student Distribution - */ -double dghst_std(const double x, const double betabar, const double nu) -{ - double *param; - param = paramghst(betabar, nu); - double beta = param[1]; - double delta = param[2]; - double mu = param[3]; - double betasqr = beta * beta; - double deltasqr = delta * delta; - double res = x - mu; - double pdf = ((1.0 - nu) / 2.0) * log(2) + nu * log(delta) + ((nu + 1.0) / 2.0) * log(fabs(beta)) + - log(bessel_k(sqrt(betasqr * (deltasqr + res * res)), (nu + 1.0)/2.0, 2.0)) - sqrt(betasqr * - (deltasqr + res * res)) + beta * res - lgammafn(nu/2.0) - log(PI)/2.0 - ((nu + 1.0) / 2.0) * - log(deltasqr + res * res) / 2.0; - free(param); - pdf = exp(pdf); - return pdf; -} - -void c_dghst(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dghst_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - -double rsghst_std(const double betabar, const double nu) -{ - // Note: for values of nu<5 it is likely that sd(r) will have a very large variance - // Existence of moment conditions (vs lower parameter bounds) are defined in the paper. - double *param; - param = paramghst(betabar, nu); - double beta = param[1]; - double delta = param[2]; - double mu = param[3]; - double y = 1.0 / rgamma(nu / 2.0, 2.0 / (delta * delta)); - double sigma = sqrt(y); - double z = rnorm(0,1); - double ans = mu + beta * sigma * sigma + sigma * z; - free(param); - return ans; -} - -void c_rghst(int *n, double *mu, double *sigma, double *skew, double *shape, double *ans) -{ - GetRNGstate(); - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + rsghst_std(skew[i], shape[i]) * sigma[i]; - } - PutRNGstate(); -} - -/* - * GH Distribution - */ -double dgh(const double x, const double alpha, const double beta, const double delta, const double mu, const double lambda) -{ - double pdf = 0.0; - if(alpha <= 0.0) { - return pdf = 0.0; - } - if(delta <= 0.0) { - return pdf = 0.0; - } - if(fabs(beta) >= alpha) { - return pdf = 0.0; - } - double arg = delta * sqrt(pow(alpha, 2.0) - pow(beta, 2.0)); - double a = (lambda / 2.0) * log(pow(alpha, 2.0) - pow(beta, 2.0)) - (log(sqrt(2.0 * PI)) + (lambda - 0.5) * log(alpha) + - lambda * log(delta) + log(bessel_k(arg, lambda, 2.0)) - arg); - double f = ((lambda - 0.5) / 2.0) * log(pow(delta, 2.0) + pow((x - mu), 2.0)); - arg = alpha * sqrt(pow(delta, 2.0) + pow((x - mu), 2.0)); - double k = log(bessel_k(arg, lambda - 0.5, 2.0)) - arg; - double e = beta * (x - mu); - pdf = exp(a + f + k + e); - return pdf; -} - -void c_dghyp(double *x, double *alpha, double *beta, double *delta, double *mu, double *lambda, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) - { - ans[i] = dgh(x[i], alpha[i], beta[i], delta[i], mu[i], lambda[i]); - if(*logr == 1) ans[i] = log(ans[i]); - } -} - - -double dgh_std(const double x, const double rho, const double zeta, const double lambda) -{ - double pdf; - double *param; - param = paramgh(rho, zeta, lambda); - double alpha = param[0]; - double beta = param[1]; - double delta = param[2]; - double mu = param[3]; - pdf = dgh(x, alpha, beta, delta, mu, lambda); - free(param); - return pdf; -} - -void c_dgh(double *x, double *mu, double *sigma, double *skew, double *shape, double *lambda, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dgh_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i], lambda[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - - -/* - * NIG Distribution - */ - -double dnig_std(const double x, const double rho, const double zeta) -{ - double pdf = 0.0; - double lambda = -0.5; - double *param; - param = paramgh(rho, zeta, lambda); - double alpha = param[0]; - double beta = param[1]; - double delta = param[2]; - double mu = param[3]; - double deltasq = delta * delta; - double res = x - mu; - pdf = -log(PI) + log(alpha) + log(delta) + log(bessel_k(alpha * sqrt(deltasq + res * res), 1, 1)) + - delta * sqrt(alpha * alpha - beta * beta) + beta * res -0.5 * log(deltasq + res * res); - pdf = exp(pdf); - free(param); - return pdf; -} - -void c_dnig(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dnig_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - - -/* - * Student Distribution - */ - -double rstd_std(const double nu) -{ - double ans = 0.0; - if(nu > 2.0) { - double s = sqrt(nu / (nu - 2.0)); - ans = rt(nu) * 1.0 / s; - } - return ans; -} - -void c_rstd(int *n, double *mu, double *sigma, double *shape, double *ans) -{ - GetRNGstate(); - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + rstd_std(shape[i]) * sigma[i]; - } - PutRNGstate(); -} - -double xdt(const double x, const double nu) -{ - double a, b, pdf; - a = gammafn((nu + 1.0)/2.0) / sqrt(PI * nu); - b = gammafn(nu / 2.0) * pow((1.0 + (x * x) / nu), ((nu + 1.0) / 2.0)); - pdf = a / b; - return pdf; -} - -double dstd_std(const double x, const double nu) -{ - double pdf, s; - if(nu <= 2) { - pdf = 999; - } else { - s = sqrt(nu / (nu - 2.0)); - pdf = s * xdt(x * s, nu); - } - return pdf; -} - -void c_dstd(double *x, double *mu, double *sigma, double *shape, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dstd_std((x[i] - mu[i]) / sigma[i], shape[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - -double pstd_std(const double q, const double mu, const double sigma, const double nu) -{ - double s = sqrt(nu / (nu - 2.0)); - double z = (q - mu) / sigma; - double p = pt(z * s, nu, 1, 0); - return p; -} - -void c_pstd(double *q, double *mu, double *sigma, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = pstd_std(q[i], mu[i], sigma[i], shape[i]); - } -} - -double qstd_std(const double p, const double mu, const double sigma, const double nu) -{ - double s = sqrt(nu / (nu - 2.0)); - double q = qt(p, nu, 1, 0) * sigma / s + mu; - return q; -} - -void c_qstd(double *p, double *mu, double *sigma, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = qstd_std(p[i], mu[i], sigma[i], shape[i]); - } -} - -/* - * Skew Student Distribution (Fernandez & Steel) - */ -double rsstd_std(const double xi, const double nu) -{ - double weight, z, rr, m1, mu, sigma, xx, ans; - ans = 0.0; - weight = xi / (xi + 1.0 / xi); - z = runif(-1.0 * weight, 1.0 - weight); - xx = (z < 0) ? 1.0 / xi : xi; - rr = -1.0 * fabs(rstd_std(nu)) / xx * sign(z); - m1 = 2.0 * sqrt(nu - 2.0) / (nu - 1.0) / beta(0.5, 0.5 * nu); - mu = m1 * (xi - 1.0 / xi); - sigma = sqrt((1.0 - (m1 * m1)) * ((xi * xi) + 1.0 / (xi * xi)) + 2 * (m1 * m1) - 1.0); - ans = (rr - mu) / sigma; - return ans; -} - -void c_rsstd(int *n, double *mu, double *sigma, double *skew, double *shape, double *ans) -{ - GetRNGstate(); - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + rsstd_std(skew[i], shape[i]) * sigma[i]; - } - PutRNGstate(); -} - -double dsstd_std(const double x, const double xi, const double nu) -{ - double mu, m1,beta, sigma, z, g,pdf,a,b, xxi; - xxi = xi; - a = 0.5; - b = nu / 2.0; - beta = (gammafn(a) / gammafn(a + b)) * gammafn(b); - m1 = 2.0 * sqrt(nu - 2.0) / (nu - 1.0) / beta; - mu = m1 * (xi - 1.0 / xi); - sigma = sqrt((1.0 - pow(m1, 2.0)) * (pow(xi, 2.0) + 1.0 / (pow(xi, 2.0))) + 2.0 * pow(m1, 2.0) - 1.0); - z = x * sigma + mu; - if(z == 0) { - xxi = 1.0; - } - if(z < 0) { - xxi = 1.0 / xi; - } - g = 2.0 / (xi + 1.0 / xi); - pdf = g * dstd_std(z / xxi, nu) * sigma; - return pdf; -} - -void c_dsstd(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dsstd_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - -double psstd_std(const double q, const double mu, const double sigma, const double xi, const double nu) -{ - double qx = (q - mu) / sigma; - double m1 = 2.0 * sqrt(nu - 2.0) / (nu - 1.0) / beta(0.5, nu / 2.0); - double mux = m1 * (xi - 1.0 / xi); - double sig = sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); - double z = qx * sig + mux; - double Xi = (z < 0) ? 1.0 / xi : xi; - double g = 2.0 / (xi + 1.0 / xi); - double p = heaviside(z, 0) - signum(z) * g * Xi * pstd_std(-fabs(z) / Xi, 0.0, 1.0, nu); - return p; -} - -void c_psstd(double *q, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = psstd_std(q[i], mu[i], sigma[i], skew[i], shape[i]); - } -} - -double qsstd_std(const double p, const double xi, const double nu) -{ - double m1 = 2.0 * sqrt(nu - 2.0) / (nu - 1.0) / beta(0.5, nu / 2.0); - double mu = m1 * (xi - 1.0 / xi); - double sigma = sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); - double g = 2.0 / (xi + 1.0 / xi); - double z = p - (1.0/(1.0 + xi * xi)); - double Xi = pow(xi, signum(z)); - double tmp = (heaviside(z, 0) - signum(z) * p) / (g * Xi); - double q = (-signum(z) * qstd_std(tmp, 0.0, 1.0, nu) * Xi - mu) / sigma; - return q; -} - -void c_qsstd(double *p, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = qsstd_std(p[i], skew[i], shape[i]) * sigma[i] + mu[i]; - } -} - -/* - * Johnson's SU Distribution - */ -//nu = skew, tau = shape (!) -double djsu_std(const double x, const double nu, const double tau) -{ - double w, z, r, omega, c, pdf = 0.0; - double rtau = 1.0 / tau; - if(rtau < 0.0000001) { - w = 1.0; - } else { - w = exp(rtau * rtau); - } - omega= -nu * rtau; - c = sqrt(1.0 / (0.5 * (w - 1.0) * (w * cosh(2.0 * omega) + 1.0))); - z = (x - (c * sqrt(w) * sinh(omega))) / c; - r = -nu + asinh(z) / rtau; - pdf = -log(c) - log(rtau) - 0.5 * log(z * z + 1.0) - 0.5 * log(2.0 * PI) - 0.5 * r * r; - pdf = exp(pdf); - return pdf; -} - -void c_djsu(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = djsu_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - -double qjsu_std(const double p, const double nu, const double tau) -{ - double rtau, rr, z, w, omega, cc, ans; - ans = 0.0; - rtau = 1.0 / tau; - rr = qnorm(p, 0.0, 1.0, 1, 0); - z = sinh(rtau * (rr + nu)); - w = (rtau < 0.0000001) ? 1 : exp(rtau * rtau); - omega = -1.0 * nu * rtau; - cc = sqrt(1.0 / (0.5 * (w - 1.0) * (w * cosh(2.0 * omega) + 1.0))); - ans = (cc * sqrt(w) * sinh(omega)) + cc * z; - return ans; -} - -void c_qjsu(double *p, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + qjsu_std(p[i], skew[i], shape[i]) * sigma[i]; - } -} - -double rjsu_std(const double nu, const double tau) -{ - double x, ans; - ans = 0.0; - x = runif(0, 1); - ans = qjsu_std(x, nu, tau); - return ans; -} - -void c_rjsu(int *n, double *mu, double *sigma, double *skew, double *shape, double *ans) -{ - GetRNGstate(); - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + rjsu_std(skew[i], shape[i]) * sigma[i]; - } - PutRNGstate(); -} - -double pjsu_std(const double q, const double mu, const double sigma, const double nu, const double tau) -{ - double rtau = 1.0 / tau; - double w = (rtau < 0.0000001) ? 1.0 : exp(rtau * rtau); - double omega = -1.0 * nu * rtau; - double c = 1.0 / sqrt(0.5 * (w - 1.0) * (w * cosh(2.0 * omega) + 1.0)); - double z = (q - (mu + c * sigma * sqrt(w) * sinh(omega))) / (c * sigma); - double r = -1.0 * nu + asinh(z) / rtau; - double p = pnorm(r, 0, 1, 1, 0); - return p; -} - -void c_pjsu(double *q, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = pjsu_std(q[i], mu[i], sigma[i], skew[i], shape[i]); - } -} - -/* - * Skew Normal Distribution - */ -double rsnorm_std(const double xi) -{ - double weight, z, rr, m1, mu, sigma, xx, ans; - weight = xi / (xi + 1.0/xi); - z = runif(-weight, 1.0 - weight); - xx = (z < 0) ? 1.0 / xi : xi; - rr = -1.0 * fabs(rnorm(0, 1)) / xx * sign(z); - m1 = 2.0 / sqrt(2.0 * PI); - mu = m1 * (xi - 1.0 / xi); - sigma = sqrt((1.0 - (m1 * m1)) * ((xi * xi) + 1.0 / (xi * xi) ) + 2 * (m1 * m1) - 1.0); - ans = (rr - mu ) / sigma; - return ans; -} - -void c_rsnorm(int *n, double *mu, double *sigma, double *skew, double *ans) -{ - GetRNGstate(); - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + rsnorm_std(skew[i]) * sigma[i]; - } - PutRNGstate(); -} - -double dsnorm_std(const double x, const double xi) -{ - double pdf; - double mu, sigma,z, xxi, g; - double m1 = 2.0 / sqrt(2.0 * PI); - double m12 = m1 * m1; - double xi2 = xi * xi; - mu = m1 * (xi - 1.0 / xi); - sigma = sqrt((1.0 - m12) * (xi2 + 1.0 / xi2) + 2.0 * m12 - 1.0); - z = x * sigma + mu; - xxi = xi; - if(z == 0.0) { - xxi = 1.0; - } - if(z < 0.0) { - xxi = 1.0 / xi; - } - g = 2.0 / (xi + 1.0 / xi); - pdf = g * dnorm_std(z / xxi) * sigma; - return pdf; -} - -void c_dsnorm(double *x, double *mu, double *sigma, double *skew, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dsnorm_std((x[i] - mu[i]) / sigma[i], skew[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - -double psnorm_std(const double q, const double mu, const double sigma, const double xi) -{ - double qx = (q - mu) / sigma; - double m1 = 2.0 / sqrt(2.0 * PI); - double mux = m1 * (xi - 1.0 / xi); - double sig = sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); - double z = qx * sig + mux; - double Xi = (z < 0) ? 1.0 / xi : xi; - double g = 2.0 / (xi + 1.0 / xi); - double p = heaviside(z, 0) - signum(z) * g * Xi * pnorm(-fabs(z)/Xi, 0, 1, 1, 0); - return p; -} - -void c_psnorm(double *q, double *mu, double *sigma, double *skew, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = psnorm_std(q[i], mu[i], sigma[i], skew[i]); - } -} - -double qsnorm_std(const double p, const double xi) -{ - double m1 = 2.0 / sqrt(2.0 * PI); - double mu = m1 * (xi - 1.0 / xi); - double sigma = sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); - double g = 2.0 / (xi + 1.0 / xi); - double z = p - (1.0 / (1.0 + xi * xi)); - double Xi = pow(xi, signum(z)); - double tmp = (heaviside(z, 0) - signum(z) * p) / (g * Xi); - double q = (-1.0 * signum(z) * qnorm(tmp, 0, Xi, 1, 0) - mu) / sigma; - return q; -} - -void c_qsnorm(double *p, double *mu, double *sigma, double *skew, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + qsnorm_std(p[i], skew[i]) * sigma[i]; - } -} - -/* - * Generalized Error Distribution - */ -double rged_std(const double nu) -{ - double lambda, rr, ans; - ans = 0.0; - lambda = sqrt(pow(0.5, 2.0 / nu) * gammafn(1.0 / nu) / gammafn(3.0 / nu)); - rr = rgamma(1.0 / nu, 1.0); - ans = lambda * pow(2.0 * rr, 1.0 / nu) * sign(runif(0, 1) - 0.5); - return ans; -} - -void c_rged(int *n, double *mu, double *sigma, double *shape, double *ans) -{ - GetRNGstate(); - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + rged_std(shape[i]) * sigma[i]; - } - PutRNGstate(); -} - -double dged_std(const double x, const double nu) -{ - double lambda, g, pdf; - lambda = sqrt(pow(1.0 / 2.0, 2.0 / nu) * gammafn(1.0 / nu) / gammafn(3.0 / nu)); - g = nu / (lambda * (pow(2.0, 1.0 + (1.0 / nu))) * gammafn(1.0 / nu)); - pdf = g * exp(-0.5 * pow(fabs(x / lambda), nu)); - return pdf; -} - -void c_dged(double *x, double *mu, double *sigma, double *shape, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dged_std((x[i] - mu[i]) / sigma[i], shape[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - -double pged_std(const double q, const double mu, const double sigma, const double nu) -{ - double qx = (q - mu) / sigma; - double lambda = sqrt(1.0 / pow(2.0, (2.0 / nu)) * gammafn(1.0 / nu) / gammafn(3.0 / nu)); - double g = nu / (lambda * (pow(2.0, (1.0 + 1.0 / nu))) * gammafn(1.0 / nu)); - double h = pow(2.0, (1.0 / nu)) * lambda * g * gammafn(1.0 / nu) / nu; - double s = 0.5 * pow(fabs(qx) / lambda , nu); - double p = 0.5 + signum(qx) * h * pgamma(s, 1.0 / nu, 1, 1, 0); - return p; -} - -void c_pged(double *q, double *mu, double *sigma, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = pged_std(q[i], mu[i], sigma[i], shape[i]); - } -} - -double qged_std(const double p, const double shape) -{ - double y = 2.0 * p - 1.0; - double lambda = sqrt(1.0 / pow(2.0, (2.0 / shape)) * gammafn(1.0 / shape) / gammafn(3.0 / shape)); - double q = lambda * pow(2.0 * qgamma(fabs(y), 1.0 / shape, 1, 1, 0), 1.0 / shape); - q = q * signum(y); - return q; -} - -void c_qged(double *p, double *mu, double *sigma, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = qged_std(p[i], shape[i]) * sigma[i] + mu[i]; - } -} - -/* - * Skew Generalized Error Distribution (Fernandez & Steel) - */ -double rsged_std(const double xi, const double nu) -{ - double weight, lambda, z, rr, m1, mu, sigma, xx, ans; - weight = xi / (xi + 1.0/xi); - z = runif(-1.0 * weight, 1.0 - weight); - xx = (z < 0) ? 1.0 / xi : xi; - rr = -1.0 * fabs(rged_std(nu)) / xx * sign(z); - lambda = sqrt (pow(0.5, 2.0/nu) * gammafn(1.0/nu) / gammafn(3.0/nu)); - m1 = pow(2.0, 1.0 / nu) * lambda * gammafn(2.0 / nu) / gammafn(1.0 / nu); - mu = m1 * (xi - 1.0 / xi); - sigma = sqrt((1.0 - (m1 * m1)) * ((xi * xi) + 1.0 / (xi* xi)) + 2.0 * (m1 * m1) - 1.0); - ans = (rr - mu) / sigma; - return ans; -} - -void c_rsged(int *n, double *mu, double *sigma, double *skew, double *shape, double *ans) -{ - GetRNGstate(); - int i; - for(i = 0; i<*n; i++) { - ans[i] = mu[i] + rsged_std(skew[i], shape[i]) * sigma[i]; - } - PutRNGstate(); -} - -double dsged_std(const double x, const double xi, const double nu) -{ - double lambda, m1, mu, sigma, z, g, pdf, xxi; - xxi = xi; - lambda = sqrt(pow(1.0 / 2.0, 2.0 / nu) * gammafn(1.0 / nu) / gammafn(3.0 / nu)); - m1 = pow(2.0, 1.0 / nu) * lambda * gammafn(2.0 / nu) / gammafn(1.0 / nu); - mu = m1 * (xi - 1.0 / xi); - sigma = (1.0 - pow(m1, 2.0)) * (pow(xi, 2.0) + 1.0 / (pow(xi, 2.0))) + 2.0 * (pow(m1, 2.0)) - 1.0; - sigma = sqrt(sigma); - z = x * sigma + mu; - if(z == 0) { - xxi = 1.0; - } - if(z < 0) { - xxi = 1 / xi; - } - g = 2.0 / (xi + 1.0 / xi); - pdf = g * dged_std(z / xxi, nu) * sigma; - return pdf; -} - -void c_dsged(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dsged_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} - -double psged_std(const double q, const double mu, const double sigma, const double xi, const double nu) -{ - double qx = (q - mu) / sigma; - double lambda = sqrt (1.0/pow(2.0, 2.0 / nu) * gammafn(1.0 / nu) / gammafn(3.0 / nu)); - double m1 = pow(2.0, 1.0 / nu) * lambda * gammafn(2.0 / nu) / gammafn(1.0 / nu); - double mux = m1 * (xi - 1.0 / xi); - double sig = sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); - double z = qx * sig + mux; - double Xi = (z < 0) ? 1.0 / xi : xi; - double g = 2.0 / (xi + 1.0 / xi); - double p = heaviside(z, 0) - signum(z) * g * Xi * pged_std(-fabs(z) / Xi, 0.0, 1.0, nu); - return p; -} - -void c_psged(double *q, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = psged_std(q[i], mu[i], sigma[i], skew[i], shape[i]); - } -} - -double qsged_std(const double p, const double xi, const double nu) -{ - double lambda = sqrt (1.0 / pow(2.0, 2.0 / nu) * gammafn(1.0 / nu) / gammafn(3.0 / nu)); - double m1 = pow(2.0, 1.0 / nu) * lambda * gammafn(2.0 / nu) / gammafn(1.0 / nu); - double mu = m1 * (xi - 1.0 / xi); - double sigma = sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); - double g = 2.0 / (xi + 1.0 / xi); - double z = p - (1.0 / (1.0 + xi * xi)); - double Xi = pow(xi, signum(z)); - double q = (heaviside(z, 0) - signum(z) * p) / (g * Xi); - q = (-signum(z) * qged_std(q, nu) * Xi - mu) / sigma; - return q; -} - -void c_qsged(double *p, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = qsged_std(p[i], skew[i], shape[i]) * sigma[i] + mu[i]; - } -} - -/* - * ----------------------------------------- - * Hypebolic Distribution - * ----------------------------------------- - */ -double dhyp(const double x, const double alpha, const double beta, const double delta, const double mu) -{ - double pdf = 0.0; - if(alpha <= 0.0) { - return pdf; - } - if(delta <= 0) { - return pdf; - } - if(fabs(beta) >= alpha) { - return pdf; - } - double g = alpha * alpha - beta * beta; - double e = x - mu; - pdf = 0.5 * log(g) - log(2.0 * alpha * delta * bessel_k(delta * sqrt(g), 1, 2)) - alpha * sqrt(delta * delta + e * e) + beta * e; - pdf = exp(pdf); - return pdf; -} - -double dhyp_std(const double x, const double rho, const double zeta) -{ - double pdf; - double *param; - param = paramgh(rho, zeta, 1.0); - double alpha = param[0]; - double beta = param[1]; - double delta = param[2]; - double mu = param[3]; - pdf = dhyp(x, alpha, beta, delta, mu); - free(param); - return pdf; -} - -void c_dhyp(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr) -{ - int i; - for(i = 0; i<*n; i++) { - ans[i] = dhyp_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; - if(*logr == 1) ans[i] = log(ans[i]); - } -} diff --git a/src/distributions.cpp b/src/distributions.cpp new file mode 100644 index 0000000..c49ccb2 --- /dev/null +++ b/src/distributions.cpp @@ -0,0 +1,1009 @@ +# include "distributions.h" +using namespace Rcpp; + +// Standard normal density function +double dnorm_std(const double x) { + double pdf = std::exp(-0.5 * x * x) / std::sqrt(2.0 * M_PI); + if (pdf == 0.0) pdf = std::numeric_limits::min(); // Smallest positive double + return pdf; +} + +// Sign function +double signum(const double x) { + return -(x < 0.0) + (x > 0.0); +} + +// Return a huge value (infinity) +double dhuge(void) { + return R_PosInf; +} + +// Heaviside step function +double heaviside(const double x, const double a) { + return (signum(x - a) + 1.0) / 2.0; +} + +// Machine epsilon +double depsilon(void) { + double r = 1.0; + while (1.0 < (1.0 + r)) { + r /= 2.0; + } + return 2.0 * r; +} + +// Kappa function for Generalized Hyperbolic Distribution +double kappagh(const double x, const double lambda) { + double kappa = 0.0; + if (lambda == -0.5) { + kappa = 1.0 / x; + } else { + kappa = (R::bessel_k(x, lambda + 1.0, 2.0) / R::bessel_k(x, lambda, 2.0)) / x; + } + return kappa; +} + +// Delta Kappa function +double deltakappagh(const double x, const double lambda) { + return kappagh(x, lambda + 1.0) - kappagh(x, lambda); +} + +// Parametrization for Generalized Hyperbolic Distribution +NumericVector paramgh(const double rho, const double zeta, const double lambda) { + double rho2 = 1.0 - std::pow(rho, 2.0); + double alpha = std::pow(zeta, 2.0) * kappagh(zeta, lambda) / rho2; + alpha *= (1.0 + std::pow(rho, 2.0) * std::pow(zeta, 2.0) * deltakappagh(zeta, lambda) / rho2); + alpha = std::sqrt(alpha); + double beta = alpha * rho; + double delta = zeta / (alpha * std::sqrt(rho2)); + double mu = -beta * std::pow(delta, 2.0) * kappagh(zeta, lambda); + + return NumericVector::create(alpha, beta, delta, mu); +} + +NumericVector paramghst(const double betabar, const double nu) { + double delta = std::sqrt(1.0 / (((2.0 * betabar * betabar) / ((nu - 2.0) * (nu - 2.0) * (nu - 4.0))) + (1.0 / (nu - 2.0)))); + double beta = betabar / delta; + double mu = -beta * (delta * delta) / (nu - 2.0); + + return NumericVector::create(nu, beta, delta, mu); +} + + +double dghst_std(const double x, const double betabar, const double nu) { + NumericVector param = paramghst(betabar, nu); + double beta = param[1]; + double delta = param[2]; + double mu = param[3]; + + double betasqr = beta * beta; + double deltasqr = delta * delta; + double res = x - mu; + + double pdf = ((1.0 - nu) / 2.0) * std::log(2) + nu * std::log(delta) + + ((nu + 1.0) / 2.0) * std::log(std::abs(beta)) + + std::log(R::bessel_k(std::sqrt(betasqr * (deltasqr + res * res)), (nu + 1.0) / 2.0, 2.0)) - + std::sqrt(betasqr * (deltasqr + res * res)) + beta * res - + R::lgammafn(nu / 2.0) - std::log(M_PI) / 2.0 - ((nu + 1.0) / 2.0) * + std::log(deltasqr + res * res) / 2.0; + + pdf = std::exp(pdf); + + return pdf; +} + +// [[Rcpp::export]] +NumericVector c_dghst(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr) { + int n = x.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = dghst_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; + if(logr == 1) ans[i] = std::log(ans[i]); + } + + return ans; +} + +double rsghst_std(const double betabar, const double nu) { + NumericVector param = paramghst(betabar, nu); + double beta = param[1]; + double delta = param[2]; + double mu = param[3]; + + // Sample y from the inverse gamma distribution + double y = 1.0 / R::rgamma(nu / 2.0, 2.0 / (delta * delta)); + double sigma = std::sqrt(y); + + // Sample z from standard normal distribution + double z = R::rnorm(0, 1); + + // Calculate the final result + double ans = mu + beta * sigma * sigma + sigma * z; + + return ans; +} + +// [[Rcpp::export]] +NumericVector c_rghst(int n, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = mu[i] + rsghst_std(skew[i], shape[i]) * sigma[i]; + } + + return ans; +} + +double dgh(const double x, const double alpha, const double beta, const double delta, const double mu, const double lambda) { + if (alpha <= 0.0 || delta <= 0.0 || std::abs(beta) >= alpha) { + return 0.0; + } + + double arg = delta * std::sqrt(std::pow(alpha, 2.0) - std::pow(beta, 2.0)); + double a = (lambda / 2.0) * std::log(std::pow(alpha, 2.0) - std::pow(beta, 2.0)) - + (std::log(std::sqrt(2.0 * M_PI)) + (lambda - 0.5) * std::log(alpha) + + lambda * std::log(delta) + std::log(R::bessel_k(arg, lambda, 2.0)) - arg); + + double f = ((lambda - 0.5) / 2.0) * std::log(std::pow(delta, 2.0) + std::pow((x - mu), 2.0)); + arg = alpha * std::sqrt(std::pow(delta, 2.0) + std::pow((x - mu), 2.0)); + double k = std::log(R::bessel_k(arg, lambda - 0.5, 2.0)) - arg; + double e = beta * (x - mu); + + return std::exp(a + f + k + e); +} + +// [[Rcpp::export]] +NumericVector c_dghyp(NumericVector x, NumericVector alpha, NumericVector beta, NumericVector delta, NumericVector mu, NumericVector lambda, int logr) { + int n = x.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = dgh(x[i], alpha[i], beta[i], delta[i], mu[i], lambda[i]); + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double dgh_std(const double x, const double rho, const double zeta, const double lambda) { + NumericVector param = paramgh(rho, zeta, lambda); + double alpha = param[0]; + double beta = param[1]; + double delta = param[2]; + double mu = param[3]; + + return dgh(x, alpha, beta, delta, mu, lambda); +} + +// [[Rcpp::export]] +NumericVector c_dgh(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, NumericVector lambda, int logr) { + int n = x.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = dgh_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i], lambda[i]) / sigma[i]; + if(logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double dnig_std(const double x, const double rho, const double zeta) { + double lambda = -0.5; + NumericVector param = paramgh(rho, zeta, lambda); + double alpha = param[0]; + double beta = param[1]; + double delta = param[2]; + double mu = param[3]; + + double deltasq = delta * delta; + double res = x - mu; + + // Calculate the PDF using the Normal Inverse Gaussian formula + double pdf = -std::log(M_PI) + std::log(alpha) + std::log(delta) + + std::log(R::bessel_k(alpha * std::sqrt(deltasq + res * res), 1.0, 1)) + + delta * std::sqrt(alpha * alpha - beta * beta) + + beta * res - 0.5 * std::log(deltasq + res * res); + + pdf = std::exp(pdf); + + return pdf; +} + +// [[Rcpp::export]] +NumericVector c_dnig(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr) { + int n = x.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = dnig_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double rstd_std(const double nu) { + if (nu > 2.0) { + double s = std::sqrt(nu / (nu - 2.0)); + return R::rt(nu) * 1.0 / s; + } + return 0.0; // Default to 0 if nu <= 2.0 (undefined) +} + +// [[Rcpp::export]] +NumericVector c_rstd(int n, NumericVector mu, NumericVector sigma, NumericVector shape) { + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = mu[i] + rstd_std(shape[i]) * sigma[i]; + } + + return ans; +} + +double xdt(const double x, const double nu) { + double a = R::gammafn((nu + 1.0) / 2.0) / std::sqrt(M_PI * nu); + double b = R::gammafn(nu / 2.0) * std::pow((1.0 + (x * x) / nu), (nu + 1.0) / 2.0); + return a / b; +} + +double dstd_std(const double x, const double nu) { + if (nu <= 2.0) { + return 999.0; // Undefined behavior for nu <= 2.0 + } else { + double s = std::sqrt(nu / (nu - 2.0)); + return s * xdt(x * s, nu); + } +} + +// [[Rcpp::export]] +NumericVector c_dstd(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector shape, int logr) { + int n = x.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = dstd_std((x[i] - mu[i]) / sigma[i], shape[i]) / sigma[i]; + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double pstd_std(const double q, const double mu, const double sigma, const double nu) { + double s = std::sqrt(nu / (nu - 2.0)); + double z = (q - mu) / sigma; + double p = R::pt(z * s, nu, 1, 0); // 1 for lower.tail = TRUE, 0 for log.p = FALSE + return p; +} + +// [[Rcpp::export]] +NumericVector c_pstd(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector shape) { + int n = q.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = pstd_std(q[i], mu[i], sigma[i], shape[i]); + } + + return ans; +} + +double qstd_std(const double p, const double mu, const double sigma, const double nu) { + double s = std::sqrt(nu / (nu - 2.0)); + double q = R::qt(p, nu, 1, 0) * sigma / s + mu; // 1 for lower.tail = TRUE, 0 for log.p = FALSE + return q; +} + +// [[Rcpp::export]] +NumericVector c_qstd(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector shape) { + int n = p.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = qstd_std(p[i], mu[i], sigma[i], shape[i]); + } + + return ans; +} + +double rsstd_std(const double xi, const double nu) { + double weight = xi / (xi + 1.0 / xi); + double z = R::runif(-1.0 * weight, 1.0 - weight); + double xx = (z < 0) ? 1.0 / xi : xi; + double rr = -1.0 * std::fabs(rstd_std(nu)) / xx * ((z < 0) ? -1.0 : 1.0); + + double a = 0.5, b = nu / 2.0; + double beta_val = (R::gammafn(a) / R::gammafn(a + b)) * R::gammafn(b); + + double m1 = 2.0 * std::sqrt(nu - 2.0) / (nu - 1.0) / beta_val; + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + + return (rr - mu) / sigma; +} + +// [[Rcpp::export]] +NumericVector c_rsstd(int n, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + NumericVector ans(n); + for(int i = 0; i < n; i++) { + ans[i] = mu[i] + rsstd_std(skew[i], shape[i]) * sigma[i]; + } + return ans; +} + +double dsstd_std(const double x, const double xi, const double nu) { + double a = 0.5, b = nu / 2.0; + double beta_val = (R::gammafn(a) / R::gammafn(a + b)) * R::gammafn(b); + + double m1 = 2.0 * std::sqrt(nu - 2.0) / (nu - 1.0) / beta_val; + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + + double z = x * sigma + mu; + double xxi = (z < 0) ? 1.0 / xi : xi; + + double g = 2.0 / (xi + 1.0 / xi); + double pdf = g * dstd_std(z / xxi, nu) * sigma; + + return pdf; +} + +// [[Rcpp::export]] +NumericVector c_dsstd(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr) { + int n = x.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = dsstd_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double psstd_std(const double q, const double mu, const double sigma, const double xi, const double nu) { + double qx = (q - mu) / sigma; + + // Calculate m1, mu, and sigma based on xi and nu + double m1 = 2.0 * std::sqrt(nu - 2.0) / (nu - 1.0) / (R::beta(0.5, nu / 2.0)); + double mux = m1 * (xi - 1.0 / xi); + double sig = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + + double z = qx * sig + mux; + double Xi = (z < 0) ? 1.0 / xi : xi; + double g = 2.0 / (xi + 1.0 / xi); + + // Compute the CDF using the heaviside and signum functions + double p = heaviside(z, 0) - signum(z) * g * Xi * pstd_std(-std::fabs(z) / Xi, 0.0, 1.0, nu); + return p; +} + +// [[Rcpp::export]] +NumericVector c_psstd(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + int n = q.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = psstd_std(q[i], mu[i], sigma[i], skew[i], shape[i]); + } + + return ans; +} + +double qsstd_std(const double p, const double xi, const double nu) { + double m1 = 2.0 * std::sqrt(nu - 2.0) / (nu - 1.0) / (R::beta(0.5, nu / 2.0)); + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + double g = 2.0 / (xi + 1.0 / xi); + + double z = p - (1.0 / (1.0 + xi * xi)); + double Xi = std::pow(xi, signum(z)); + + double tmp = (heaviside(z, 0) - signum(z) * p) / (g * Xi); + double q = (-signum(z) * qstd_std(tmp, 0.0, 1.0, nu) * Xi - mu) / sigma; + + return q; +} + +// [[Rcpp::export]] +NumericVector c_qsstd(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + int n = p.size(); + NumericVector ans(n); + + for(int i = 0; i < n; i++) { + ans[i] = qsstd_std(p[i], skew[i], shape[i]) * sigma[i] + mu[i]; + } + + return ans; +} + +double djsu_std(const double x, const double nu, const double tau) { + double w, z, r, omega, c, pdf; + double rtau = 1.0 / tau; + + if (rtau < 1e-7) { + w = 1.0; + } else { + w = std::exp(rtau * rtau); + } + + omega = -nu * rtau; + c = std::sqrt(1.0 / (0.5 * (w - 1.0) * (w * std::cosh(2.0 * omega) + 1.0))); + z = (x - (c * std::sqrt(w) * std::sinh(omega))) / c; + r = -nu + std::asinh(z) / rtau; + + pdf = -std::log(c) - std::log(rtau) - 0.5 * std::log(z * z + 1.0) - 0.5 * std::log(2.0 * M_PI) - 0.5 * r * r; + return std::exp(pdf); +} + +// [[Rcpp::export]] +NumericVector c_djsu(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr) { + int n = x.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = djsu_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double qjsu_std(const double p, const double nu, const double tau) { + double rtau = 1.0 / tau; + double rr = R::qnorm(p, 0.0, 1.0, 1, 0); + double z = std::sinh(rtau * (rr + nu)); + + double w = (rtau < 1e-7) ? 1.0 : std::exp(rtau * rtau); + double omega = -nu * rtau; + double cc = std::sqrt(1.0 / (0.5 * (w - 1.0) * (w * std::cosh(2.0 * omega) + 1.0))); + + return (cc * std::sqrt(w) * std::sinh(omega)) + cc * z; +} + +// [[Rcpp::export]] +NumericVector c_qjsu(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + int n = p.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = mu[i] + qjsu_std(p[i], skew[i], shape[i]) * sigma[i]; + } + + return ans; +} + +double pjsu_std(const double q, const double mu, const double sigma, const double nu, const double tau) { + double rtau = 1.0 / tau; + double w = (rtau < 1e-7) ? 1.0 : std::exp(rtau * rtau); + double omega = -nu * rtau; + double c = 1.0 / std::sqrt(0.5 * (w - 1.0) * (w * std::cosh(2.0 * omega) + 1.0)); + + double z = (q - (mu + c * sigma * std::sqrt(w) * std::sinh(omega))) / (c * sigma); + double r = -nu + std::asinh(z) / rtau; + + return R::pnorm(r, 0, 1, 1, 0); // pnorm(r, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) +} + +// [[Rcpp::export]] +NumericVector c_pjsu(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + int n = q.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = pjsu_std(q[i], mu[i], sigma[i], skew[i], shape[i]); + } + + return ans; +} + + +double rjsu_std(const double nu, const double tau) { + double x = R::runif(0.0, 1.0); + return qjsu_std(x, nu, tau); +} + +// [[Rcpp::export]] +NumericVector c_rjsu(int n, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + NumericVector ans(n); + for (int i = 0; i < n; i++) { + ans[i] = mu[i] + rjsu_std(skew[i], shape[i]) * sigma[i]; + } + return ans; +} + +double rsnorm_std(const double xi) { + double weight = xi / (xi + 1.0 / xi); + double z = R::runif(-weight, 1.0 - weight); + double xx = (z < 0) ? 1.0 / xi : xi; + double rr = -1.0 * std::fabs(R::rnorm(0.0, 1.0)) / xx * ((z < 0) ? -1.0 : 1.0); + + double m1 = 2.0 / std::sqrt(2.0 * M_PI); + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m1 * m1) * ((xi * xi) + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + + return (rr - mu) / sigma; +} + +// [[Rcpp::export]] +NumericVector c_rsnorm(int n, NumericVector mu, NumericVector sigma, NumericVector skew) { + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = mu[i] + rsnorm_std(skew[i]) * sigma[i]; + } + + return ans; +} + +double dsnorm_std(const double x, const double xi) { + double m1 = 2.0 / std::sqrt(2.0 * M_PI); + double m12 = m1 * m1; + double xi2 = xi * xi; + + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m12) * (xi2 + 1.0 / xi2) + 2.0 * m12 - 1.0); + + double z = x * sigma + mu; + double xxi = (z == 0.0) ? 1.0 : ((z < 0.0) ? 1.0 / xi : xi); + + double g = 2.0 / (xi + 1.0 / xi); + double pdf = g * dnorm_std(z / xxi) * sigma; + + return pdf; +} + +// [[Rcpp::export]] +NumericVector c_dsnorm(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, int logr) { + int n = x.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = dsnorm_std((x[i] - mu[i]) / sigma[i], skew[i]) / sigma[i]; + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double psnorm_std(const double q, const double mu, const double sigma, const double xi) { + double qx = (q - mu) / sigma; + double m1 = 2.0 / std::sqrt(2.0 * M_PI); + double mux = m1 * (xi - 1.0 / xi); + double sig = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + + double z = qx * sig + mux; + double Xi = (z < 0) ? 1.0 / xi : xi; + double g = 2.0 / (xi + 1.0 / xi); + + double p = heaviside(z, 0) - signum(z) * g * Xi * R::pnorm(-std::fabs(z) / Xi, 0.0, 1.0, 1, 0); + return p; +} + +// [[Rcpp::export]] +NumericVector c_psnorm(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector skew) { + int n = q.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = psnorm_std(q[i], mu[i], sigma[i], skew[i]); + } + + return ans; +} + +double qsnorm_std(const double p, const double xi) { + double m1 = 2.0 / std::sqrt(2.0 * M_PI); + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + + double g = 2.0 / (xi + 1.0 / xi); + double z = p - (1.0 / (1.0 + xi * xi)); + double Xi = std::pow(xi, signum(z)); + + double tmp = (heaviside(z, 0) - signum(z) * p) / (g * Xi); + double q = (-signum(z) * R::qnorm(tmp, 0.0, Xi, 1, 0) - mu) / sigma; + + return q; +} + +// [[Rcpp::export]] +NumericVector c_qsnorm(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector skew) { + int n = p.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = mu[i] + qsnorm_std(p[i], skew[i]) * sigma[i]; + } + + return ans; +} + +double rged_std(const double nu) { + double lambda = std::sqrt(std::pow(0.5, 2.0 / nu) * R::gammafn(1.0 / nu) / R::gammafn(3.0 / nu)); + double rr = R::rgamma(1.0 / nu, 1.0); + double ans = lambda * std::pow(2.0 * rr, 1.0 / nu) * ((R::runif(0.0, 1.0) < 0.5) ? -1.0 : 1.0); + return ans; +} + +// [[Rcpp::export]] +NumericVector c_rged(int n, NumericVector mu, NumericVector sigma, NumericVector shape) { + NumericVector ans(n); + for (int i = 0; i < n; i++) { + ans[i] = mu[i] + rged_std(shape[i]) * sigma[i]; + } + return ans; +} + +double dged_std(const double x, const double nu) { + double lambda = std::sqrt(std::pow(1.0 / 2.0, 2.0 / nu) * R::gammafn(1.0 / nu) / R::gammafn(3.0 / nu)); + double g = nu / (lambda * (std::pow(2.0, 1.0 + (1.0 / nu))) * R::gammafn(1.0 / nu)); + double pdf = g * std::exp(-0.5 * std::pow(std::fabs(x / lambda), nu)); + return pdf; +} + +// [[Rcpp::export]] +NumericVector c_dged(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector shape, int logr) { + int n = x.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = dged_std((x[i] - mu[i]) / sigma[i], shape[i]) / sigma[i]; + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double pged_std(const double q, const double mu, const double sigma, const double nu) { + double qx = (q - mu) / sigma; + double lambda = std::sqrt(1.0 / std::pow(2.0, 2.0 / nu) * R::gammafn(1.0 / nu) / R::gammafn(3.0 / nu)); + double g = nu / (lambda * (std::pow(2.0, 1.0 + 1.0 / nu)) * R::gammafn(1.0 / nu)); + double h = std::pow(2.0, 1.0 / nu) * lambda * g * R::gammafn(1.0 / nu) / nu; + double s = 0.5 * std::pow(std::fabs(qx) / lambda, nu); + double p = 0.5 + ((qx >= 0.0) ? 1.0 : -1.0) * h * R::pgamma(s, 1.0 / nu, 1.0, 1, 0); + return p; +} + +// [[Rcpp::export]] +NumericVector c_pged(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector shape) { + int n = q.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = pged_std(q[i], mu[i], sigma[i], shape[i]); + } + + return ans; +} + +double qged_std(const double p, const double shape) { + double y = 2.0 * p - 1.0; + double lambda = std::sqrt(std::pow(0.5, 2.0 / shape) * R::gammafn(1.0 / shape) / R::gammafn(3.0 / shape)); + double q = lambda * std::pow(2.0 * R::qgamma(std::fabs(y), 1.0 / shape, 1.0, 1, 0), 1.0 / shape); + q = q * ((y >= 0) ? 1.0 : -1.0); + return q; +} + +// [[Rcpp::export]] +NumericVector c_qged(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector shape) { + int n = p.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = qged_std(p[i], shape[i]) * sigma[i] + mu[i]; + } + + return ans; +} + +double rsged_std(const double xi, const double nu) { + double weight = xi / (xi + 1.0 / xi); + double z = R::runif(-weight, 1.0 - weight); + double xx = (z < 0) ? 1.0 / xi : xi; + double rr = -1.0 * std::fabs(rged_std(nu)) / xx * ((z < 0) ? -1.0 : 1.0); + + double lambda = std::sqrt(std::pow(0.5, 2.0 / nu) * R::gammafn(1.0 / nu) / R::gammafn(3.0 / nu)); + double m1 = std::pow(2.0, 1.0 / nu) * lambda * R::gammafn(2.0 / nu) / R::gammafn(1.0 / nu); + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + + return (rr - mu) / sigma; +} + +// [[Rcpp::export]] +NumericVector c_rsged(int n, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = mu[i] + rsged_std(skew[i], shape[i]) * sigma[i]; + } + + return ans; +} + +double dsged_std(const double x, const double xi, const double nu) { + double lambda = std::sqrt(std::pow(1.0 / 2.0, 2.0 / nu) * R::gammafn(1.0 / nu) / R::gammafn(3.0 / nu)); + double m1 = std::pow(2.0, 1.0 / nu) * lambda * R::gammafn(2.0 / nu) / R::gammafn(1.0 / nu); + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + + double z = x * sigma + mu; + double xxi = (z == 0) ? 1.0 : ((z < 0) ? 1.0 / xi : xi); + + double g = 2.0 / (xi + 1.0 / xi); + double pdf = g * dged_std(z / xxi, nu) * sigma; + + return pdf; +} + +// [[Rcpp::export]] +NumericVector c_dsged(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr) { + int n = x.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = dsged_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + + +double psged_std(const double q, const double mu, const double sigma, const double xi, const double nu) { + double qx = (q - mu) / sigma; + double lambda = std::sqrt(1.0 / std::pow(2.0, 2.0 / nu) * R::gammafn(1.0 / nu) / R::gammafn(3.0 / nu)); + double m1 = std::pow(2.0, 1.0 / nu) * lambda * R::gammafn(2.0 / nu) / R::gammafn(1.0 / nu); + double mux = m1 * (xi - 1.0 / xi); + double sig = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + double z = qx * sig + mux; + double Xi = (z < 0) ? 1.0 / xi : xi; + double g = 2.0 / (xi + 1.0 / xi); + double p = heaviside(z, 0) - signum(z) * g * Xi * pged_std(-std::fabs(z) / Xi, 0.0, 1.0, nu); + return p; +} + +// [[Rcpp::export]] +NumericVector c_psged(NumericVector q, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + int n = q.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = psged_std(q[i], mu[i], sigma[i], skew[i], shape[i]); + } + + return ans; +} + +double qsged_std(const double p, const double xi, const double nu) { + double lambda = std::sqrt(1.0 / std::pow(2.0, 2.0 / nu) * R::gammafn(1.0 / nu) / R::gammafn(3.0 / nu)); + double m1 = std::pow(2.0, 1.0 / nu) * lambda * R::gammafn(2.0 / nu) / R::gammafn(1.0 / nu); + double mu = m1 * (xi - 1.0 / xi); + double sigma = std::sqrt((1.0 - m1 * m1) * (xi * xi + 1.0 / (xi * xi)) + 2.0 * m1 * m1 - 1.0); + double g = 2.0 / (xi + 1.0 / xi); + double z = p - (1.0 / (1.0 + xi * xi)); + double Xi = std::pow(xi, signum(z)); + double q = (heaviside(z, 0) - signum(z) * p) / (g * Xi); + q = (-signum(z) * qged_std(q, nu) * Xi - mu) / sigma; + + return q; +} + +// [[Rcpp::export]] +NumericVector c_qsged(NumericVector p, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape) { + int n = p.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = qsged_std(p[i], skew[i], shape[i]) * sigma[i] + mu[i]; + } + + return ans; +} + + +double dhyp(const double x, const double alpha, const double beta, const double delta, const double mu) { + double pdf = 0.0; + if (alpha <= 0.0 || delta <= 0 || std::fabs(beta) >= alpha) { + return pdf; + } + + double g = alpha * alpha - beta * beta; + double e = x - mu; + double besselK = R::bessel_k(delta * std::sqrt(g), 1.0, 2); + + if (besselK == 0.0) return 0.0; + + pdf = 0.5 * std::log(g) - std::log(2.0 * alpha * delta * besselK) - alpha * std::sqrt(delta * delta + e * e) + beta * e; + return std::exp(pdf); +} + +double dhyp_std(const double x, const double rho, const double zeta) { + NumericVector param = paramgh(rho, zeta, 1.0); + double alpha = param[0]; + double beta = param[1]; + double delta = param[2]; + double mu = param[3]; + + return dhyp(x, alpha, beta, delta, mu); +} + +// [[Rcpp::export]] +NumericVector c_dhyp(NumericVector x, NumericVector mu, NumericVector sigma, NumericVector skew, NumericVector shape, int logr) { + int n = x.size(); + NumericVector ans(n); + + for (int i = 0; i < n; i++) { + ans[i] = dhyp_std((x[i] - mu[i]) / sigma[i], skew[i], shape[i]) / sigma[i]; + if (logr == 1) { + ans[i] = std::log(ans[i]); + } + } + + return ans; +} + +double g_function(double y, double beta, double m, double lambda) { + return 0.5 * beta * std::pow(y, 3) - std::pow(y, 2) * (0.5 * beta * m + lambda + 1) + y * (-0.5 * beta) + 0.5 * beta * m; +} + +double f_function(double y, double beta, double m, double lambda) { + return 0.5 * beta * std::pow(y, 3) - std::pow(y, 2) * (0.5 * beta * m + lambda + 1) + + y * ((lambda - 1) * m - 0.5 * beta) + 0.5 * beta * m; +} + +NumericVector rgig1(int n, NumericVector param) { + double chi_val = param[0]; + double psi_val = param[1]; + double lambda = 1.0; + + double alpha = std::sqrt(psi_val / chi_val); + double beta = std::sqrt(psi_val * chi_val); + double m = std::abs(beta) / beta; + + Rcpp::Function uniroot("uniroot"); + + List yM_result = uniroot(Named("f") = Rcpp::InternalFunction(&g_function), + Named("interval") = NumericVector::create(0.0, m), + Named("beta") = beta, + Named("m") = m, + Named("lambda") = lambda); + + double lower2 = m; + double upper2 = 2 * m; + + while (g_function(lower2, beta, m, lambda) * g_function(upper2, beta, m, lambda) > 0) { + upper2 *= 2; + if (upper2 > 1000) { + stop("No root found for yP within a reasonable interval."); + } + } + + List yP_result = uniroot(Named("f") = Rcpp::InternalFunction(&g_function), + Named("interval") = NumericVector::create(lower2, upper2), + Named("beta") = beta, + Named("m") = m, + Named("lambda") = lambda); + + double yM = as(yM_result["root"]); + double yP = as(yP_result["root"]); + + double a = (yP - m) * std::exp(-0.25 * beta * (yP + 1 / yP - m - 1 / m)); + double b = (yM - m) * std::exp(-0.25 * beta * (yM + 1 / yM - m - 1 / m)); + double c = -0.25 * beta * (m + 1 / m); + + NumericVector output(n); + + for (int i = 0; i < n; i++) { + bool needValue = true; + + while (needValue) { + double R1 = R::runif(0.0, 1.0); + double R2 = R::runif(0.0, 1.0); + double Y = m + a * R2 / R1 + b * (1 - R2) / R1; + + if (Y > 0) { + if (-std::log(R1) >= 0.25 * beta * (Y + 1 / Y) + c) { + needValue = false; + output[i] = Y; + } + } + } + } + + return output / alpha; +} + + +NumericVector rgig(int n, NumericVector param) { + if (param.size() != 3) { + stop("Parameter vector must have exactly 3 elements."); + } + double chi = param[0]; + double psi = param[1]; + double lambda = param[2]; + double alpha = std::sqrt(psi / chi); + double beta = std::sqrt(psi * chi); + double m = (lambda - 1 + std::sqrt((lambda - 1) * (lambda - 1) + beta * beta)) / beta; + Rcpp::Function uniroot("uniroot"); + double upper = m; + while (f_function(upper, beta, m, lambda) <= 0) { + upper *= 2; + } + List yM_result = uniroot(Named("f") = Rcpp::InternalFunction(&f_function), + Named("interval") = NumericVector::create(0.0, m), + Named("beta") = beta, + Named("m") = m, + Named("lambda") = lambda); + + List yP_result = uniroot(Named("f") = Rcpp::InternalFunction(&f_function), + Named("interval") = NumericVector::create(m, upper), + Named("beta") = beta, + Named("m") = m, + Named("lambda") = lambda); + + double yM = as(yM_result["root"]); + double yP = as(yP_result["root"]); + + double a = (yP - m) * std::pow(yP / m, 0.5 * (lambda - 1)) * std::exp(-0.25 * beta * (yP + 1 / yP - m - 1 / m)); + double b = (yM - m) * std::pow(yM / m, 0.5 * (lambda - 1)) * std::exp(-0.25 * beta * (yM + 1 / yM - m - 1 / m)); + double c = -0.25 * beta * (m + 1 / m) + 0.5 * (lambda - 1) * std::log(m); + + NumericVector output(n); + + for (int i = 0; i < n; i++) { + bool needValue = true; + + while (needValue) { + double R1 = R::runif(0.0, 1.0); + double R2 = R::runif(0.0, 1.0); + double Y = m + a * R2 / R1 + b * (1 - R2) / R1; + + if (Y > 0) { + if (-std::log(R1) >= -0.5 * (lambda - 1) * std::log(Y) + 0.25 * beta * (Y + 1 / Y) + c) { + needValue = false; + output[i] = Y; + } + } + } + } + + return output / alpha; +} + +// [[Rcpp::export]] +NumericVector c_rghyp(int n, double mu = 0, double delta = 1, double alpha = 1, double beta = 0, double lambda = 1) { + double chi = delta * delta; + double psi = alpha * alpha - beta * beta; + NumericVector X; + if (lambda == 1) { + X = rgig1(n, NumericVector::create(chi, psi)); + } else { + X = rgig(n, NumericVector::create(chi, psi, lambda)); + } + NumericVector sigma = sqrt(X); + NumericVector Z = rnorm(n); + NumericVector Y = mu + beta * pow(sigma, 2) + sigma * Z; + return Y; +} diff --git a/src/distributions.h b/src/distributions.h index 18eab33..9cbf79b 100644 --- a/src/distributions.h +++ b/src/distributions.h @@ -1,131 +1,84 @@ #ifndef DISTRIBUTIONS_H #define DISTRIBUTIONS_H -/* - * ----------------------------------------- - * Key Functions - * ----------------------------------------- - */ -double dnorm_std(const double ); -double signum(const double ); + +#include +#include +#include + +// Declare all the functions +double dnorm_std(const double x); +double signum(const double x); double dhuge(void); -double Heaviside(const double , const double ); +double heaviside(const double x, const double a); double depsilon(void); -double kappagh(const double , const double ); -double deltakappagh(const double , const double ); -double* paramgh(const double , const double , const double ); -double* paramghskt(const double , const double ); -/* - * ----------------------------------------- - * GH Skew Student Distribution - * ----------------------------------------- - */ -double dghst_std(const double , const double , const double ); -void c_dghst(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int* logr); -double rsghst_std(const double , const double ); -void c_rghst(int *n, double *mu, double *sigma, double *skew, double *shape, double *ans); -/* - * ----------------------------------------- - * GH Distribution - * ----------------------------------------- - */ -// alpha-beta-delta-mu parameterization -double dgh(const double , const double , const double , const double , const double , const double); -void c_dghyp(double *, double *, double *, double *, double *, double *, double *, int *, int *); -// rho-zeta parameterization -double dgh_std(const double , const double , const double , const double); -void c_dgh(double *x, double *mu, double *sigma, double *skew, double *shape, double *lambda, double *ans, int *n, int *logr); -/* - * ----------------------------------------- - * NIG Distribution - * ----------------------------------------- - */ -double dnig_std(const double , const double , const double); -void c_dnig(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr); -/* - * ----------------------------------------- - * Student Distribution - * ----------------------------------------- - */ -double rstd_std(const double ); -void c_rstd(int *n, double *mu, double *sigma, double *shape, double *ans); -double xdt(const double , const double ); -double dstd_std(const double , const double ); -void c_dstd(double *x, double *mu, double *sigma, double *shape, double *ans, int *n, int *logr); -double pstd_std(const double , const double , const double , const double ); -void c_pstd(double *q, double *mu, double *sigma, double *shape, double *ans, int *n); -double qstd_std(const double , const double , const double , const double ); -void c_qstd(double *p, double *mu, double *sigma, double *shape, double *ans, int *n); -/* - * ----------------------------------------- - * Skew Student Distribution (Fernandez & Steel) - * ----------------------------------------- - */ -double rsstd_std(const double , const double ); -void c_rsstd(int *n, double *mu, double *sigma, double *skew, double *shape, double *ans); -double dsstd_std(const double , const double , const double ); -void c_dsstd(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr); -double psstd_std(const double , const double , const double , const double , const double ); -void c_psstd(double *q, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n); -double qsstd_std(const double , const double , const double ); -void c_qsstd(double *p, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n); -/* - * ----------------------------------------- - * Johnson's SU Distribution - * ----------------------------------------- - */ -double djsu_std(const double , const double , const double ); -void c_djsu(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr); -double qjsu_std(const double , const double , const double ); -void c_qjsu(double *p, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n); -double rjsu_std(const double n, const double ); -void c_rjsu(int *n, double *mu, double *sigma, double *skew, double *shape, double *ans); -double pjsu_std(const double , const double , const double , const double , const double ); -void c_pjsu(double *q, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n); -/* - * ----------------------------------------- - * Skew Normal Distribution - * ----------------------------------------- - */ -double rsnorm_std(const double ); -void c_rsnorm(int *n, double *mu, double *sigma, double *skew, double *ans); -double dsnorm_std(const double , const double ); -void c_dsnorm(double *x, double *mu, double *sigma, double *skew, double *ans, int *n, int *logr); -double psnorm_std(const double , const double , const double , const double ); -void c_psnorm(double *q, double *mu, double *sigma, double *skew, double *ans, int *n); -double qsnorm_std(const double , const double ); -void c_qsnorm(double *p, double *mu, double *sigma, double *skew, double *ans, int *n); -/* - * ----------------------------------------- - * Generalized Error Distribution - * ----------------------------------------- - */ -double rged_std(const double ); -void c_rged(int *n, double *mu, double *sigma, double *shape, double *ans); -double dged_std(const double , const double ); -void c_dged(double *x, double *mu, double *sigma, double *shape, double *ans, int *n, int *logr); -double pged_std(const double , const double , const double , const double ); -void c_pged(double *q, double *mu, double *sigma, double *shape, double *ans, int *n); -double qged_std(const double , const double ); -void c_qged(double *p, double *mu, double *sigma, double *shape, double *ans, int *n); -/* - * ----------------------------------------- - * Skew Generalized Error Distribution (Fernandez & Steel) - * ----------------------------------------- - */ -double rsged_std(const double , const double ); -void c_rsged(int *n, double *mu, double *sigma, double *skew, double *shape, double *ans); -double dsged_std(const double , const double , const double ); -void c_dsged(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr); -double psged_std(const double , const double , const double , const double , const double ); -void c_psged(double *q, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n); -double qsged_std(const double , const double, const double ); -void c_qsged(double *p, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n); -/* - * ----------------------------------------- - * Hypebolic Distribution - * ----------------------------------------- - */ -double dhyp(const double , const double , const double , const double , const double); -double dhyp_std(const double , const double , const double); -void c_dhyp(double *x, double *mu, double *sigma, double *skew, double *shape, double *ans, int *n, int *logr); -#endif /* DISTRIBUTIONS_H */ +double kappagh(const double x, const double lambda); +double deltakappagh(const double x, const double lambda); +Rcpp::NumericVector paramgh(const double rho, const double zeta, const double lambda); +Rcpp::NumericVector paramghst(const double betabar, const double nu); +double dghst_std(const double x, const double betabar, const double nu); +Rcpp::NumericVector c_dghst(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape, int logr); +double rsghst_std(const double betabar, const double nu); +Rcpp::NumericVector c_rghst(int n, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double dgh(const double x, const double alpha, const double beta, const double delta, const double mu, const double lambda); +Rcpp::NumericVector c_dghyp(Rcpp::NumericVector x, Rcpp::NumericVector alpha, Rcpp::NumericVector beta, Rcpp::NumericVector delta, Rcpp::NumericVector mu, Rcpp::NumericVector lambda, int logr); +double dgh_std(const double x, const double rho, const double zeta, const double lambda); +Rcpp::NumericVector c_dgh(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape, Rcpp::NumericVector lambda, int logr); +double dnig_std(const double x, const double rho, const double zeta); +Rcpp::NumericVector c_dnig(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape, int logr); +double rstd_std(const double nu); +Rcpp::NumericVector c_rstd(int n, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector shape); +double xdt(const double x, const double nu); +double dstd_std(const double x, const double nu); +Rcpp::NumericVector c_dstd(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector shape, int logr); +double pstd_std(const double q, const double mu, const double sigma, const double nu); +Rcpp::NumericVector c_pstd(Rcpp::NumericVector q, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector shape); +double qstd_std(const double p, const double mu, const double sigma, const double nu); +Rcpp::NumericVector c_qstd(Rcpp::NumericVector p, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector shape); +double rsstd_std(const double xi, const double nu); +Rcpp::NumericVector c_rsstd(int n, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double dsstd_std(const double x, const double xi, const double nu); +Rcpp::NumericVector c_dsstd(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape, int logr); +double psstd_std(const double q, const double mu, const double sigma, const double xi, const double nu); +Rcpp::NumericVector c_psstd(Rcpp::NumericVector q, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double qsstd_std(const double p, const double xi, const double nu); +Rcpp::NumericVector c_qsstd(Rcpp::NumericVector p, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double djsu_std(const double x, const double nu, const double tau); +Rcpp::NumericVector c_djsu(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape, int logr); +double qjsu_std(const double p, const double nu, const double tau); +Rcpp::NumericVector c_qjsu(Rcpp::NumericVector p, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double pjsu_std(const double q, const double mu, const double sigma, const double nu, const double tau); +Rcpp::NumericVector c_pjsu(Rcpp::NumericVector q, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double rsnorm_std(const double xi); +Rcpp::NumericVector c_rsnorm(int n, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew); +double dsnorm_std(const double x, const double xi); +Rcpp::NumericVector c_dsnorm(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, int logr); +double psnorm_std(const double q, const double mu, const double sigma, const double xi); +Rcpp::NumericVector c_psnorm(Rcpp::NumericVector q, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew); +double qsnorm_std(const double p, const double xi); +Rcpp::NumericVector c_qsnorm(Rcpp::NumericVector p, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew); +double rged_std(const double nu); +Rcpp::NumericVector c_rged(int n, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector shape); +double dged_std(const double x, const double nu); +Rcpp::NumericVector c_dged(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector shape, int logr); +double pged_std(const double q, const double mu, const double sigma, const double nu); +Rcpp::NumericVector c_pged(Rcpp::NumericVector q, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector shape); +double qged_std(const double p, const double shape); +Rcpp::NumericVector c_qged(Rcpp::NumericVector p, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector shape); +double rsged_std(const double xi, const double nu); +Rcpp::NumericVector c_rsged(int n, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double dsged_std(const double x, const double xi, const double nu); +Rcpp::NumericVector c_dsged(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape, int logr); +double psged_std(const double q, const double mu, const double sigma, const double xi, const double nu); +Rcpp::NumericVector c_psged(Rcpp::NumericVector q, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double qsged_std(const double p, const double xi, const double nu); +Rcpp::NumericVector c_qsged(Rcpp::NumericVector p, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape); +double dhyp(const double x, const double alpha, const double beta, const double delta, const double mu); +double dhyp_std(const double x, const double rho, const double zeta); +Rcpp::NumericVector c_dhyp(Rcpp::NumericVector x, Rcpp::NumericVector mu, Rcpp::NumericVector sigma, Rcpp::NumericVector skew, Rcpp::NumericVector shape, int logr); +double g_function(double y, double beta, double m, double lambda); +double f_function(double y, double beta, double m, double lambda); +Rcpp::NumericVector rgig1(int n, Rcpp::NumericVector param); +Rcpp::NumericVector rgig(int n, Rcpp::NumericVector param); +Rcpp::NumericVector c_rghyp(int n, double mu, double delta, double alpha, double beta, double lambda); + +#endif diff --git a/src/tsdistributions_init.c b/src/tsdistributions_init.c deleted file mode 100644 index 8ae1b78..0000000 --- a/src/tsdistributions_init.c +++ /dev/null @@ -1,79 +0,0 @@ -// Dummy file required so that useDynLib(tsdistributions, .registration=TRUE) doesn't fail on empty 'src' -#include -#include -#include -#include // for NULL -#include -/* .C calls */ -extern void c_dghst(double *, double *, double *, double *, double *, double *, int *, int*); -extern void c_rghst(int *, double *, double *, double *, double *, double *); -extern void c_dghyp(double *, double *, double *, double *, double *, double *, double *, int *, int *); -extern void c_dgh(double *, double *, double *, double *, double *, double *, double *, int *, int *); -extern void c_rgh(int *, double *, double *, double *, double *, double *, double *); -extern void c_dnig(double *, double *, double *, double *, double *, double *, int *, int *); -extern void c_rnig(int *, double *, double *, double *, double *, double *); -extern void c_rstd(int *, double *, double *, double *, double *); -extern void c_dstd(double *, double *, double *, double *, double *, int *, int *); -extern void c_pstd(double *, double *, double *, double *, double *, int *); -extern void c_qstd(double *, double *, double *, double *, double *, int *); -extern void c_rsstd(int *, double *, double *, double *, double *, double *); -extern void c_dsstd(double *, double *, double *, double *, double *, double *, int *, int *); -extern void c_psstd(double *, double *, double *, double *, double *, double *, int *); -extern void c_qsstd(double *, double *, double *, double *, double *, double *, int *); -extern void c_djsu(double *, double *, double *, double *, double *, double *, int *, int *); -extern void c_qjsu(double *, double *, double *, double *, double *, double *, int *); -extern void c_rjsu(int *, double *, double *, double *, double *, double *); -extern void c_pjsu(double *, double *, double *, double *, double *, double *, int *); -extern void c_rsnorm(int *, double *, double *, double *, double *); -extern void c_dsnorm(double *, double *, double *, double *, double *, int *, int *); -extern void c_psnorm(double *, double *, double *, double *, double *, int *); -extern void c_qsnorm(double *, double *, double *, double *, double *, int *); -extern void c_rged(int *, double *, double *, double *, double *); -extern void c_dged(double *, double *, double *, double *, double *, int *, int *); -extern void c_pged(double *, double *, double *, double *, double *, int *); -extern void c_qged(double *, double *, double *, double *, double *, int *); -extern void c_rsged(int *, double *, double *, double *, double *, double *); -extern void c_dsged(double *, double *, double *, double *, double *, double *, int *, int *); -extern void c_psged(double *, double *, double *, double *, double *, double *, int *); -extern void c_qsged(double *, double *, double *, double *, double *, double *, int *); -extern void c_dhyp(double *, double *, double *, double *, double *, double *, int *, int *); - - -static const R_CMethodDef CEntries[] = { - {"c_dghst", (DL_FUNC) &c_dghst, 8}, - {"c_rghst", (DL_FUNC) &c_rghst, 6}, - {"c_dghyp", (DL_FUNC) &c_dghyp, 9}, - {"c_dgh", (DL_FUNC) &c_dgh, 9}, - {"c_dnig", (DL_FUNC) &c_dnig, 8}, - {"c_rstd", (DL_FUNC) &c_rstd, 5}, - {"c_dstd", (DL_FUNC) &c_dstd, 7}, - {"c_pstd", (DL_FUNC) &c_pstd, 6}, - {"c_qstd", (DL_FUNC) &c_qstd, 6}, - {"c_rsstd", (DL_FUNC) &c_rsstd, 6}, - {"c_dsstd", (DL_FUNC) &c_dsstd, 8}, - {"c_psstd", (DL_FUNC) &c_psstd, 7}, - {"c_qsstd", (DL_FUNC) &c_qsstd, 7}, - {"c_djsu", (DL_FUNC) &c_djsu, 8}, - {"c_qjsu", (DL_FUNC) &c_qjsu, 7}, - {"c_rjsu", (DL_FUNC) &c_rjsu, 6}, - {"c_pjsu", (DL_FUNC) &c_pjsu, 7}, - {"c_rsnorm", (DL_FUNC) &c_rsnorm, 5}, - {"c_dsnorm", (DL_FUNC) &c_dsnorm, 7}, - {"c_psnorm", (DL_FUNC) &c_psnorm, 6}, - {"c_qsnorm", (DL_FUNC) &c_qsnorm, 6}, - {"c_rged", (DL_FUNC) &c_rged, 5}, - {"c_dged", (DL_FUNC) &c_dged, 7}, - {"c_pged", (DL_FUNC) &c_pged, 6}, - {"c_qged", (DL_FUNC) &c_qged, 6}, - {"c_rsged", (DL_FUNC) &c_rsged, 6}, - {"c_dsged", (DL_FUNC) &c_dsged, 8}, - {"c_psged", (DL_FUNC) &c_psged, 7}, - {"c_qsged", (DL_FUNC) &c_qsged, 7}, - {"c_dhyp", (DL_FUNC) &c_dhyp, 8}, - {NULL, NULL, 0} -}; - -void attribute_visible R_init_tsdistributions(DllInfo *dll) { - R_registerRoutines(dll, CEntries, NULL, NULL, NULL); - R_useDynamicSymbols(dll, FALSE); -} diff --git a/tests/testthat.R b/tests/testthat.R old mode 100644 new mode 100755 diff --git a/tests/testthat/test-liktest.R b/tests/testthat/test-liktest.R old mode 100644 new mode 100755 diff --git a/tests/testthat/test-quantile.R b/tests/testthat/test-quantile.R old mode 100644 new mode 100755 diff --git a/vignettes/custom.css b/vignettes/custom.css old mode 100644 new mode 100755 diff --git a/vignettes/estimation_demo.Rmd b/vignettes/estimation_demo.Rmd old mode 100644 new mode 100755 diff --git a/vignettes/location_scale_distributions.Rmd b/vignettes/location_scale_distributions.Rmd old mode 100644 new mode 100755 diff --git a/vignettes/profile_demo.Rmd b/vignettes/profile_demo.Rmd old mode 100644 new mode 100755 diff --git a/vignettes/references.bib b/vignettes/references.bib old mode 100644 new mode 100755 diff --git a/vignettes/spd_demo.Rmd b/vignettes/spd_demo.Rmd old mode 100644 new mode 100755