Generieren Sie eine Zufallsvariable mit einer definierten Korrelation zu einer oder mehreren vorhandenen Variablen.

Für eine Simulationsstudie muss ich Zufallsvariablen generieren, die eine vorab festgelegte (Populations-) Korrelation zu einer vorhandenen Variablen .$Y$

Ich sah in die RPakete copulaund CDVineder Zufall multivariate Verteilungen mit einer bestimmten Abhängigkeitsstruktur erzeugen kann. Es ist jedoch nicht möglich, eine der resultierenden Variablen an eine vorhandene Variable zu binden.

Anregungen und Links zu bestehenden Funktionen sind willkommen!

Schlussfolgerung: Es wurden zwei gültige Antworten mit unterschiedlichen Lösungen gefunden:

1. Ein R Skript von caracal, das eine Zufallsvariable mit einer exakten (Beispiel-) Korrelation zu einer vordefinierten Variablen berechnet
2. Eine R Funktion, die ich selbst gefunden habe und die eine Zufallsvariable mit einer definierten Populationskorrelation zu einer vordefinierten Variablen berechnet

[@ttnphns 'Zusatz: Ich habe mir die Freiheit genommen, den Fragentitel von einem einzelnen Fall fester Variablen auf eine beliebige Anzahl fester Variablen zu erweitern; dh wie man eine Variable mit vordefinierten Korretationen mit einer oder mehreren festen, existierenden Variablen erzeugt]

Hier ist eine andere: Bei Vektoren mit dem Mittelwert 0 entspricht ihre Korrelation dem Kosinus ihres Winkels. Ein Weg, um einen Vektor mit genau der gewünschten Korrelation , die einem Winkel :r $x$$r$$\theta$

1. erhalten Sie den festen Vektor und einen zufälligen Vektorx $x_1$$x_2$
2. Zentrieren Sie beide Vektoren (Mittelwert 0) und geben Sie die Vektoren , ≤ x$\dot{x}_{1}$$\dot{x}_{2}$
3. mache orthogonal zu (Projektion auf den orthogonalen Unterraum) und gebe ˙ x 1 ˙ x ⊥ $\dot{x}_{2}$$\dot{x}_{1}$$\dot{x}_{2}^{\perp}$
4. Skalieren Sie und auf die Länge 1, und geben Sie und ˙ x ⊥ 2 ˉ x 1 ˉ x ⊥ $\dot{x}_{1}$$\dot{x}_{2}^{\perp}$$\bar{x}_{1}$$\bar{x}_{2}^{\perp}$
5. ˉ x 1θ ˉ x 1Rx$\bar{x}_{2}^{\perp} + (1/\tan(\theta)) \cdot \bar{x}_{1}$ ist der Vektor, dessen Winkel zu ist , und dessen Korrelation mit also . Dies ist auch die Korrelation zu da lineare Transformationen die Korrelation unverändert lassen.$\bar{x}_{1}$$\theta$$\bar{x}_{1}$$r$$x_1$

Hier ist der Code:

n     <- 20                    # length of vector
rho   <- 0.6                   # desired correlation = cos(angle)
theta <- acos(rho)             # corresponding angle
x1    <- rnorm(n, 1, 1)        # fixed given data
x2    <- rnorm(n, 2, 0.5)      # new random data
X     <- cbind(x1, x2)         # matrix
Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)

Id   <- diag(n)                               # identity matrix
Q    <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))      # QR-decomposition, just matrix Q
P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
x2o  <- (Id-P) %*% Xctr[ , 2]                 # x2ctr made orthogonal to x1ctr
Xc2  <- cbind(Xctr[ , 1], x2o)                # bind to matrix
Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1

x <- Y[ , 2] + (1 / tan(theta)) * Y[ , 1]     # final new vector
cor(x1, x)                                    # check correlation = rho

Für die orthogonale Projektion ich die Zerlegung verwendet, um die numerische Stabilität zu verbessern, da dann einfach .Q R P = Q Q $P$$QR$$P = Q Q'$

Ich werde die allgemeinste mögliche Lösung beschreiben. Die Lösung des Problems in dieser Allgemeinheit ermöglicht es uns, eine bemerkenswert kompakte Softwareimplementierung zu erzielen: Nur zwei kurze RCodezeilen reichen aus.

Wähle einen Vektor , der die gleiche Länge wie , nach einem der Verteilung Sie mögen. Lassen die Residuen der Regression der kleinsten Quadrate der seine gegen : Diese extrahiert die - Komponente von . Indem wir ein geeignetes Vielfaches von zu , können wir einen Vektor mit jeder gewünschten Korrelation mit erzeugen . Bis zu einer beliebigen additiven Konstante und positiven multiplikativen Konstante - die Sie nach Belieben wählen können - ist die LösungY Y ⊥ X Y Y X Y Y ⊥ & rgr; $X$$Y$$Y^\perp$$X$$Y$$Y$$X$$Y$$Y^\perp$$\rho$$Y$

(" " steht für eine Berechnung, die proportional zu einer Standardabweichung ist.)$\operatorname{SD}$

Hier ist RArbeitscode. Wenn Sie kein angeben, der Code seine Werte aus der multivariaten Standardnormalverteilung.$X$

complement <- function(y, rho, x) {
  if (missing(x)) x <- rnorm(length(y)) # Optional: supply a default if `x` is not given
  y.perp <- residuals(lm(x ~ y))
  rho * sd(y.perp) * y + y.perp * sd(y) * sqrt(1 - rho^2)
}

Zur Veranschaulichung habe ich ein zufälliges mit 50 Komponenten erzeugt und X Y erzeugt ; ρ mit verschiedenen spezifizierten Korrelationen mit diesem Y . Sie wurden alle mit demselben Startvektor X = ( 1 , 2 , … , 50 ) erstellt . Hier sind ihre Streudiagramme. Die "Rugplots" am unteren Rand jedes Panels zeigen den gemeinsamen Y- Vektor.$Y$$50$$X_{Y;\rho}$$Y$$X=(1,2,\ldots, 50)$$Y$

Es gibt eine bemerkenswerte Ähnlichkeit zwischen den Handlungen, nicht wahr :-).

Wenn Sie experimentieren möchten, finden Sie hier den Code, der diese Daten erzeugt hat, und die Abbildung. (Ich habe mich nicht darum gekümmert, die Freiheit zu nutzen, die Ergebnisse zu verschieben und zu skalieren. Das sind einfache Operationen.)

y <- rnorm(50, sd=10)
x <- 1:50 # Optional
rho <- seq(0, 1, length.out=6) * rep(c(-1,1), 3)
X <- data.frame(z=as.vector(sapply(rho, function(rho) complement(y, rho, x))),
                rho=ordered(rep(signif(rho, 2), each=length(y))),
                y=rep(y, length(rho)))

library(ggplot2)
ggplot(X, aes(y,z, group=rho)) + 
  geom_smooth(method="lm", color="Black") + 
  geom_rug(sides="b") + 
  geom_point(aes(fill=rho), alpha=1/2, shape=21) +
  facet_wrap(~ rho, scales="free")

$Y$$X_{Y_1,Y_2,\ldots,Y_k;\rho_1,\rho_2,\ldots,\rho_k}$$Y_i$$Y_i$$X$$Y_i$$Y$$Y$

R$Y_i$y

y <- scale(y)             # Makes computations simpler
e <- residuals(lm(x ~ y)) # Take out the columns of matrix `y`
y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
return(y.dual %*% rho + sqrt(sigma2)*e)

Das Folgende ist eine vollständigere Implementierung für diejenigen, die experimentieren möchten.

complement <- function(y, rho, x) {
  #
  # Process the arguments.
  #
  if(!is.matrix(y)) y <- matrix(y, ncol=1)
  if (missing(x)) x <- rnorm(n)
  d <- ncol(y)
  n <- nrow(y)
  y <- scale(y) # Makes computations simpler
  #
  # Remove the effects of `y` on `x`.
  #
  e <- residuals(lm(x ~ y))
  #
  # Calculate the coefficient `sigma` of `e` so that the correlation of
  # `y` with the linear combination y.dual %*% rho + sigma*e is the desired
  # vector.
  #
  y.dual <- with(svd(y), (n-1)*u %*% diag(ifelse(d > 0, 1/d, 0)) %*% t(v))
  sigma2 <- c((1 - rho %*% cov(y.dual) %*% rho) / var(e))
  #
  # Return this linear combination.
  #
  if (sigma2 >= 0) {
    sigma <- sqrt(sigma2) 
    z <- y.dual %*% rho + sigma*e
  } else {
    warning("Correlations are impossible.")
    z <- rep(0, n)
  }
  return(z)
}
#
# Set up the problem.
#
d <- 3           # Number of given variables
n <- 50          # Dimension of all vectors
x <- 1:n         # Optionally: specify `x` or draw from any distribution
y <- matrix(rnorm(d*n), ncol=d) # Create `d` original variables in any way
rho <- c(0.5, -0.5, 0)          # Specify the correlations
#
# Verify the results.
#
z <- complement(y, rho, x)
cbind('Actual correlations' = cor(cbind(z, y))[1,-1],
      'Target correlations' = rho)
#
# Display them.
#
colnames(y) <- paste0("y.", 1:d)
colnames(z) <- "z"
pairs(cbind(z, y))

Hier ist ein weiterer rechnerischer Ansatz (die Lösung stammt aus einem Forumsbeitrag von Enrico Schumann). Laut Wolfgang (siehe Kommentare) ist dies rechnerisch identisch mit der von ttnphns vorgeschlagenen Lösung.

$\rho$$\rho$

$\rho$x

# returns a data frame of two variables which correlate with a population correlation of rho
# If desired, one of both variables can be fixed to an existing variable by specifying x
getBiCop <- function(n, rho, mar.fun=rnorm, x = NULL, ...) {
     if (!is.null(x)) {X1 <- x} else {X1 <- mar.fun(n, ...)}
     if (!is.null(x) & length(x) != n) warning("Variable x does not have the same length as n!")

     C <- matrix(rho, nrow = 2, ncol = 2)
     diag(C) <- 1

     C <- chol(C)

     X2 <- mar.fun(n)
     X <- cbind(X1,X2)

     # induce correlation (does not change X1)
     df <- X %*% C

     ## if desired: check results
     #all.equal(X1,X[,1])
     #cor(X)

     return(df)
}

Die Funktion kann auch nicht normale Randverteilungen verwenden, indem Parameter angepasst werden mar.fun. Beachten Sie jedoch, dass das Fixieren einer Variablen nur mit einer normalverteilten Variablen zu funktionieren scheint x! (was sich auf Macros Kommentar beziehen könnte).

Beachten Sie auch, dass der "kleine Korrekturfaktor" aus dem ursprünglichen Beitrag entfernt wurde, da er die resultierenden Korrelationen zu verzerren scheint, zumindest im Fall von Gaußschen Verteilungen und Pearson-Korrelationen (siehe auch Kommentare).

$X$$Y$$X$$r$$X$$r$$Y= rX+E$$E$$0$$\text{sd}=\sqrt{1-r^2}$$X$$Y$$r$$X$$Y$$X$$\rho=r$

$r$$E$$X$$E$$X$$Y$$X_1, X_2, X_3,...$

$X$$r$$Y$$Y$$r$$Y$

Update 11. November 2017. Ich bin heute auf diesen alten Thread gestoßen und habe beschlossen, meine Antwort zu erweitern, indem ich den Algorithmus der iterativen Anpassung zeige, über die ich ursprünglich gesprochen habe.

$Y$ $X$

Disclamer: Dieser iterative Lösung , die ich schlechter als die ausgezeichnete eine gefunden haben , basierend auf der Erkenntnis , duale Basis und vorgeschlagen von @whuber in diesem Thread heute. Die Lösung von @ whuber ist nicht iterativ und scheint, was für mich wichtiger ist, die Werte der Eingangsvariablen "pig" etwas weniger zu beeinflussen als "my" -Algorithmus (es wäre dann von Vorteil, wenn die Aufgabe darin besteht, "zu korrigieren". die vorhandene Variable und nicht von Grund auf zufällig zu generieren). Trotzdem veröffentliche ich meine aus Neugier und weil es funktioniert (siehe auch Fußnote).

$X_1, X_2,...,X_m$$Y$$Y$$r_1, r_2,...,r_m$$X$

$Y$$X$$Y$$Y$

1. $r$$\text{df}=n-1$$S_j=r_j \text{df}$$j$$X$

2. $\text{df}$$Y$$X$$\text{df}$

3. $Y$$X$$r$$\bf b=(X'X)^{-1} S$

4. $Y$$\hat{Y}=\bf Xb$

5. $E=Y-\hat{Y}$

6. $SS_S=\text{df}-SS_{\hat {Y}}$

7. $E$$X_j$$C_j= \sum_{i=1}^n E_i X_{ij}$

8. $E$$C$$0$$i$

   $$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m C_j X_{ij}} {n\sum_{j=1}^m X_{ij}^2}$$

   (Der Nenner ändert sich bei Iterationen nicht. Berechnen Sie ihn im Voraus.)

   $E$$0$ $E$$C$

   $$E_i[\text{corrected}]=E_i-\frac{\sum_{j=1}^m \frac{C_j X_{ij}^3}{\sum_{i=1}^n X_{ij}^2}} {\sum_{j=1}^m X_{ij}^2}$$

   $^1$

9. $SS_E$$E_i[\text{corrected}]=E_i \sqrt{SS_S/SS_E}$

   $m$$r$$SS_S$$n$

10. $C$$E$$r$$Y$$Y[\text{corrected}]=\hat{Y}+E$

11. $Y$

12. $Y$$r$

$Y$$r$$Y$

$^1$$Y$$X$

Ich wollte etwas programmieren, also habe ich @ Adams gelöschte Antwort genommen und mich dazu entschlossen, eine nette Implementierung in R zu schreiben. Ich konzentriere mich auf die Verwendung eines funktional ausgerichteten Stils (dh Looping mit lappigem Stil). Die allgemeine Idee ist, zwei Vektoren zu nehmen, die zufällig einen der Vektoren permutieren, bis eine bestimmte Korrelation zwischen ihnen erreicht ist. Dieser Ansatz ist sehr brachial, aber einfach zu implementieren.

Zuerst erstellen wir eine Funktion, die den Eingabevektor zufällig permutiert:

randomly_permute = function(vec) vec[sample.int(length(vec))]
randomly_permute(1:100)
  [1]  71  34   8  98   3  86  28  37   5  47  88  35  43 100  68  58  67  82
 [19]  13   9  61  10  94  29  81  63  14  48  76   6  78  91  74  69  18  12
 [37]   1  97  49  66  44  40  65  59  31  54  90  36  41  93  24  11  77  85
 [55]  32  79  84  15  89  45  53  22  17  16  92  55  83  42  96  72  21  95
 [73]  33  20  87  60  38   7   4  52  27   2  80  99  26  70  50  75  57  19
 [91]  73  62  23  25  64  51  30  46  56  39

... und einige Beispieldaten erstellen

vec1 = runif(100)
vec2 = runif(100)

... schreibe eine Funktion, die den Eingabevektor permutiert und ihn mit einem Referenzvektor korreliert:

permute_and_correlate = function(vec, reference_vec) {
    perm_vec = randomly_permute(vec)
    cor_value = cor(perm_vec, reference_vec)
    return(list(vec = perm_vec, cor = cor_value))
  }
permute_and_correlate(vec2, vec1)
$vec
  [1] 0.79072381 0.23440845 0.35554970 0.95114398 0.77785348 0.74418811
  [7] 0.47871491 0.55981826 0.08801319 0.35698405 0.52140366 0.73996913
 [13] 0.67369873 0.85240338 0.57461506 0.14830718 0.40796732 0.67532970
 [19] 0.71901990 0.52031017 0.41357545 0.91780357 0.82437619 0.89799621
 [25] 0.07077250 0.12056045 0.46456652 0.21050067 0.30868672 0.55623242
 [31] 0.84776853 0.57217746 0.08626022 0.71740151 0.87959539 0.82931652
 [37] 0.93903143 0.74439384 0.25931398 0.99006038 0.08939812 0.69356590
 [43] 0.29254936 0.02674156 0.77182339 0.30047034 0.91790830 0.45862163
 [49] 0.27077191 0.74445997 0.34622648 0.58727094 0.92285322 0.83244284
 [55] 0.61397396 0.40616274 0.32203732 0.84003379 0.81109473 0.50573325
 [61] 0.86719899 0.45393971 0.19701975 0.63877904 0.11796154 0.26986325
 [67] 0.01581969 0.52571331 0.27087693 0.33821824 0.52590383 0.11261002
 [73] 0.89840404 0.82685046 0.83349287 0.46724807 0.15345334 0.60854785
 [79] 0.78854984 0.95770015 0.89193212 0.18885955 0.34303707 0.87332019
 [85] 0.08890968 0.22376395 0.02641979 0.43377516 0.58667068 0.22736077
 [91] 0.75948043 0.49734797 0.25235660 0.40125309 0.72147500 0.92423638
 [97] 0.27980561 0.71627101 0.07729027 0.05244047

$cor
[1] 0.1037542

... und tausendmal durchlaufen:

n_iterations = lapply(1:1000, function(x) permute_and_correlate(vec2, vec1))

Beachten Sie, dass die Bereichsregeln von R sicherstellen, dass vec1und vec2in der globalen Umgebung außerhalb der oben verwendeten anonymen Funktion gefunden werden. Die Permutationen sind also alle relativ zu den ursprünglichen Testdatensätzen, die wir generiert haben.

Als nächstes finden wir die maximale Korrelation:

cor_values = sapply(n_iterations, '[[', 'cor')
n_iterations[[which.max(cor_values)]]
$vec
  [1] 0.89799621 0.67532970 0.46456652 0.75948043 0.30868672 0.83244284
  [7] 0.86719899 0.55623242 0.63877904 0.73996913 0.71901990 0.85240338
 [13] 0.81109473 0.52571331 0.82931652 0.60854785 0.19701975 0.26986325
 [19] 0.58667068 0.52140366 0.40796732 0.22736077 0.74445997 0.40125309
 [25] 0.89193212 0.52031017 0.92285322 0.91790830 0.91780357 0.49734797
 [31] 0.07729027 0.11796154 0.69356590 0.95770015 0.74418811 0.43377516
 [37] 0.55981826 0.93903143 0.30047034 0.84776853 0.32203732 0.25235660
 [43] 0.79072381 0.58727094 0.99006038 0.01581969 0.41357545 0.52590383
 [49] 0.27980561 0.50573325 0.92423638 0.11261002 0.89840404 0.15345334
 [55] 0.61397396 0.27077191 0.12056045 0.45862163 0.18885955 0.77785348
 [61] 0.23440845 0.05244047 0.25931398 0.57217746 0.35554970 0.34622648
 [67] 0.21050067 0.08890968 0.84003379 0.95114398 0.83349287 0.82437619
 [73] 0.46724807 0.02641979 0.71740151 0.74439384 0.14830718 0.82685046
 [79] 0.33821824 0.71627101 0.77182339 0.72147500 0.08801319 0.08626022
 [85] 0.87332019 0.34303707 0.45393971 0.47871491 0.29254936 0.08939812
 [91] 0.35698405 0.67369873 0.27087693 0.78854984 0.87959539 0.22376395
 [97] 0.02674156 0.07077250 0.57461506 0.40616274

$cor
[1] 0.3166681

... oder finden Sie den nächsten Wert zu einer Korrelation von 0,2:

n_iterations[[which.min(abs(cor_values - 0.2))]]
$vec
  [1] 0.02641979 0.49734797 0.32203732 0.95770015 0.82931652 0.52571331
  [7] 0.25931398 0.30047034 0.55981826 0.08801319 0.29254936 0.23440845
 [13] 0.12056045 0.89799621 0.57461506 0.99006038 0.27077191 0.08626022
 [19] 0.14830718 0.45393971 0.22376395 0.89840404 0.08890968 0.15345334
 [25] 0.87332019 0.92285322 0.50573325 0.40796732 0.91780357 0.57217746
 [31] 0.52590383 0.84003379 0.52031017 0.67532970 0.83244284 0.95114398
 [37] 0.81109473 0.35554970 0.92423638 0.83349287 0.34622648 0.18885955
 [43] 0.61397396 0.89193212 0.74445997 0.46724807 0.72147500 0.33821824
 [49] 0.71740151 0.75948043 0.52140366 0.69356590 0.41357545 0.21050067
 [55] 0.87959539 0.11796154 0.73996913 0.30868672 0.47871491 0.63877904
 [61] 0.22736077 0.40125309 0.02674156 0.26986325 0.43377516 0.07077250
 [67] 0.79072381 0.08939812 0.86719899 0.55623242 0.60854785 0.71627101
 [73] 0.40616274 0.35698405 0.67369873 0.82437619 0.27980561 0.77182339
 [79] 0.19701975 0.82685046 0.74418811 0.58667068 0.93903143 0.74439384
 [85] 0.46456652 0.85240338 0.34303707 0.45862163 0.91790830 0.84776853
 [91] 0.78854984 0.05244047 0.58727094 0.77785348 0.01581969 0.27087693
 [97] 0.07729027 0.71901990 0.25235660 0.11261002

$cor
[1] 0.2000199

Um eine höhere Korrelation zu erhalten, müssen Sie die Anzahl der Iterationen erhöhen.

$Y_1$$Y_2,\dots,Y_n$$R$

Lösung:

1. $CC^T=R$
2. $X_2,\dots,X_n$$Y_1$
3. $Y_1$
4. $Y=CX$$Y_i$$Y_1$

Python-Code:

import numpy as np
import math
from scipy.linalg import toeplitz, cholesky
from statsmodels.stats.moment_helpers import cov2corr

# create the large correlation matrix R
p = 4
h = 2/p
v = np.linspace(1,-1+h,p)
R = cov2corr(toeplitz(v))

# create the first variable
T = 1000;
y = np.random.randn(T)

# generate p-1 correlated randoms
X = np.random.randn(T,p)
X[:,0] = y
C = cholesky(R)
Y = np.matmul(X,C)

# check that Y didn't change
print(np.max(np.abs(Y[:,0]-y)))

# check the correlation matrix
print(R)
print(np.corrcoef(np.transpose(Y)))

Testausgang:

0.0
[[ 1.   0.5  0.  -0.5]
 [ 0.5  1.   0.5  0. ]
 [ 0.   0.5  1.   0.5]
 [-0.5  0.   0.5  1. ]]
[[ 1.          0.50261766  0.02553882 -0.46259665]
 [ 0.50261766  1.          0.51162821  0.05748082]
 [ 0.02553882  0.51162821  1.          0.51403266]
 [-0.46259665  0.05748082  0.51403266  1.        ]]

Generieren Sie normale Variablen mit der angegebenen SAMPLING-Kovarianzmatrix

covsam <- function(nobs,covm, seed=1237) {; 
          library (expm);
          # nons=number of observations, covm = given covariance matrix ; 
          nvar <- ncol(covm); 
          tot <- nvar*nobs;
          dat <- matrix(rnorm(tot), ncol=nvar); 
          covmat <- cov(dat); 
          a2 <- sqrtm(solve(covmat)); 
          m2 <- sqrtm(covm);
          dat2 <- dat %*% a2 %*% m2 ; 
          rc <- cov(dat2);};
          cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3);
          cm; 
          res <- covsam(10,cm)  ;
          res;

Generieren Sie normale Variablen mit der angegebenen POPULATION-Kovarianzmatrix

covpop <- function(nobs,covm, seed=1237) {; 
          library (expm); 
          # nons=number of observations, covm = given covariance matrix;
          nvar <- ncol(covm); 
          tot <- nvar*nobs;  
          dat <- matrix(rnorm(tot), ncol=nvar); 
          m2 <- sqrtm(covm);
          dat2 <- dat %*% m2;  
          rc <- cov(dat2); }; 
          cm <- matrix(c(1,0.5,0.1,0.5,1,0.5,0.1,0.5,1),ncol=3);
          cm; 
          res <- covpop(10,cm); 
          res

Erstellen Sie einfach einen zufälligen Vektor und sortieren Sie, bis Sie das gewünschte r erhalten.