Choliy.NumericaMethods.2011Dec01

9. Статистика

У цьому роздiлi розглядаються окремi питання прикладної статистики: генерацiя випадкових чисел, моменти, порiвняння середнього та дисперсiї за статистичними критерiями. Не для всiх приведених результатiв та алгоритмiв виконаний строгий та повний аналiз. В основному вони приведенi так, щоб ними було можна користуватися, а доведення необхiдних фактiв залишено для курсу теорiї iмовiрностi та статистики.

Математична статистика - галузь математики, що займається методами збору та iнтерпретацiї даних. Її головною задачею є отримання осмислених висновкiв на основi даних, якi не завжди узгоджуються мiж собою. Ця наука має свiй апарат i методи. Один iз її методiв уже розглядався як метод найменших квадратiв. Область використання методiв статистики надзвичайно широка.

Основним поняттям статистики є поняття випадкового числа, єдиною властивiстю якого є його розподiл. Сукупнiсть випадкових чисел утворює вибiрку. Властивостi вибiрки називаються статистиками.

Генерацiя випадкових чисел комп’ютером, що є повнiстю детермiнованим автоматом - справа неможлива. Насправдi комп’ютери генерують псевдовипадковi числа - такi числа, що за своїми властивостями нагадують випадковi. Справжнiм генератором випадкових чисел мiг би бути лiчильник Гейгера з iнтерфейсом для вводу його вiдлiкiв у комп’ютер.

Рiвномiрно розподiленi випадковi числа зазвичай генеруються лiнiйним конгруентним методом, за яким наступне число в послiдовностi залежить вiд попереднього:

Ij+1 = aIj + b mod m,

де mod m означає ”за модулем m“, тобто остачу вiд дiлення на m. Iншi величини - a, b - параметри, як правило, цiлочисельнi сталi. Дуже вiдомi в колах програмiстiв стандартнi функцiї srand та rand теж використовують цей метод, як видно з наведеного нижче шматка коду. Цей алгоритм не може розпочати роботу самостiйно, для генерацiї першого випадкового числа йому потрiбно передати на вхiд випадкове число.

81

#define RAND_MAX 32768 unsigned long next = 1; int rand(void)

{

next = next * 1103515245 + 12345;

return (unsigned int)(next/65536)%RAND_MAX;

}

void srand(unsigned int seed)

{

next = seed;

}

Якщо привести згенерованi цим алгоритмом випадковi числа до iнтервалу, наприклад, [0, 1] i вiдкласти на графiку залежнiсть випадкового числа вiд попереднього, то легко переконатися, що числа не бажають рiвномiрно заповнювати одиничний квадрат, а менш-бiльш регулярно шикуються вздовж прямих лiнiй. Iнший момент полягає в тому, що навiть у найкращому випадку через RAND_MAX викликiв числа почнуть повторюватись. Це може спричинити дуже цiкавi фiзичнi результати: якщо будуючи фiзичну модель ми викликаємо цю функцiю, наприклад, 100*RAND_MAX разiв, то це аж нiяк не означає, що ми отримаємо 100*RAND_MAX рiзних випадкових чисел, насправдi це буде 100 разiв одних i тих же RAND_MAX чисел, якi утворюють перiодичну послiдовнiсть. Результати такого фiзичного моделювання можна серйозно критикувати.

Лiнiйний конгруентний метод неодноразово пiддавався глибокому аналiзу. Виявилось, що старшi бiти генерованих чисел значно випадковiшi, нiж молодшi. Тому нiяк не можна рекомендавати таку процедуру для отримання випадкового числа в iнтервалi [p, q]:

j = p + (rand()%(q - p));

значно краще робити це так:

j = p + (int)((double)(q - p)*rand()/(RAND_MAX + 1.));

другий алгоритм використовує усi бiти згенерованого випадкового числа, а попереднiй - лише кiлька останнiх.

Незважаючи на критику, метод залишається широко розповсюдженим. На його основi будуються iншi методи. Алгоритм Левiса використовує ще простiшу версiю формули:

Ij+1 = aIj mod m,

82

з такими, наприклад, параметрами: a = 16807 = 75, m = 231−1 = 2147483647. Однак, пряма реалiзацiя цього способу на комп’ютерi неможлива, оскiльки добуток у правiй частинi бiльший, нiж найбiльше з допустимих цiлих чисел. Якщо ж m попередньо розбити m = aq + r, то можна показати, що

Ij+1 = a(Ij mod q) − r[Ij/q],

де [] - операцiя взяття цiлої частини, i якщо отримане число Ij+1 виходить вiд’ємним, то до нього треба додати m. Для нашого першого прикладу q = 127773, r = 2836. Перiод такого генератора = RAND_MAX*RAND_MAX. Якщо не дуже замислюватися про якiсть отриманої випадковостi, то осмислений вибiр a, b, m - це все, що потрiбно для швидкого (i неякiсного) генератора випадкових чисел.

Лiнiйний конгруентний метод, на жаль, генерує числа, якi погано задовольняють критерiю випадковостi, крiм того, вони досить швидко повторюються. Криптографiя, однак, потребує значно кращих випадкових чисел, що продовжує стимулювати роботи в цьому напрямку. Наведем тут для прикладу два генератори Дж.Марсальї [20], якi мають величезнi перiоди повторення, порядку 2125.

unsigned long mwc()

{

static unsigned long x = 123456789L,

y = 362436069L, z = 77465321L, c = 13579L.

unsigned long long t; // unsigned __int64 для Windows t = 916905990LL * x + c;

x = y; y = z; c = (t >> 32);

return z = (t & 0xffffffff);

}

unsigned long x128()

{

static unsigned long x = 123456789L,

y = 362436069L, z = 521288629L, w = 88675123L. unsigned long t;

t = (x ^ (x << 11)); x = y; y = z; z = w;

return (w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)));

}

83

Ми не будемо зупинятися тут на поясненнi того як вони працюють i чому результат достатньо ”гарний“, однак для потреб нашого посiбника цих генераторiв цiлком достатньо. Ще кращi алгоритми слiд шукати в лiтературi, для початку в [21], [2].

Наведенi приклади за кожним викликом функцiї генерують випадкове число у iнтервалi [0, 232 −1]. Щоб привести його до iнтервалу [0, 1], перетворимо результат до типу double i подiлимо на 232 − 1 = 2147483647.

Єдиною характеристикою випадкового числа є його розподiл. Якщо p(x) - розподiл числа x, то p(x)dx - iмовiрнiсть того, що число x лежить в iнтервалi [x, x + dx]. Iнтеграл вiд функцiї розподiлу по всiй областi її визначення має бути строго рiвен одиницi. Випадковi числа з довiльним розподiлом можна згенерувати, прийнявши до уваги те, що якщо величина x розподiлена за p(x), то y(x) буде розподiлена як

Якщо, наприклад, y(x) = − ln x, то p(y)dy = exp (−x)dx, а y(x) має експоненцiйний розподiл, а величина y/λ матиме розподiл Пуассона:

p(y)dy = λ exp (−λx)dx.

Випадковi числа з нормальним розподiлом можна отримати таким методом: вiзьмемо двi випадковi величини x1, x2, розподiленi рiвномiрно у iнтервалi [0, 1] i перетворимо їх до y1 та y2:

p p

y1 = −2 ln x1 cos (2πx2), y2 = −2 ln x1 sin (2πx2),

якi будуть розподiленi нормально.

Ще один спосiб генерацiї нормально розподiлених на [0, 1] випадкових чисел використовує той факт, що середнє арифметичне значення достатньо ве-

ликої кiлькостi n випадкових чисел розподiлених рiвномiрно на [0, 1] буде

√

розподiлене нормально з середнiм значенням 1/2 та дисперсiєю n. Моменти випадкових величин. Це найбiльш вiдомi статистики. Сере-

днє значення (перший момент):

84

варiацiя (другий момент):

Середнє арифметичне значення є найбiльш iмовiрною оцiнкою справжнього значення за результатами окремих спостережень. Воно також мiнiмiзує середньоквадратичну похибку (дисперсiю). Асиметрiя показує наскiльки розподiл несиметричний вiдносно середнього арифметичного, значення бiльше нуля вказує на те, що хвiст розподiлу сильнiше виражений при x > x¯. Крутiсть описує наскiльки розподiл гострий бiля x¯: додатнє значення говорить про те, що розподiл крутiший, нiж нормальний.

Перевiрка гiпотез. Це статистична процедура, що дозволяє вияснити iмовiрнiсть статистичного результату. Для перевiрки, чи спiвпадають середнi значення двох вибiрок xi, i = 1, . . . , Nx та yi, i = 1, . . . , Ny за результатами спостережень обчислюємо

i

t = x¯ − y¯. SD

Ця величина t має розподiл Стьюдента з Nx + Ny − 2 ступенями вiльностi. Сам розподiл Стьюдента виглядає так:

−t

85

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

Для перевiрки того, чи вiдрiзняються дисперсiї двох сукупностей за результатами двох вибiрок, користуються розподiлом Фiшера з ν1 та ν2 ступенями вiльностi:

Якщо обидва значення ν1 та ν2 цiлi, орiєнтуємося по меншому з них. Величиною x тут виступає вiдношення вибiркових дисперсiй.

Методи Монте-Карло. Пiд такою назвою вiдомi методи, що суттєво використовують випадковi числа. З цiєї великої областi науки про алгоритми розглянемо лише метод знаходження iнтегралiв.

86

Нехай потрiбно знайти iнтеграл

b

Z

f(x)dx.

a

Найбiльше та найменше значення функцiї f(x) на iнтегралi [a, b] - вiдомi. Згенеруємо велику кiлькiсть випадкових точок у квадратi [a, b] × [min, max]. Якщо точки достатньо випадковi, то кiлькiсть точок, що попаде пiд криву f(x) буде пропорцiйна шуканому iнтегралу, а сам iнтеграл можна оцiнити

так:

b

Z

f(x)dx = кiлькiсть точок пiд кривою.

загальна кiлькiсть точок

a

Для того, щоб точка (x, y) виявилась пiд кривою, досить перевiрити умову y < f(x). Цей метод у аналогiчнiй постановцi можна використовувати також для знаходження площ складних фiгур чи об’ємiв складних тiл.

87

Задачi. У задачах цього роздiлу, що мають парнi номери, потрiбно знайти площу поверхнi чи об’єм фiгури, заданих вказаними рiвняннями методом Монте-Карло.

1. Згенеруйте 10000 рiвномiрно розподiлених випадкових чисел на iнтервалi [0, 1] методом mvc, знайдiть частоти попадання випадкового числа у кожен з iнтервалiв [ai, ai+1], де ai = 0.1i, i = 0, . . . , 10.

2.

y2 = x3, y2 = 8(6 − x)2.

3. Згенеруйте двi вибiрки нормально розподiлених випадкових чисел з рiзними середнiми значеннями та однаковими дисперсiями. Знайдiть для них

A(t|ν).

4.

x2 + y2 = 4, x + y + z = 4, x = 0.

5. Згенеруйте 10000 рiвномiрно розподiлених випадкових чисел на iнтервалi [0, 1] методом x128, знайдiть частоти попадання випадкового числа у кожен з iнтервалiв [ai, ai+1], де ai = 0.1i, i = 0, . . . , 10.

6.

3x2 = 25y, 5y2 = 9x.

7. Змоделюйте 10000 кидань пари кубикiв. Вияснiть, з якою частотою випадає кожна з можливих комбiнацiй.

8.

x + y + z = 2, 3x + y = 2, 3x + 2y = 4, y = 0, z = 0.

9. Згенеруйте 10000 рiвномiрно розподiлених випадкових чисел на iнтервалi [0, 1] методом mvc, i побудуйте на їх основi 10000 випадкових чисел з нормальним розподiлом. Знайдiть частоти попадання випадкового числа у кожен з iнтервалiв [ai, ai+1], де ai = 0.1i, i = −10, 9, . . . , 0, 1, . . . , 10.

10.

y = exp x, y = exp(2x), x = 1.

11. Побудуйте послiдовнiсть нулiв та одиниць, у якiй одиничка випадає з частотою 1/4, а нулик - з частотою 3/4.

12.

13. Згенеруйте 10000 рiвномiрно розподiлених випадкових чисел на iнтервалi [0, 1] методом x128, i побудуйте на їх основi 10000 випадкових чисел з

88

нормальним розподiлом. Знайдiть частоти попадання випадкового числа у кожен з iнтервалiв [ai, ai+1], де ai = 0.1i, i = −10, 9, . . . , 0, 1, . . . , 10.

14.

ρ = 4 sin ϕ, ρ = 2 sin ϕ.

15. Побудуйте послiдовнiсть нулiв, одиниць та двiйок, у якiй нулик випадає з частотою 1/8, одиничка - з частотою 1/4, двiйка - у всiх iнших випадках.

16.

17. Згенеруйте 10000 рiвномiрно розподiлених випадкових чисел на iнтервалi [0, 1] методом mvc, знайдiть чотири першi моменти цiєї вибiрки.

18.

xy = 4, x + y = 5.

19. Згенеруйте двi вибiрки нормально розподiлених випадкових чисел з рiзними середнiми значеннями та однаковими дисперсiями. Знайдiть для них

F (x).

20.

x2 + y + 2 + z2 = 8, x2 + y2 = 4.

21. Згенеруйте 10000 рiвномiрно розподiлених випадкових чисел на iнтервалi [0, 1] методом x128, знайдiть чотири першi моменти цiєї вибiрки.

22.

ρ= 2 cos 2ϕ.

23.Змоделюйте викидання трiйки монет. Яка комбiнацiя випадає найчастiше.

24.

y = z = 2, y − z = 2, x2 + y2 = 1.

25. Згенеруйте 10000 рiвномiрно розподiлених випадкових чисел на iнтервалi [0, 1] методом mvc, i побудуйте на їх основi 10000 випадкових чисел з нормальним розподiлом. Знайдiть першi чотири моменти цiєї вибiрки.

26.

ρ= 2 sin 3ϕ.

27.Змоделюйте кидання монети. З якою частотою випадають комбiнацiї "решка-орел"та "орел-решка"i "решка-решка-орел"та "орел-орел-решка".

28.

x2 + y2 = z2, x2 + y2 = 2y, z = 0.

89

29. Згенеруйте 10000 рiвномiрно розподiлених випадкових чисел на iнтервалi [0, 1] методом x128, i побудуйте на їх основi 10000 випадкових чисел з нормальним розподiлом. Знайдiть першi чотири моменти цiєї вибiрки.

30.

(x2 + y2)2 = 2y3.

31. Змоделюйте викидання пари неправильних кубикiв, у яких шiстка випадає вдвiчi рiдше, нiж будь-яка з iнших граней, кожна з яких випадає з однаковою частотою. Яка комбiнацiя випадатиме в такому випадку найчастiше.

32.

(x2 + y2)3 = 4(x4 + y4).

90