Алгоритм

Даний алгоритм належить до множин алгоритмів послідовностей на основі обчислення класів відрахувань. Властивості таких послідовностей описані в [1]. У розглянутому тут окремому випадку цілі числа формуються за рекурентною формулою

I_n+1=(J I_n+1)mod M, n=1, 2, …, M-1..... (A.1)

Потім послідовність цілих чисел нормується для одержання значень в інтервалі від 0 до 1. Нижче розглядається початкове значення I₀.

Очевидно, що період послідовності (А.1) не може бути більше М. Відповідно до [1] період дорівнює М, зокрема тоді, коли

J=4K+1 і M=2^L, (А. 2)

де К и L – такі цілі числа, що М>J. Найбільше значення I_n у (А. 1) дорівнює М-1, тому для формування по (А. 1) чисел, тільки менших 2³¹, необхідно, щоб

J(M-1)+1<2³¹, чи (4K+2³¹)/(4K+1)>2^L, (A. 3)

Таким чином, (А. 3) установлює співвідношення між J і M у (А. 1). Відзначимо, що при малому початковому значенні I₀ (А. 1) сегмент, що починається з нього, є послідовністю монотонно зростаючих чисел доти, поки J I_n+1 не стані більше М. Наприклад, у найпростішому випадку при J=1 і I₀=0 послідовність має вид 0, 1, 2, …, М-1. Таким чином, послідовність стає випадково. при великих значеннях J.

Нехай мінімальний період послідовності дорівнює 10⁶, тоді L=20 і з (А. 3) найбільше можливе значення K=511. Звідси для (А. 1) маємо

M=2²⁰=1048576, J=4(511)+1=2045 (A. 4)

Крім цього, виберемо (довільно) I₀=12357, (A. 5)

тоді перші 70 цілих чисел

Відзначимо, що всі цілі числа лежати в інтервалі від 0 до 1048576 і виявляється, що ця коротка послідовність, принаймні, не має небажаних властивостей. Оскільки в звичайній програмі побудова випадкових чисел числа формуються в інтервалі від 0 до 1, можна задати n-й відлік формулою

R_n=(I_n+1)/(M+1) (A. 6)

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

На основі співвідношень (А. 1), (А. 4) – (А. 6) можна написати підпрограму формування I₁, I₂, …, I_n..... При практичній реалізації цієї програми виникає проблема забезпечення автоматичної ініціалізації. Схема ініціалізації обов'язково залежить від типу ЕОМ. Однак для більшості ЕОМ працездатною є версія 1:

FUNCTION RANDOM (X)

IF (INIT.EQ.12357) GO TO 1

INIT= 12357

I=INIT

1 I=2045*I+1

I=I-(I/1048576)*1048576

RANDOM=FLOAT(I+1)/1048577.0

RETURN

END

Відзначимо, що під час завантаження даної програми внутрішня змінна INIT не винна дорівнювати 12357, а в інтервалах між викликами програми як INIT, так і I не повинні змінюватися. Початкова установа I=INIT забезпечує, як було описано вище, початок випадкової послідовності. У даному варіанті алгоритму Х – фіктивна змінна.

Якщо не потрібна автоматична ініціалізація чи вона винна бути змінною, те становить інтерес версія 2:

FUNCTION RANDOM (I)

I=2045*I+1

I=I-(I/1048576)*1048576

RANDOM=FLOAT(I+1)/1048577.0

RETURN

END

У цій версії І змінюється тільки при ініціалізації, коли встановлюється її початкове значення. Наприклад, для установи початкового значення I=12357 і формування перших двох випадкових чисел R1 і R2 можна використовувати припущення

K=12357

R1=RANDOM (K).

R2=RANDOM (K).

При звертанні до підпрограми RANDOM у версії 2 аргумент (К), природньо, не повинний змінюватися при її виконанні. Переваги іншого варіанта в порівнянні з першим полягають у тому, що:

1) немає необхідності приймати припущення про ті, що проміжні перемінні між звертаннями до підпрограми залишаються незмінними;

2) ініціалізацію програм можна здійснювати в основній програмі довільним числом в інтервалі від 0 до 1048575.

До недоліків відноситься те, що користувачу необхідно пам'ятати про правильну ініціалізацію I (чи еквівалента I в основній програмі) в інтервалі від 0 до 1048575 і про те, що між послідовними викликами програм I не повинно змінюватися.

У прикладі програми, що приводитися нижче, версії 2 розраховані перші 70 випадкових чисел відповідно з визначеними вище даними. Оскільки в даному прикладі початкове значення I=12357, у результаті отримані ті ж числа, що й у приведеному вище прикладі версії 1:

PROGRAM VERS2(INPUT, OUTPUT)

C – PRINT 70 RANDOM NUMBERS.

DIMENSION R(70)

K=12357

DO 1 J=1, 70

1 R (J) =RANDOM (K)

PRINT 2, R

2 FORMAT(7F8.6)

STOP

END

C

·FUNCTION RANDOM (I)

I=2045*I+1

I=I-(I/1048576) *1048576

·RANDOM=FLOAT(I+1)/1048577.0

RETURN

END

/LGO

.099414 .199705 .122502 .857597 .136420 .852693 .591484 .090877 .069850 .152838

.299419 .395030 .515310 .784634 .977466 .756376 .584924 .842606 .841235 .551257

.310731 .834445 .807589 .575788 .909412 .789466 .169059 .128435 .325905 .318828

.442812 .440026 .518642 .485529 .748327 .457776 .723296 .647902 .475414 .127320

.548520 .851800 .621811 .906695 .327288 .150114 .140660 .958640 .216346 .368044

.725532 .930929 .603319 .192025 .303295 .981302 .648336 .419003 .425186 .647826

.712078 .750182 .758750 .689554 .236604 .763126 .846010 .860181 .503743 .802619

Відзначимо ще раз, що в іншому варіанті незалежно від початкового значення I період випадкової послідовності дорівнює 1048576.

Щодо інших датчиків випадкових величин, те стандартним методом моделювання непрерервної випадкової величини x з функцією розподілу F(x), коли існує зворотна до неї F^-1(x), є використання алгоритму виду

x=F^-1(a). (1.2.1)

Цей алгоритм можна використовувати в тих випадках, коли існує аналітичне вираз для F^-1(x) або в математичному забезпеченні є стандартна процедура, що обчислює цю функцію.

Для дискретної випадкової величини x із законом розподілу p_k=P{x=x_k}, k=0,1, 2,…,універсальний алгоритм реалізує <<метод віднімання>> (рис. 1.8). Реалізація цього алгоритму припускає наявність у пам'яті ЕОМ усього набору чисел p_k, k=0,1, 2,…

Рис.1.8. Алгоритм віднімання

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

Завдання для виконання

1.Нормальний розподіл

p(x)=(2p)^1/2exp(-x ^2/2), xÎ(-¥, ¥). (1.2.2)

Алгоритм А1.

0. Зарезервовано константу с =2p.

1. R=

2. j=ca₂.

3. U1=R cosj, U2=R sinj.

Дата добавления: 2015-08-13