Como gerar números aleatórios correlacionados (médias, variações e grau de correlação)?

53

Sinto muito se isso parece um pouco básico, mas acho que estou apenas tentando confirmar o entendimento aqui. Tenho a sensação de que eu teria que fazer isso em duas etapas e comecei a tentar criar matrizes de correlação, mas isso está começando a parecer realmente envolvido. Estou procurando uma explicação concisa (idealmente com dicas para uma solução de pseudocódigo) de uma maneira boa e idealmente rápida de gerar números aleatórios correlacionados.

Dadas duas variáveis ​​pseudo-aleatórias, altura e peso, com médias e variações conhecidas e uma determinada correlação, acho que estou basicamente tentando entender como deve ser a segunda etapa:

   height = gaussianPdf(height.mean, height.variance)
   weight = gaussianPdf(correlated_mean(height.mean, correlation_coefficient), 
                        correlated_variance(height.variance, 
                        correlation_coefficient))
  • Como calculo a média e variância correlacionadas? Mas quero confirmar que esse é realmente o problema relevante aqui.
  • Preciso recorrer à manipulação de matrizes? Ou tenho outra coisa muito errada na minha abordagem básica para esse problema?
Joseph Weissman
fonte
11
Não sei se o entendi corretamente, mas você não precisa calcular "média e variância correlacionadas". Se você está assumindo que as variáveis ​​são bivariadas normais, deve ser suficiente especificar as médias e variações individuais e a correlação. Existe algum software específico que você deseja usar para isso?
mark999

Respostas:

44

Para responder à sua pergunta sobre "uma maneira boa e idealmente rápida de gerar números aleatórios correlacionados": Dada a matriz variância-covariância desejada que é definida por definição positiva positiva, a decomposição de Cholesky é: C = L L T ; L sendo a matriz triangular inferior.CCLLTL

Se você agora usar essa matriz para projetar um vetor variável aleatório não correlacionado X , a projeção resultante Y = L X será a das variáveis ​​aleatórias correlacionadas.LXY=LX

Você pode encontrar uma explicação concisa sobre por que isso acontece aqui .

usεr11852 diz Reinstate Monic
fonte
Obrigado! Isso foi extremamente útil. Eu acho que pelo menos tenho uma noção melhor do que eu preciso olhar a seguir.
Joseph Weissman
7
Esse método se aplica apenas a distribuições gaussianas (conforme especificado na pergunta) ou pode ser usado para gerar variáveis ​​correlatas que seguem outras distribuições? Caso contrário, você conhece um método que poderia ser usado nesse caso?
user000001
11
@ Michael: Sim. Dito isto, dado que é uma matriz de covariância válida, a decomposição de Cholesky é o caminho mais rápido. Você também pode obter a matriz X de raiz quadrada (simétrica) de C usando SVD (então C = X X = X X T , em que X = U S 0,5 V T de C = U S V T ), mas isso seria mais caro também. CXCC=XX=XXTX=US0.5VTC=USVT
usεr11852 diz Reinstate Monic
11
@ Michael: Claro. Sua covariância será (aproximadamente) a mesma, não os números em si.
usεr11852 diz Reinstate Monic 03/04
11
@ Sid: Qualquer distribuição contínua não suportada em toda a linha real falhará imediatamente. Por exemplo, se usarmos um uniforme , não podemos garantir que os "números correlatos" estarão em [ 0 , 1 ] ; da mesma forma para um Poisson, acabaremos com números não discretos. Além disso, qualquer distribuição em que a soma das distribuições ainda não seja a mesma (por exemplo, somar distribuição t não resulta em distribuições t ) também falhará. Em todos os casos mencionados, os números produzidos serão correlacionados de acordo com a CU[0,1][0,1]ttCmas eles não corresponderão à distribuição que começamos.
usεr11852 diz Reinstate Monic
36

+1 a @ user11852 e @ jem77bfp, são boas respostas. Deixe-me abordar isso de uma perspectiva diferente, não porque eu acho que é necessariamente melhor na prática , mas porque eu acho que é instrutivo. Aqui estão alguns fatos relevantes que já sabemos:

  1. é a inclinação da reta de regressão quando X e Y sãopadronizados, isto é, N ( 0 , 1 ) , rXYN(0,1)
  2. é a proporção da variação em Y atribuível à variação em X , r2YX



    (também, das regras para variações ):

  3. a variação de uma variável aleatória multiplicada por uma constante é a constante ao quadrado vezes a variação original:
    Var[aX]=a2Var[X]
  4. as variações adicionam , ou seja, a variação da soma de duas variáveis ​​aleatórias (supondo que sejam independentes) é a soma das duas variações:
    Var[X+ε]=Var[X]+Var[ε]

rρXN(0,1)aveYN(0,a2+ve) . (Observe que | a | deve ser1 para que isso funcione e que, além disso, a = r .) Assim, você começa com o r que deseja; esse é o seu coeficiente, a . Então você descobre a variação de erro que precisará, é 1 - r 2 . (Se o seu software exigir que você use o desvio padrão, pegue a raiz quadrada desse valor.) Finalmente, para cada variável pseudo-aleatória, x i , que você gerou, gere uma variável pseudo-aleatória, e ia2+ve=1|a| 1a=rra1r2xiei, com a variação de erro apropriada , e calcule a variável pseudo-aleatória correlacionada, y i , multiplicando e adicionando. veyi

Se você quiser fazer isso no R, o seguinte código pode funcionar para você:

correlatedValue = function(x, r){
  r2 = r**2
  ve = 1-r2
  SD = sqrt(ve)
  e  = rnorm(length(x), mean=0, sd=SD)
  y  = r*x + e
  return(y)
}

set.seed(5)
x = rnorm(10000)
y = correlatedValue(x=x, r=.5)

cor(x,y)
[1] 0.4945964

0

Novamente, isso, em sua forma mais simples, permite gerar apenas um par de variáveis ​​correlacionadas (isso pode ser ampliado, mas fica muito rápido) e certamente não é a maneira mais conveniente de realizar o trabalho. Em R, você desejaria usar ? Mvrnorm no pacote MASS , tanto porque é mais fácil quanto porque você pode gerar muitas variáveis ​​com uma matriz de correlação populacional. No entanto, acho que vale a pena ter percorrido esse processo para ver como alguns princípios básicos funcionam de maneira simples.

- Reinstate Monica
fonte
Essa abordagem essencialmente regressiva é particularmente agradável, permitindo gerar um Y aleatório correlacionado com qualquer número de "preditores" X existentes . Estou certo em tal entendimento?
ttnphns
Depende exatamente de qual padrão de correlação entre as variáveis ​​que você deseja, @ttnphns. Você pode iterar isso um após o outro, mas seria tedioso. Para criar muitas variáveis ​​correlacionadas com um determinado padrão, é melhor usar a decomposição de Cholesky.
gung - Restabelece Monica
você sabe como usar Cholesky para gerar um Y correlacionado (aproximadamente, como no seu método) de acordo com um vetor de correlações com vários Xs existentes (não simulados)?
ttnphns
@ttnphns, você deseja gerar um único Y com uma determinada correlação populacional com um conjunto de Xs, não um conjunto de variáveis ​​p com todas as correlações populacionais pré-especificadas? Uma maneira simples seria escrever uma equação de regressão para gerar um único Y-hat a partir dos seus Xs e, em seguida, use o método acima para gerar Y como um correlato do seu Y-hat. Você pode fazer uma nova pergunta, se quiser.
gung - Restabelece Monica
11
Isso é o que eu quis dizer no meu comentário inicial: esse método seria uma extensão direta do que você fala na sua resposta: essencialmente um método regressivo (Hat).
ttnphns
16

Em geral, isso não é uma coisa simples de se fazer, mas acredito que existem pacotes para geração de variáveis normais multivariadas (pelo menos em R, veja mvrnormno MASSpacote), onde você apenas insere uma matriz de covariância e um vetor médio.

(X1,X2)F(x1,x2)Fx2

FX1(x1)=F(x1,x2)dx2.
FX11FX1ξ1[0,1]x^1=FX11(ξ)

F(x1,x2)x1=x^1

F(x2|X1=x^1)=F(x^1,x2)fX1(x^1),
fX1X1FX1(x1)=fX1(x1)

ξ2[0,1]ξ1F(x2|X1=x^1)x^2=(F(x2|X1=x^1))1(ξ)x^2F(x^2|X1=x^1)=ξ

Se você não entender o significado de conectar uma variável uniforme a uma função de distribuição de probabilidade inversa, tente fazer um esboço do caso univariado e lembre-se de qual é a interpretação geométrica da função inversa.

jem77bfp
fonte
Idéia inteligente! Tem apelo intuitivo simples. Mas sim parece caro computacionalmente.
MichaelChirico 28/09
(+1) ponto muito bom. Seria melhor no começo dizerfX,Y(x,y)=fX(x)fY|X(y)
1

Se você estiver pronto para desistir da eficiência, poderá usar um alogoritmo descartável. Sua vantagem é que ela permite qualquer tipo de distribuição (não apenas gaussiana).

{xi}i=1N{yi}i=1NC

cold=corr({xi},{yi})

n1n2:1n1,2N

xn1xn2

cnew=corr({xi},{yi})

|Ccnew|<|Ccold|

|Cc|<ϵ

xi

Boa sorte!

F. Jatpil
fonte
xicorr(xi,yi)
xi{xi}ycorr(xi,yi)corr({xi},{yi})=(1/N)Σi=1N(xix¯)(yyy¯)
Entendo, faz todo sentido. Ignorei o " " em c o r r ( { x i }{}corr({xi},{yi})