Générer une variable aléatoire avec une corrélation définie avec une ou plusieurs variables existantes

71

Pour une étude de simulation , je dois générer des variables aléatoires qui montrent une corrélation prefined (population) à une variable existante .Y

J'ai examiné les Rpackages copulaet ceux CDVinequi peuvent produire des distributions multivariées aléatoires avec une structure de dépendance donnée. Cependant, il n'est pas possible de fixer l'une des variables résultantes à une variable existante.

Toutes les idées et les liens vers des fonctions existantes sont appréciés!


Conclusion: Deux réponses valables ont été trouvées, avec des solutions différentes:

  1. Un R script de caracal, qui calcule une variable aléatoire avec une corrélation exacte (échantillon) avec une variable prédéfinie
  2. Une R fonction que je me suis trouvée, qui calcule une variable aléatoire avec une corrélation de population définie à une variable prédéfinie

[Ajout de @ttnphns: j’ai pris la liberté d’élargir le titre de la question d’un cas à une variable fixe à un nombre arbitraire de variables fixes; comment générer une variable ayant une ou plusieurs corrections prédéfinies avec une ou plusieurs variables fixes, existantes]

Felix S
la source
2
Consultez cette question connexe stats.stackexchange.com/questions/13382/… qui aborde directement votre question (au moins son aspect théorique).
Macro
Le Q suivant est également étroitement lié et sera d’intérêt: Comment générer des nombres aléatoires corrélés (variances moyennes et degré de corrélation donnés) .
gung - Rétablir Monica

Réponses:

56

En voici un autre: pour les vecteurs de moyenne 0, leur corrélation est égale au cosinus de leur angle. Donc, une façon de trouver un vecteur avec exactement la corrélation souhaitée r , correspondant à un angle θ :xrθ

  1. obtenir un vecteur fixe et un vecteur aléatoire x 2x1x2
  2. centrer les deux vecteurs (moyenne 0), donnant les vecteurs , ˙ x 2x˙1x˙2
  3. faire orthogonal à ˙ x 1 (projection sur orthogonal), donnant ˙ x 2x˙2x˙1x˙2
  4. échelle et ˙ x 2 à la longueur 1, donnant ˉ x 1 et ˉ x 2x˙1x˙2x¯1x¯2
  5. est le vecteur dontangle d' ° x 1estθ, et dontcorrélation avec ˉ x 1 est doncr. C’est aussi la corrélation avecx1puisque les transformations linéaires laissent la corrélation inchangée.x¯2+(1/tan(θ))x¯1x¯1θx¯1rx1

Voici le 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

entrez la description de l'image ici

Pour la projection orthogonale , j'ai utilisé la décomposition Q R pour améliorer la stabilité numérique, puis simplement P = Q Q .PQRP=QQ

caracal
la source
J'essayais de réécrire le code dans la syntaxe SPSS. Je trébuche sur votre décomposition QR qui renvoie 20x1 colonne. Dans SPSS, j'ai une orthonormalisation de Gram-Schmidt (qui est également une décomposition QR) mais je ne peux pas répliquer votre colonne Q résultante. Pouvez-vous mâcher votre action QR pour moi s'il vous plaît. Ou indiquez une solution pour obtenir la projection. Merci.
ttnphns
@caracal, P <- X %*% solve(t(X) %*% X) %*% t(X)ne produit pas r = 0,6, ce n'est donc pas la solution de rechange . Je suis encore confus. (Je serais heureux d'imiter votre expression Q <- qr.Q(qr(Xctr[ , 1, drop=FALSE]))dans SPSS mais je ne sais pas comment.)
ttnphns
@ttnphns Désolé pour la confusion, mon commentaire était pour le cas général. Application à la situation dans l'exemple: Obtenir la matrice de projection via une décomposition QR sert uniquement à la stabilité numérique. Vous pouvez obtenir la matrice de projection si le sous - espace est engendré par les colonnes de la matrice X . Dans R, vous pouvez écrire ici car le sous-espace est recouvert par la première colonne de . La matrice pour la projection sur le complément orthogonal est alors IP. P=X(XX)1XXXctr[ , 1] %*% solve(t(Xctr[ , 1]) %*% Xctr[ , 1]) %*% t(Xctr[ , 1])Xctr
Caracal
4
Quelqu'un pourrait-il expliquer comment réaliser quelque chose de similaire pour plus de deux échantillons? Dites, si je voulais 3 échantillons corrélés par paire par rho, comment puis-je transformer cette solution pour y parvenir?
Andre Terra
pour le cas limite rho=1je l' ai trouvé utile de faire quelque chose comme ceci: if (isTRUE(all.equal(rho, 1))) rho <- 1-10*.Machine$double.eps, sinon je devenais NaNs
PatrickT
19

Je décrirai la solution la plus générale possible. Résoudre le problème dans cette généralité nous permet de réaliser une implémentation logicielle remarquablement compacte: deux lignes de Rcode suffisent.

Choisissez un vecteur , de la même longueur que Y , selon la distribution de votre choix . Soit Y soit les résidus de la régression des moindres carrés de X contre Y : cet extrait le Y composant de X . En ajoutant de nouveau un multiple approprié de Y à Y , on peut produire un vecteur ayant une corrélation souhaitée ρ avec Y . Jusqu'à une constante additive arbitraire et une constante multiplicative positive - que vous êtes libre de choisir, de toute façon que ce soit - la solution est:XYYXYYXYYρY

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

SD


RX

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

Y50XY;ρYX=(1,2,,50)Y

Figure

Il y a une similitude remarquable entre les parcelles, n'est-ce pas :-).


Si vous souhaitez expérimenter, voici le code qui a généré ces données et la figure. (Je n'ai pas pris la peine d'utiliser la liberté de modifier et d'ajuster les résultats, ce qui est une opération facile.)

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)

Ce qui suit est une implémentation plus complète pour ceux qui souhaitent expérimenter.

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
la source
YBTW, this method readily generalizes to more... Just use ordinary least squares... and form a suitable linear combination
1
@ttnphns je l'ai fait.
whuber
1
Merci beaucoup! Je vois, et j'ai codé votre approche aujourd'hui dans SPSS pour moi-même. Vraiment super proposition de votre part. Je n'ai jamais pensé à la notion de double base comme applicable pour résoudre la tâche.
ttnphns
Est-il possible d'utiliser une approche similaire pour créer un vecteur uniformément distribué? C’est-à-dire que j’ai un vecteur existant xet que je veux générer un nouveau vecteur en ycorrélation avec xle yvecteur mais que je veuille aussi que le vecteur soit uniformément distribué.
Skumin
@Skumin Pensez à utiliser une copule pour pouvoir contrôler la relation entre les deux vecteurs.
whuber
6

Voici une autre approche informatique (la solution est adaptée d’un post de forum d’Enrico Schumann). Selon Wolfgang (voir commentaires), il s’agit d’un calcul identique à la solution proposée par tnphns.

ρρ

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

La fonction peut également utiliser des distributions marginales non normales en ajustant le paramètre mar.fun. Notez cependant que fixer une variable ne semble fonctionner qu'avec une variable normalement distribuée x! (qui pourrait se rapporter au commentaire de Macro).

Notez également que le "petit facteur de correction" de la publication d'origine a été supprimé car il semble biaiser les corrélations résultantes, du moins dans le cas des distributions gaussiennes et des corrélations de Pearson (voir également les commentaires).

Felix S
la source
ρ
1
Il est facile de montrer que, à l’exception de la "petite correction à apporter à rho" (dont le but dans ce contexte m’échappe), c’est exactement la même chose que ce que nous suggéraient plus tôt. La méthode est simplement basée sur la décomposition de Choleski de la matrice de corrélation pour obtenir la matrice de transformation souhaitée. Voir, par exemple: en.wikipedia.org/wiki/… . Et oui, cela ne vous donnera que deux vecteurs dont la corrélation de population est égale à rho.
Wolfgang
La "petite correction à rho" était dans le post original et est décrite ici . En fait, je ne le comprends pas vraiment. mais une étude de 50000 corrélations simulées avec rho = .3 montre que, sans la "petite correction", une moyenne de rs de .299 est produite, tandis qu'avec la correction, une moyenne de .312 (qui est la valeur du rho corrigé) est produit. Par conséquent, j'ai supprimé cette partie de la fonction.
Felix S
Je sais que cela est ancien, mais je tiens également à noter que cette méthode ne fonctionnera pas pour les matrices de corrélation définies non positives. Par exemple, une corrélation de -1.
vendredi
1
Merci; J'ai remarqué que si x1 est pas normalisée moyenne = 0, sd = 1, et que vous préférez ne pas redimensionnez, vous devrez modifier la ligne: X2 <- mar.fun(n)pour X2 <- mar.fun(n,mean(x),sd(x))obtenir la corrélation souhaitée entre x1 et x2
Dave M
6

XYXrXrY=rX+EE0sd=1r2XYrXYXρ=r

rEXEXYX1,X2,X3,...

XrYYrY


Mise à jour du 11 novembre 2017. J'ai rencontré cet ancien fil de discussion aujourd'hui et décidé d'élargir ma réponse en montrant l'algorithme de l'ajustement itératif dont je parlais initialement.

Y X

Disclamer: Cette solution itérative que j'ai trouvée est inférieure à l'excellente solution basée sur la recherche de la double base et proposée par @whuber dans ce fil de discussion aujourd'hui. La solution de @Wuber n'est pas itérative et, plus important encore pour moi, elle semble affecter les valeurs de la variable d'entrée "pig" un peu moins que "mon" algorithme (ce serait alors un atout si la tâche est de "corriger" la variable existante et ne pas générer une variable aléatoire à partir de zéro). Néanmoins, je publie le mien par curiosité et parce que cela fonctionne (voir aussi la note en bas de page).

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

    (le dénominateur ne change pas lors des itérations, calculez-le à l'avance)

    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

tnphns
la source
1
Merci pour votre réponse. C’est une solution empirique / itérative à laquelle je pensais également. Pour mes simulations, cependant, il me faut une solution plus analytique sans procédure d’ajustement coûteuse. Heureusement, je viens de trouver une solution que je publierai bientôt ...
Felix S
Cela fonctionne pour générer des normales à deux variables mais ne fonctionne pas pour une distribution arbitraire (ou une distribution non 'additive')
Macro
1
Je ne vois pas pourquoi vous proposez l'itération quand vous pouvez produire directement le cône entier de solutions. Y a-t-il un but particulier à cette approche?
whuber
1
Y
1
@ Whuber, votre commentaire est ce que j'attendais; En fait, ma réponse (à propos de l'hétéroscédasticité, à laquelle je renvoie) a été conçue pour vous interpeller: c'est peut-être une invitation à publier votre solution - aussi complète et brillante que vous le faites habituellement.
ttnphns
4

J'avais envie de faire de la programmation, alors j'ai pris la réponse supprimée de @ Adam et j'ai décidé d'écrire une belle implémentation en R. Je me suis concentré sur un style orienté fonctionnellement (c'est-à-dire une boucle de style lapply). L'idée générale est de prendre deux vecteurs, permuter de façon aléatoire l'un des vecteurs jusqu'à ce qu'une certaine corrélation soit atteinte entre eux. Cette approche est très brutale, mais simple à mettre en œuvre.

Tout d'abord, nous créons une fonction qui permute de manière aléatoire le vecteur d'entrée:

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

... et créer des exemples de données

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

... écrivez une fonction qui permute le vecteur d'entrée et le corrèle à un vecteur de référence:

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

... et itérer mille fois:

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

Notez que les règles de portée de R garantissent que vec1et vec2se trouvent dans l'environnement global, en dehors de la fonction anonyme utilisée ci-dessus. Ainsi, les permutations sont toutes relatives aux jeux de données de test d'origine que nous avons générés.

Ensuite, nous trouvons la corrélation maximale:

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

... ou trouvez la valeur la plus proche d'une corrélation de 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

Pour obtenir une corrélation plus élevée, vous devez augmenter le nombre d'itérations.

Paul Hiemstra
la source
2

Y1Y2,,YnR

Solution:

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

Code Python:

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

Test de sortie:

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
la source
Y1
@whuber c'était une faute de frappe
Aksakal
0

Générer des variables normales avec la matrice de covariance SAMPLING telle que donnée

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;

Générer des variables normales avec la matrice de covariance de POPULATION telle que donnée

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
utilisateur3635627
la source
2
Vous devez apprendre à formater le code dans la réponse! Il existe une option spécifique pour marquer le texte en tant que fragments de code, utilisez-le!
kjetil b halvorsen
-6

Il suffit de créer un vecteur aléatoire et de trier jusqu'à obtenir le résultat souhaité r.

Adam
la source
Dans quelles situations cela serait-il préférable aux solutions ci-dessus?
Andy W
Une situation où un utilisateur veut une réponse simple. J'ai lu une question similaire sur le forum, et c'est la réponse qui a été donnée.
Adam
3
r
3
Si cette réponse a été donnée sur le forum de r-help, je suppose que c'était soit (a) ironique (c'est-à-dire destiné à être une blague), soit (b) proposé par quelqu'un qui n'est pas très sophistiqué sur le plan statistique. Pour le dire plus succinctement, c’est une mauvaise réponse à la question. -1
gung - Rétablir Monica