# corr01lib.R : Correlation bounds for various distributions
#               Tools for computation and simulation
#               Version 1.0
#               (c) Bahjat F. Qaqish (2017)
#
# This code is free for personal research use. 
# Please cite the following sources upon which the code is based:
#
#   Leonov, S. & Qaqish, B. (2017).
#   Correlated endpoints: simulation, modeling, and extreme correlations.
#   Statistical Papers. https://doi.org/10.1007/s00362-017-0960-2
#
#   Preisser, J. S. & Qaqish, B. F. (2014).
#   A comparison of methods for simulating correlated binary
#   variables with specified marginal means and correlations.
#   Journal of Statistical Computation and Simulation,
#   volume 84, issue 11,  pages 2441-2452.
#
# For commercial use, please contact the author.
#
# Function collection.
# Function names are a_b_max or a_b_min or a_b_min_max
# where a and b are distribution names.
# a_b_min returns the minimum (largest negative)
# a_b_max returns the maximum
# a_b_min_max returns both in a list(min=, max=)
# The R names (norm, pois, binom, etc) are used when applicable.
# a and b are sorted alphabetically, so there is norm_pois_max()
# but no pois_norm_max()
# All parameters are scalars except where noted.
# bdd_ is used for a generic bivariate discrete distribution
#
#**************************************************************
# General tools and building blocks 
#**************************************************************

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

  # min amd max correlation by simulation
  # Qx and Qy are two quantile functions
  # nsim = number of (X,Y) points to simulate
  # Can be used on continuous, discrete and mixed types

  u = runif(nsim)
  x = Qx(u)
  y = Qy(u)
  rmax = cor(x,y)
  y = Qy(1-u)
  rmin = cor(x,y)
  return(list(min=rmin, max=rmax))
}

#**************************************************************
ctrp2 <- function(y1, y2, p) {
  # random picking out of two
  # y1 w.p. p
  # y2 w.p. 1-p

  if (p<0 | p>1) {return(NA)}
  if (length(y1) != length(y2))  {return(NA)}
  u = (runif(length(y1)) <= p)
  return(u*y1 + (1-u)*y2)
}

ctrp3 <- function(y1, y2, y3, p1, p2) {
  # random picking out of three
  # y1 w.p. p1
  # y2 w.p. p2
  # y3 w.p. 1-(p1+p2)
  if (p1<0 | p1>1) {return(NA)}
  if (p2<0 | p2>1) {return(NA)}
  if (p1+p2 > 1)   {return(NA)}
  if (length(y1) != length(y2))  {return(NA)}
  if (length(y1) != length(y3))  {return(NA)}
  u = runif(length(y1))
  u = (u > p1+p2) + (u > p1) + 1
  z = (u==1) * y1 + (u==2) * y2 + (u==3) * y3
  return(z)
}

#**************************************************************
# Bivariate discrete distributions (bdd).
# i.e. two correlated multinomials
#**************************************************************

bdd_max_matrix <- function(px, py) {
  # Both  arguments are vectors
  # px = pmf(X), py=pmf(Y)
  # Returns the joint pmf(X,Y) in a matrix pxy such
  # that corr(X,Y) is maximized.
  # pxy[i,j] = P(X=i, Y=j), will be mostly zeros,
  # See below for a "sparse matrix" version
  # Algorithm: Walk from the upper left corner to 
  # the lower right corner,
  # moving right or down or right-and-down, 
  # tracing a staircase path through the table.
  m = length(px)
  n = length(py)
  i = 1; j = 1;
  a = px[1]; b=py[1]
  px = c(px, 0) # we'll access past the last element
  py = c(py, 0) # same
  pxy = matrix(0, m, n)
  while (i <= m & j <= n) {
    if (a < b) { pxy[i,j] = a; b = b - a; i=i+1; a=px[i]; }
    else 
    if (b < a) { pxy[i,j] = b; a = a - b; j=j+1; b=py[j]; }
    else # a == b
               { pxy[i,j] = a; i=i+1; j=j+1; a=px[i]; b=py[j];}
  }  
  return(pxy)
}

bdd_min_matrix <- function(px, py) {
  # Both  arguments are vectors
  # px = pmf(X), py=pmf(Y)
  # Returns the joint pmf(X,Y) in a matrix pxy such
  # that corr(X,Y) is minimized.
  A   = bdd_max_matrix(rev(px), py)
  pxy = apply(A, 2, rev)
  return(pxy)
}
#**************************************************************

bdd_max_list <- function(px, py) {
  # Both  arguments are vectors
  # px = pmf(X), py=pmf(Y), vectors
  # Returns a list of:
  # the joint pmf(X,Y) in a 1-dim array pxy, 
  # indices in i and j, such
  # that corr(X,Y) is maximized.
  # pxy[k] = P(X=i[k], Y=j[k])
  # Algorithm: Walk from the upper left corner to 
  # the lower right corner,
  # moving right or down or right-and-down, 
  # tracing a staircase path through the table.
  # Array pxy will contain only the non-zero cells
  m = length(px)
  n = length(py)
  i = 1; j = 1;
  a = px[1]; b=py[1]
  px = c(px, 0) # we'll access past the last element
  py = c(py, 0) # same
  pxy = double (m + n - 1)
  ii  = integer(m + n - 1)
  jj  = integer(m + n - 1)
  k = 0
  while (i <= m & j <= n) {
    k = k + 1
    ii[k]=i; jj[k]=j;
    if (a < b) { pxy[k] = a; b = b - a; i=i+1; a=px[i]; }
    else
    if (b < a) { pxy[k] = b; a = a - b; j=j+1; b=py[j]; }
    else # b == a
               { pxy[k] = a; i=i+1; j=j+1; a=px[i]; b=py[j];}
  } 
  if (k < m + n - 1) { # happens if we made any diagonal moves
    pxy = pxy[1:k]
    ii  = ii[1:k]
    jj  = jj[1:k]
  }     
  return(list(pxy=pxy, i=ii, j=jj))
}
#**************************************************************

bdd_min_list <- function(px, py) {
  # Both  arguments are vectors
  # px = pmf(X), py=pmf(Y), vectors
  # Returns a list of:
  # the joint pmf(X,Y) in a 1-dim array pxy, 
  # indices in i and j, such
  # that corr(X,Y) is minimized.
  # pxy[k] = P(X=i[k], Y=j[k])
  ml = bdd_max_list(rev(px), py)  # maximize with px reversed
  ml$i = length(ml$i) - ml$i      # flip the x values
  return(ml)
}
#**************************************************************

bdd_corr_list <- function(xvalues, yvalues, joint) {
  # Correlation in a bivariate discrete distribution
  # xvalues, yvalues: vectors
  # joint: list(i, j, pxy)
  EX = xvalues[joint$i] %*% joint$pxy 
  EY = yvalues[joint$j] %*% joint$pxy
  dx = xvalues[joint$i] - as.vector(EX)
  dy = yvalues[joint$j] - as.vector(EY)
  VXX = dx^2 %*% joint$pxy 
  VYY = dy^2 %*% joint$pxy 
  VXY = (dx * dy) %*% joint$pxy 
  r   = VXY / sqrt(VXX * VYY)
  return(r)
}
#**************************************************************

bdd_corr_matrix <- function(xvalues, yvalues, joint) {
  # Correlation in a bivariate discrete distribution
  # xvalues, yvalues: vectors
  # joint: matrix form of pxy
  # *** under construction ***
  px = apply(joint, 1, sum)
  py = apply(joint, 2, sum)
  EX = xvalues %*% px 
  EY = yvalues %*% py
  dx = xvalues - c(EX)
  dy = yvalues - c(EY)
  VXX = dx^2 %*% px
  VYY = dy^2 %*% py 
  VXY = sum(outer(dx, dy) * joint)
  r   = VXY / sqrt(VXX * VYY)
  return(r)
}
#**************************************************************

bdd_max <- function(xvalues, px, yvalues, py) {
  # All 4 arguments are vectors of the same length
  # xvalues and yvalues must be sorted in ascending order
  # px and py are the corresponding pmfs, each sums to 1
  joint = bdd_max_list(px, py)
  r = bdd_corr_list(xvalues, yvalues, joint)
  return(r)
}

bdd_min <- function(xvalues, px, yvalues, py) {
  # All 4 arguments are vectors of the same length
  # xvalues and yvalues must be sorted in ascending order
  # px and py are the corresponding pmfs, each sums to 1
  joint = bdd_max_list(rev(px), py)
  r = bdd_corr_list(rev(xvalues), yvalues, joint)
  return(r)
}
#**************************************************************

bdd_max_spearman <- function(px, py)  {

  # Max Spearman correlation in the bivariate discrete
  # px, py : marginal pmf's

  return(bdd_max(cumsum(px), px, cumsum(py), py))

}
  
bdd_min_spearman <- function(px, py)  {

  # Min Spearman correlation in the bivariate discrete
  # px, py : marginal pmf's

  return(bdd_min(cumsum(px), px, cumsum(py), py))

}
#**************************************************************
# Bivariate continuous-and-multinomial distributions
#**************************************************************

contin_multi_max_spearman <- function(p) {
  # p[] is a pmf, sum=1  
  # returns the max Spearmen corr(Continuous, X),
  # X ~ discrete with pmf p (ordered by the values X takes)

  multi_unif_max (cumsum(p))
}

contin_multi_min_spearman <- function(p) {
  # p[] is a pmf, sum=1  
  # returns the max Spearmen corr(Continuous, X),
  # X ~ discrete with pmf p (ordered by the values X takes)

  multi_unif_min (cumsum(p))
}
#**************************************************************

multi_unif_max <- function(xvals, p) {
  # p[]: the pmf of X, sum=1  
  # xvals[]: X values
  # returns the max corr(U, X)
  # U ~ uniform(0,1)
  # X ~ discrete with pmf p 
  
  EX   = p %*% xvals
  varX = p %*% (xvals - as.vector(EX))^2
  EU = 0.5 # uniform
  varU = 1/12  # uniform
  G = cumsum(p)  # cdf of X
  left = c(0, G[-length(G)])
  mid = (left + G) / 2
  EUX = p %*% (mid * xvals)
  covUX = EUX - EU * EX
  corUX = covUX / sqrt(varU * varX)
  return(corUX)
}

multi_unif_min <- function(xvals, p) {
  - multi_unif_max (xvals, p)
}
#**************************************************************
# Specific bivariate continuous distributions
#**************************************************************
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))}
#**************************************************************
# Specific bivariate mixed type distributions
#**************************************************************
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) {
  a = bern_norm_max(p)
  list(min = -a, max = a)
}
#**************************************************************
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) {
  a = bern_unif_max(p)
  list(min = -a, max = a)
}
#**************************************************************

exp_pois_max <- function(lambda) {
  # X ~ gamma(1, 1) = exp(1)
  # Y ~ Poisson(lambda)
  # returns max corr(X,Y)
  eps = 1e-6
  m = qpois(eps, lambda, lower.tail=FALSE)
  y = 0:m
  G = ppois(y, lambda)
  a = -log(1 - G)
  s = (1 + a) * (1 - G)    # S() = 1 - F(), gamma(2, 1)
  s[m+1] = 0    # G[m+1] = 1
  h = sum(s)    # mean = integral of S()
  corrxy = (h - lambda) / sqrt(lambda)
  return(corrxy)
}
#**************************************************************

exp_pois_min <- function(lambda) {
  # X ~ gamma(1, 1) = exp(1)
  # Y ~ Poisson(lambda) 
  # returns min corr(X,Y)
  eps = 1e-6
  m = qpois(eps, lambda, lower.tail=FALSE)
  y = 0:m   
  G = ppois(y, lambda, lower.tail=FALSE)
  a = -log(1 - G)
  s = 1 - (1 + a) * (1 - G)  # S() = 1 - F(), gamma(2, 1)
  s[m+1] = 0    # G[m+1] = 0
  h = sum(s)
  corrxy = (h - lambda) / sqrt(lambda)
  return(corrxy)
}
#**************************************************************
exp_pois_min_max <- function(mu) {
  list(min = exp_pois_min(mu),  
       max = exp_pois_max(mu))
}
#**************************************************************

gamma_multi_max <- function(alpha, pmfy, yval) {
  # X ~ gamma(alpha, 1)
  # Y ~ multinomial (pmfy, yvalues)
  # yval[] must be in ascending order
  # returns max corr(X,Y)
  G = cumsum(pmfy)
  a = qgamma(G, alpha)
  F = pgamma(a, alpha + 1)
  F[length(F)] = 1  
  p = c(F[1], diff(F))
  EYs = yval %*% p     # E[Y*], E[XY]=E[X] E[Y*]
  EY  = yval %*% pmfy  # E[Y]
  delta = EYs - EY     # cov(X,Y) = E[X] delta
  VY  = (yval - c(EY))^2 %*% pmfy   # var(Y)
  corrxy = sqrt(alpha / VY) * delta
  return(corrxy)
}
#**************************************************************

gamma_multi_min <- function(alpha, pmfy, yval) {
  # X ~ gamma(alpha, 1)
  # Y ~ multinomial (pmfy, yval)
  # yval[] must be in ascending order
  # returns min corr(X,Y)
  return(gamma_multi_max(alpha, rev(pmfy), rev(yval)))
}
#**************************************************************
gamma_multi_min_max <- function(alpha, pmfy, yval) {
  # X ~ gamma(alpha, 1)
  # Y ~ multinomial (pmfy, yval)
  # yval[] must be in ascending order
  # returns min and max corr(X,Y)
  list(min = gamma_multi_min(alpha, pmfy, yval),
       max = gamma_multi_max(alpha, pmfy, yval))
}
#**************************************************************

gamma_pois_max <- function(alpha, lambda) {
  # X ~ gamma(alpha, 1)
  # Y ~ Poisson(lambda)
  # returns max corr(X,Y)
  eps = 1e-6
  m = qpois(eps, lambda, lower.tail=FALSE)
  y = 0:m
  G = ppois(y, lambda)
  a = qgamma(G, alpha)
  s = pgamma(a, alpha + 1, lower.tail=FALSE)  # S() = 1 - F()
  s[m+1] = 0    # G[m+1] = 1
  h = sum(s)    # mean = integral of S()
  corrxy = sqrt(alpha / lambda) * (h - lambda)
  return(corrxy)
}
#**************************************************************

gamma_pois_min <- function(alpha, lambda) {
  # X ~ gamma(alpha, 1)
  # Y ~ Poisson(lambda)
  # returns min corr(X,Y)
  eps = 1e-6
  m = qpois(eps, lambda, lower.tail=FALSE)
  y = 0:m   
  G = ppois(y, lambda, lower.tail=FALSE)
  a = qgamma(G, alpha)
  s = pgamma(a, alpha + 1)
  s[m+1] = 0    # G[m+1] = 0
  h = sum(s)
  corrxy = sqrt(alpha / lambda) * (h - lambda)
  return(corrxy)
}
#**************************************************************
gamma_pois_min_max <- function(alpha, lambda) {
  list(min = gamma_pois_min(alpha, lambda),
       max = gamma_pois_max(alpha, lambda))
}
#**************************************************************

norm_pois_max <- function(mu) {
  # X ~ normal(0,1)
  # Y ~ Poisson(lambda)
  # returns max corr(X,Y)
  eps = 1e-6
  m = qpois(eps, mu, lower.tail=FALSE)
  y = 0:m
  G = ppois(y, mu)
  a = qnorm(G)
  sum(dnorm(a)) / sqrt(mu)
}
#**************************************************************

norm_pois_min <- function(mu) {
  # X ~ normal(0,1)
  # Y ~ Poisson(mu)
  # returns min corr(X,Y)
  - norm_pois_max(mu)
}
#**************************************************************

norm_pois_min_max <- function(mu) {
  # X ~ normal(0,1)
  # Y ~ Poisson(mu)
  # returns min corr(X,Y)
  r = norm_pois_max(mu)
  return(list(min = -r, max = r))
}
#**************************************************************
# Specific bivariate discrete distributions.
#**************************************************************  

bern_bern_max <- function(p1, p2) {
  # X ~ Bern(p1), Y ~ Bern(p2)
  # return max corr(X,Y)
  a = sqrt( (p1*(1-p2)) / ((1-p1)*p2) )
  min(a, 1 / a)
}

bern_bern_min <- function(p1, p2) {
  # X ~ Bern(p1), Y ~ Bern(p2)
  # return min corr(X,Y)
  a = sqrt( (p1*p2) / ((1-p1)*(1-p2)) )
  - min(a, 1 / a)
}

bern_bern_min_max <- function(p1, p2) {
  # X ~ Bern(p1), Y ~ Bern(p2)
  # return min and max corr(X,Y)
  a = sqrt( p1 / (1 - p1) )
  b = sqrt( p2 / (1 - p2) )
  ab = a * b
  return(list(min = -min(ab, 1 / ab),
              max =  min(a / b, b / a)))
}

#**************************************************************

binom_binom_max <- function(m1, p1, m2, p2) {
  # X ~ Binomial(m1,p1), Y ~ Binomial(m2,p2)
  # return max corr(X,Y)
  eps = 1e-6
  x = 0:m1
  y = 0:m2
  px=dbinom(x, m1, p1)
  py=dbinom(y, m2, p2)
  return(bdd_max (x, px, y, py))
}

binom_binom_min <- function(m1, p1, m2, p2) {
  # X ~ Binomial(m1,p1), Y ~ Binomial(m2,p2)
  # return min corr(X,Y)
  eps = 1e-6
  x = 0:m1
  y = 0:m2
  px=dbinom(x, m1, p1)
  py=dbinom(y, m2, p2)
  return(bdd_min (x, px, y, py))
}

binom_binom_min_max <- function(m1, p1, m2, p2) {
  # X ~ Binomial(m1,p1), Y ~ Binomial(m2,p2)
  # return min and max corr(X,Y)
  eps = 1e-6
  x = 0:m1
  y = 0:m2
  px=dbinom(x, m1, p1)
  py=dbinom(y, m2, p2)
  return(list(min= bdd_min (x, px, y, py),
              max= bdd_max (x, px, y, py) ))
}

#**************************************************************

binom_nbinom_max <- function(m1, p1, m2, p2) {
  # X ~ Binomial(m1, p1), Y ~ NegBin(m2, p2)
  # return max corr(X,Y)
  if (m1 == 1) return(bern_nbinom_max(p1, m2, p2))
  eps = 1e-6
  x = 0: m1
  y = 0: qnbinom(eps, m2, p2, lower.tail=FALSE)
  px=dbinom (x, m1, p1)
  py=dnbinom(y, m2, p2)
  py = py / sum(py)
  return(bdd_max (x, px, y, py))
}

binom_nbinom_min <- function(m1, p1, m2, p2) {
  # X ~ Binomial(m1, p1), Y ~ NegBin(m2, p2)
  # return min corr(X,Y)
  if (m1 == 1) return(bern_nbinom_min(p1, m2, p2))
  eps = 1e-6
  x = 0: m1
  y = 0: qnbinom(eps, m2, p2, lower.tail=FALSE)
  px=dbinom (x, m1, p1)
  py=dnbinom(y, m2, p2)
  py = py / sum(py)
  return(bdd_min (x, px, y, py))
}

binom_nbinom_min_max <- function(m1, p1, m2, p2) {
  # X ~ Binomial(m1, p1), Y ~ NegBin(m2, p2)
  # return list(min, max) corr(X,Y)
  # Use if you need both min and max
  if (m1 == 1) return(bern_nbinom_min_max(p1, m2, p2))
  eps = 1e-6
  x = 0: m1
  y = 0: qnbinom(eps, m2, p2, lower.tail=FALSE)
  px=dbinom (x, m1, p1)
  py=dnbinom(y, m2, p2)
  py = py / sum(py)
  return(list(min= bdd_min (x, px, y, py),
              max= bdd_max (x, px, y, py) ))
}

#**************************************************************
bern_nbinom_max <- function(p, m, a) {
  # X ~ Bern(p), Y ~ NegBin(m, a)
  # return max corr(X,Y)
  q = 1 - p
  k = qnbinom(q, m, a)
  delta = pnbinom(k, m, a) - q  # spell over to X=1
  alpha = dnbinom(0:k, m, a)    # k can be 0
  alpha[1+k] = alpha[1+k] - delta
  ey = m * (1 - a) / a
  exy = ey - (0:k) %*% alpha  # E[XY] = E[Y] - E[(1-X)Y]
  cxy = exy - p * ey
  vy = ey / a
  rxy = cxy / sqrt(p*q*vy)
  return(rxy)
}

bern_nbinom_min <- function(p, m, a) {
  # X ~ Bern(p), Y ~ NegBin(m, a)
  # return min corr(X,Y)
  return(- bern_nbinom_max(1-p, m, a))
}

bern_nbinom_min_max <- function(p, m, a) {
  # X ~ Bern(p), Y ~ NegBin(m, a)
  # return min and max corr(X,Y)
  return(list(min = -bern_nbinom_max(1-p, m, a),
              max =  bern_nbinom_max(p, m, a) ))
}
#**************************************************************
bern_pois_max <- function(p, mu) {
  # X ~ Bern(p), Y ~ Poisson(mu)
  # return max corr(X,Y)
  q = 1 - p
  k = qpois(q, mu)
  delta = ppois(k, mu) - q  # spell over to X=1
  alpha = dpois(0:k, mu)    # k can be 0
  alpha[1+k] = alpha[1+k] - delta
  cxy = mu * q - (0:k) %*% alpha
  rxy = cxy / sqrt(mu*p*q)
  return(rxy)
}

bern_pois_min <- function(p, mu) {
  # X ~ Bern(p), Y ~ Poisson(mu)
  # return min corr(X,Y)
  return(- bern_pois_max(1-p, mu))
}

bern_pois_min_max <- function(p, mu) {
  # X ~ Bern(p), Y ~ Poisson(mu)
  # return min and max corr(X,Y)
  return(list(min = -bern_pois_max(1-p, mu),
              max =  bern_pois_max(p, mu) ))
}
#**************************************************************

binom_pois_max <- function(m, p, mu) {
  # X ~ Binomial(m,p), Y ~ Poisson(mu)
  # return max corr(X,Y)
  if (m == 1) return(bern_pois_max(p, mu))
  eps = 1e-6
  x = 0: m
  y = 0: qpois(eps, mu, lower.tail=FALSE)
  px=dbinom(x, m, p)
  py=dpois(y, mu)
  py=py/sum(py)
  return(bdd_max (x, px, y, py))
}

binom_pois_min <- function(m, p, mu) {
  # X ~ Binomial(m,p), Y ~ Poisson(mu)
  # return min corr(X,Y)
  if (m == 1) return(bern_pois_min(p, mu))
  eps = 1e-6
  x = 0: m
  y = 0: qpois(eps, mu, lower.tail=FALSE)
  px=dbinom(x, m, p)
  py=dpois(y, mu)
  py=py/sum(py)
  return(bdd_min (x, px, y, py))
}

binom_pois_min_max <- function(m, p, mu) {
  # X ~ Binomial(m,p), Y ~ Poisson(mu)
  # return list(min, max) corr(X,Y)
  # Use if you need both min and max
# if (m == 1) return(bern_pois_min_max(p, mu))
  eps = 1e-6
  x = 0: m
  y = 0: qpois(eps, mu, lower.tail=FALSE)
  px=dbinom(x, m, p)
  py=dpois(y, mu)
  py=py/sum(py)
  return(list(min= bdd_min (x, px, y, py),
              max= bdd_max (x, px, y, py) ))
}
#**************************************************************

nbinom_nbinom_max <- function(m1, p1, m2, p2) {
  # X ~ NegBin(m1, p1), Y ~ NegBin(m2, p2)
  # return max corr(X,Y)
  eps = 1e-6
  x = 0: qnbinom(eps, m1, p1, lower.tail=FALSE)
  y = 0: qnbinom(eps, m2, p2, lower.tail=FALSE)
  px=dnbinom(x, m1, p1)
  py=dnbinom(y, m2, p2)
  px = px / sum(px)
  py = py / sum(py)
  return(bdd_max (x, px, y, py))
}

nbinom_nbinom_min <- function(m1, p1, m2, p2) {
  # X ~ NegBin(m1, p1), Y ~ NegBin(m2, p2)
  # return min corr(X,Y)
  eps = 1e-6
  x = 0: qnbinom(eps, m1, p1, lower.tail=FALSE)
  y = 0: qnbinom(eps, m2, p2, lower.tail=FALSE)
  px=dnbinom(x, m1, p1)
  py=dnbinom(y, m2, p2)
  px = px / sum(px)
  py = py / sum(py)
  return(bdd_min (x, px, y, py))
}

nbinom_nbinom_min_max <- function(m1, p1, m2, p2) {
  # X ~ NegBin(m1, p1), Y ~ NegBin(m2, p2)
  # return list(min, max) corr(X,Y)
  # Use if you need both min and max
  eps = 1e-6
  x = 0: qnbinom(eps, m1, p1, lower.tail=FALSE)
  y = 0: qnbinom(eps, m2, p2, lower.tail=FALSE)
  px=dnbinom(x, m1, p1)
  py=dnbinom(y, m2, p2)
  px = px / sum(px)
  py = py / sum(py)
  return(list(min= bdd_min (x, px, y, py),
              max= bdd_max (x, px, y, py) ))
}

#**************************************************************

nbinom_pois_max <- function(m, p, mu) {
  # X ~ Binomial(m,p), Y ~ Poisson(mu)
  # return max corr(X,Y)
  eps = 1e-6
  x = 0: qnbinom(eps, m, p, lower.tail=FALSE)
  y = 0: qpois(eps, mu, lower.tail=FALSE)
  px=dnbinom(x, m, p)
  py=dpois(y, mu)
  px=px/sum(px)
  py=py/sum(py)
  return(bdd_max (x, px, y, py))
}

nbinom_pois_min <- function(m, p, mu) {
  # X ~ NegativeBinomial(m,p), Y ~ Poisson(mu)
  # return min corr(X,Y)
  eps = 1e-6
  x = 0: qnbinom(eps, m, p, lower.tail=FALSE)
  y = 0: qpois(eps, mu, lower.tail=FALSE)
  px=dnbinom(x, m, p)
  py=dpois(y, mu)
  px=px/sum(px)
  py=py/sum(py)
  return(bdd_min (x, px, y, py))
}

nbinom_pois_min_max <- function(m, p, mu) {
  # X ~ Binomial(m,p), Y ~ Poisson(mu)
  # return list(min, max) corr(X,Y)
  # Use if you need both min and max
  eps = 1e-6
  x = 0: qnbinom(eps, m, p, lower.tail=FALSE)
  y = 0: qpois(eps, mu, lower.tail=FALSE)
  px=dnbinom(x, m, p)
  py=dpois(y, mu)
  px=px/sum(px)
  py=py/sum(py)
  return(list(min= bdd_min (x, px, y, py),
              max= bdd_max (x, px, y, py) ))
}

#**************************************************************

pois_pois_max <- function(mu_x, mu_y) {
  # X ~ Poisson(mu_x), Y ~ Poisson(mu_y)
  # return max corr(X,Y)
  eps = 1e-6
  x = 0: qpois(eps, mu_x, lower.tail=FALSE)
  y = 0: qpois(eps, mu_y, lower.tail=FALSE)
  px=dpois(x, mu_x)
  py=dpois(y, mu_y)
  px=px/sum(px)
  py=py/sum(py)
  return(bdd_max (x, px, y, py))
}
#**************************************************************

pois_pois_min <- function(mu_x, mu_y) {
  # X ~ Poisson(mu_x), Y ~ Poisson(mu_y)
  # return min corr(X,Y)
  eps = 1e-6
  x = 0: qpois(eps, mu_x, lower.tail=FALSE)
  y = 0: qpois(eps, mu_y, lower.tail=FALSE)
  px=dpois(x, mu_x)
  py=dpois(y, mu_y)
  px=px/sum(px)
  py=py/sum(py)
  return(bdd_min (x, px, y, py))
}
#**************************************************************

pois_pois_min_max <- function(mu_x, mu_y) {
  # X ~ Poisson(mu_x), Y ~ Poisson(mu_y)
  # return list(min, max) corr(X,Y)
  # Use if you need both min and max
  eps = 1e-6
  x = 0: qpois(eps, mu_x, lower.tail=FALSE)
  y = 0: qpois(eps, mu_y, lower.tail=FALSE)
  px=dpois(x, mu_x)
  py=dpois(y, mu_y)
  px=px/sum(px)
  py=py/sum(py)
  return(list(min= bdd_min (x, px, y, py),
              max= bdd_max (x, px, y, py) ))
}

#**************************************************************
# pmf's with extremal correlations
#**************************************************************

dpois_pois_max_matrix <- function(mu_x, mu_y) {
  # X ~ Poisson(mu_x), Y ~ Poisson(mu_y)
  # return distribution with max corr(X,Y)
  eps = 1e-6
  x = 0: qpois(eps, mu_x, lower.tail=FALSE)
  y = 0: qpois(eps, mu_y, lower.tail=FALSE)
  px=dpois(x, mu_x)
  py=dpois(y, mu_y)
  px=px/sum(px)
  py=py/sum(py)
  return(bdd_max_matrix (px, py))
}

#**************************************************************

dpois_pois_min_matrix <- function(mu_x, mu_y) {
  # X ~ Poisson(mu_x), Y ~ Poisson(mu_y)
  # return distribution with min corr(X,Y)
  eps = 1e-6
  x = 0: qpois(eps, mu_x, lower.tail=FALSE)
  y = 0: qpois(eps, mu_y, lower.tail=FALSE)
  px=dpois(x, mu_x)
  py=dpois(y, mu_y)
  px=px/sum(px)
  py=py/sum(py)
  return(bdd_min_matrix (px, py))
}

#**************************************************************

dpois_pois_max_list   <- function(mu_x, mu_y) {
  # X ~ Poisson(mu_x), Y ~ Poisson(mu_y)
  # return distribution with max corr(X,Y)
  eps = 1e-6
  x = 0: qpois(eps, mu_x, lower.tail=FALSE)
  y = 0: qpois(eps, mu_y, lower.tail=FALSE)
  px=dpois(x, mu_x)
  py=dpois(y, mu_y)
  px=px/sum(px)
  py=py/sum(py)
  return(bdd_max_list   (px, py))
}

#**************************************************************

dpois_pois_min_list   <- function(mu_x, mu_y) {
  # X ~ Poisson(mu_x), Y ~ Poisson(mu_y)
  # return distribution with min corr(X,Y)
  eps = 1e-6
  x = 0: qpois(eps, mu_x, lower.tail=FALSE)
  y = 0: qpois(eps, mu_y, lower.tail=FALSE)
  px=dpois(x, mu_x)
  py=dpois(y, mu_y)
  px=px/sum(px)
  py=py/sum(py)
  return(bdd_min_list   (px, py))
}

#**************************************************************