CC BY
89
УДК 519.2
Вестник СПбГУ. Сер. 1, 2006, вып. 2
Т. М. Товстик
СРАВНЕНИЕ НЕКОТОРЫХ СТАТИСТИЧЕСКИХ СВОЙСТВ КВАЗИСЛУЧАЙНЫХ И ПСЕВДО-СЛУЧАЙНЫХ ПОСЛЕДОВАТЕЛЬНОСТЕЙ*
Исследуются статистические свойства квазислучайных последовательностей Холтона и Варнока, а также квазислучайной и псевдослучайной последовательностей, предложенных автором. Сравниваются корреляционные свойства всех четырех последовательностей, в частности, на примере моделирования стационарных процессов.
1. Введение. Широко известны квазислучайные последовательности LnT [1], Холтона [2, 3], Варнока [4]. Это детерминированные зависимые последовательности одномерные или многомерные с компонентами в диапазоне Kn = [0,1]n. Корреляционные свойства рядов LnT подробно рассмотрены в [5, 6], поэтому в данной статье остановимся на рядах Холтона и Варнока. Кроме этих последовательностей рассматриваются также две другие последовательности, одна из которых является квазислучайной, а другая — псевдослучайной. Алгоритмы используют идею Холтона о представлении чисел в системах, отличных от десятичной. Сравниваются корреляционные свойства всех четырех последовательностей, в частности, на примере моделирования стационарных процессов.
2. Алгоритм Холтона. Каждое простое число m > 2 может создать со значениями в области K = [0,1] одномерную последовательность Холтона [2]
hm(i), i = No,No + 1,..., 0 < hm(i) < 1. (2.1)
Здесь целое No > 0. В частности, при No = 0 последовательность начинается с точки hm(0) = 0. Для получения hm(i) используется следующий алгоритм: если в m-ичной системе счисления целое число i записывается в виде
i = ajaj—i . . . ai, 0 < ak < m — 1, k < j, (2.2)
то соответствующее число Холтана равно
hm(i)=0,aia2 . ..aj. (2.3)
Связь чисел i и hm(i) в десятичной и m-ичной системах такова:
jj i = 53 akmk—1, hm(i) = ^2 akm~k. (2.4)
k = 1 k=1
Для получения чисел hm(i) c точностью s знаков после запятой достаточно найти ak при 1 < k < v = s log 10/ log m + 2.
* Работа выполнена при финансовой поддержке РФФИ (грант №05.01.00865-а).
© Т.М.Товстик, 2006
При фиксированном т и N = 1 числа /т(г), г = 1, 2,..., сначала возрастают и достигают максимума при г = т — 1. Далее, при г = т идет спад и монотонный подъем до очередного максимума с длиной серии равной т.
Например, при т = 199 имеем /199(1) = 0.00503, /199(120) = 0.60302, /199(198) =
0.99497, /¿199(199) = 0.00003. Если промоделировать 120 точек, то в интервал (0.603, 1) не попадет ни одной точки.
Для заполнения интервала [0,1] число моделируемых точек N должно быть не менее т. Свойства последовательности лучше проявляются, если N > т2.
Последовательность /2(г), г = 1, 2,..., сильно коррелирована. При моделировании ее в пределах 1 < г < 3000 коэффициенты корреляций удовлетворяли неравенству \гк \ > 0.55 при всех к. Последовательные корреляции постоянно меняют знак.
о к 30 200 к 230
Рис. 1. Коэффициенты корреляций ряда h 2 (i).
На рис. 1 приведены коэффициенты корреляций ряда h2(i), i = 1, 2,..., 1000 в пределах 0 < к < 30 и 200 < к < 230.
Коэффициенты корреляций Tk рядов hm(i), m < 2, имеют максимумы, близкие к единице, через каждые m значений.
Хотелось бы внести некоторую «случайность» в строгую закономерность последовательности Холтона, особенно это важно (как указывалось выше) при больших m и небольших объемах выборки. Эта цель достигается в алгоритме Варнока.
3. Алгоритм Варнока. Варнок [3] предложил использовать алгоритм Холтона, но образовывать последовательность не от подряд идущих величин i, как указано в (2.1), а в виде
wm (i) = hm (Yi ), m > 1, Yi = i Л(т), i = Nq,Nq + 1,... (3.1)
Здесь m — простое число, Л(2) = 1, Л(3) = 2, а при m > 5 числа A(m) вычисляются по следующей схеме.
Обозначим через [х] и {ж} целую и дробную части х. Пусть Y = А = [У],
B = A + 1, тогда A < Y < B. Представим A/m и B/m в виде конечных цепных дробей
A 1 B 1
m п , 1 ’ m . 1 '
h + ---------------- 9i + —-------------------
f2 + ... 1 92 + ... 1
fk 9j
Если выполнено хотя бы одно из трех условий
к 1=1 3 <£; 1=1
к 1=1 3 = 1=1
3 =3 1=1
О = < я а е
где г и = шах1<<к В частности, Л(5) = 2, Л(7) = 5, Л(17) = 3.
4. Квази-случайная последовательность. Рассмотрим также более простую последовательность со свойствами близкими к тем, которыми обладает последовательность Варнока:
гт(г) = 1гт(7*), то > 1, 7* = *[а/то], г = М0, М0 + 1, .. . (4.1)
5. Корреляционные свойства квази-случайных последовательностей. Под
£т(г) будем подразумевать одну из последовательностей кт(г), и)т(г) или ¿т(г). При моделировании последовательностей £т (г) будем считать, что индекс г пробегает подряд идущие номера, иначе среднее и дисперсия ряда не будут такими, как у равномерного распределения. Если номер г нам не важен, то будем его опускать.
Заметим, что
Л-2 = №2 = ¿2 = £2, Л-з = ¿з, №5 = ¿5, №11 = ¿ц. (5.1)
Моделируем выборку £т(1), £т(2),..., £тN) объема N, находим выборочные ковариации
1 N-к 1 N
^к = N -к + ^ ~г))(^(Я {'=—^2£т{]) (5.2)
3=1 3=1
и коэффициенты корреляций
Гк = &к/%. (5.3)
Пусть V = тп,п > 1,« > 13 и объем N выборки достаточно велик, тогда коэффи-
циенты корреляций Гк обладают следующими свойствами:
г к « 1, к = «з, 3 > 0, (5.4)
Гу 3+г ~ гу ж, 3, I > 0, (5.5)
Гу 3+1 « Гу ь-1, 3,1 > 0. (5.6)
Для подтверждения справедливости соотношений (5.5) приведем максимумы расхождений между указанными там коэффициентами корреляций (аналогичные результаты имеют место между коэффициентами корреляций (5.6)). Для процессов Нт(г), №т(г) и ¿т(г) в столбцах IV-VI таблицы 1, соответственно, приводятся величины
шах шах \Гуз+^ - Гу + \ . (5.7)
0<г<у—1 0<3,1<3
3
В таблице 2 столбцы IV-VI объединены ввиду (5.1).
Таблица 1
т V N IV V VI
49 5000 0.00158 0.00071 0.00058
15000 0.00023 0.00055 0.00035
61 61 5000 0.00196 0.00033 0.00187
15000 0.00066 0.00031 0.00031
103 103 5000 0.00210 0.00164 0.00094
15000 0.00050 0.00039 0.00036
Таблица 2
т v{m) N IV-VI т v{m) N IV-VI
2 32 1000 0.00255 2 64 5000 0.00061
2 32 5000 0.00165 2 64 15000 0.00039
Замечание. Приближенные соотношения (5.4)-(5.6) для выборочных корреляций не влекут соответствующих равенств для рк — теоретических коэффициентов корреляций. Действительно, тогда для случайных величин
Ст(к) = ^¿£т(0-1)г;-Н)
при V = тп дисперсии совпадали бы с ковариациями, однако для соответствующих выборочных величин последнее совпадение не имеет места.
6. Моделирование последовательности со свойствами псевдо-случайных чисел. Образуем последовательность вт (г) с помощью алгоритма
%(*) = Ьт(тг), п = [г'/1Л}, г = И0, N0 + 1,. •. (6.1)
где Ь > 1 —целое число, Нт(т) —числа Холтона. Последовательность вт(г) обладает очень слабой корреляционной зависимостью, т. е. имеет свойства псевдослучайных чисел.
1
rk О
-1
Рис. 2. Выборочные коэффициенты корреляций рядов (a — h\7, b — W17, c — £17, d — S17).
На рис. 2 приведены выборочные коэффициенты корреляций rk (к = 0,1,..., 50) рядов h17(i), w17(i) (Л(17) = 3), t17(i), s17(i) при N = 1000. У квазислучайных последовательностей (рис. 2 (a)-(c)) ярко выражен период длины m = 17.
7. Проверка выборки на нормальность и равномерность. С помощью критерия х2 проверялась принадлежность выборок каждой из последовательностей к двум законам распределения: нормальному и равномерному. Проверка производилась при
2 < т < 107, N = 4000 и N = 15000. Если не делать каких-либо группировок, то все последовательности оказывались равномерно распределенными. К нормальному распределению может приводить группировка.
Введем обозначение
1
ст(к) = ^ У2^({к-1)т + г) к=1,...,Ыт, ЛГт = [ЛГ/т].
Случай
п > 1, т > 15
(7.1)
(7.2)
является особым, так как при его выполнении ряд (7.1), как правило, подчиняется равномерному распределению. При т = 2 для т = 16, 32, 64 процесс ^т имеет равномерное распределение. При т = 3 для т = 27 и = кт процесс £т имел равномерное распределение, а при процесс £т не согласовался ни с нормальным, ни с
равномерным распределением.
Если соотношение (7.2) не выполнено, то процесс Ст(к), в основе которого лежат ряды и Ьт, чаще приводит к нормальному распределению.
8. Зависимость дискрепанса от объема выборки и параметра т. Равномерность выборочного ряда X1,...,XN проверяется дискрепансом
ВN = 8Ир
еК
1 *
1-х- -х)
к = 1
который в одномерном случае можно вычислить по формуле
В* = тах тах 1< к<Ы
1 - хк —
N - к + 1
N
1 - хк —
N — к
N
(8.1)
(8.2)
Здесь г* <,...,< г* — вариационный ряд исходного. В данном случае дискрепанс является максимальным уклонением выборочной функции распределения от равномерной. Большие значения В* свидетельствуют о кучности величин XI в некоторой области.
Включим в вариационный ряд две точки х0 = 0,х* +* = 1 и у полученного ряда найдем максимальное и минимальное расстояния между соседними точками
Ртах^) = тах (х*+1 - X*i), рmin(N) = тт (х*+* - х*).
0<г<N +1 0<г<N +1
Будем моделировать последовательность Холтона, начиная с N0 = 1 при фиксированном т. Если объем выборки равен целой степени т, т. е. N = тп, то
ВN = Ртах^) = 1^ = 1/тП, ртш^) = 1/тП+1 = Вт»+1,
(8.3)
т. е. максимальное расстояние между соседними точками совпадает с дискрепансом выборки, а минимальное — с дискрепансом при N = тп+1. При N = тп дискрепансы последовательностей ,ш(г) и Ь(г) также равны 1/N = 1/тп, а рmax(N) и ртт^), как
правило, не совпадают со значениями у Холтона.
п
т
Рис. 3. Величины дискрепансов трех квази-случайных рядов (а — Нш, Ь — , с — ).
На рис. 3 приводятся величины дискрепансов трех квазислучайных рядов, полученных при N = 3500, 2 < т <, 541. На графике 3а масштаб в 10 раз меньше, чем на двух других графиках. На рис. 3 обозначены линии 1п N N (такой аситмптотический порядок при N имеет дискрепанс [2]).
За исключением вариантов (8.3), дискрепанс у ряда Холтона, как правило, больше, чем у двух других рядов. Особенно это свойство проявляется при N < т2, и тем более, при N < т. Как видно из рис. 3, дискрепанс у ряда Варнока, как правило, меньше, чем у двух других рядов.
9. Многомерные квазислучайные последовательности. В п-мерном единичном кубе Кп с помощью различных простых чисел т\, т2,..., тп можно генерировать последовательность Холтона размерности п
Н(п) (г) = (Нт1 (г),...,Нтп (г)), г = 1, 2,..., 0 < Нт (г) < 1. (9.1)
Последовательности Нтк (г) и Нтз (г) при любых к = 5 статистически независимы.
Аналогично образуются многомерные точки, в основе которых лежат последовательности (3.1), (4.1) и (6.1). Они будут обозначаться W(n) (г), Т(п) (г) и Б(п) (г). Если нанести на график пары Н(2)(г), то проявляется внутренняя зависимость каждой компоненты, которая особенно при больших т выражается в регулярных полосах (см. рис. 4(1)).
При т\ = 307, т2 = 643 были промоделированы 3500 пар точек всех типов, начиная с N0 = 1001. У последовательности (3.1) Л(307) = 161, Л(643) = 230.
На рис. 4 приведены соответственно пары точек (Н307(г),Ь643(г)), ^307(г),гш643(г)), (^307(г), ¿643(г)) и (в3от(г),вб43(г)).
Таблица 3
1711 ГП2 £>(Л) -О(ги) т £>(*)
1. 307 0.0150 0.0012 0.0019 0.0068
2. 643 0.0453 0.0010 0.0013 0.0068
3. 307 643 0.0669 0.0031 0.0027 0.0104
4. 8 0.0034 0.0010 0.0007 0.0187
5. 9 0.0038 0.0009 0.0011 0.0186
6. 8 9 0.0064 0.0036 0.0035 0.0187
В таблице 3 выписаны дискрепансы всех четырех рядов, причем в третьей и шестой строках для двумерных рядов, а в остальных — для одномерных.
10. Моделирование нормальных случайных величин. Наиболее употребительные формулы, используемые для формирования стандартных нормальных слу-
Рис. 4. Расположение точек на плоскости (1- H(2) (i), 2- W(2) (i), 3- T(2) (i), 4- ^(2) (i)).
чайных величин, это: _____________
is (к) = \/—2. log(ai(fc)) cos 2тга2 (/с), (10.1)
12
n(к “j (к) - 6 (10-2)
j=1
В них при всех i величины “i — независимые равномерно распределенные на (0,1). Были промоделированы последовательности с элементами (10.1) и (10.2), в которых “i заменялись на hm (i), wm (ji), tm (ji), sm (ji). В дальнейшем под “i будем подразумевать один из этих вариантов.
Если пользоваться формулой (10.1), то “1 (k) и “2(к) следует вычислять по разным простым числам m1 и m2. Получающиеся последовательности распределены нормально, они значимо коррелированы, если образуются с помощью квазислучайных последовательностей, и статистически независимы, если основаны на рядах sm.
Аналогично, в формуле (10.2) при фиксированном к величины “j (к) при всех j следует моделировать с помощью разных (например, идущих подряд) простых чисел, начиная с m0. Наиболее слабую статистическую зависимость для ц(к) дает алгоритм Варнока, а алгоритм на основе чисел Холтона при небольших m0 имеет значимый первый коэффициент корреляции, а при m0 > 23 для всех к < 4, коэффициент корреляции гк > 0.55.
На рис. 5 приведены коэффициенты корреляции рядов (10.2) при N = 1000, m0 = 2. При других mo характер кривых на рис. 5 меняется.
11. Моделирование гауссовского стационарных процессов. Рассмотрим моделирование стационарного процесса авторегрессии с остатками скользящего среднего АРСС(3,2) [6]
3 2
С(t) = — ^ aiС(t — i) + ^ Ькц(t — к). (11.1)
i=1 к=0
Рис. 5. Коэффициенты корреляции рядов (10.2) (а Нт 5 Ь Мт 5 с tm 5 & 5т).
Пусть коэффициенты в (11.1) равны ах = —2.483, а2 = 2.299, а3 = —0.795, Ь0 = 0.467, Ьх = —0.643, Ь2 = 0.291.
Рис. 6. Теоретические и выборочные коэффициенты корреляции (а Нт5 Ь Мт5 с tm5 & 5т).
Если в (11.1) случайные величины п (к) являются независимыми стандартными нормальными, то С(к) —гауссовский стационарный процесс с ¡1 = ЕС(к) = 0. Случайные числа п(к) получались с помощью формулы (10.2), но в ней кроме ряда Нт подставлялись также тт, Ьт и вт. Во всех случаях т0 = 2, N0 = 1, N = 1000.
На рис. 6 для 0 < к < 50 приведены теоретические (рь) и выборочные (гь) коэффициенты корреляций процесса С (к), полученные, соответственно, на основе каждого из рядов Нт, и!т кт и вт. Налицо расхождение между теоретическими и выборочными корреляциями. При теоретическом значении а = 1.02 выборочные среднеквадратические отклонения таковы: ан = 1.34, = 0.73.
Следовательно, прохождение через фильтр (11.1) квазислучайных последовательностей приводит к искажению среднеквадратических отклонений (<г) и коэффициентов корреляций.
Автор благодарит С. М. Ермакова за постановку задачи и обсуждение результатов. Summary
T. M. Tovstik. Comparison of some statistical properties of quasi-random and pseudo-random sequences.
Some statistical properties of quasi-random Halton and Warnock sequences and statistical properties of quasi-random and pseudo-random sequences offered by the author are investigated. The correlation properties of these four sequences are compared. As an example the simulation of stationary processes based on these sequences is studied.
Литература
1. Соболь И. М. Многомерные квадртурные формулы и функции Хаара. М.: Наука. 1969. 288 с.
2. Halton J. H. On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals // Numer. Math. 1960. Vol. 2. P. 84-90.
3. Halton J. H. Quasi-probability // Monte Carlo Methods and Appl. 2005. Vol. 11. N 3. P. 203350.
4. Warnock T. T. Effective error estlmates for quasi-Monte-Corlo computations // Int. Conf. on Monte Carlo and Quasi-Monte-Carlo Methods. 1998. Claremont. CA.
5. Ермаков С. М., Товстик Т. М. Квазислучайные последовательности в алгоритмах моделирования случайных процессов // Вестник Санкт-Петерб. ун-та. 1999. Сер. 1. Вып. 4. C. 26-33.
6. Ermakov S. M., Tovstik T. M. On the quasi-random sequences in the random processes modeling algorithms // J. Applied Statist. Science (JASS), 2002. V. 11. N1. P. 45-56.
Статья поступила в редакцию 14 февраля 2006 г.
УДК 519.246.8+519.254
Вестник СПбГУ. Сер. 1, 2006, вып. 2
Ф. И. Александров
ВЫДЕЛЕНИЕ АДДИТИВНЫХ КОМПОНЕНТ ВРЕМЕННОГО РЯДА ПРИ ПАКЕТНОЙ ОБРАБОТКЕ МЕТОДОМ «ГУСЕНИЦА»-SSA
1. Постановка задачи, история, актуальность
Временным рядом называется последовательность N наблюдений некоторого процесса Ен = (/о, /1,..., /н_х), /п € К, N > 2. Важной задачей является задача выделения аддитивной составляющей ряда, когда требуется найти (или аппроксимировать) Е^ по наблюдаемому ряду Ен = Е^ + Е^ (поэлементно). Существует много вариантов данной задачи: выделение тренда, сезонности, сигнала, — а также подходов для их решения: параметрические и авторегрессионные модели, линейные фильтры, преобразование Фурье, вейвлеты. Среди них можно выделить группу методов, основанных на представлении ряда в некотором базисе. В использованном нами методе «Гусеница»-ЭЭА базис порождается самой структурой ряда, что делает метод гибким и позволяет обрабатывать широкий спектр рядов.
Метод «Гусеница»-88А [1, 2] зародился в 1970-1980-х годах для поиска аттракторов динамической системы, улучшения соотношения сигнал/шум, идентификации сигнала в шуме. В его основе лежит преобразование ряда в матрицу и её сингулярное разложение, используемое в анализе главных компонент. Данный метод позволяет выделять различные аддитивные составляющие ряда, причём выбор нужных для выделения сингулярных векторов делается визуально.
К настоящему времени метод «Гусеница»-ЭЭА хорошо себя зарекомендовал для решения задач анализа и прогноза в различных областях: геофизике, метеорологии, эконометрике, социальных науках (см. [1, 2]). В связи с обработкой значительных объемов данных возникла проблема автоматизации визуальной части метода. В работе [3] были предложены методы автоматического выбора собственных векторов при выделении тренда и периодической составляющей, но не были изучены вопросы выбора их параметров, а также методики их применения. Нами были разработаны подходы к выбору параметров, предложен более устойчивый метод для выделения тренда, а также выработана методика для применения этих методов при обработке серии рядов.
Подобная постановка задачи, в частности, была мотивирована необходимостью применения метода к данными по экспрессии гена, которые представляют собой естественные классы (задаваемые типом гена и возрастом эмбриона), каждый из которых содержит множество незначительно отличающихся рядов (соответствующих конкретному эмбриону).
2. Базовый алгоритм метода «Гусеница»-SSA
Опишем вкратце алгоритм метода «Гусеница»-88А для выделения аддитивной составляющей ЕД1) ряда Ен = Е(1) + Е^ (подробное описание см. в [1]).
На первом этапе алгоритма мы выбираем параметр метода — длину окна Ь, 1 < Ь < N, и строим траекторную матрицу ряда X € Кь хК , К = N — Ь +1, столбцами которой являются вектора X^ = (/j_l,..., /j+L_2)T,] = 1,...,К. Затем вычисляется матрица ХХТ, её собственные числа {А^}^=1, Аj ^ Аj+l, собственные вектора {Uj}^=х, соответ-
© Ф. И. Александров, 2006
