Как создать псевдослучайную положительно определенную матрицу с ограничениями на недиагональные элементы - PullRequest
4 голосов
/ 24 июня 2009

Пользователь хочет наложить уникальную, нетривиальную верхнюю / нижнюю границу для корреляции между каждой парой переменной в матрице var / covar.

Например: я хочу матрицу отклонений, в которой все переменные имеют 0,9> | rho (x_i, x_j) | > 0.6, где rho (x_i, x_j) - корреляция между переменными x_i и x_j.

Спасибо.


Хорошо, было найдено какое-то быстрое и грязное решение, но если кто-то знает более точный точный способ добраться туда, это будет приветствоваться.


Я потерял свой первоначальный логин, поэтому я пересылаю вопрос под новым логином. Предыдущая итерация получила следующий ответ

* вы имеете в виду псевдослучайный, это правильная терминология для полу случайных - Роберт Гулд

* Хороший вопрос, но я думаю, что он имел в виду полупсевдослучайный (псевдослучайность подразумевается, когда речь идет о компьютерной случайности :-p) - fortran

* Под "корреляцией" вы подразумеваете "ковариацию"? - Сванте

* нет, я действительно имею в виду корреляцию. Я хочу создать положительно определенную матрицу, такую, чтобы все корреляции имели более жесткие, чем тривиальные границы. - вак

* Смотрите мой ответ. Вы настаиваете, чтобы выборочные корреляции лежали в указанных пределах, или только корреляции совокупности, которые генерируют выборку? Я предлагаю идею, которая может сработать, если ваша проблема - первая. - щепа

* woodship: нет, боюсь, ваше решение не сработает, пожалуйста, смотрите мой ответ в оригинальной угрозе (ссылка выше). Спасибо.

Ответы [ 4 ]

2 голосов
/ 11 октября 2013

Вы можете создать набор из N случайных векторов размера M и единичной дисперсии. И добавить к ним случайный вектор (размер N и единичная дисперсия), умноженный на определенное число k. Затем вы берете корреляцию между всеми этими векторами, которая будет положительно определенной матрицей. Если M очень большое, тогда не будет никакой дисперсии в распределении корреляции, и корреляция будет: k ^ 2 / (1 + k ^ 2) Чем меньше М, тем шире распределение недиагональных элементов. В качестве альтернативы, вы можете позволить M быть очень большим и умножить «общий вектор» на разные k каждый. Вы можете получить более жесткий контроль, если будете правильно играть с этими параметрами. Вот код Matlab для этого:

clear all;
vecLarg=10;
theDim=1000;
corrDist=0*randn(theDim,1);
Baux=randn(vecLarg,theDim)+  (corrDist*randn(1,vecLarg))'+(k*ones(theDim,1)*randn(1,vecLarg))'  ;
A=corrcoef(Baux);
hist(A(:),100);
2 голосов
/ 25 июня 2009

Вот ваш ответ на мой ответ в оригинальной теме:

"Давай люди, должно быть что-то попроще"

Извините, но нет. Желания выиграть в лотерею недостаточно. Требовать, чтобы новички выиграли серию, недостаточно. Также вы не можете просто требовать решения математической задачи и вдруг обнаружите, что это легко.

Проблема генерирования псевдослучайных отклонений с параметрами выборки в указанном диапазоне является нетривиальной, по крайней мере, если отклонения должны быть действительно псевдослучайными в любом смысле. В зависимости от диапазона, одному может повезти. Я предложил схему отказа, но также заявил, что это вряд ли будет хорошим решением. Если в корреляциях много измерений и узких диапазонов, то вероятность успеха мала. Также важен размер выборки, так как он будет определять дисперсию выборки результирующих корреляций.

Если вы действительно хотите найти решение, вам нужно сесть и указать свою цель, четко и точно. Вам нужна случайная выборка с номинальной заданной структурой корреляции, но строгими границами для корреляций? Будет ли удовлетворительной любая выборочная корреляционная матрица, которая удовлетворяет ограничению целей? Различия также даны?

1 голос
/ 26 июня 2009

ОК, фантастический Грегг: мы куда-то добираемся. Комбинируя вашу идею с древесной щепой, вы получаете этот альтернативный подход. Математически очень грязно, но, похоже, работает:

library(MCMCpack)
library(MASS)
p<-10
lb<-.6
ub<-.8
zupa<-function(theta){
    ac<-matrix(theta,p,p)
    fe<-rwish(100*p,ac%*%t(ac))
    det(fe)
}
ba<-optim(runif(p^2,-10,-5),zupa,control=list(maxit=10))
ac<-matrix(ba$par,p,p)
fe<-rwish(100*p,ac%*%t(ac))
me<-mvrnorm(p+1,rep(0,p),fe)
A<-cor(me)
bofi<-sqrt(diag(var(me)))%*%t(sqrt((diag(var(me)))))
va<-A[lower.tri(A)]
l1=100
while(l1>0){
    r1<-which(va>ub)
    l1<-length(r1)
    va[r1]<-va[r1]*.9
}
A[lower.tri(A)]<-va
A[upper.tri(A)]<-va
vari<-bofi*A
mk<-mvrnorm(10*p,rep(0,p),vari)
pc<-sign(runif(p,-1,1))
mf<-sweep(mk,2,pc,"*")
B<-cor(mf)
summary(abs(B[lower.tri(B)]))

По сути, это идея (скажем, верхняя граница = .8 и нижняя граница = .6), у нее достаточно хороший показатель приемлемости, который не равен 100%, но она подойдет на данном этапе проекта. ,

1 голос
/ 26 июня 2009

Возможно, этот ответ поможет реализовать его на практике:

Одним классом матриц, обладающим этим свойством неотрицательной определенности, является Wishart Distribution . А сэмплы из ~ W () такие, что все недиагональные записи находятся между границами [l, u], будут соответствовать вашему вопросу. Однако я не считаю, что это то же самое, что распределение всех положительно определенных матриц с недиагоналями в [l, u].

На странице википедии есть алгоритм для вычисления из ~ W ().

Более простое, хакерское решение (возможно, приближающееся к этому):

(учитывая, что u> l и l> 0)

  1. взять из многовариантной нормали, где сигма = среднее (l, u).
  2. Затем берём выборку, вычисляем её матрицу корреляции => C
  3. Эта матрица будет иметь некоторую случайность (нечеткость), но математика того, сколько она будет иметь, немного не в моей области. Значения недиагностиков в этой C-матрице ограничены [-1,1] со средним значением (l, u). По глазному яблоку я предполагаю какую-то бета / экспоненциальную. В любом случае, это непрерывное распределение отключенных диаграмм в C гарантирует, что оно не будет вести себя и лежать внутри границ (l, u), если (l, u) = [-1,1].
  4. Вы можете отрегулировать количество «пуха», увеличив / уменьшив длину выборки на шаге 1. Я бы поспорил (бездоказательно), что величина дисперсии в нечетных диагонали C пропорциональна квадратному корню из количество образцов.

Так что, кажется, нетривиально, чтобы действительно ответить!

Как и предлагали другие постеры, вы можете сгенерировать из Wishart, а затем оставить те, где свойство, которое вы хотите, соответствует действительности, но вы можете пробовать в течение длительного времени! Если вы исключите тех, кто является 0-определенным (это слово?), То это должно хорошо работать для генерации хороших матриц. Однако это не является истинным распределением всех матриц pos-def, чьи недиагностики находятся в [l, u].

Код (в R) для предложенной выше схемы немой выборки

sigma1 <- function(n,sigma) {
    out <- matrix(sigma,n,n)
    diag(out) <- 1
    return (out)
}

library(mvtnorm)
sample_around_sigma <- function(size, upper,lower, tight=500) {
    #  size:  size of matrix
    #  upper, lower:  bounds on the corr, should be > 0
    #  tight:  number of samples to use.  ideally this
    #     would be calcuated such that the odd-diags will
    #     be "pretty likely" to fall in [lower,upper]
    sigma <- sigma1(size,mean(c(upper,lower)))
    means <- 0*1:size
    samples <- rmvnorm(n=tight, mean=means,sigma=sigma)
    return (cor(samples))
}

> A <- sample_around_sigma(5, .3,.5)
> A
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.3806354 0.3878336 0.3926565 0.4080125
[2,] 0.3806354 1.0000000 0.4028188 0.4366342 0.3801593
[3,] 0.3878336 0.4028188 1.0000000 0.4085453 0.3814716
[4,] 0.3926565 0.4366342 0.4085453 1.0000000 0.3677547
[5,] 0.4080125 0.3801593 0.3814716 0.3677547 1.0000000
> 
> summary(A[lower.tri(A)]); var(A[lower.tri(A)])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3678  0.3808  0.3902  0.3947  0.4067  0.4366 
[1] 0.0003949876
...