Метод Санделиуса для отримання випадкових перестановок

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

Однак у додатках поряд з послідовностями випадкових і псевдовипадкових чисел потрібно отримувати перестановки чисел, що мають рівномірний розподіл. Наприклад, потреба в таких перестановки періодично з'являється в криптографічних додатках.

Метод описаний нижче запропоновано Санделиусом (М. Sandelius) ще в 1962 р. у роботі [1].



Описаний алгоритм дозволяє отримувати випадкові перестановки з n елементів, які мають рівномірний розподіл. Останнє означає, що ймовірність отримати одну з n! можливих підстановок, використовуючи даний метод дорівнює 1/n!..

Алгоритм отримання перестановок

Метод Санделиуса простіше описати рекурсивно.

На кожному кроці обробляється масив P. Для масиву P генеруються випадкові біти в кількості, що дорівнює числу елементів у масиві P. i-ий біт послідовності зіставляється i-того елементу масиву Масив P. P ділиться на два масиву P0 і P1 за принципом: всі елементи, яким зіставлені нулі заносяться в масив P0, інші в масив P1. Для кожного масиву P0 і P1 виконується перемішування (рекурсія). Перемішані масиви об'єднуються в один. Спочатку йдуть елементи з масиву P0, потім з масиву P1.

Процедура перемішування Sandelius(P):
1. n := |P|; - потужність (кількість елементів) P
2. Якщо n=1, то повернути P;
3. g:=[gi, i=1..n]; – отримання послідовності випадкових біт gi
4. P0:=[]; P1:=[]; - порожні масиви
5. i0:=0; i1:=0; - індекси для запису в масиви
6. k:= 0;
7. Якщо g[k]=0, то
7.1. P0[i0] := P[k];
7.2. i0 := i0+1;
8. Якщо g[k]=1, то
8.1. P1[i1] := P[k];
8.2. i1 := i1+1;
9. k := k+1;
10. Якщо k < n, то переходимо до кроку 7.
11. P0 := Sandelius(P0); 
12. P1 := Sandelius(P1);
13. P := P0||P1 – об'єднуємо масиви в один великий масив.
14. Повертаємо P.

Перестановка за методом Санделиуса виходить з послідовності чисел від 1 до n:
P=[1,2, ...,n];
Sandelius(P);


Алгоритм, досить легко програмується.
Ось так, наприклад, можна запрограмувати в Maple
Sandelius:=proc(p)
local A,m,i,p1,p2;
m:=nops(p);
A:=[seq(getNextRndBit(), i=1..m)];
p1:=[]; p2:=[];
for i from 1 to m do
if A[i]=0 then p1:=[m(p1),p[i]];
else p2:=[m(p2),p[i]];
fi;
od;
if nops(p1)>1 then p1:=Sandelius(p1); fi;
if nops(p2)>1 then p2:=Sandelius(p2); fi;
return [op(p1),op(p2)]; 
end proc:



Ось приклад моєї реалізації алгоритму на C++, який я використовував в одному своєму дослідженні
unsigned __int8 *bits,*tmp_perm;
Sandelius(unsigned __int8 *perm,int n)
{
tmp_perm = (unsigned __int8 *)malloc(n);
bits = (unsigned __int8 *)malloc(n);
for(int i=0;i < n;i++) perm[i]=i;
recursiveSandelius(perm,n);
free(bits);
free(tmp_perm);
}

recursiveSandelius(unsigned __int8 *perm,int n)
{
if (n < =1) return;
for(int i=0;i < n;i++)
bits[i]=getNextRndBit();
k=0;
for(int i=0;i < n;i++)
if(!bits[i])
tmp_perm[k++]=perm[i];
zeros=k;
for(int i=0;i < n;i++)
if(bits[i])
tmp_perm[k++]=perm[i];
memcpy(perm,tmp_perm,n);
recursiveSandelius (perm,zeros);
recursiveSandelius (&perm[zeros],n-zeros);
}



Особливості

Сподіваюся, алгоритм генерації перестановки описаний мною досить зрозуміло. Тепер хочеться трохи обговорити, особливості його роботи.

Для роботи алгоритм вимагає послідовність випадкових біт. Головна вимога до цієї послідовності – біти повинні бути незалежними. В цьому випадку алгоритм генерує рівномірно розподілені перестановки навіть.

Слід враховувати, що випадкові біти можуть мати нерівномірне розподіл, але вони повинні бути незалежними, інакше розподіл підстановки не буде рівномірним.

До недоліків можна віднести той факт, що число випадкових біт, які потрібні для генерації перестановки не визначено заздалегідь.

Практична перевірка

Відразу скажу, що факти зазначені в попередньому пункті доводяться аналітично, але все одно захотілося перевірити їх на практиці.

Для перевірки рівномірності розподілу одержуваних підстановки написана проста програма в математичному пакеті Maple. У програмі я генерував велике число перестановок, підраховував кількість підстановок кожного виду (щось на зразок гістограми). Для отриманого масиву перевірялася гіпотеза про рівномірності за критерієм Пірсона. Крім того, підраховували розподіл кількості біт, необхідних для генерації одного підстановки.

Приводити вихідний текст програми тут не бачу сенсу, але якщо є бажання порахувати самим, то файли з вихідними кодами можна знайти тут.

Досліджувалися перестановки довжини n=7. Генерувалося N=n!*1000 перестановок. Випадкові біти генерувалися так: 0 з ймовірністю 0,5+d, 1 з ймовірністю 0,5-d. d дорівнює 0 для рівномірно розподілених біт. Для отримання залежного біта генерувався випадковий біт і складався з попереднім бітом.

Число n=7 взято з міркувань розумного часу виконання (у мене це 10-20 хвилин).

Результати моделювання:
Варіант Рівноймовірні підстановки (критерій Пірсона)? Гістограма Середнє число випадкових біт/середньоквадратичне відхилення
Незалежні біти, d=0 Так Незалежні, d=0 28.24/28.26
Незалежні біти, d=0.05 Так Незалежні, d=0.05 28.50/28,54
Незалежні біти, d=0.4 Так Незалежні, d=0.4 71.47/74.77
Залежні біти, d=0.05 Немає Залежні, d=0.05 30.15/30.32


На малюнках по горизонталі зазначені номери перестановок, а по вертикалі частоти їх появи. Не озброєним оком видно, що в останньому випадку розподіл сильно відрізняється від рівномірного.

Так само видно, що не залежно від перекосу ймовірності в незалежних випадкових біт алгоритм генерує рівномірно розподілені перестановки. Однак чим більше перекіс у ймовірності, тим більше випадкових біт вимагає алгоритм.

Якщо ж біти мають залежність, то генеруються підстановки мають розподіл відмінне від рівномірного розподілу.

Кількість підстановок, кожного типу повинно мати розподіл близький до нормального з середнім N/n! і дисперсією N/n!(1-1/n!).

У перших трьох випадках гістограма виглядає приблизно так:

Гістограма (випадок 2)

В останньому випадку видно, розподіл далеко від очікуваного:

Гістограма (4 випадок)

Література

1. М. Sandelius. A Simple Randomisation Procedure, J. of the Royal Statistical Society. Ser. Ст., V. 24, № 2, 1962.

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.