Вычисление дискрепанса конечного числа точек в n-мерном единичном Кубе Текст научной статьи по специальности «Математика»

УДК 519.6 Вестник СПбГУ. Сер. 1, 2007, вып. 3

Т. М. Товстик

ВЫЧИСЛЕНИЕ ДИСКРЕПАНСА КОНЕЧНОГО ЧИСЛА ТОЧЕК В п-МЕРНОМ ЕДИНИЧНОМ КУБЕ*

1. Введение. При обращении к методам Монте-Карло часто возникает необходимость подсчета дискрепанса конечного числа точек, который служит мерой равномерности их распределения. Например, такая задача может возникнуть при проверке свойств псевдо- и квазислучайных последовательностей, при оценке погрешности вычисления интегралов и т.д. [1, 2].

В данной работе приводится алгоритм вычисления дискрепанса конечного числа точек в двумерном и п-мерном единичном кубе [0,1)”. Сравниваются некоторые свойства и время подсчета дискрепансов с помощью данного алгорита и алгоритмов, предложенных П. Бундшук и Я. Жу [3, 4]. При 2 < п < 4 время вычисления дискрепанса по алгоритму, предложенному автором, значительно меньше. При больших п выбор более быстродействующего алгоритма зависит от количества проверяемых точек.

Пусть п — размерность компонент ряда Х1,..., :

,(1)

>)^ є кп = [0,1)п, 3 =

N.

(1)

Равномерность ряда (1) проверяется дискрепансом

В

у,п = вир ХеКп

^1 — х(1) ^ • • • ^1 — х(п)^

1

N

3=1

(п)

••• X (х3 - х

(п)

. (2)

Здесь X —произвольная точка в К” (X = (х(1),..., х(п)) € К”). Такое выражение для дискрепанса впервые использовал Бергстром [5]. Однако ранее Колмогоров для проверки равномерности распределения одномерных случайных величин [6] вывел функцию распределения правой части (2).

Будем предполагать, что последовательность ранжирована по первой переменной

(3)

х11) < • • • < х((1)

^1 ^ xN .

Ранжируем остальные компоненты и обозначим их, соответственно

к = 2,..., п.

(к)* ^ ^ (к)* х1 < • • • < ху ,

(4)

Для общности записи будем использовать обозначение х^1)* = х31). При каждом к =

2,..., п и з = 1,..., N обозначим через пк (з) порядковый номер х^к)* в ряду (1).

Величина дискрепанса не изменится, если к ряду (1) добавить точки Хо = (0,..., 0) и Ху +1 = (1,..., 1), при этом Введем обозначение

(к)*

= 0, х^к+*1 = 1, к = 1, .. ., п.

V

N

х(1)*

31

х(п)* 3п ,

х(1)*

х3і + 1 •

х(п)*

3п+1

* Работа выполнена при финансовой поддержке РФФИ (грант №05-01-00865-а). © Т.М.Товстик, 2007

3

тогда дискрепанс можно вычислить по формуле

DN,n = S(ji, j2, • • • , jn; )• (6)

c1<jk <N

1 < k<n

Здесь vj1 ... jn/N — величина выборочной функции распределения в n-мерном кубике,

/ (1)* (n)*\

минимальная и максимальная координаты которого, соответственно (xji , • • •, xjn ),

(xj1+*i, • • •, xjn+*i). Заметим, что при вычислении дискрепанса по формуле (6), не обязательно подсчитывать функцию S во всех перечисленных точках.

Укажем алгоритм выбора нужных точек сначала в двухмерном, а затем n-мерном случае.

2. Двухмерный случай. При n = 2 обозначим выборку

(х^ш),^ (xN ,VN), x1 <•••< xN, (7)

ранжируем ряд вторых компонент

y* <•••< yN, (8)

и пусть при каждом j = 1, • • • N, n2(j) —номер y* в ряду (7). Соотношение (6) примет вид

Dn,2 = max S(i,j,; ), (9)

1<i,j<N

S(i, j; v) = max(v/N - xjyj, xi+iyj+i - v/N)• (10)

Укажем, при каких i и j нужно вычислять S(i, j; v) и как найти v^

Пусть ¿2 = max(x1, y*) При каждом j (j = 1, • • •, N) счетчик v суммирует число выполнений неравенств

П2 (i) < j (11)

при изменении i от 1 до N, а счетчик ц — число выполнений обратных неравенств между соседними скачками v, т. е. изменениями v на единицу. Если при данном j и очередном i выполнено неравенство (11), то величина v увеличивается на единицу и вычисляются si = S(i,j; v), ¿2 = max(d2,si). Если в этот момент ^ > 1, то вычисляются также S2 = S(i — 1,j; v — 1), ¿2 = max(d2,S2). Указанные операции выполняются при всех j = 1, • • •, N и i = 1, • • •, N, и после их завершения получаем дискрепанс Dn,2 = Благодаря использованию счетчика ^ удается сократить число операций.

3. Общий случай. Пусть размерность компонент Xj равна n и пусть nk(j) —номер xjk) * в ряду (1). Полагаем dn = max(xi1 ), • • •, xin )). Для любого набора j2, • • •, jn

(1 < j2 , • • • , jn < N) счетчик v при изменении j1 от единицы до N накапливает число одновременно выполненных неравенств

П2 (j 1) < j2, Пз(ji) < j3, • • • П2(j 1) < jn• (12)

Между каждыми двумя последовательными выполнениями неравенств (12) счетчик ^ суммирует число последовательных проверок, в которых не выполнены эти неравенства. Если при очередном j1 неравенства (12) справедливы, то к счетчику v добавляется единица и с этим v вычисляется si = S(ji, j2, • • •, jn; v) и dn = max(dn, si). Если в

этот момент ^ > 1, то находим в2 = S(31 — 1,32, • • •, 3п; V — 1) и ¿п = шах(^„, в2). При завершении вычислений по всем переменным величина ¿п будет равна дискрепансу:

DN,n = ^п-

4. Примеры. Для сравнения алгоритмов фиксировалось время вычисления дис-крепансов при размерностях вектора п, 2 < п < 5 и различных объемах выборки N. (Разумеется, время счета зависит от быстродействия РС, поэтому имеет смысл лишь отношение приводимых ниже величин. В таблице приведено время в секундах для машины, используемой автором.) Программы были написаны на языке ЕОКТИ,АМ-4. Алгоритмы [3, 4] и автора для дискрепанса дают одинаковый результат. Время, относящееся к алгоритмам автора и [3, 4], будем обозначать, соответственно, tN,n и TN,n. В таблице приведены величины дискрепансов для рядов Холтона [1, 7] при п = 2, 3 и Варнока [1, 8] при п = 4, 5. Ряды Холтона и Варнока размерности п моделировались на основе первых п из следующих простых чисел 13, 17, 19, 3, 11.

Таблица

п N ,71 ТМ,п ,71 п N ,71 ТМ ,п ,71

2 500 0.15 5.7 0.0262 3 100 0.37 1.1 0.1300

2 1000 0.57 43 0.0160 3 200 2.85 17.6 0.0881

2 3000 8 1140 0.0067 3 500 38 960 0.0390

2 4000 13 2520 0.0053 3 2000 2400 - 0.0131

4 100 28 28 0.0852 5 50 75 15 0.1629

4 150 120 191 0.0689 5 100 2460 490 0.1078

4 200 200 1110 0.0539 5 150 15900 5340 0.0783

Как видно из таблицы при 2 < п < 4 подсчет дискрепанса по алгоритму, описанному в данной статье занимает меньше времени, чем по алгоритмам [3, 4] при всех N и, чем больше N, тем эта разница существеннее. При п = 5, наоборот, tN,5 > т^б-

Обсудим результаты, полученные при численных экспериментах. Рассмотрим отношение роста времен tN,5 и т^б в зависимости от параметров N и п.

При п = 5, N0 = 50 и N = ^-N0 = 100 из таблицы имеем приближенные соотношения (которые справедливы и при других п):

tN,n ~ tNo,n • кп, TN,n ~ TN0,n • кп. (13)

При увеличении N время для алгоритмов [3, 4] растет быстрее, чем заложено в формуле (13), а для алгоритма автора, наоборот. А именно, Т1бо,б ~ 1.5 • 35 тбо,б, ^50,5 ~ 0.9 • 35 тбо,5 •

Если выбирать алгоритм, основываясь на минимизации времени счета, то можно найти N0 такое, что при N < N0 следует использовать алгоритм [4], а при N > N0 — алгоритм, предлагаемый нами.

Описанный алгоритм прост для программирования и требует дополнительно п — 1 массив целых чисел размерности N + 2, а для алгоритмов [3, 4] нужно дополнительно п(п — 1)/2 массивов реальных чисел размерности N + 2.

T. M. Tovstik. The algorithm of the discrepancy calculation for the finite number of points in the unit n-dimensional cube [0,1)n is given.

This algorithm is simple for programming. For 2 < n < 4 the time of the discrepancy calculation by this algorithm is essentially lass than for the P. Bundschuh and Y. Zhu’s algorithm. If n is larger then the choice of the quick algorithm depends of the number of points.

Литература

1. Товстик Т. М. Сравнение некоторых статистических свойств квази-случайных и псевдослучайных чисел // Вестн. С.-Петерб. ун-та. Сер. 1. 2006. Вып. 1. С. 62-70.

2. Niederreiter H. Random number generation and quasi-Monte-Carlo methods // Soc. Indust. and Appl. Math. Philadelphia (Pennsylvania), 1992. 243 p.

3. Bundschuh P., Zhu Y. The formulas of exact calculation of the disgrepancyof low-dimensinal finite point sets (I) // Chinese Sci. Bul. 1993. Vol. 38. N15. P. 1318-1319.

4. Bundschuh P., Zhu Y. The formulas of exact calculation of the disgrepancyof low-dimensinal finite point sets (II) // Chinese Sci. Bul. 1995. Vol. 40. N7. P. 610-612.

5. Bergstrom V. Einige Bemerkunger zur theorie der diophantischen approximationen // Fysiogr. Salsk. Lund. Forh. 6. 1936. N13. P. 1-19.

6. Kolmogoroff A. N. Sulla determinazione empirica di una legge di distribuzione // G. Ist. ital. attuar. 1933. Vol. 4. N1. P. 83-91.

7. 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.

8. Warnock T. Effective error estlmates for quasi-Monte-Carlo computations // Int. Conf. on Monte Carlo and Quasi-Monte-Carlo Methods. 1998. Claremont. CA.

Статья поступила в редакцию 15 марта 2007 г.