Генерация коррелированных чисел

Question

Вот забавный пример: мне нужно сгенерировать случайные пары х / у, которые коррелируются с заданным значением коэффициента корреляции моментов произведения Пирсона или Пирсона r . Вы можете представить это как два массива, массив X и массив Y, где значения массива X и массива Y должны быть заново сгенерированы, переупорядочены или преобразованы, пока они не будут коррелированы друг с другом на данном уровне Пирсона r. Вот кикер: Array X и Array Y должны иметь равномерное распределение.

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

Есть идеи?

-

Вот код AS3 для случайных коррелированных гауссиан, чтобы заставить колеса вращаться:

public static function nextCorrelatedGaussians(r:Number):Array{             
         var d1:Number;
         var d2:Number;
         var n1:Number;
         var n2:Number;
         var lambda:Number;
         var r:Number;
         var arr:Array = new Array();
         var isNeg:Boolean; 

        if (r<0){
            r *= -1;
              isNeg=true;
        }            
        lambda= (   (r*r)  -  Math.sqrt(  (r*r) - (r*r*r*r)  )     )   /   ((  2*r*r ) - 1  );

        n1 = nextGaussian();
        n2 = nextGaussian();           
        d1 = n1;            
        d2 = ((lambda*n1) + ((1-lambda)*n2)) / Math.sqrt( (lambda*lambda) + (1-lambda)*(1-lambda));

        if (isNeg) {d2*= -1}           
        arr.push(d1);
        arr.push(d2);
        return arr;
    }

Answer 1

В итоге я написал короткую статью об этом

Он не включает ваш метод сортировки (хотя на практике я думаю, что он похож на мой первый метод, окольным путем), но описывает два способа, которые не требуют итерации.

Answer 2

Странная проблема требует странного решения - вот как я ее решил.

-Генеративный массив X

-Клон массива X для создания массива Y

-Сортировка массива X (вы можете использовать любой метод сортировки массива X - быстрая сортировка, heapsort что-нибудь стабильное.)

-Измеряет начальный уровень Пирсона R с отсортированным массивом X и массивом Y.

WHILE the correlation is outside of the range you are hoping for

   IF the correlation is to low
         run one iteration of CombSort11 on array Y then recheck correlation
   ELSE IF the correlation is too high
         randomly swap two values and recheck correlation

И это все! Комбинация - это настоящий ключ, она медленно и неуклонно увеличивает корреляцию. Посмотрите демо Джейсона Харрисона , чтобы понять, что я имею в виду. Чтобы получить отрицательную корреляцию, вы можете инвертировать сортировку или инвертировать один из массивов после завершения всего процесса.

Вот моя реализация в AS3:

public static function nextReliableCorrelatedUniforms(r:Number, size:int, error:Number):Array {
        var yValues:Array = new Array;
        var xValues:Array = new Array;
        var coVar:Number = 0;
        for (var e:int=0; e < size; e++) { //create x values            
            xValues.push(Math.random());
    }
    yValues = xValues.concat();
    if(r != 1.0){
        xValues.sort(Array.NUMERIC);
    }
    var trueR:Number = Util.getPearson(xValues, yValues);

        while(Math.abs(trueR-r)>error){
            if (trueR < r-error){   // combsort11 for y     
                var gap:int = yValues.length;
                var swapped:Boolean = true; 
                while (trueR <= r-error) {
                    if (gap > 1) {
                        gap = Math.round(gap / 1.3);
                    }
                    var i:int = 0;
                    swapped = false;
                    while (i + gap < yValues.length  &&  trueR <= r-error) {
                        if (yValues[i] > yValues[i + gap]) {
                            var t:Number = yValues[i];
                            yValues[i] = yValues[i + gap];
                            yValues[i + gap] = t;
                            trueR = Util.getPearson(xValues, yValues)
                            swapped = true;
                        }
                        i++;
                    }
                }
            }

            else { // decorrelate
                while (trueR >= r+error) {
                    var a:int = Random.randomUniformIntegerBetween(0, size-1);
                    var b:int = Random.randomUniformIntegerBetween(0, size-1);
                    var temp:Number = yValues[a];
                    yValues[a] = yValues[b];
                    yValues[b] = temp;
                    trueR = Util.getPearson(xValues, yValues)
               }                
            }
        }
        var correlates:Array = new Array;
        for (var h:int=0; h < size; h++) {
            var pair:Array = new Array(xValues[h], yValues[h]);
            correlates.push(pair);}
        return correlates;
    }

Answer 3

Чтобы получить корреляцию 1, X и Y должны быть одинаковыми, поэтому скопируйте X в Y, и вы получите корреляцию 1. Чтобы получить корреляцию -1, сделайте Y = 1 - X. (при условии, что значения X [0,1])

Answer 4

Этот, по-видимому, простой вопрос запутался в моей памяти со вчерашнего вечера! Я искал тему моделирования распределений с зависимостью, и лучшее, что я нашел, это: моделирование зависимых случайных величин . Суть в том, что вы можете легко смоделировать 2 нормали с заданной корреляцией, и они наметят метод преобразования этих независимых переменных, но это не сохранит корреляцию. Корреляция преобразования будет, так сказать, коррелированной, но не идентичной. Смотрите параграф «Ранговые коэффициенты корреляции».

Редактировать: из того, что я понял из второй части статьи, метод копулы позволил бы вам моделировать / генерировать случайные величины с ранговой корреляцией.

Answer 5

Вот реализация алгоритма twolfe18, написанного в Actionscript 3:

for (var j:int=0; j < size; j++) {
 xValues[i]=Math.random());
}
var varX:Number = Util.variance(xValues);
var varianceE:Number = 1/(r*varX) - varX;

for (var i:int=0; i < size; i++) {
 yValues[i] = xValues[i] + boxMuller(0, Math.sqrt(varianceE));
}

boxMuller - это просто метод, который генерирует случайный гауссиан с аргументами (mean, stdDev). size - размер дистрибутива.

Пример вывода

Target p: 0.8
Generated p: 0.04846346291280387
variance of x distribution: 0.0707786253165176
varianceE: 17.589920412141158

Как видите, мне еще далеко. Есть предложения?

Answer 6

начать с модели y = x + e, где e - ошибка (нормальная случайная величина). e должно иметь среднее значение 0 и дисперсию k.

Короче говоря, вы можете написать формулу для ожидаемого значения Пирсона в терминах k и решить для k. обратите внимание, что вы не можете случайным образом генерировать данные с Пирсоном, точно равным конкретному значению, только с ожидаемым Пирсоном определенного значения.

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

РЕДАКТИРОВАТЬ: хорошо, у меня есть волнистое решение, которое, вероятно, правильно (но для подтверждения потребуется тестирование). пока предположим, что желаемый Pearson = p > 0 (вы можете выяснить случай p < 0). Как я уже упоминал ранее, установите для вашей модели значение Y = X + E (X однородно, E нормально).

1. образец, чтобы получить ваши х
2. вычислить var (x)
3. дисперсия E должна быть: (1/(rsd(x)))^2 - var(x)
4. генерирует ваши y на основе ваших x и выборки из вашей обычной случайной величины E

для p < 0, установите Y = -X + E. действуйте соответственно.

В основном это следует из определения Пирсона: cov (x, y) / var (x) * var (y). когда вы добавляете шум к x (Y = X + E), ожидаемая ковариационная cov (x, y) не должна изменяться по сравнению с отсутствием шума. var (x) не изменяется. var (y) является суммой var (x) и var (e), отсюда и мое решение.

ВТОРОЕ РЕДАКТИРОВАНИЕ: хорошо, мне нужно лучше прочитать определения. определение Пирсона: cov (x, y) / (sd (x) sd (y)). Исходя из этого, я думаю, что истинное значение var (E) должно быть (1 / (r sd (x))) ^ 2 - var (x). посмотрите, работает ли это.