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

71

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]

Felix S
quelle
2
Sehen Sie sich diese verwandte Frage stats.stackexchange.com/questions/13382/… an, die Ihre Frage direkt anspricht (zumindest die theoretische Seite).
Makro
Das folgende Q ist ebenfalls stark verwandt und wird von Interesse sein: Wie man korrelierte Zufallszahlen erzeugt (gegeben bedeutet Varianzen und Grad der Korrelation) .
gung - Wiedereinsetzung von Monica

Antworten:

56

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 θxrθ

  1. erhalten Sie den festen Vektor und einen zufälligen Vektorx 2x1x2
  2. Zentrieren Sie beide Vektoren (Mittelwert 0) und geben Sie die Vektoren ,x 2x˙1x˙2
  3. mache orthogonal zu (Projektion auf den orthogonalen Unterraum) und gebe ˙ x 1 ˙ x 2x˙2x˙1x˙2
  4. Skalieren Sie und auf die Länge 1, und geben Sie und ˙ x 2 ˉ x 1 ˉ x 2x˙1x˙2x¯1x¯2
  5. ˉ x 1θ ˉ x 1Rx1x¯2+(1/tan(θ))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.x¯1θx¯1rx1

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

Bildbeschreibung hier eingeben

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

caracal
quelle
Ich habe versucht, den Code in die SPSS-Syntax umzuschreiben. Ich stolpere über deine QR-Zerlegung, die eine 20x1-Spalte zurückgibt. In SPSS habe ich eine Gram-Schmidt-Orthonormalisierung (die auch eine QR-Zerlegung ist), kann aber Ihre resultierende Q-Spalte nicht replizieren. Kannst du mir bitte deine QR-Aktion durchkauen? Oder geben Sie eine Umgehungsmöglichkeit an, um die Projektion zu erhalten. Vielen Dank.
TTNPHNS
@caracal, P <- X %*% solve(t(X) %*% X) %*% t(X)erzeugt kein r = 0.6, das ist also nicht die Umgehung . Ich bin immer noch verwirrt. (Ich würde gerne Ihren Ausdruck Q <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))in SPSS imitieren , weiß aber nicht wie.)
ttnphns
@ttnphns Entschuldigung für die Verwirrung, mein Kommentar war für den allgemeinen Fall. Anwendung auf die Situation im Beispiel: Das Abrufen der Projektionsmatrix über QR-Zerlegung dient nur der numerischen Stabilität. Sie können die Projektionsmatrix als wenn der Unterraum von den Spalten der Matrix überspannt wird . In R können Sie hier schreiben, da der Unterraum von der ersten Spalte von überspannt wird . Die Matrix für die Projektion auf das orthogonale Komplement ist dann IP. XP=X(XX)1XXXctr[ , 1] %*% solve(t(Xctr[ , 1]) %*% Xctr[ , 1]) %*% t(Xctr[ , 1])Xctr
Karakal
4
Könnte jemand erklären, wie man etwas Ähnliches für mehr als nur zwei Proben ausführt? Angenommen, ich wollte 3 Samples, die durch Rho paarweise korreliert werden. Wie kann ich diese Lösung transformieren, um dies zu erreichen?
Andre Terra
Für den Grenzfall rho=1fand ich es nützlich, so etwas zu machen: if (isTRUE(all.equal(rho, 1))) rho <- 1-10*.Machine$double.epsAnsonsten bekam ich NaNs
PatrickT
19

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; YXYYXYYXYYρY

XY;ρ=ρSD(Y)Y+1ρ2SD(Y)Y.

(" " steht für eine Berechnung, die proportional zu einer Standardabweichung ist.)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.Y50XY;ρYX=(1,2,,50)Y

Zahl

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")

YXY1,Y2,,Yk;ρ1,ρ2,,ρkYiYiXYiYY

RYiy

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))
whuber
quelle
YBTW, this method readily generalizes to more... Just use ordinary least squares... and form a suitable linear combination
1
@ttnphns Ich habe es getan.
whuber
1
Vielen Dank! Ich verstehe, und ich habe Ihren Ansatz heute in SPSS für mich selbst codiert. Ein wirklich großartiger Vorschlag von dir. Ich hätte nie gedacht, dass der Begriff der doppelten Basis für die Lösung der Aufgabe anwendbar ist.
ttnphns
Ist es möglich, mit einem ähnlichen Ansatz einen gleichmäßig verteilten Vektor zu erhalten? Das heißt, ich habe einen vorhandenen Vektor xund möchte einen neuen Vektor erzeugen, der ymit diesem korreliert ist x, möchte aber auch, dass der yVektor gleichmäßig verteilt ist.
Skumin
@Skumin Verwenden Sie dazu eine Kopula, damit Sie die Beziehung zwischen den beiden Vektoren steuern können.
whuber
6

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.

ρρ

ρ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).

Felix S
quelle
ρ
1
Es ist leicht zu zeigen, dass dies bis auf die "kleine Korrektur von Rho" (deren Zweck sich mir in diesem Zusammenhang entzieht) genau das ist , was ttnphns zuvor vorgeschlagen hat. Das Verfahren basiert einfach auf der Choleski-Zerlegung der Korrelationsmatrix, um die gewünschte Transformationsmatrix zu erhalten. Siehe zum Beispiel: en.wikipedia.org/wiki/… . Und ja, dies ergibt nur zwei Vektoren, deren Populationskorrelation gleich ist rho.
Wolfgang
Die "kleine Korrektur zu Rho" war im ursprünglichen Beitrag und wird hier beschrieben . Eigentlich verstehe ich es nicht wirklich; Eine Untersuchung von 50000 simulierten Korrelationen mit rho = .3 zeigt jedoch, dass ohne die "kleine Korrektur" ein Durchschnitt von rs von .299 erzeugt wird, während mit der Korrektur ein Durchschnitt von .312 (das ist der Wert des korrigierten rho) ist produziert. Deshalb habe ich diesen Teil aus der Funktion entfernt.
Felix S
Ich weiß, dass dies alt ist, aber ich möchte auch darauf hinweisen, dass diese Methode für nicht positiv definierte Korrelationsmatrizen nicht funktioniert. ZB - eine Korrelation von -1.
zzk
1
Vielen Dank; Ich habe bemerkt , dass , wenn x1 nicht Mittelwert normiert = 0, sd = 1, und Sie würden lieber nicht neu skaliert, werden Sie die Zeile ändern müssen: X2 <- mar.fun(n)auf X2 <- mar.fun(n,mean(x),sd(x))die gewünschte Korrelation zwischen x1 und x2 zu erhalten
Dave M
6

XYXrXrY=rX+EE0sd=1r2XYrXYXρ=r

rEXEXYX1,X2,X3,...

XrYYrY


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).

X1,X2,...,XmYYr1,r2,...,rmX

YXYY

  1. rdf=n1Sj=rjdfjX

  2. dfYXdf

  3. YXrb=(XX)1S

  4. YY^=Xb

  5. E=YY^

  6. SSS=dfSSY^

  7. EXjCj=i=1nEiXij

  8. EC0i

    Ei[corrected]=Eij=1mCjXijnj=1mXij2

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

    E0 EC

    Ei[corrected]=Eij=1mCjXij3i=1nXij2j=1mXij2

    1

  9. SSEEi[corrected]=EiSSS/SSE

    mrSSSn

  10. CErYY[corrected]=Y^+E

  11. Y

  12. Yr

YrY


1YX

ttnphns
quelle
1
Danke für deine Antwort. Das ist eine empirische / iterative Lösung, über die ich auch nachgedacht habe. Für meine Simulationen benötige ich jedoch eine analytischere Lösung ohne kostspielige Anpassungsprozedur. Zum Glück habe ich gerade eine Lösung gefunden, die ich in Kürze veröffentlichen werde ...
Felix S
Dies funktioniert zum Generieren von bivariaten Normalen, funktioniert jedoch nicht für eine beliebige Verteilung (oder eine beliebige nicht 'additive' Verteilung)
Makro
1
Ich verstehe nicht, warum Sie eine Iteration vorschlagen, wenn Sie den gesamten Konus von Lösungen direkt produzieren können. Hat dieser Ansatz einen besonderen Zweck?
whuber
1
Y
1
@whuber, auf deinen Kommentar habe ich gewartet; Tatsächlich war meine Antwort (zu Heteroskedastizität, auf die ich verweise) als Herausforderung für Sie gedacht: Vielleicht ist es eine Einladung, Ihre Lösung zu veröffentlichen - so gründlich und brillant wie gewöhnlich.
TTNPHNS
4

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.

Paul Hiemstra
quelle
2

Y1Y2,,YnR

Lösung:

  1. CCT=R
  2. X2,,XnY1
  3. Y1
  4. Y=CXYiY1

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.        ]]
Aksakal
quelle
Y1
@whuber es war ein Tippfehler
Aksakal
0

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
user3635627
quelle
2
Sie müssen lernen, den Code in der Antwort zu formatieren! Es gibt eine spezielle Option, um Text als Codefragmente zu markieren. Verwenden Sie diese Option!
kjetil b halvorsen
-6

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

Adam
quelle
In welchen Situationen wäre dies den obigen Lösungen vorzuziehen?
Andy W
Eine Situation, in der ein Benutzer eine einfache Antwort wünscht. Ich habe eine ähnliche Frage im Forum gelesen, und es ist die Antwort, die gegeben wurde.
Adam
3
r
3
Wenn diese Antwort im r-help-Forum gegeben wurde, war sie vermutlich entweder (a) ironisch (dh als Scherz gedacht) oder (b) von jemandem angeboten, der statistisch nicht sehr ausgefeilt ist. Um es kurz zu machen, dies ist eine schlechte Antwort auf die Frage. -1
Gung - Monica Wiedereinstellung