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

Вот реализация 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). размер - это размер распределения.

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

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

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

start with the model y = x + e where e is the error (a normal random variable). e should have a mean of 0 and variance k.

long story short, you can write a formula for the expected value of the Pearson in terms of k, and solve for k. note, you cannot randomly generate data with the Pearson exactly equal to a specific value, only with the expected Pearson of a specific value.

i'll try to come back and edit this post to include a closed form solution when i have access to some paper.

EDIT: ok, i have a hand-wavy solution that is probably correct (but will require testing to confirm). for now, assume desired Pearson = p > 0 (you can figure out the p < 0 case). like i mentioned earlier, set your model for Y = X + E (X is uniform, E is normal).

1. sample to get your x's
2. compute var(x)
3. the variance of E should be: (1/(rsd(x)))^2 - var(x)
4. generate your y's based on your x's and sample from your normal random variable E

for p < 0, set Y = -X + E. proceed accordingly.

basically, this follows from the definition of Pearson: cov(x,y)/var(x)*var(y). when you add noise to the x's (Y = X + E), the expected covariance cov(x,y) should not change from that with no noise. the var(x) does not change. the var(y) is the sum of var(x) and var(e), hence my solution.

SECOND EDIT: ok, i need to read definitions better. the definition of Pearson is cov(x, y)/(sd(x)sd(y)). from that, i think the true value of var(E) should be (1/(rsd(x)))^2 - var(x). see if that works.

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

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

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

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

-Сгенерировать массив X

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

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

-Измерьте начальный уровень 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

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

Вот моя реализация в 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;
    }