This function calculates a k_nb x k_nb intermediate matrix of correlations for the Negative Binomial variables by
    extending the method of Yahav & Shmueli (2012, 10.1002/asmb.901). The intermediate correlation between Z1 and Z2 (the
    standard normal variables used to generate the Negative Binomial variables Y1 and Y2 via the inverse cdf method) is
    calculated using a logarithmic transformation of the target correlation.  First, the upper and lower Frechet-Hoeffding bounds
    (mincor, maxcor) on \(\rho_{y1,y2}\) are simulated.  Then the intermediate correlation is found as follows:
    $$\rho_{z1,z2} = (1/b) * log((\rho_{y1,y2} - c)/a)$$, where \(a = -(maxcor * mincor)/(maxcor + mincor)\),
    \(b = log((maxcor + a)/a)\), and \(c = -a\).  The function adapts code from Amatya & Demirtas' (2016) package
    PoisNor-package by:
1) allowing specifications for the number of random variates and the seed for reproducibility
2) providing the following checks: if \(\rho_{z1,z2}\) >= 1, \(\rho_{z1,z2}\) is set to 0.99; if \(\rho_{z1,z2}\) <= -1, \(\rho_{z1,z2}\) is set to -0.99
3) simulating Negative Binomial variables.
The function is used in findintercorr and rcorrvar.
    This function would not ordinarily be called by the user.
findintercorr_nb(rho_nb, size, prob, mu = NULL, nrand = 100000,
  seed = 1234)a k_nb x k_nb matrix of target correlations
a vector of size parameters for the Negative Binomial variables (see NegBinomial)
a vector of success probability parameters
a vector of mean parameters (*Note: either prob or mu should be supplied for all Negative Binomial variables,
not a mixture; default = NULL)
the number of random numbers to generate in calculating the bound (default = 10000)
the seed used in random number generation (default = 1234)
the k_nb x k_nb intermediate correlation matrix for the Negative Binomial variables
Please see references for findintercorr_pois.
PoisNor-package, findintercorr_pois,
    findintercorr_pois_nb,
    findintercorr, rcorrvar