sim_min_max <- function(Qx, Qy, nsim) {

ctrp2 <- function(y1, y2, p) {
ctrp3 <- function(y1, y2, y3, p1, p2) {

bdd_max_matrix <- function(px, py) {
bdd_min_matrix <- function(px, py) {

bdd_max_list <- function(px, py) {
bdd_min_list <- function(px, py) {

bdd_corr_list <- function(xvalues, yvalues, joint) {
bdd_corr_matrix <- function(xvalues, yvalues, joint) {

bdd_max <- function(xvalues, px, yvalues, py) {
bdd_min <- function(xvalues, px, yvalues, py) {

bdd_max_spearman <- function(px, py)  {
bdd_min_spearman <- function(px, py)  {

contin_multi_max_spearman <- function(p) {
contin_multi_min_spearman <- function(p) {

multi_unif_max <- function(xvals, p) {
multi_unif_min <- function(xvals, p) {

exp_exp_max <- function() {1}
exp_exp_min <- function() {1 - pi^2 / 6}
exp_exp_min_max <- function() {list(min = 1 - pi^2 / 6, max = 1)}

exp_norm_max <- function() { 0.903 }
exp_norm_min <- function() {-0.903 }
exp_norm_min_max <- function() {list(min = -0.903, max = 0.903)}

exp_unif_max <- function() { sqrt(3)/2 }
exp_unif_min <- function() {-sqrt(3)/2 }
exp_unif_min_max <- function() {list(min = -sqrt(3)/2, max = sqrt(3)/2)}

norm_unif_max <- function() { sqrt(3/pi) }
norm_unif_min <- function() {-sqrt(3/pi) }
norm_unif_min_max <- function() {list(min = -sqrt(3/pi), max = sqrt(3/pi))}

bern_norm_max <- function(p) {dnorm(qnorm(1-p))/sqrt(p*(1-p))}
bern_norm_min <- function(p) {- bern_norm_max(p)}
bern_norm_min_max <- function(p) {

bern_unif_max <- function(p) {sqrt(3*p*(1-p))}
bern_unif_min <- function(p) {- bern_unif_max(p)}
bern_unif_min_max <- function(p) {

exp_pois_max <- function(lambda) {
exp_pois_min <- function(lambda) {
exp_pois_min_max <- function(mu) {

gamma_multi_max <- function(alpha, pmfy, yval) {
gamma_multi_min <- function(alpha, pmfy, yval) {
gamma_multi_min_max <- function(alpha, pmfy, yval) {

gamma_pois_max <- function(alpha, lambda) {
gamma_pois_min <- function(alpha, lambda) {
gamma_pois_min_max <- function(alpha, lambda) {

norm_pois_max <- function(mu) {
norm_pois_min <- function(mu) {
norm_pois_min_max <- function(mu) {

bern_bern_max <- function(p1, p2) {
bern_bern_min <- function(p1, p2) {
bern_bern_min_max <- function(p1, p2) {

binom_binom_max <- function(m1, p1, m2, p2) {
binom_binom_min <- function(m1, p1, m2, p2) {
binom_binom_min_max <- function(m1, p1, m2, p2) {

binom_nbinom_max <- function(m1, p1, m2, p2) {
binom_nbinom_min <- function(m1, p1, m2, p2) {
binom_nbinom_min_max <- function(m1, p1, m2, p2) {

bern_nbinom_max <- function(p, m, a) {
bern_nbinom_min <- function(p, m, a) {
bern_nbinom_min_max <- function(p, m, a) {

bern_pois_max <- function(p, mu) {
bern_pois_min <- function(p, mu) {
bern_pois_min_max <- function(p, mu) {

binom_pois_max <- function(m, p, mu) {
binom_pois_min <- function(m, p, mu) {
binom_pois_min_max <- function(m, p, mu) {

nbinom_nbinom_max <- function(m1, p1, m2, p2) {
nbinom_nbinom_min <- function(m1, p1, m2, p2) {
nbinom_nbinom_min_max <- function(m1, p1, m2, p2) {

nbinom_pois_max <- function(m, p, mu) {
nbinom_pois_min <- function(m, p, mu) {
nbinom_pois_min_max <- function(m, p, mu) {

pois_pois_max <- function(mu_x, mu_y) {
pois_pois_min <- function(mu_x, mu_y) {
pois_pois_min_max <- function(mu_x, mu_y) {

dpois_pois_max_matrix <- function(mu_x, mu_y) {
dpois_pois_min_matrix <- function(mu_x, mu_y) {

dpois_pois_max_list   <- function(mu_x, mu_y) {
dpois_pois_min_list   <- function(mu_x, mu_y) {