为了创建一个数据集来测试统计计算包，我希望能够生成一个与具有给定希尔曼系数的参考样本相关的样本。
我成功地为皮尔逊系数做到了这一点，但至少可以说，斯皮尔曼在排名上工作的事实使得它相当棘手！
作为示例，使用生成 Spearman 和 Pearson 相关样本的代码：
from statistics import correlation
import numpy as np
from scipy import stats
def generate_pearson_correlated_to_sample(x, correlation):
"""
Generate a variable with a specified Pearson correlation coefficient to a given sample.
Parameters
----------
x : array-like
The fixed first sample.
correlation : float
Desired Pearson correlation coefficient (-1 <= correlation <= 1).
Returns
-------
array-like
The generated sample y with specified Pearson correlation to x.
"""
# standardize the pixel data
x_std = (x - x.mean()) / x.std()
# generate independent standard normal data
z = np.random.normal(loc=0, scale=1, size=x_std.shape)
# Create correlated variable (standardized)
y_std = correlation * x_std + np.sqrt(1 - correlation ** 2) * z
# Scale to desired standard deviation and add mean
mean = x.mean()
std = x.std()
y = y_std * std + mean
return y
def generate_spearman_correlated_to_sample(x, correlation):
"""
Generate a variable with a specified Spearman correlation coefficient to a given sample.
Parameters
----------
x : array-like
The fixed first sample.
correlation : float
Desired Spearman correlation coefficient (-1 <= correlation <= 1).
Returns
-------
array-like
The generated sample y with specified Spearman correlation to x.
"""
n_samples = len(x)
# Convert x to ranks (normalized between 0 and 1)
x_ranks = stats.rankdata(x) / (n_samples + 1)
# Convert ranks to normal distribution
normal_x = stats.norm.ppf(x_ranks)
# Generate correlated normal variable
normal_y = correlation * normal_x + np.sqrt(1 - correlation ** 2) * np.random.normal(0, 1, n_samples)
# Convert back to uniform distribution
y_uniform = stats.norm.cdf(normal_y)
# Convert uniform to same marginal distribution as x using empirical CDF
x_sorted = np.sort(x)
y = np.interp(y_uniform, np.linspace(0, 1, n_samples), x_sorted)
return y
def verify_correlations(x, y):
"""
Calculate both Spearman and Pearson correlations between two variables.
Parameters:
-----------
x, y : array-like
The two variables to check
Returns:
--------
tuple
(spearman_correlation, pearson_correlation)
"""
spearman_corr = stats.spearmanr(x, y)[0]
pearson_corr = stats.pearsonr(x, y)[0]
return spearman_corr, pearson_corr
# Example usage
if __name__ == "__main__":
# Set random seed for reproducibility
np.random.seed(42)
# Create different types of example data
x_normal = np.random.normal(0, 1, 10000) # Normal distribution
x_exp = np.random.exponential(2, 10000) # Exponential distribution
x_bimodal = np.concatenate([np.random.normal(-2, 0.5, 5000),
np.random.normal(2, 0.5, 5000)]) # Bimodal distribution
# Test with different distributions and correlations
test_cases = [
(x_normal, 0.7, "Normal Distribution"),
(x_exp, 0.5, "Exponential Distribution"),
(x_bimodal, 0.8, "Bimodal Distribution")
]
# Run examples
for x, target_corr, title in test_cases:
print(f"\nTesting with {title}")
print(f"Target correlation: {target_corr}")
# Generate correlated sample
y_spearman = generate_spearman_correlated_to_sample(x, correlation=target_corr)
y_pearson = generate_pearson_correlated_to_sample(x, correlation=target_corr)
# Calculate actual correlations
spearman_corr, _= verify_correlations(x, y_spearman)
_, pearson_corr = verify_correlations(x, y_pearson)
print(f"Achieved Spearman correlation: {spearman_corr:.4f}")
print(f"Achieved Pearson correlation: {pearson_corr:.4f}")
使用上述代码，由于代码的随机性，生成的皮尔逊系数当然不完全等于目标值。 但我发现 Spearman 的系统偏差要大得多，这让我怀疑我的代码中存在问题。
我使用 Python 工作，但感谢任何帮助！
一个简单的方法是将其放入最小二乘法作为抛光方法，并用您的初始方法作为估计：
def generate_pearson_correlated_to_sample(
x: np.ndarray,
correlation: float,
) -> np.ndarray:
xerror = x - x.mean()
lhs = correlation*np.sqrt(xerror.dot(xerror))
# standardize the pixel data
mean = x.mean()
std = x.std()
x_std = (x - mean)/std
# generate independent standard normal data
z = np.random.normal(loc=0, scale=1, size=x_std.shape)
# Create correlated variable (standardized)
y_std = correlation*x_std + np.sqrt(1 - correlation**2)*z
# Scale to desired standard deviation and add mean
y0 = y_std*std + mean
def error_fun(y: np.ndarray) -> float:
yerror = y - y.mean()
corr_error = xerror.dot(yerror) - lhs*np.sqrt(yerror.dot(yerror))
return corr_error
result = least_squares(
fun=error_fun, x0=y0,
method='dogbox', jac='2-point', x_scale='jac', loss='linear')
assert result.success
return result.x
仅与 Pearson 一起显示；斯皮尔曼也会一样：
Testing with Normal Distribution
Target correlation: 0.7
Achieved Pearson correlation: 0.7000
Testing with Exponential Distribution
Target correlation: 0.5
Achieved Pearson correlation: 0.5000
Testing with Bimodal Distribution
Target correlation: 0.8
Achieved Pearson correlation: 0.8000