Identify position in each list element when one value becomes greater than another

Question

I have a list, lprobs. Each list element is a simple vector of probabilities (P[Z≥z]) under different parameter values (defined below in kvec).

When comparing these probabilities, the lower parameter value initially always has a lower probability when z>=2. Then there is a "crossover" at some value of z where the lower parameter value begins yielding a higher probability. I need to determine the the z value and probability when these crossovers occur for all combinations of k values.

The following plot may help visualize this:

The black (k = 0.01) is the easiest example to point out: I need to find the z value and P(Z=z) when P(Z=z) for k = 0.01 crosses over with k = 0.25, 0.5, 0.75, 1.0 (somewhere after z > 300). If there is not a crossover by the maximum z value (defined in the below code by cutoff), it should return NA (i.e. when comparing k values of 0.01 and 0.25).

I have tried several for loops and manually finding these values, but I have not found an elegant solution to systematically record the z value and probability for each crossover. The actual data are much larger.

The desired result would be as below (though open to suggestions, doesn't have to be exactly this), where z_val and prob_z identify the specific z value and probability of each crossover event:

## putting 1234 as a dummy z value and 0.567 as dummy probability 

want <- data.frame(expand.grid(kvec, kvec), z_val = 1234, prob_z = 0.567)
want <- want[want[,1] != want[,2],c(2:1,3:4)]
colnames(want)[1:2] <- c("k1", "k2")

#    k1   k2    z_val prob_z
# 2  0.01 0.25  1234  0.567
# 3  0.01 0.50  1234  0.567
# 4  0.01 0.75  1234  0.567
# 5  0.01 1.00  1234  0.567
#...
# 21 1.00 0.01  1234  0.567
# 22 1.00 0.25  1234  0.567
# 23 1.00 0.50  1234  0.567
# 24 1.00 0.75  1234  0.567

The data to recreate these probabilities and the figure is here:


## Functions to calculate probabilities
clprob <- function(R, k, z, n = 1){
  lnumer = log(n) + lgamma(k*z+z-n) + (z-n)*(log(R)-log(k))  
  ldenom = log(z) + lgamma(k*z) + lgamma(z-n+1) + (k*z+z-n)*(log(k+R)-log(k))  
  return(exp(lnumer - ldenom))
}
clcdf <- function(R, k, z, n = 1){
  1-sum(clprob(R,k,1:(z-1),n))
}

# Define parameters
kvec <- c(0.01, 0.25, 0.50, 0.75, 1.0) 
cutoff <- 500
clsize <- 2:cutoff 

# List to store probabilities
lprobs <- lapply(kvec, function(x) sapply(clsize, clcdf, R = 0.90, k = x))

# Plot
plot(NA, xlim = c(2, cutoff), ylim = c(1e-3, 1), log = 'y', 
     ylab = expression("P[Z">="z] (log scale)"), xlab = expression("z (">="2)"))
for (i in 1:length(lprobs)){points(lprobs[[i]], pch = 19, col = i, cex = 0.25)}
legend("topright", legend = paste("k =", kvec), col=1:length(kvec), pch=16, bty='n')

user_4447540 · Answer

I believe you can do something like this:

names(lprobs) = kvec
result = as.data.frame(cbind(
  t(combn(kvec,2)),
  t(sapply(combn(kvec,2, simplify = F), function(x) {
      zval = which(lprobs[[as.character(x[1])]] >= lprobs[[as.character(x[2])]])[1]
      prob = lprobs[[as.character(x[1])]][zval]
      return(c(zval,prob))}))
))
names(result) = c("k1", "k2", "z_val", "prob_z")

Output:

     k1   k2 z_val      prob_z
1  0.01 0.25    NA          NA
2  0.01 0.50   426 0.002858335
3  0.01 0.75   362 0.003171162
4  0.01 1.00   328 0.003373438
5  0.25 0.50   144 0.013438024
6  0.25 0.75   128 0.015014513
7  0.25 1.00   119 0.016052089
8  0.50 0.75   104 0.019002732
9  0.50 1.00    98 0.020167960
10 0.75 1.00    88 0.022700673

As the above shows, the black line (0.01) never crosses the red line (0.25). The black line does cross the light blue line (1.00) at Z=328.

Identify position in each list element when one value becomes greater than another

1 Answers