如何用R生成固定均值标准差的80×80对称联合概率矩阵？

Got it, let's work through how to build that 80x80 symmetric joint probability matrix in R. You've got specific constraints on marginal probabilities (the diagonals) and the joint probabilities (off-diagonals), plus fixed mean and standard deviation for those marginals. Here's a step-by-step solution that covers all your requirements, including validation to make sure the matrix is statistically valid:

First: Define Parameters & Generate Valid Marginal Probabilities #

First, set your target mean (mu_marginal) and standard deviation (sigma_marginal) for the marginal probabilities (the diagonal elements of your matrix). Since these probabilities have to stay between 0 and 1, using a truncated normal distribution is the easiest way to generate values that hit your mean/std dev while respecting the 0-1 bound.

We'll use the truncnorm package for this—if you don't have it installed, the code below will handle that:

# Install/load the truncated norm package if needed
if (!require(truncnorm)) {
  install.packages("truncnorm")
  library(truncnorm)
}

# Core parameters
n <- 80  # Matrix dimension
mu_marginal <- 0.5  # Replace with your desired marginal mean
sigma_marginal <- 0.15  # Replace with your desired marginal std dev

# Generate marginal probabilities constrained to [0, 1]
marginal_probs <- rtruncnorm(n, a = 0, b = 1, mean = mu_marginal, sd = sigma_marginal)

# Quick check to confirm we're close to target stats (randomness will mean slight differences)
cat("Marginal Probability Mean:", round(mean(marginal_probs), 4), "\n")
cat("Marginal Probability Std Dev:", round(sd(marginal_probs), 4), "\n")

If you don't want to use external packages, you can generate a normal sample and manually truncate it, but note this might push your mean/std dev slightly off target:

# No-package alternative: Generate normal, then clamp to 0-1
marginal_probs <- rnorm(n, mean = mu_marginal, sd = sigma_marginal)
marginal_probs <- pmax(pmin(marginal_probs, 1), 0)

Second: Initialize the Symmetric Matrix #

Create an empty 80x80 matrix and populate the diagonal with your marginal probabilities (since the diagonal is p(Ai), the probability that variable Ai is 1):

# Initialize empty symmetric matrix
joint_matrix <- matrix(0, nrow = n, ncol = n)
# Fill diagonal with marginal probabilities
diag(joint_matrix) <- marginal_probs

Third: Fill Off-Diagonal Joint Probabilities #

For the off-diagonal elements (p(AiAj)), we need to generate values that satisfy:
max(0, p(Ai) + p(Aj) - 1) ≤ p(AiAj) ≤ min(p(Ai), p(Aj))

Since the matrix is symmetric, we only need to fill the upper triangle and mirror it to the lower triangle to save computation:

# Fill upper triangle and mirror to lower triangle
for (i in 1:(n-1)) {
  for (j in (i+1):n) {
    p_i <- marginal_probs[i]
    p_j <- marginal_probs[j]
    # Calculate valid bounds for the joint probability
    lower_bound <- max(0, p_i + p_j - 1)
    upper_bound <- min(p_i, p_j)
    # Generate a random value within the bounds
    joint_matrix[i, j] <- runif(1, min = lower_bound, max = upper_bound)
    # Mirror to lower triangle for symmetry
    joint_matrix[j, i] <- joint_matrix[i, j]
  }
}

Fourth: Validate the Matrix (Critical Step!) #

Just satisfying the element-wise constraints isn't enough—your joint probability matrix corresponds to a covariance matrix that must be positive semi-definite (all eigenvalues ≥ 0) to be statistically valid (meaning there actually exists a set of binary variables with these probabilities).

We'll check this, and if needed, fix the matrix using the Matrix package's nearPD function (which finds the closest positive semi-definite matrix):

# Calculate covariance matrix from joint probabilities
cov_matrix <- joint_matrix - outer(marginal_probs, marginal_probs)

# Check if covariance matrix is positive semi-definite (allow tiny numerical errors)
eigen_vals <- eigen(cov_matrix, only.values = TRUE)$values
is_valid <- all(eigen_vals >= -1e-8)

if (is_valid) {
  cat("✅ Joint probability matrix is statistically valid!\n")
} else {
  cat("⚠️ Warning: Covariance matrix is not positive semi-definite. Fixing it...\n")
  # Install/load Matrix package if needed
  if (!require(Matrix)) {
    install.packages("Matrix")
    library(Matrix)
  }
  # Get the closest positive semi-definite covariance matrix
  corrected_cov <- nearPD(cov_matrix, corr = FALSE)$mat
  # Reconstruct the joint probability matrix from the corrected covariance
  joint_matrix <- corrected_cov + outer(marginal_probs, marginal_probs)
  # Ensure all values stay within valid [0,1] bounds after correction
  joint_matrix <- pmax(pmin(joint_matrix, 1), 0)
  # Restore original marginal probabilities (correction might tweak them slightly)
  diag(joint_matrix) <- marginal_probs
  cat("✅ Matrix corrected to be statistically valid.\n")
}

Final Checks #

Let's verify all constraints are still satisfied after any corrections:

# Check marginal probabilities are still in [0,1]
marginal_ok <- all(diag(joint_matrix) >= 0 & diag(joint_matrix) <= 1)
# Check all off-diagonal elements meet their constraints
joint_ok <- TRUE
for (i in 1:n) {
  for (j in 1:n) {
    if (i != j) {
      p_i <- diag(joint_matrix)[i]
      p_j <- diag(joint_matrix)[j]
      lower <- max(0, p_i + p_j - 1)
      upper <- min(p_i, p_j)
      if (joint_matrix[i,j] < lower - 1e-8 || joint_matrix[i,j] > upper + 1e-8) {
        joint_ok <- FALSE
        break
      }
    }
  }
  if (!joint_ok) break
}

if (marginal_ok && joint_ok) {
  cat("✅ All constraints are satisfied! Your matrix is ready to use.\n")
} else {
  cat("❌ Some constraints are violated. Double-check the matrix.\n")
}

Quick Notes #

• If you need exact mean and standard deviation for your marginal probabilities (not just close), you'll need to tweak the generated vector. One way is to iteratively adjust values until you hit the target stats, while keeping them within 0-1.
• The nearPD correction is a robust way to fix invalid covariance matrices, but it will make small adjustments to your joint probabilities. If you need minimal changes, this is the best approach.

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