如何用R生成满足约束条件的对称联合概率矩阵？

Alright, let's figure out how to generate that symmetric joint probability matrix in R while meeting all three constraints you listed. I'll walk you through a couple of approaches, starting with a straightforward method and then moving to a more efficient stepwise version.

First, let's recap the matrix structure to make sure we're aligned:

• It's an n×n symmetric matrix where:
  • Diagonal elements: p(Ai) (probability that binary variable Ai is 1), satisfying 0 ≤ p(Ai) ≤ 1 (constraint 1)
  • Off-diagonal elements: p(AiAj) (probability both Ai and Aj are 1), which must fit between max(0, p(Ai)+p(Aj)−1) and min(p(Ai), p(Aj)) (constraint 2)
  • Every triplet of variables must satisfy p(Ai)+p(Aj)+p(Ak)−p(AiAj)−p(AiAk)−p(AjAk) ≤ 1 (constraint 3, ensuring the joint distribution is valid)

This approach generates random values that satisfy constraints 1 and 2 first, then checks constraint 3. If invalid, it retries until a valid matrix is found.

generate_joint_prob_matrix <- function(n, max_tries = 1000) {
  for (try in 1:max_tries) {
    # Step 1: Generate valid single-variable probabilities (constraint 1)
    p_single <- runif(n, 0, 1)
    
    # Step 2: Initialize symmetric matrix and fill diagonal
    prob_mat <- matrix(0, nrow = n, ncol = n)
    diag(prob_mat) <- p_single
    
    # Step 3: Fill off-diagonal elements (constraint 2)
    for (i in 1:(n-1)) {
      for (j in (i+1):n) {
        lower_bound <- max(0, p_single[i] + p_single[j] - 1)
        upper_bound <- min(p_single[i], p_single[j])
        prob_mat[i,j] <- runif(1, lower_bound, upper_bound)
        prob_mat[j,i] <- prob_mat[i,j]  # Enforce symmetry
      }
    }
    
    # Step 4: Validate constraint 3 (add small tolerance for floating point errors)
    is_valid <- TRUE
    for (i in 1:(n-2)) {
      for (j in (i+1):(n-1)) {
        for (k in (j+1):n) {
          lhs <- p_single[i] + p_single[j] + p_single[k] -
            prob_mat[i,j] - prob_mat[i,k] - prob_mat[j,k]
          if (lhs > 1 + 1e-8) {
            is_valid <- FALSE
            break
          }
        }
        if (!is_valid) break
      }
      if (!is_valid) break
    }
    
    if (is_valid) {
      cat(paste("Valid matrix generated after", try, "attempts.\n"))
      return(prob_mat)
    }
  }
  stop(paste("Failed to generate a valid matrix after", max_tries, "attempts."))
}

How to Use This Function #

# Set seed for reproducibility
set.seed(456)
# Generate a 4×4 joint probability matrix
my_prob_matrix <- generate_joint_prob_matrix(4)
print(my_prob_matrix)

# Verify all constraints manually (optional)
# Check constraint 1
all(diag(my_prob_matrix) >= 0 & diag(my_prob_matrix) <= 1)
# Check constraint 2
n <- nrow(my_prob_matrix)
constraint2_ok <- TRUE
for (i in 1:(n-1)) {
  for (j in (i+1):n) {
    p_i <- diag(my_prob_matrix)[i]
    p_j <- diag(my_prob_matrix)[j]
    p_ij <- my_prob_matrix[i,j]
    if (p_ij < max(0, p_i + p_j -1) - 1e-8 || p_ij > min(p_i, p_j) + 1e-8) {
      constraint2_ok <- FALSE
      break
    }
  }
  if (!constraint2_ok) break
}
constraint2_ok
# Check constraint 3
constraint3_ok <- TRUE
for (i in 1:(n-2)) {
  for (j in (i+1):(n-1)) {
    for (k in (j+1):n) {
      lhs <- diag(my_prob_matrix)[i] + diag(my_prob_matrix)[j] + diag(my_prob_matrix)[k] -
        my_prob_matrix[i,j] - my_prob_matrix[i,k] - my_prob_matrix[j,k]
      if (lhs > 1 + 1e-8) {
        constraint3_ok <- FALSE
        break
      }
    }
    if (!constraint3_ok) break
  }
  if (!constraint3_ok) break
}
constraint3_ok

The first method can waste retries when n is big, since constraint 3 is only checked after filling the entire matrix. This stepwise approach builds the matrix one variable at a time, ensuring all constraints are satisfied as we go.

generate_joint_prob_stepwise <- function(n) {
  # Initialize empty matrix
  prob_mat <- matrix(0, nrow = n, ncol = n)
  
  # First variable
  prob_mat[1,1] <- runif(1, 0, 1)
  
  # Add variables one by one
  for (k in 2:n) {
    # Generate p(Ak) (constraint 1)
    p_k <- runif(1, 0, 1)
    prob_mat[k,k] <- p_k
    
    # Generate p(AiAk) for all existing i < k
    for (i in 1:(k-1)) {
      # Base bounds from constraint 2
      lower <- max(0, prob_mat[i,i] + p_k - 1)
      upper <- min(prob_mat[i,i], p_k)
      
      # Add extra bounds from constraint 3 with all previous j < i
      if (i > 1) {
        extra_lowers <- sapply(1:(i-1), function(j) {
          prob_mat[i,i] + p_k + prob_mat[j,j] - 1 - prob_mat[i,j] - prob_mat[j,k]
        })
        lower <- max(lower, max(extra_lowers))
      }
      
      # If bounds are invalid, restart the entire process (simplified backtracking)
      if (lower > upper) {
        return(generate_joint_prob_stepwise(n))
      }
      
      prob_mat[i,k] <- runif(1, lower, upper)
      prob_mat[k,i] <- prob_mat[i,k]
    }
  }
  
  return(prob_mat)
}

Why This Works #

By adding variables incrementally, we ensure that every new joint probability p(AiAk) satisfies both the pairwise constraint (2) and all triplet constraints (3) with the variables we've already added. This drastically reduces the chance of needing full retries, especially for larger n.

内容的提问来源于stack exchange，提问作者BS.Mira