生成10000组(X,Y)样本，手动估算统计量与相关系数（禁用R专用函数）

Hey there! Let's work through this problem step by step, building everything from scratch in R as you requested. No fancy built-in correlation functions here—just good old manual calculations based on sample expectations.

First, we need a function to generate bivariate (X,Y) samples. I'll use a multivariate normal distribution as an example (you can swap this out with your own custom function if you have one). This function takes the number of samples, expectation vector, and covariance matrix as inputs:

generate_bivariate_samples <- function(n, mu, cov_matrix) {
  # n: Number of (X,Y) pairs to generate
  # mu: Vector of true expectations (c(E[X], E[Y]))
  # cov_matrix: 2x2 covariance matrix for X and Y
  samples <- MASS::mvrnorm(n = n, mu = mu, Sigma = cov_matrix)
  colnames(samples) <- c("X", "Y")
  return(as.data.frame(samples))
}

Note: We use MASS::mvrnorm here to generate multivariate normal samples, but you can replace this with your own custom generation logic if needed.

Let's pick some true parameters for demonstration (you can adjust these to match your specific expectations/covariances):

# Set true population parameters
true_mu <- c(2, 5)          # E[X] = 2, E[Y] = 5
true_cov <- matrix(c(4, 3, 3, 9), nrow = 2)  # Cov(X,X)=4, Cov(X,Y)=3, Cov(Y,Y)=9

# Generate 10,000 samples (set seed for reproducibility)
set.seed(123)
sample_data <- generate_bivariate_samples(n = 10000, mu = true_mu, cov_matrix = true_cov)

We'll use sample means to estimate each expectation—since sample means are unbiased estimators of population expectations:

# Calculate sample-based estimates
est_E_X <- mean(sample_data$X)
est_E_X2 <- mean(sample_data$X^2)
est_E_Y <- mean(sample_data$Y)
est_E_Y2 <- mean(sample_data$Y^2)
est_E_XY <- mean(sample_data$X * sample_data$Y)

# Print the results
cat("Estimated E(X):", round(est_E_X, 4), "\n")
cat("Estimated E(X²):", round(est_E_X2, 4), "\n")
cat("Estimated E(Y):", round(est_E_Y, 4), "\n")
cat("Estimated E(Y²):", round(est_E_Y2, 4), "\n")
cat("Estimated E(XY):", round(est_E_XY, 4), "\n")

Recall the formula for the Pearson correlation coefficient:

ρ = Cov(X,Y) / (σ_X * σ_Y)
Where:

• Cov(X,Y) = E[XY] - E[X]E[Y]
• σ_X = sqrt(E[X²] - (E[X])²) (standard deviation of X)
• σ_Y = sqrt(E[Y²] - (E[Y])²) (standard deviation of Y)

Let's implement this using our estimated expectations:

# Estimate covariance between X and Y
est_cov_XY <- est_E_XY - est_E_X * est_E_Y

# Estimate standard deviations of X and Y
est_sd_X <- sqrt(est_E_X2 - est_E_X^2)
est_sd_Y <- sqrt(est_E_Y2 - est_E_Y^2)

# Calculate the correlation coefficient
est_correlation <- est_cov_XY / (est_sd_X * est_sd_Y)

# Print the final result
cat("Estimated Correlation Coefficient:", round(est_correlation, 4), "\n")

Optional Verification (For Your Reference) #

If you want to confirm our manual calculation is correct, you can compare it to R's built-in cor function (though we're not using this for the actual task):

cat("Built-in Correlation Check:", round(cor(sample_data$X, sample_data$Y), 4), "\n")

All steps here avoid using R's dedicated correlation functions—we built every part from basic sample mean calculations, which meets your requirements perfectly.

内容的提问来源于stack exchange，提问作者Han