twolodzko / extradistr

Additional univariate and multivariate distributions

extraDistr

Density, distribution function, quantile function and random generation for a number of univariate and multivariate distributions.

This package follows naming convention that is consistent with base R, where density (or probability mass) functions, distribution functions, quantile functions and random generation functions names are followed by d*, p*, q*, and r* prefixes.

Behaviour of the functions mimics the base R, where for invalid parameters NaN's are returned, while for values beyond function support 0's are returned (e.g. for non-integers in discrete distributions, or for negative values in functions with non-negative support).

All the functions vectorized and coded in C++11 using Rcpp.

extradistr's Issues

Citing?

How do you want us to cite extraDistr?

Feature Request: Implement Continuous Bernoulli and Continuous Categorical distributions

I'm a big fan of your package 👍. I was curious if you might consider adding two new exponential family distributions for percent or compositional data: the Continuous Bernoulli and Continuous Categorical which have been recently added into tensor flow. The paper for Continuous Categorical is here.

pbnbinom and pgpois are slow for some combinations of parameters

pbnbinom and pgpois (and potentially bbbinom) are slow for some combinations of parameters. This is caused by the fact that they are computed as F(x) = sum(f(1:x)) and it is inefficient (no matter of optimizations) for large x or for some combinations of parameters (e.g. for beta negative binomial small alpha and large beta).

Numerical instability in rwald

Great package, thanks for publishing. I think I have found problem in the rng_wald for very large lambda. I am getting zero and negative values from rwald, for example try rwald(n=100, lambda = 1, mu = 1e20).

Line 72 and 73 of wald-distribution.cpp should be using something like Horner's method for numerical stability. I think the following should be acceptable:

x = mu * (1 + (mu/(2.0*lambda)) * ( y - sqrt(4.0*lambda*y/mu+(y*y))));

versus the original

x = mu + (mu*mu*y)/(2.0*lambda) - mu/(2.0*lambda) * sqrt(4.0*mu*lambda*y+(mu*mu)*(y*y));

I can issue a pull request late next week if that helps. Might need to factor out lambda too.

zinb gives NAs when size is non-integer

rzinb samples NAs when the size parameter is non-integer.
In the documentation, it states that, how you expect from rnbinom, the size should be > 0 but need not be integer.

sz <- seq(0.5,5.5,0.5)
rzinb(n = sz, size = sz, prob = 0.5, pi = 0)
## [1] NA  0 NA  0 NA  3 NA  2 NA  6 NA

With the same paramegers, the rnbinom works fine

rnbinom(n = sz, size = sz, prob = 0.5)
## [1] 1 0 1 2 3 6 0 3 8 5 3

It seems that within rng_zinb it throws a warning if the size parameter is non-integer, and automatically returns NA_REAL, while the R::rnbinom function would otherwise work fine also with non-integer size parameter.

Issue when only 1 color is present

Running rmvhyper(nn = 1, n = 2, k = 1) returns 0 whereas I think it should return 1, since we're making 1 draw from an urn with 2 balls of a single color. Am I misunderstanding something?

Potential speedup in `pbbinom` exploiting the symmetry of the distribution

Hi,
due to the sequential nature of the algorithm in pbbinom (which is still wicked fast compared to competition!), it takes long time to evaluate with large q.

A simple improvement would be to take advantage of the symmetry of the distribution, i.e.:

extraDistr::pbbinom(q, size, alpha = alpha, beta = beta) ==
  1 - extraDistr::pbbinom(size - q - 1, size, alpha = beta, beta = alpha)

so we can choose whichever results in lower number of steps and always compute at most size / 2 steps.

A quick benchmark shows predictable speedup for large q:

N <- 100
q <- sample(550:950, size = N)
size <- 1000
alpha <- rexp(N)
beta <- rexp(N)

microbenchmark::microbenchmark(
  extraDistr::pbbinom(q, size, alpha = alpha, beta = beta),
  1 - extraDistr::pbbinom(size - q - 1, size, alpha = beta, beta = alpha),
  check = "equal"
)

# Unit: milliseconds
#                                                                     expr       min      lq      mean   median        uq       max neval
#                 extraDistr::pbbinom(q, size, alpha = alpha, beta = beta) 10.012301 10.4201 10.921302 10.70250 10.853801 18.183300   100
#  1 - extraDistr::pbbinom(size - q - 1, size, alpha = beta, beta = alpha)  4.763001  4.9289  5.197128  5.12875  5.244451  7.730002   100

Happy to provide a pull request if you agree this is worth implementing.

Thanks for all the hard work on the package.

Only NaNs produced

I'm using rmvhyper and my vector n have a length of 8614 and is consisted of one entry with 5000073631and the remaining are just one... When I sample 5e+06 I get
Warning message:
In cpp_rmvhyper(nn, n, k) : NaNs produced

rmnom returns NaNs for small probabilities

For very small values of prob parameter, rmnom returns NaNs, e.g.

p <- c(0, 0, 1, 0, 2.77053929981958e-18)
rmnom(1, 100, p)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0  100  NaN  NaN

Same with rdirmnom:

rdirmnom(1, 100, p+1e-5)
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0  100  NaN  NaN

Feature request: making extraDistr functions callable from C++

Hello,

Would it be possible to export the functions in the extraDistr library to make them callable from C++?

Really love this package!

Thanks,
Jon

Dependencies for unit tests not declared in description

Issue raised by CRAN in auto-email:

$extraDistr
'::' or ':::' imports not declared from:
‘Compositional’ ‘VGAM’ ‘evd’ ‘hoa’ ‘skellam’ ‘triangle’

doubt in ddnorm

how to draw the curve of ddnorm (discrete normal) in R

Rcpp implementations may not always be faster then pure-R

With simple functions like dlaplace the Rcpp implementation can be slighty slower then pure-R benchmarks:

library(rbenchmark)

dlaplaceR <- function(x, mu, sigma) {
  z <- (x-mu)/sigma
  1/(2*sigma) * exp(-abs(z))
}

dlaplaceRlog <- function(x, mu, sigma) {
  LOG_2F <- 0.6931471805599452862268
  z <- abs(x-mu)/sigma
  exp(-z - LOG_2F - log(sigma))
}

x <- seq(-100, 100, by = 0.001)

benchmark(
  dlaplace(x, 0, 1),
  dlaplaceR(x, 0, 1),
  dlaplaceRlog(x, 0, 1)
)
##                    test replications elapsed relative user.self sys.self user.child sys.child
## 1     dlaplace(x, 0, 1)          100    2.78    2.376      2.75     0.04         NA        NA
## 2    dlaplaceR(x, 0, 1)          100    1.17    1.000      1.14     0.03         NA        NA
## 3 dlaplaceRlog(x, 0, 1)          100    1.25    1.068      1.20     0.04         NA        NA

x <- -10:10

benchmark(
  dlaplace(x, 0, 1),
  dlaplaceR(x, 0, 1),
  dlaplaceRlog(x, 0, 1), 
  replications = 5000
)
##                    test replications elapsed relative user.self sys.self user.child sys.child
## 1     dlaplace(x, 0, 1)         5000    0.05        5      0.05        0         NA        NA
## 2    dlaplaceR(x, 0, 1)         5000    0.01        1      0.02        0         NA        NA
## 3 dlaplaceRlog(x, 0, 1)         5000    0.02        2      0.02        0         NA        NA

Some speedup could be achieved by parallelizing the functions using OpenMP, i.e. adding

#ifdef _OPENMP
#include <omp.h>
#endif

and

#ifdef _OPENMP
#pragma omp parallel for reduction(|| : throw_warning)
#endif

before each of the loops to parallelize (see this). This would influence only the outer for loops.

The problem is that while with large problems it can lead to x2 improvement in speed, with small examples it can lead even to x13 decreased speed as compared to pure-R benchmarks, so the overhead can be significant.

pbbinom(): vectorised evaluation fails for certain combinations of values

The following reproducible example shows that I can evaluate the cumulative probability for individual values, but putting the same values into vectors fail.

library(extraDistr)

test_dat <- 
  data.frame(
    q = c(1, 10, 10, 100),
    size = c(10, 100, 1000, 1000),
    alpha = 1, 
    beta = 1
  )

## Cumulative probabilites computed individually
test_dat$cumprobs <- 
  apply(test_dat, 1, function(x) pbbinom(x[1], x[2], x[3], x[4]))

## Cumulative probs computed with vectors

## First three rows: ok
with(
  head(test_dat, 3),
  identical(
    pbbinom(q, size, 1, 1),
    cumprobs
  )
)
#> [1] TRUE

## Last three rows: ok
with(
  tail(test_dat, 3),
  identical(
    pbbinom(q, size, 1, 1),
    cumprobs
  )
)
#> [1] TRUE

## But all four rows not ok
with(
  test_dat,
  pbbinom(q, size, 1, 1)
) # inadmissible values
#> Error in cpp_pbbinom(q, size, alpha, beta, lower.tail[1L], log.p[1L]): inadmissible values

## Not a matter of number of rows, but of combination of values:
with(
  test_dat[c(1, 4),],
  pbbinom(q, size, 1, 1)
) # inadmissible values
#> Error in cpp_pbbinom(q, size, alpha, beta, lower.tail[1L], log.p[1L]): inadmissible values

^{Created on 2018-10-22 by the reprex package (v0.2.1)}

The problem seems to arise when the range of values is large, which makes me think of some numerical issue. But I was using the function for evaluating MCMC samples, and the ranges of values is naturally large.

When possible, use more numerically stable implementations

Issues like #9 can happen also in other cases, more numerically stable implementations need to be used.

rgpd returns negative samples

From e-mail:

I just noticed that the command
extraDistr:rgpd(10)
return negative values.
That doesn't seem right.
Interestingly,
rgpd(10,xi=0.000001)
and
rgpd(10,xi=-0.000001)
both only return positive values.
So it's just at xi=0 that it seems broken.

Running pgev, pgpd, pbbinom functions with numeric(0) result in R crashing

Found this through using the fitdistrplus package when estimating a generalized pareto distribution, but realized it does the same effect on a number of p distribution functions (tested pgpd, pbbinom, pgev all fail, but seems okay with some like pnorm and pdgamma).

Basically, if I run pgpd(q = numeric(0), mu = 0, sigma= 1, xi = 0, lower.tail = TRUE, log.p = FALSE), R crashes. Similarly for other functions mentioned (not sure if more, didn't check all of your p functions).

Reason why this matters: I'm estimating a distribution using censored data via MLE using fitdistcens from the package fitidistrplus and it builds the likelihood function as if we have left/right/interval censored data (meaning, if I actually don't have all three, part of the function ends up calling the distribution function from extraDistr with empty numeric(0) quantile argument. Normally that's fine as, for example, ppois just returns numeric(0) which is exactly what one would expect (and not a crash).

Guessing the underlying C++ should have a length check or something? I already bothered the fitdistrplus authors to include checks if there's actually data for each type of censoring when building these function calls, but figured you'd like to know.

Cannot submit to CRAN

I cannot submit v 1.10.0 to CRAN because their checks throw:

Undocumented code objects:
'dbbinom' 'dbern' 'dbetapr' 'dbhatt' 'dbnbinom' 'dbvnorm' 'dbvpois'
'dcat' 'ddgamma' 'ddirichlet' 'ddirmnom' 'ddlaplace' 'ddnorm'
'ddunif' 'ddweibull' 'dfatigue' 'dfrechet' 'dgev' 'dgompertz' 'dgpd'
'dgpois' 'dgumbel' 'dhcauchy' 'dhnorm' 'dht' 'dhuber' 'dinvchisq'
'dinvgamma' 'dkumar' 'dlaplace' 'dlgser' 'dlomax' 'dlst' 'dmixnorm'
'dmixpois' 'dmnom' 'dmvhyper' 'dnhyper' 'dnsbeta' 'dpareto' 'dpower'
'dprop' 'drayleigh' 'dsgomp' 'dskellam' 'dslash' 'dtbinom' 'dtnorm'
'dtpois' 'dtriang' 'dwald' 'dzib' 'dzinb' 'dzip' 'pbbinom' 'pbern'
'pbetapr' 'pbhatt' 'pbnbinom' 'pcat' 'pdgamma' 'pdlaplace' 'pdnorm'
'pdunif' 'pdweibull' 'pfatigue' 'pfrechet' 'pgev' 'pgompertz' 'pgpd'
'pgpois' 'pgumbel' 'phcauchy' 'phnorm' 'pht' 'phuber' 'pinvchisq'
'pinvgamma' 'pkumar' 'plaplace' 'plgser' 'plomax' 'plst' 'pmixnorm'
'pmixpois' 'pnhyper' 'pnsbeta' 'ppareto' 'ppower' 'pprop' 'prayleigh'
'psgomp' 'pslash' 'ptbinom' 'ptnorm' 'ptpois' 'ptriang' 'pwald'
'pzib' 'pzinb' 'pzip' 'qbern' 'qbetapr' 'qcat' 'qdunif' 'qdweibull'
'qfatigue' 'qfrechet' 'qgev' 'qgompertz' 'qgpd' 'qgumbel' 'qhcauchy'
'qhnorm' 'qht' 'qhuber' 'qinvchisq' 'qinvgamma' 'qkumar' 'qlaplace'
'qlgser' 'qlomax' 'qlst' 'qnhyper' 'qnsbeta' 'qpareto' 'qpower'
'qprop' 'qrayleigh' 'qtbinom' 'qtlambda' 'qtnorm' 'qtpois' 'qtriang'
'qzib' 'qzinb' 'qzip' 'rbbinom' 'rbern' 'rbetapr' 'rbhatt' 'rbnbinom'
'rbvnorm' 'rbvpois' 'rcat' 'rcatlp' 'rdgamma' 'rdirichlet' 'rdirmnom'
'rdlaplace' 'rdnorm' 'rdunif' 'rdweibull' 'rfatigue' 'rfrechet'
'rgev' 'rgompertz' 'rgpd' 'rgpois' 'rgumbel' 'rhcauchy' 'rhnorm'
'rht' 'rhuber' 'rinvchisq' 'rinvgamma' 'rkumar' 'rlaplace' 'rlgser'
'rlomax' 'rlst' 'rmixnorm' 'rmixpois' 'rmnom' 'rmvhyper' 'rnhyper'
'rnsbeta' 'rpareto' 'rpower' 'rprop' 'rrayleigh' 'rsgomp' 'rsign'
'rskellam' 'rslash' 'rtbinom' 'rtlambda' 'rtnorm' 'rtpois' 'rtriang'
'rwald' 'rzib' 'rzinb' 'rzip'

Which doesn't seem to be true, but I have no idea how to "fix" it for now.

Use _PACKAGE for docs

From CRAN:

Dear maintainer,

You have file 'extraDistr/man/extraDistr.Rd' with \docType{package},
likely intended as a package overview help file, but without the
appropriate PKGNAME-package \alias as per "Documenting packages" in
R-exts.

This seems to be the consequence of the breaking change

Using @doctype package no longer automatically adds a -package alias.
Instead document _PACKAGE to get all the defaults for package
documentation.

in roxygen2 7.0.0 (2019-11-12) having gone unnoticed, see
r-lib/roxygen2#1491.

As explained in the issue, to get the desired PKGNAME-package \alias
back, you should either change to the new approach and document the new
special sentinel

"_PACKAGE"

or manually add

@Aliases extraDistr-package

if remaining with the old approach.

Please fix in your master sources as appropriate, and submit a fixed
version of your package within the next few months.

Best,
-k

dgpd and pgpd return wrong value with xi->0

dgpd and pgpd are not handling xi=0 special case correctly

dgpd(1)
[1] 0
dgpd(1, xi=0)
[1] 0
dgpd(1, xi=-1e-10)
[1] 0.3678794
dgpd(1, xi=1e-10)
[1] 0.3678794
pgpd(1)
[1] 0
pgpd(1, xi=0)
[1] 0
pgpd(1, xi=-1e-10)
[1] 0.6321206
pgpd(1, xi=1e-10)
[1] 0.6321206

Also cases with xi close to 0 return non-smooth results

dgpd(1, xi=1e-16)
[1] 1
dgpd(1, xi=1e-15)
[1] 0.3294855
dgpd(1, xi=1e-14)
[1] 0.3681736
dgpd(1, xi=1e-13)
[1] 0.3681736
dgpd(1, xi=1e-12)
[1] 0.3678467
dgpd(1, xi=1e-11)
[1] 0.3678794
dgpd(1, xi=1e-10)
[1] 0.3678794

The correct results with xi->0 can be obtained from dexp and pexp

dexp(1)
[1] 0.3678794
pexp(1)
[1] 0.6321206

Error while compiling in Linux

Found issues compiling while trying to use both the CRAN as well as the github latest version. I was wondering if you had any suggestions?

• installing source package ‘extraDistr’ ...
  ** package ‘extraDistr’ successfully unpacked and MD5 sums checked
  ** libs
  g++ -std=c++0x -I/usr/local/thirdparty/R/R-3.2.3/lib64/R/include -DNDEBUG -I/usr/local/include -I"/opt/home/u479871/R/x86_64-pc-linux-gnu-library/3.2/Rcpp/include" -fpic -g -O2 -c RcppExports.cpp -o RcppExports.o
  g++ -std=c++0x -I/usr/local/thirdparty/R/R-3.2.3/lib64/R/include -DNDEBUG -I/usr/local/include -I"/opt/home/u479871/R/x86_64-pc-linux-gnu-library/3.2/Rcpp/include" -fpic -g -O2 -c bernoulli-distribution.cpp -o bernoulli-distribution.o
  g++ -std=c++0x -I/usr/local/thirdparty/R/R-3.2.3/lib64/R/include -DNDEBUG -I/usr/local/include -I"/opt/home/u479871/R/x86_64-pc-linux-gnu-library/3.2/Rcpp/include" -fpic -g -O2 -c beta-binomial-distribution.cpp -o beta-binomial-distribution.o
  In file included from /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/functional:58,
  from /opt/home/u479871/R/x86_64-pc-linux-gnu-library/3.2/Rcpp/include/RcppCommon.h:61,
  from /opt/home/u479871/R/x86_64-pc-linux-gnu-library/3.2/Rcpp/include/Rcpp.h:27,
  from beta-binomial-distribution.cpp:1:
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple: In constructor ‘std::_Head_base<_Idx, _Head, false>::_Head_base(_UHead&&) [with _UHead = std::_Head_base<0ul, long int, false>, long unsigned int _Idx = 0ul, _Head = int]’:
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:179: instantiated from ‘std::_Tuple_impl<_Idx, _Head, _Tail ...>::_Tuple_impl(std::_Tuple_impl<_Idx, _UElements ...>&&) [with _UElements = long int, long int, long int, long unsigned int _Idx = 0ul, _Head = int, _Tail = int, int]’
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:256: instantiated from ‘std::tuple<_Elements>::tuple(std::tuple<_UElements ...>&&) [with _UElements = long int, long int, long int, _Elements = int, int, int]’
  beta-binomial-distribution.cpp:196: instantiated from here
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:94: error: cannot convert ‘std::_Head_base<0ul, long int, false>’ to ‘int’ in initialization
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple: In constructor ‘std::_Head_base<_Idx, _Head, false>::_Head_base(_UHead&&) [with _UHead = std::_Head_base<1ul, long int, false>, long unsigned int _Idx = 1ul, _Head = int]’:
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:179: instantiated from ‘std::_Tuple_impl<_Idx, _Head, _Tail ...>::_Tuple_impl(std::_Tuple_impl<_Idx, _UElements ...>&&) [with _UElements = long int, long int, long unsigned int _Idx = 1ul, _Head = int, _Tail = int]’
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:179: instantiated from ‘std::_Tuple_impl<_Idx, _Head, _Tail ...>::_Tuple_impl(std::_Tuple_impl<_Idx, _UElements ...>&&) [with _UElements = long int, long int, long int, long unsigned int _Idx = 0ul, _Head = int, _Tail = int, int]’
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:256: instantiated from ‘std::tuple<_Elements>::tuple(std::tuple<_UElements ...>&&) [with _UElements = long int, long int, long int, _Elements = int, int, int]’
  beta-binomial-distribution.cpp:196: instantiated from here
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:94: error: cannot convert ‘std::_Head_base<1ul, long int, false>’ to ‘int’ in initialization
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple: In constructor ‘std::_Head_base<_Idx, _Head, false>::_Head_base(_UHead&&) [with _UHead = std::_Head_base<2ul, long int, false>, long unsigned int _Idx = 2ul, _Head = int]’:
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:179: instantiated from ‘std::_Tuple_impl<_Idx, _Head, _Tail ...>::_Tuple_impl(std::_Tuple_impl<_Idx, _UElements ...>&&) [with _UElements = long int, long unsigned int _Idx = 2ul, _Head = int, _Tail = ]’
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:179: instantiated from ‘std::_Tuple_impl<_Idx, _Head, _Tail ...>::_Tuple_impl(std::_Tuple_impl<_Idx, _UElements ...>&&) [with _UElements = long int, long int, long unsigned int _Idx = 1ul, _Head = int, _Tail = int]’
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:179: instantiated from ‘std::_Tuple_impl<_Idx, _Head, _Tail ...>::_Tuple_impl(std::_Tuple_impl<_Idx, _UElements ...>&&) [with _UElements = long int, long int, long int, long unsigned int _Idx = 0ul, _Head = int, _Tail = int, int]’
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:256: instantiated from ‘std::tuple<_Elements>::tuple(std::tuple<_UElements ...>&&) [with _UElements = long int, long int, long int, _Elements = int, int, int]’
  beta-binomial-distribution.cpp:196: instantiated from here
  /usr/lib/gcc/x86_64-redhat-linux/4.4.7/../../../../include/c++/4.4.7/tuple:94: error: cannot convert ‘std::_Head_base<2ul, long int, false>’ to ‘int’ in initialization
  make: *** [beta-binomial-distribution.o] Error 1
  ERROR: compilation failed for package ‘extraDistr’
  removing ‘/opt/home/u479871/R/x86_64-pc-linux-gnu-library/3.2/extraDistr’

 sessionInfo()

R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8
[8] LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] Rcpp_0.12.9 httr_1.2.1 devtools_1.10.0
loaded via a namespace (and not attached):
[1] R6_2.2.0 tools_3.2.3 withr_1.0.1 curl_2.3 memoise_1.0.0 git2r_0.13.1 digest_0.6.9

CRAN warnings

e-mail from CRAN:

Dear maintainer,

Please see the problems shown on
https://cran.r-project.org/web/checks/check_results_extraDistr.html.

Please correct before 2023-12-12 to safely retain your package on CRAN.

consider adding the lkj correlation distribution

Hi, I'm a big fan of the package, maybe you want to consider adding the lkj correlation distribution

The code in R for dlkj and rlkj provided in https://www.rdocumentation.org/packages/rethinking/versions/1.59/topics/dlkjcorr is quite straightforward. But I don't know anything about C++.

rdirmnom produces NaNs

alpha <- c(1.480592e+00, 1.394943e-03, 4.529932e-06, 3.263573e+00, 4.554952e-06)
rdirmnom(10, 100, alpha)

More often than not, I'll get NaNs in the result.

ZIP probability density function definition inconsistent

The documentation states that the pdf is:

f(x) = [if x = 0:] (1-π)+π * exp(-λ) [else:] (1-π) * dpois(x, lambda)

Which is not self-inconsistent, because in the first case the exponential has a factor π and in the second case it has a factor (1-π). I assume it should have a factor (1-π) in both cases and it got swapped by accident with the π term in the first case. So I think the pdf should be:

f(x) = [if x = 0:] π + (1-π) * exp(-λ) [else:] (1-π) * dpois(x, lambda)

I think this is just wrong in the documentation and not in the code, as I am getting sensible numbers from the dnorm.

Remove implicit integer constraint on size parameter for beta-negative binomial

When trying to compute the PMF of beta-negative binomial with non-integer size ("r") parameter, a value of NaN is returned. For example:

dbnbinom(0,size=2.001,alpha=1.5,beta=2)

Looking at the source code, the NaN is produced programmatically (not due to numerical overflow or underflow) whenever r is not an integer. Since the r parameter can take on any positive real value, this restricts evaluation of the PMF function to a subset of the true parameter space.

I suggest either changing the documentation to reflect that the function only works for integer "r" (aka "size"), or remove the !isInteger(r, false) checking from the logpmf_bnbinom function.