четвер, 24 травня 2018 р.

Розв'язки до вправ з "An Introduction to the Analysis of Algorithms" Роберта Седжвіка. Розділ II.

Вправа 2.2
Зовнішній цикл робить $N_{max}$ ітерацій. Внутрішній $\frac{N_{max}(N_{max}-1)}{2}$. На кожній ітерації внутрішнього циклу виконується два додавання і одне ділення. В підсумку $3\frac{N_{max}^2(N_{max}-1)}{2}$.
Вправа 2.3
int Cn(int n) {
    if (n == 1) return 2;
    return (1 + 1. / n) * Cn(n - 1) + 2;
}
На кожному рівні рекурсії ми маємо дві операції з $C_n$ - множення і додавання. Всього таких рівнів $N-1$, отже $2(n-1)$.
Вправа 2.6
\[ a_n = q \operatorname{Fib}(n-1) + p \operatorname{Fib}(n-2) \]
Вправа 2.7
\[ a_n = q \operatorname{Fib}(n-1) + p \operatorname{Fib}(n-2) + r \operatorname{Fib}(n-2) \] p, q, p + q + r, p + 2q + r, 2p + 3q + 2r, 3p + 5q + 3r, 5p + 8q + 5r.
Вправа 2.8
\[ f^*(a_{n-1}, a_{n-2}) = f(a_{n-1}, a_{n-2}) - f(0,0) \] Маємо, що $f^*$ однорідна. Нехай \begin{align*} u_n &= f^*(u_{n-1}, u_{n-2}),\quad \text{для}\ n > 1, u_0 = 1, u_1 = 0,\\ v_n &= f^*(v_{n-1}, v_{n-2}),\quad \text{для}\ n > 1, v_0 = 0, v_1 = 1. \end{align*} Тоді \[ a_n = a_0 u_n + a_1 v_n + f(0,0)\operatorname{Fib}(n-2). \]
Вправа 2.9
\[ a_n = \frac{n}{n+2}a_{n-1} = \prod_{k=1}^n \frac{k}{k+2} = \frac{1\cdot 2\cdot \cdots n}{3\cdot 4 \cdots (n+2)} = \frac{1\cdot 2}{(n+1)(n+2)}. \]
Вправа 2.10
\[ a_n = 1 + \underbrace{(-1) + 2}_1 + \underbrace{(-3) + 4}_1 \cdots, \] \[ a_n= \begin{cases} 1 + \frac{n}{2} & n\ \text{парне}\\ 1 - \frac{n}{2} & n\ \text{непарне} \end{cases} \]
Вправа 2.11
\begin{align*} n(n-1)a_n &= (n-1)(n-2)a_{n-1} + 2(n-1)\\ (n-1)(n-2)a_{n-1} &= (n-2)(n-3)a_{n-4} + 2(n-3)\\ \dots&\\ 2\cdot 1\cdot a_2 &= 0\cdot(-1)\cdot 1 + 2\cdot(n-(n-2)-1) \end{align*} Тепер складемо всі рівності. \[ n(n-1)a_n = 2\sum_{k=1}^{n-1}k = n(n-1) \] Отримуємо: \[ a_n = 1 \]
Вправа 2.12
\begin{align*} 2^{-n}a_n &= 2^{-(n-1)}a_{n-1} + 2^{-n}\\ \dots &\\ 2^{-2}a_2 &= 2^{-1}\cdot 1 + 2^{-2}\\ \end{align*} Тепер складемо всі рівності. \[ 2^{-n}a_n = 2^{-1}\cdot 1 + \sum_{k=2}^{n}2^{-k} \] Отримуємо: \[ a_n = 2^{n-1} + \sum_{k=2}^n2^{n-k} = 2^{n-1} + \sum_{k=0}^{n-2}2^k = 2^n - 1. \]
Вправа 2.13
Тут нам треба звернути увагу на те, що на відміну від теореми 2.1, нулю дорівнює не $a_1$, а $a_0$. Тому сума починається з $j = 0$: \begin{align*} a_n &= 1 + \sum_{j=0}^{n-1}\frac{j+1}{j+2}\cdots\frac{n}{n+1}\\ &= 1 + \sum_{j=0}^{n-1}\frac{j+1}{n+1} = 1 + \frac{1}{n+1}\sum_{j=0}^{n-1} (j+1)\\ &= \frac{1}{n+1}\sum_{j=1}^{n+1}j = \frac{n+2}{2}. \end{align*}
Вправа 2.14
Тут сума починається з $t+1$, також присутній член $a_t$, бо не сказано, що він дорівнює нулю. \begin{equation} a_n = y_n + \sum_{j=t+1}^{n-1} y_{j} x_{j+1} \dots x_n + a_tx_{t+1}\dots x_n.\label{eq:rr_at} \end{equation}
Вправа 2.15
\begin{align*} a_n &= 2 \sum_{j=1}^{n-1} \frac{j+2}{j+1} \dots \frac{n+1}{n} + 2\\ &= 2 \sum_{j=1}^{n-1}\frac{n+1}{j+1} + 2\\ &= 2(n+1)\left(\sum_{j=1}^{n-1}\frac{1}{j+1}+\frac{1}{n+1}\right)\\ &= 2(n+1)(H_{n+1}-1). \end{align*} Під час перевірки можна скористатись тим, що $nH_n - n = \sum_{k=1}^{n-1}H_k$. \begin{align*} 2n\sum_{k=1}^nH_k &= 2(n+1)\sum_{k=1}^{n-1}H_k + 2n\\ nH_n &= \sum_{k=1}^{n-1}H_k + n \end{align*}
Вправа 2.16
Скориставшись теоремою 2.1: \begin{align*} a_n = 12 H_n + 12 \sum_{j=5}^{n-1} H_j \frac{(j-3)(j-2)(j-1)j}{n(n-1)(n-2)(n-3)}, \end{align*} або ми можемо використати сумувальний множник $(n-1)(n-2)(n-3)$: \begin{align*} n(n-1)(n-2)(n-3)a_n &= (n-1)(n-2)(n-3)(n-4)a_{n-1}\\ &\phantom{=}+ 12 n(n-1)(n-2)(n-3) H_n,\\ (5\cdot 4\cdot 3\cdot 2)a_5 &= 12 (5\cdot 4\cdot 3\cdot 2) H_5. \end{align*} Сумуючи отримуємо той самий вислід.
Вправа 2.17
Зауважимо, що кожен 2-вузол містить одну вершину, а кожен 3-вузол містить дві вершини, тобто розпадається на два 2-вузли. \begin{align*} A_N &= A_{N-1}\left(\frac{N-6}{N}\right) +2 \\ &= 2 + 2 \sum_{j=1}^{N-1}\left(\frac{j-5}{j+1}\right)\dots\left(\frac{N-7}{N-1}\right)\\ \end{align*} Тут і у попередній вправі $x_n$ дуже схожі. І виконуючи теорему 2.1 нам може захотітись дати таку відповідь: \[ A_N = 2 + 2 \sum_{j=1}^{N-1} \frac{(j-5)(j-4)(j-3)(j-2)(j-1)j}{N(N-1)(N-2)(N-3)(N-4)(N-5)}, \] але це не правильно, бо ми не можемо ділити на нуль. Ми мусимо використати \eqref{eq:rr_at}. Отже, відповіддю буде \begin{align*} A_N &= 2 + 2 \sum_{j=6}^{N-1} \frac{(j-5)(j-4)(j-3)(j-2)(j-1)j}{N(N-1)(N-2)(N-3)(N-4)(N-5)}\\ &\phantom{=} + A_5\frac{0\cdot1\cdot2\cdot3\cdot4\cdot5}{N(N-1)(N-2)(N-3)(N-4)(N-5)},\quad \text{для }\ N > 5. \end{align*} Нам навіть не потрібно вручну обчислювати $A_5$, бо його потрібно множити на $0$.
Щоб краще уявити собі кількість 2-вузлів і перевірити формулу, можна запустити таку програму:
float tree23(int n) {
    if (n == 0) return 0;
    float prev = tree23(n - 1);
    return prev - 2 * prev / n + 2 * (1 - 2 * prev / n);
}

float tree23_sum(int n) {
    float sum = 0.;
    for (int j = 6; j < n; ++j) {
        float num = 1.f, denom = 1.f;
        for (int k = 0; k < 6; ++k) {
            num *= j - k;
            denom *= n - k;
        }
        sum += num / denom;
    }
    return 2 + 2 * sum;
}

int main() {
    for (int i = 6; i < 30; ++i) {
        cout << setw(6) << setprecision(3) << fixed << tree23(i) << " " << tree23_sum(i) << endl;
    }
    return 0;
}

 2.000 2.000
 2.286 2.286
 2.571 2.571
 2.857 2.857
 3.143 3.143
 3.429 3.429
 3.714 3.714
 4.000 4.000
 4.286 4.286
 4.571 4.571
 4.857 4.857
 5.143 5.143
 5.429 5.429
 5.714 5.714
 6.000 6.000
 6.286 6.286
 6.571 6.571
 6.857 6.857
 7.143 7.143
 7.429 7.429
 7.714 7.714
 8.000 8.000
 8.286 8.286
 8.571 8.571
Вправа 2.18
\begin{align*} a_1 &= \frac{1 + 0\cdot a_0}{1+1\cdot a_0} = \frac{\operatorname{Fib}_1 + \operatorname{Fib}_0\cdot a_0}{\operatorname{Fib}_2 + \operatorname{Fib}_1\cdot a_0}\\ a_2 &= \frac{1 + 1\cdot a_0}{2+1\cdot a_0} = \frac{\operatorname{Fib}_2 + \operatorname{Fib}_1\cdot a_0}{\operatorname{Fib}_3 + \operatorname{Fib}_2\cdot a_0}\\ \dots&\\ a_n &= \frac{\operatorname{Fib}_n + \operatorname{Fib}_{n-1}\cdot a_0}{\operatorname{Fib}_{n+1} + \operatorname{Fib}_n\cdot a_0}\\ \end{align*} Зі збільшенням $n$ перші доданки починають домінувати і ми отримуємо $\frac{\operatorname{Fib}_{n}}{\operatorname{Fib}_{n+1}}$, тоді як ми знаємо, що $\lim_{n\to \infty}\frac{\operatorname{Fib}_{n+1}}{\operatorname{Fib}_{n}} = \varphi$ - золотий перетин, а $\frac{1}{\varphi}$ і є потрібним результатом, з яким $\lim_{n\to \infty} b_n = 0$.
Вправа 2.19
На рисунку праворуч показани графіки функцій $ y = \cos(x)$ і $x = \cos(y)$ або інакше $y = \arccos(x)$. Точка перетину має абсцису $0.73908\dots$.
Вправа 2.20
Щоб зрозуміти чому формула працює можна розглянути такий її варіант: $$2a_na_{n-1} = a_{n-1}^2+b$$ Якщо $b$ від'ємне, то нема причин очікуватись збіжність $a_n$ до $b$.
Вправа 2.21
  1. Спершу доведемо, що $a_n \in \Theta(1/n)$. \begin{align*} \frac{1}{a_n} &= \frac{1}{a_{n-1}}\frac{1}{1-a_{n-1}}&\\ \frac{1}{a_n} &< \frac{1}{a_{n-1}}\frac{n-1}{n-2}&\text{використали}\ a_n < \frac{1}{n}\\ a_n &> \frac{1}{n-1}a_3&\text{де}\ a_3 = 39/256. \end{align*} Отже, для достатньо великих $n$, $a_n > \frac{1}{n-1}$.
  2. Тепер оцінимо $c$:
        float a = 0.5f;
        for (int i = 1; i < 1000; ++i) {
            a = a * (1 - a);
            cout << setw(15) << a << " " << a * i << endl;
        }
    
    Така програма видає нам:
    ...
        0.000996334 0.991353
        0.000995342 0.99136
        0.000994351 0.991368
        0.000993362 0.991376
        0.000992376 0.991383
    
    Отже, схоже що $c$ дорівнює $1$.
  3. З пункту 1 і з книги ми знаємо, що $$\frac{1}{n-1} < a_n < \frac{1}{n}$$ отже $na_n \to 1$.
Вправа 2.22
  1. На початку $\sin(x)$ зростає повільніше ніж $x$, тобто $x > \sin(x)$ окрім як в $0$. Отже послідовність $a_n$ обмежена спадна функція, тобто вона має границю. $\sin(x)$ - неперервна, отже ця границя в стаціонарній точці, єдина така точка в $[0,1]$ це $0$.
  2. Тепер доведемо, що $a_n = O(1/\sqrt{n})$.
Вправа 2.23
Розглянемо $$a_n = \alpha + \varepsilon_{n} = f(a_{n-1}) = f(\alpha + \varepsilon_{n-1}).$$ Використаємо розклад в ряд Тейлора: $$f(\alpha + \varepsilon_{n-1}) = f(\alpha) + f'(\alpha) * \varepsilon_{n-1} + f''(\alpha) * \varepsilon_{n-1} / 2!.$$ Зауважимо, що $f(\alpha) = \alpha$. \begin{equation} \varepsilon_{n} = f'(\alpha) * \varepsilon_{n-1} + f''(\alpha) * \varepsilon_{n-1} / 2!\label{eq:cnvrg_rate_drv} \end{equation} Отже, якщо $f'(\alpha) > 1$, то ряд розбігається.
Вправа 2.24
Використовуючи \eqref{eq:cnvrg_rate_drv} маємо квадратну збіжність якщо $f'(\alpha)$ рівне $0$ з коефіцієнтом $f''(\alpha)/2$. Аналогічно можна отримати коефіцієнти для $f'(\alpha) \ne 0$.
Вправа 2.25
Для того, щоб збагнути як утворити необхідне рекурентне відношення розглянемо таке рекурентне відношення: \begin{equation} a_n - 3*a_{n-1} = 4*(a_{n-1} - 3*a_{n-2}),\quad a_0 = 0, a_1 = 1.\label{eq:rr_25_43} \end{equation} Якщо пройти рекурсивно до базового випадку, маємо: \begin{align*} a_n - 3 a_{n-1} &= 4^{n-1}\\ 3^1(a_{n-1} - 3 a_{n-2}) &= 3^1 4^{n-2}\\ \dots &\\ 3^i(a_{n-i} - 3 a_{n-i-1}) &= 3^i 4^{n-i-1}\\ \dots &\\ 3^{n-1}(a_1 - 3 a_0) &= 3^{n-1} \end{align*} Складаємо і скорочуємо: \begin{align*} a_n = \sum_{i=1}^n 4^{n-i}3^{i-1} \end{align*} Зауважимо, що $4^n = 1\dot 4^{n-1} + 3\cdot 4^{n-1}$ або ж $3^{n-1}\dot 4^0 + 3^n\cdot 4^0 = 4\cdot 3^n$. Застосовуючи останню формулу маємо для $n=5$: $$ \underbrace{4^4\cdot 3^0 + \underbrace{4^3\cdot3 + \underbrace{4^2\cdot3^2+ \underbrace{4\cdot3^3 + \underbrace{3^4\cdot 4^0 + 3^5\cdot 4^0 }_{3^4\cdot 4} }_{3^3\cdot 4^2} }_{3^2\cdot 4^3} }_{3^1\cdot 4^4} }_{4^5} $$ Зверніть увагу, що тут ми додали $3^5$, отже формула така: $$\sum_{i=1}^n 4^{n-i}3^{i-1} = 4^n - 3^n.$$ Тобто розв'язок до \refeq{eq:rr_25_43} такий $a_n = 4^n - 3^n$. Гм.. це маже те, що нам треба! А ось і наша відповідь: $$a_n - 3*a_{n-1} = 4*(a_{n-1} - 3*a_{n-2})+2^{n-1}.$$
long a(int n) { 
    switch (n) {
    case 0: return 1;
    case 1: return 3;
    default: return powl(4, n) - powl(3, n) + powl(2,n);
    }
}

    for (int n = 2; n < 25; ++n) {
        long lhs = a(n) - 3 * a(n - 1);
        long rhs = 4 * (a(n - 1) - 3 * a(n - 2)) + powl(2,n-1);
        cout << setw(2) << n << ": " << setw(10) << lhs << " " << rhs << endl;
    }

 2:          2 2
 3:         12 12
 4:         56 56
 5:        240 240
 6:        992 992
 7:       4032 4032
...
Вправа 2.26
  1. Ввести $b_n = a_n - a_{n-1}$, це дасть нам однорідне рекурентне відношення.
  2. Розписати в стовпчик усі $a_i - a_{i-1}$ і скласти, завдяки телескопуванню, ми зможемо отримати формулу для $a_n$, спрощення цієї формули може бути досить складним.
Вправа 2.27
Характеристичний многочлен такий: $$\beta^2 - 5\beta + 6 = 0.$$ Коренями є $2, 3$. Значить маємо розв'язки у вигляді $c_12^n + c_23^n$. Звідси бачимо, що для розв'язку $2^n$ $c_1 = 1, c_2 = 0$, а розв'язок $2^n - 1$ отримати не можливо.
Вправа 2.28
Корені $2, -i, i$.
  1. $a_0 = a_1 = a_2 = 0$.
  2. $a_0 = a_1 = a_2 = 1$.
  3. $a_0 = i, a_1 = -1, a_2 = i$.
Вправа 2.29
Коренями є $1+\sqrt{5}, 1-\sqrt{5}$. Відповідь $\frac{(1+\sqrt{5})^n - (1-\sqrt{5})^n}{2\sqrt{5}}$.
Вправа 2.30
Перепишемо рекурентне відношення у вигляді $$a_n - a_{n-1} = a_{n-1} - a_{n-2}.$$ Ми бачимо, що зі зростанням $n$ на одиницю $a$ зростає на сталу. Звідси нам легко знайти відповіді.
  1. $a_n = n$.
  2. $a_n = 1$.
Вправа 2.31
Тут корені це $\frac{1\pm i\sqrt{3}}{2}$. Отже, маємо два рівняння: \begin{cases} c_1 + c_2 &= 0\\ c_1\frac{1 + i\sqrt{3}}{2} + c_2\frac{1 - i\sqrt{3}}{2} &= 1 \end{cases} Що дає $c_1 = -c_2 = \frac{1}{i\sqrt{3}}$. І, власне, розв'язок: $$a_n = \frac{1}{i\sqrt{3}}\left(\frac{1 + i\sqrt{3}}{2}\right)^n -\frac{1}{i\sqrt{3}}\left(\frac{1 - i\sqrt{3}}{2}\right)^n$$
Вправа 2.32
Корені тут $\frac{1}{2}, \frac{1\pm\sqrt{3}}{2}$. \begin{cases} c_1+c_2+c_3 &= 0\\ c_1\frac{1}{2}+c_2 \frac{1+\sqrt{3}}{2}+c_3 \frac{1-\sqrt{3}}{2} &= 1\\ c_1\left(\frac{1}{2}\right)^2+c_2 \left(\frac{1+\sqrt{3}}{2}\right)^2+c_3 \left(\frac{1-\sqrt{3}}{2}\right)^2 &= 2 \end{cases}
Вправа 2.33
$$a_n = 4(a_{n-2}-a_{n-4}).$$ При $a_0 = 1, a_1 = -1, a_2 = 2, a_3 = -2$. Для парних маємо $1,2,4,8,\dots$ і для непарних те саме тільки з мінусом.
Вправа 2.34
>> p = [1 -1 -1 -1];
>> r = roots(p)
r =

   1.83929 + 0.00000i
  -0.41964 + 0.60629i
  -0.41964 - 0.60629i

>> A = [r.^0'; r.^1'; r.^2'];
>> b = [0 0 1]';
>> c = A\b
c =

   0.182804 + 0.000000i
  -0.091402 - 0.340547i
  -0.091402 + 0.340547i

>> F20_exact = r.^19' * c
F20_exact = 1.9513e+004 + 6.5683e-013i
>> F20_approx = r(1)^19' * c(1)
F20_approx = 1.9513e+004 + 6.5683e-013i
Вправа 2.35
Домножити обидва боки на $(n-2)!$. Маємо Фібоначчі: $$n!a_n = (n-1)!a_{n-1} + (n-2)!a_{n-2}.$$ Отже, $a_n = F_n / n!$.
Вправа 2.36
Входження $t_i$ означає відсутність $a_i$ і $a_{i+1}$. Всі одночлени - унікальні.
Вправа 2.37
Характеристичний многочлен для послідовності $b_n$ такий: $x^2 - \frac{1}{2}x - \frac{1}{2} = 0$. Корені - $1, -\frac{1}{2}$. \begin{cases} c_1+c_2 &= 0\\ c_1 + c_2\left(-\frac{1}{2}\right) & = 1 \end{cases} Звідси, $c_1 = -c_2 = \frac{2}{3}$. І відповідь $$b_n = \frac{2}{3} - \frac{2}{3}\left(-\frac{1}{2}\right)^n.$$ $$a_n = 2^{\frac{2}{3} - \frac{2}{3}\left(-\frac{1}{2}\right)^n}.$$
Вправа 2.38
Введемо заміну змінних $b_n = a_n^2$. Маємо: $$b_n = b_{n-1} + 1.$$ Це відношення має розв'язок $b_n = n$. Відповідно, $a_n = \sqrt{n}$.
Вправа 2.39
Для отримання формули в для $a_0 = 3, a_0 = 4$ необхідно підставити конкретні значення в формулу наведену в книзі. У випадку коли $a_0 = \frac{3}{2}$ послідовність "хаотично" рухається в $[-2, 2 ]$.
Вправа 2.40
Для великих $\epsilon$: \begin{align*} a_n &\approx \left(\frac{1}{2}(2+\epsilon + \sqrt{(2+\epsilon)^2 - 4}\right)^{2^n}\\ &= \left(\frac{1}{2}(2+\epsilon+\sqrt{4+4\epsilon+\epsilon^2 -4})\right)^{2^n}&\quad\text{в підкорінневому виразі домінує}\ \epsilon^2\\ &\approx \left(1 + \epsilon\right)^{2^n}\\ &\approx \epsilon^{2^n}. \end{align*}
Вправа 2.41
Якщо $b_n = f\cdot a_n + g$, тоді $a_n = \frac{b_n-g}{f}$. \begin{align*} \frac{b_n - g}{f} &= \alpha\left(\frac{b_{n-1} - g}{f}\right)^2 + \beta\frac{b_{n-1} - g}{f} + \gamma\\ b_n &= \frac{\alpha}{f}(b_{n-1} - 2b_{n-1}g + g^2)+\beta b_{n-1}-\beta g + \gamma f + g\\ b_n &= \underbrace{\frac{\alpha}{f}}_{=1}b_{n-1}^2 - \underbrace{\left(\frac{2\alpha g}{f}-\beta\right)}_{=0}b_{n-1} + \underbrace{\frac{\alpha g^2}{f} - \beta g + \gamma f + g}_{=-2} \end{align*} З цього ми миожемо зробити такі висновки: $f = \alpha$, $g = \frac{\beta}{2}$. А те, що $\frac{\alpha g^2}{f} - \beta g + \gamma f + g = -2$ дозволяє нам обрати значення для $\alpha, \beta, \gamma$, ми маємо два ступені свободи тут.
Наприклад, $\beta = 2, \alpha = 1, \gamma = -2$. $$a_n = a_{n-1}^2 + 2 a_{n-1} - 2$$ Отже, $f = 1$, $g = 1$ і $b_n = a_b + 1$. Перевіримо: \begin{align*} a_n + 1 &= (a_{n-1} + 1)^2 - 2\\ & = a_{n-1}^2 + 2 a_{n-1} + 1 - 2 \end{align*} Що дає $$a_n = a_{n-1}^2 + 2 a_{n-1} -2.$$
Вправа 2.42
\begin{align*} a_n &= 2 a_{n-1} \sqrt{1-a_{n-1}^2}&\quad n>0, a_0=\frac{1}{2}\\ a_n^2 &= 4 a_{n-1}^2 (1-a_{n-1}^2)\\ b_n &= 4b_{n-1}-4b_{n-1}^2&\quad b_n=a_n^2 \end{align*} Тут ми можемо використати результат попередньої вправи з $\alpha = -4$, $\beta = -4$, $\gamma = 0$. Тоді $c_n = -4b_n+2$. Маємо $$c_n = c_{n-1}^2 - 2.$$ А це вже алгоритм розподілення регістрів з книги.
Вправа 2.43
Ми можемо переписати нашу рекуренцію так: $$a_n = \underbrace{\frac{\alpha a_{n-1}}{\gamma a_{n-1} + \delta}}_{a'_{n}} + \underbrace{\frac{\beta}{\gamma a_{n-1} + \delta}}_{a''_{n}}.$$ Тепер ми можемо працювати з двома рекуренціями і скласти їх результати. \begin{align*} a'_n &= \frac{\alpha a'_{n-1}}{\gamma a'_{n-1} + \delta}\\ \frac{1}{a'_n} &= \frac{\gamma}{\alpha} + \frac{\delta}{\alpha}\frac{1}{a'_{n-1}}&\quad b'_n = \frac{1}{a'_n}\\ b'_n &= \frac{\gamma}{\alpha} + \frac{\delta}{\alpha}b'_{n-1}\\ \end{align*} Тепер подібно до теореми 2.1: $$b'_n = \frac{\gamma}{\alpha}\sum_{i=0}^{n-1}\left(\frac{\delta}{\alpha}\right)^i + \left(\frac{\delta}{\alpha}\right)^nb'_0.$$ Друга рекуренція майже як в прикладі з книги: \begin{align*} a''_n &= \frac{\beta}{\gamma a''_{n-1} + \delta}&\quad a''_n = \frac{b''_{n-1}}{b''_n}\\ b''_n &= \frac{\delta}{\beta}b''_{n-1} + \frac{\gamma}{\beta}b''_{n-2}\\ \end{align*} Тут ми знайдемо два корені і нам буде потрібно вибрати значення коефіцієнтів залежно початкових умов. Тобто, початкова умова $a_0 = a'_0 + a''_0 = 1$ має виконуватись.
Вправа 2.44
Розпишемо перші члени послідовності $(a_n)$: \begin{align*} a_0 &= 1\\ a_1 &= \frac{1}{s_1+t_1}\\ a_2 &= \frac{s_1+t_1}{s_2s_1+s_2t_1 + t_2}\\ a_3 &= \frac{s_2s_1+s_2t_1 + t_2}{s_3s_2s_1+s_3s_2t_1 + s_3t_2 + t3s_1 + t_3t_1}\\ a_4 &= \frac{s_3s_2s_1+s_3s_2t_1 + s_3t_2 + t_3s_1 + t_3t_1}{s_4s_3s_2s_1+s_4s_3s_2t_1 + s_4s_3t_2 + s_3t_2s_1 + s_3t_2t_1 + t_4s_2s_1+t_4s_2t_1 + t_4t_2}\\ \end{align*} Тут ми вже бачимо взірець за яким можна будувати наступні члени. Ми можемо записати послідовність трошки інакше: $$a_n = \frac{s_{n-1}\operatorname{den}{a_{n-2}} + t_{n-1}\operatorname{nom}{a_{n-2}}}{s_{n}\operatorname{den}{a_{n-1}} + t_{n}\operatorname{nom}{a_{n-1}}}$$ Зверніть увагу на послідовність опрацьовану в вправі 2.36, вона дуже подібна. Якщо $a_n = \frac{b_{n}}{b_{n+1}}$, то $$b_n = s_{n-1}b_{n-1}+t_{n-1}b_{n-2},\quad n>1,\ b_0 = b_1 = 1.$$
Вправа 2.45
Найперше ми маємо знайти форму $a_n$, яка дасть нам входження $n^3$ в $a_n - \left(2\sum_{j=0}^{n-1}\right)/n$. Для цього ми можемо вибрати $n(n-1)(n-2)$. \begin{align*} \sum_{k=0}^n k(k-1)(k-2) & = \sum_{k=0}^n (k^2 - 3k^2 + 2k)\\ &= \frac{n^2(n+1)^2}{4} - \frac{n(n+1)(2n+1)}{2} + n(n+1)\\ &= \frac{1}{4}(n^4 - 2n^3 - n^2 + 2n). \end{align*} Перевіримо: $$a_3 = (3\cdot 2 \cdot 1) = 6 = \frac{1}{4}(81 - 54 - 9 + 6).\quad\text{збіг!}$$ Отже, \begin{align*} a_n - \left(2\sum_{j=0}^{n-1}\right)/n & = n(n-1)(n-2)\\ &\phantom{=}- \frac{1}{2}\big((n-1)^4 - 2(n-1)^3 - (n-1)^2 + 2(n-1)\big)/n \end{align*} Тепер маючи компонент з $n^3$ ми можемо за допомогою нескладної, але нудної алгебри зкомпонувати його з уже відомими розв'язками, щоб отримати $f(n) = n^3$.
Вправа 2.46
$$c_n = n+1+\sum_{k=1}^n\frac{(n-k)(k-1)}{\binom{n}{3}}(c_{k-1} + c_{n-k}),\quad n>3.$$ Перепишемо: \begin{align*} c_n &= n+1+\frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}(n-(k+1))((k+1)-1))c_{k}\\ &=n+1+\frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}(nk- k^2 - k)c_{k} \end{align*} Заповнимо репертуарну таблицю:
  • Для $c_n = 1$:
\begin{align*} &\phantom{=\ }1 - \frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}(nk - k^2 - k)\\ &= 1 - \frac{2}{\binom{n}{3}}\left(\frac{n^2(n-1)}{2} - \frac{(n-1)n(2n-1)}{6} + \frac{n(n-1)}{2}\right)\\ &= 1 - 2\frac{3\cdot 2}{n(n-1)(n-2)}\left(\frac{1}{6}(n^3-3n^2+2n)\right)\\ &= -1. \end{align*}
  • Для $c_n = n$:
\begin{align*} &\phantom{=\ }n - \frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}(nk^2 - k^3 - k^2)\\ &= n - \frac{2}{\binom{n}{3}}\left((n-1)\frac{(n-1)n(2n-1)}{6} - \frac{(n-1)^2n^2}{4}\right)\\ &= n - 2\frac{3\cdot 2}{n(n-1)(n-2)}\frac{1}{12}n(n-1)^2(n-2)\\ &= 1. \end{align*}
  • Для $c_n = \sum_{k=1}^n H_n$:
\begin{align*} &\phantom{=\ }n - \frac{2}{\binom{n}{3}}\sum_{k=0}^{n-1}((n-1)k - k^2)\sum_{i=1}^k H_i\\ \end{align*} Не завершено.
Вправа 2.47
Ми бачимо, що $$\frac{2}{n+1} < a_n < \frac{2}{n}.$$ Тепер підставимо значення для $a_{n-1}$ в рекуренцію: \begin{align*} \frac{2}{n+\frac{2}{n}} &< a_n < \frac{2}{n + \frac{2}{n-1}}\\ \frac{2n}{n^2+2} &< a_n < \frac{2(n-1)}{n^2 - n + 2}\\ \end{align*} Але $\frac{2}{n+1} < \frac{2n}{n^2+2}$ і $\frac{2(n-1)}{n^2 - n + 2} < \frac{2}{n}$, тобто завдяки бутстрепу вдалось звузити діапазон можливих начень $a_n$.
Вправи з 48 і далі нерозв'язані.

вівторок, 17 квітня 2018 р.

Розв'язки до вправ з "An Introduction to the Analysis of Algorithms" Роберта Седжвіка. Розділ I.

Вправа 1.4
Зауважимо одразу, що $C_1 = 0, C_2 = 2$. \begin{align*} C_{N+1} - C_N &= C_{\lceil \frac{N+1}{2}\rceil} + C_{\lfloor \frac{N+1}{2}\rfloor} + N + 1 - \left(C_{\lceil\frac{N}{2}\rceil} + C_{\lfloor\frac{N}{2}\rfloor} + N\right)\\ &=C_{\lceil \frac{N+1}{2}\rceil} - C_{\lfloor \frac{N}{2}\rfloor} + 1\\ &=\lfloor\lg N\rfloor + 2 \end{align*} Тут $\lfloor\lg N\rfloor$ - це сума одиничок. Вона дорівнює кількості разів, що ми маємо поділити $N$ на $2$ із взяттям найближчого меншого цілого, тобто ось за такою формулою: $$ \left\lfloor\begin{array}{@{} c @{}} \frac{ \substack{ \left\lfloor\begin{array}{@{} c @{}} \frac{ \left\lfloor\begin{array}{@{} c @{}}\frac{\lfloor \frac{N}{2}\rfloor}{2}\end{array}\right\rfloor }{2} \end{array}\right\rfloor \\ \\ \\ \dots } }{2} \end{array}\right\rfloor = 1. $$ Навіть якщо припустити, що ми відкидаємо $0.5$ після кожного ділення, то зворотній процес виглядає так: \[ (\cdots(((1\cdot 2+1)\cdot2+1)\cdot2+1)\cdots) = 2^m + 2^{m-1} + \dots < 2^{m+1} \] Отже, $\lfloor\lg N\rfloor = m$. $2$ з'являється із $(C_2 - C_1)$ - остання ітерація. Тепер ми можемо записати $C_{N}$ як телескопічний ряд: \begin{equation*} C_{N} = \sum_{k=1}^{N-1} \left(C_{k+1} - C_k\right) = \sum_{k=1}^{N-1}(\lfloor\lg k\rfloor + 2). \end{equation*}
Вправа 1.5
Перевіримо нерівність для $N = 1$: \[C_1 = 1\cdot\lceil \lg 1\rceil + 1 - 2^{\lceil \lg 1\rceil} = 0.\] Тепер перевіримо, що якщо нерівність виконується для $N$, то вона виконується і для $N+1$. Тут у нас є два випадки:
  1. $\lceil \lg (N)\rceil = \lceil \lg (N+1)\rceil$, тобто $N \ne 2^n$. Тоді, \begin{align*} C_{N+1} &= (N + 1)\cdot\lceil \lg (N + 1)\rceil + N + 1 - 2^{\lceil \lg (N + 1)\rceil}\\ &= N\cdot\lceil \lg N\rceil + \lceil \lg N\rceil + N + 1 - 2^{\lceil \lg N\rceil}\\ &= C_N + \lceil \lg N\rceil + 1\\ \end{align*} Тут $\lceil \lg N\rceil + 1 = \lfloor \lg N\rfloor + 2$ і використовуючи результат вправи 1.4 рівність доведена.
  2. $N = 2^n$, \begin{align*} C_{N+1} &= (N + 1)\cdot\lceil \lg (N + 1)\rceil + N + 1 - 2^{\lceil \lg (N + 1)\rceil}\\ &= (N + 1)\cdot (\lg N + 1) + N + 1 - 2\cdot N\\ &= C_N + \lg N + 2\\ \end{align*} Знов-таки, використовуючи вправу 1.4 завершуємо доведення.
Вправа 1.7
Розглянемо варіант зі складністю $N\lg N$: \begin{align*} t_{2N} &= 2 t_N \left(1 + 1 / \lg N\right)\\ &=2N(\lg N + 1)\\ &=2N\lg (2N). \end{align*}
Вправа 1.8
void merge_sort_internal(
    const vector<unsigned>::iterator a_beg,
    const vector<unsigned>::iterator a_end,
    const vector<unsigned>::iterator b_beg,
    const vector<unsigned>::iterator b_end
) {
    assert(distance(a_beg, a_end) == distance(b_beg, b_end));
    const auto dist = distance(a_beg, a_end);
    if (dist <= 1) return;
    const auto a_mid = a_beg + (dist + 1) / 2;
    const auto b_mid = b_beg + (dist + 1) / 2;
    merge_sort_internal(a_beg, a_mid, b_beg, b_mid);
    merge_sort_internal(a_mid, a_end, b_mid, b_end);
    copy(a_beg, a_end, b_beg);
    auto a = a_beg;
    auto b1 = b_beg;
    auto b2 = b_mid;
    while (b1 < b_mid && b2 < b_end) {
        *a++ = *b1 < *b2 ? *b1++ : *b2++;
    }
    copy(b1, b_mid, a);
    copy(b2, b_end, a);
}

void merge_sort(vector<unsigned>& data) {
    vector<unsigned> buffer(data.size());
    merge_sort_internal(data.begin(), data.end(), buffer.begin(), buffer.end());
}

long long test(size_t size) {
    random_device rnd_device;
    minstd_rand lce(rnd_device());
    uniform_int_distribution<unsigned> dist(1, numeric_limits<unsigned>::max());
    vector<unsigned> v(size);
    generate(v.begin(), v.end(), [&dist, &lce]() { return dist(lce); });
    auto start = chrono::system_clock::now();
    merge_sort(v);
    auto end = chrono::system_clock::now();
    auto elapsed = chrono::duration_cast<chrono::microseconds>(end - start);
    return elapsed.count();
}

int main() {
    long long t_10_6 = test(1000 * 1000);
    long long t_10_7 = test(10 * 1000 * 1000);
    long long t_10_8 = test(100 * 1000 * 1000);
    double ratio = 10. * (1. + 1. / log2(10));
    cout << t_10_6 << ", " << t_10_7 << ", " << t_10_8 << endl;
    cout << "expected ratio:" << ratio << endl;
    cout << "actual ratio 10^7 to 10^6:" << t_10_7 / static_cast<double>(t_10_6) << endl;
    cout << "actual ratio 10^8 to 10^7:" << t_10_8 / static_cast<double>(t_10_7) << endl;
    system("pause");
    return 0;
}

161000, 1776000, 20196000
expected ratio:13.0103
actual ratio 10^7 to 10^6:11.0311
actual ratio 10^8 to 10^7:11.3716

Вправи 1.9, 1.10
Використайте сортування підрахунком.
Вправи 1.11
Тут можна використати той самий підхід, що й у теоремі 1.3. $(N-1)! - $ це кількість можливих перестановок елементів для кожного розбиття. $(N+1) - $ кількість порівнянь. \[ C'_N = \sum_{i=1}^N\left(C'_{i-1}+C'_{N-i}\right)+(N-1)!(N+1) \]
Вправи 1.12
Вибираючи елементи в праву і ліві частині розбиття ми порівнюємо елементи лише з опірним елементом $v$ і ніколи між собою. Через ми ніяк не впорядковуємо елементи в кожній з частин. Якщо замість j = hi; ми напишемо j = hi + 1;, то першим же обміном в циклі while (true) ми перенесемо елемент $v$ в ліву частину розбиття, і всі елементи перед ним менші ніж цей елемент. Наприклад, \[3,5,2,11,7,10,12,9\] розбиваємо на \[3,5,2,9,7,10\quad\mbox{і}\quad12\] і в лівій частині ми бачимо, що до першої $9$ всі елементи менше ніж $9$.
Вправа 1.13
Під час розбиття ми не робимо порівнянь між елементами лівого і правого підмасивів, таким чином ніяк не впорядковоючи їх, тому ці підмасиви також випадкові перестановки. Якщо ж ми ініціалізуємо правий вказівник j = hi + 1, то умова while (v < a[--j]) завжди хиба на першій ітерації зовнішнього циклу. Через це ми точно знаємо, що до елемента зі значенням $v$ всі елементи менші ніж він, а за ним вже можуть бути й більші. Тобто перестановка невипадкова. $$1,3,9,5,7,2,4$$ $$1,3,4,5,7,2\quad9$$ Тут $v = 4$, зауважте, що ліворуч всі елементи до $4$ менші ніж $4$.
Вправа 1.14
\begin{align*} A_N N &= N + 2\sum_{1\le j\le N} A_{j-1}\\ A_{N-1}(N-1) &= N-1 + 2\sum_{1\le j\le N-1} A_{j-1}. \end{align*} Віднімемо друге рівняння від першого. \begin{align*} A_N N - A_{N-1}(N-1) &= 1 + 2A_{N-1}\\ A_N N - A_{N-1}(N+1) &= 1 \end{align*} Розділемо на $N(N+1)$. \begin{align*} \frac{A_N}{N+1} - \frac{A_{N-1}}{N} = \frac{1}{N(N+1)} \end{align*} Тепер двічі використаємо прийом телескопічного ряду. \begin{align*} \frac{A_N}{N+1} - \frac{A_0}{1} &= \sum_1^N\frac{1}{k(k+1)}\\ \frac{A_N}{N+1} - \frac{A_0}{1} &= \sum_1^N\left(\frac{1}{k}-\frac{1}{k+1}\right)\\ \frac{A_N}{N+1} - \frac{A_0}{1} &= \frac{1}{1}-\frac{1}{N+1}\\ \end{align*} Якщо $A_0 = 0$, то \begin{equation*} A_N = N \end{equation*} Дуже цікавий результат, і його нескладно перевірити.
Вправа 1.15
Нас цікавить скільки обмінів відбулось в тілі цього циклу:
while (true)
{
    while (a[++i] < v) ;
    while (v < a[--j]) if (j == lo) break;
    if (i >= j) break;
    t = a[i]; a[i] = a[j]; a[j] = t;
}
Фактично нам потрібно знайти математичне сподівання кількості обмінів. На останню позицію в масиві рівноймовірно може потрапити будь-який елемент. Позначимо його позицію в відсортованому масиві через $p$, а сам елемент через $v$. Отже, у нас є $N-p$ елементів більше ніж $v$ і $p-1$ елементів менше ніж $v$. Імовірність натрапити на елемент менший ніж $v$ становить $\frac{p-1}{N}$ і більший - $\frac{N-p}{N}$. Відповідно, щоб зустріти елемент більший ніж $v$ ми в середньому пройдемо $\frac{N}{N-p}$, а менший - $\frac{N}{p-1}$. В цьому циклі ми маємо зробити $N$ порівнянь, останній елемент порівнюється двічі. Отже, \begin{align*} \frac{1}{N}\sum_{p=1}^N\frac{N}{\frac{N}{N-p}+\frac{N}{p-1}} &= \sum_{p=1}^N\frac{1}{\frac{N(N-1)}{(N-p)(p-1)}}\\ &=\frac{1}{N(N-1)}\sum_{p=1}^N(Np - N - P^2 +p)\\ &=\frac{1}{N-1}\left(-N + \frac{N(N+1)}{2} - \frac{1}{N}\sum_{p=1}^N p^2 + \frac{N+1}{2}\right)\\ &=\frac{1}{N-1}\left(-N + \frac{N(N+1)}{2} - \frac{(N+1)(2N+1)}{6} + \frac{N+1}{2}\right)\\ &=\frac{1}{N-1}\frac{N^2-3N+2}{6}\\ &=\frac{N-2}{6}. \end{align*}
Вправа 1.16
\begin{align*} T_N &= \frac{1}{N}\sum_{k=1}^N\left(T_{k-1}+T_{N-k}\right)\\ & = \frac{2}{N}\sum_{k=0}^{N-1}T_k \end{align*} Тепер розглянувши $NT_N - (N-1)T_{N-1}$ ми отримуємо: \begin{align*} T_N &= \frac{N+1}{N}T_{N-1}\\ &= \frac{N+1}{N}\cdots\frac{5}{4}T_3\\ &= \frac{N+1}{6} \end{align*} Базовий випадок: \[ T_3 = 1/3(1 + 1) = 2/3 \]
Вправа 1.17
\begin{align*} C_N &= N+1 + \frac{1}{N}\sum_{j=1}^{N}(C_{j-1}+C_{N-j})\\ &= N+1 + \frac{2}{N}\sum_{j=0}^{N-1}C_j\\ &= N+1 + \frac{2}{N}\left(\sum_{j=M+1}^{N-1}C_j + \frac{1}{4}\sum_{1}^M j(j-1)\right) \end{align*} \begin{align*} NC_N - (N-1)C_{N-1} &= (N+1)N - N(N-1) + 2C_{N-1}\\ NC_N - (N+1)C_{N-1} &= 2N \end{align*} Використовуємо принцип телескопа складаємо всі рівності від $\frac{C_N}{N+1} - \frac{C_{N-1}}{N}$ до $\frac{C_{M+1}}{M+2} - \frac{\frac{1}{4}M(M-1)}{M+1}$: \[ \frac{C_N}{N+1} - \frac{1}{4}\frac{M(M-1)}{M+1} = 2 (H_{N+1} - H_M) \] \[ C_N = \frac{1}{4}(N+1)\frac{M(M-1)}{M+1} + 2(N+1)(H_{N+1} - H_M) \]
Вправа 1.18
\begin{align*} C_N &\approx \frac{1}{4}N\frac{M(M-1)}{M+1} + 2N(\ln(N+1) - \ln M)\\ &\approx 2N\ln N + N\left(\frac{1}{4}\frac{M(M-1)}{M+1} - 2 \ln M\right) \end{align*} Мінімуму функція сягає коли $M = 8$.
Вправа 1.19
Кількість порівнянь, яку виконує швидке сортування $\approx 2N\ln N - 1.846N$. При переході з $17$ на $18$, $f(M)$ перестрибує $-1.846$.
Вправа 1.20
Це те саме, що запитати скільки урн міститиме більш ніж одну кулю якщо покласти $10^4$ куль в $10^6$ урн. Досить просто це можна порахувати за допомогою індикаторних змінних. Нехай $X_i=1$ якщо $i$-та урна містить більше ніж одну кулю і $0$ якщо ні. Тоді $E[X_i]=1−P_0−P_1$, де $P_J$ позначає ймовірність, що ця урна містить рівно $j$ куль. \begin{align*} P_0 &= \left(\frac{U-1}{U}\right)^B,\\ P_1 &= U\left(\frac{U-1}{U}\right)^{B-1}\left(1-\frac{U-1}{U}\right). \end{align*} Розрахунок в Matlab:
>> U = 10^6;
>> B = 10^4;
>> p = (U-1)/U;
>> (1 - p^B - B*p^(B-1)*(1-p))*U
ans =  49.663
І невеличка симуляція:
    const int U = 1000 * 1000;
    const int B = 10 * 1000;
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<> dis(0, U-1);
    double mean = 0;
    const int roundsCount = 100;
    for (int round = 0; round < roundsCount; ++round) {
        unordered_map<int, int> file;
        for (int i = 0; i < B; ++i) {
            ++file[dis(gen)];
        }
        int duplicated = 0;
        for (auto& el : file) {
            if (el.second > 1) {
                ++duplicated;
            }
        }
        mean += duplicated / static_cast<double>(roundsCount);
    }
    cout << mean << endl;
Вивід не значно відрізняється від отриманого теоретичного результату: 48.87
Вправа 1.21
Чим з меншого діапазону вибираються елементи масиву, тим швидше працює швидке сортування. В крайньому випадку, якщо всі елементи однакові, то розбиття відбувається навпіл і алгоритм працює найшвидше.
Вправа 1.22
Вправа 1.23
Реалізація наведена в цій частині завжди ділить навпіл, і виконує порівняння лише під час злиття і кількість порівнянь не залежить від входових даних. Тому дисперсія дорівнює $0$.

неділя, 11 березня 2018 р.

Міра Лебега

Означення
Сім'я множин $\mathfrak{R}$ називається кільцем якщо $A,B \in \mathfrak{R}$ тягне за собою \begin{equation} A\cup B \in \mathfrak{R},\quad A - B \in \mathfrak{R}. \end{equation} Також, якщо $\mathfrak{R}$ - це кільце, то $A\cap B \in \mathfrak{R}$, бо $A\cap B = A - (A - B)$.
Означення
Кільце $\mathfrak{R}$ називається $\sigma-$кільцем, якщо для будь-яких $A_n \in \mathfrak{R} (n = 1,2,3,\dots)$ виконується \begin{equation} \cup_{n=1}^{\infty} A_n \in \mathfrak{R}. \end{equation} Якщо $A_n \in \mathfrak{R}, i = 1,2,3, \dots$, то \begin{equation*} \cap_{n=1}^{\infty} A_n \in \mathfrak{R}, \end{equation*} тому що \begin{equation*} \cap_{n=1}^{\infty} A_n = A_1 - \cup_{n=1}^{\infty}(A_1-A_n). \end{equation*}
Означення
$\phi$ - це функція множини визначена на $\mathfrak{R}$ якщо $\phi$ співставляє кожному $A \in \mathfrak{R}$ число $\phi(A)$ з розширеної системи дійсних чисел. $\phi$ адитивна, якщо $A\cap B = 0$ тягне, що \begin{equation} \phi(A\cup B) = \phi(A) + \phi(B).\label{eq:phi_additive} \end{equation} $\phi$ зліченно адитивна, якщо $A_i\cap A_j = 0$ тягне, що \begin{equation} \phi(\cup_{n=1}^{\infty} A_n) = \sum_{n=1}^{\infty}\phi(A_n).\label{eq:phi_count_add} \end{equation}
Ми вважатимемо, що область значень $\phi$ не включає одночасно $+\infty$ і $-\infty$, інакше правий бік \eqref{eq:phi_additive} безсенсовий.

Зауважимо, що лівий бік \eqref{eq:phi_count_add} не залежить від того як впорядковані $A_n$, отже, згідно з теоремою Рімана про умовно збіжний ряд, якщо правий бік збігається, то він збігається абсолютно.

Для адитивної $\phi$ легко довести такі властивості: \begin{equation} \phi(0)=0.\label{eq:prop_add_zero} \end{equation} Якщо $A_i \cap A_j = 0$ при $i\neq j$ \begin{equation} \phi(A_1\cup\cdots A_n)=\phi(A_1)+\cdots+\phi(A_n).\label{eq:prop_add_finite} \end{equation} \begin{equation} \phi(A \cup B) + \phi(A \cap B) = \phi(A) + \phi(B).\label{eq:prop_add_cupcap} \end{equation} Якщо $\phi(A)\ge 0$ для всіх $A$ і $A_1 \subset A_2$, тоді \begin{equation} \phi(A_1)\le \phi(A_2). \label{eq:prop_add_mono} \end{equation} Завдяки \eqref{eq:prop_add_mono} невід'ємні адитивні функці множин часто називають монотонними.
Якщо $b \subset A$, і $|\phi(B)|<+\infty$, то \begin{equation} \phi(A-B)=\phi(A) - \phi(B) \end{equation}.
Теорема
Припустимо, що $\phi$ зліченно-адитивна на кільці $\mathfrak{R}$. Припустимо також, що $A_n\in \mathfrak{R}, n = 1,2,3,\dots, A_1\subset A_2\subset A_3\subset \cdots$, $A \in \mathfrak{R}$, і \begin{equation*} A = \cup_{n=1}^{\infty}A_n. \end{equation*} Тоді, як $n \to \infty$, \begin{equation*} \phi(A_n)\to\phi(A). \end{equation*}
Позначимо $C_1 = A_1$ і $C_i = A_{i} - A_{i-1}$, тоді $C_i \cap C_j = 0$ для $i\neq j$ і $A_n = \cup_{i=1}^n C_i$. Згідно з \eqref{eq:phi_count_add} маємо, що $$ \phi(A_n) = \sum_{i=1}^{n} \phi(C_i) \quad \mbox{i}\quad \phi(A) = \sum_{i=1}^{\infty} \phi(C_i) $$

Побудова міри Лебега

Означення
Нехай $R^p$ позначає $p-$вимірний евклідовий простір. Тоді інтервал в $R^p$ - це множина точок $\mathbf{x}=(x_1,\dots,x_p)$ таких, що \begin{equation} a_i\le x_i\le b_i \quad i=1,\dots,p,\label{eq:interval} \end{equation} або множина точок охарактеризована \eqref{eq:interval} з деякими чи усіма $\le$ замінинеми на $<$. Порожня множина такоє інтервал.

Якщо $A$ - це об'єднання скінченної кількості інтервалів, то $A$ називаємо елементарною множиною.

Якщо $I$ - це інтервал, ми визначимо $$m(I) = \prod_{i=1}^p(b_i-a_i).$$ Якщо $A=I_1\cup\cdots \cup I_n$ і якщо ці інтервали попарно неперетинні, то $$m(A) = m(I_1)+\cdots+m(I_n).$$ $\mathcal{E}$ позначає сім'ю всіх елементарних підмножин $R^p$.

Означення
Невід'ємна додатня функція множини $\phi$ визначена на $\mathcal{E}$ регулярна, якщо: для кожного $A \in \mathcal{E}$ і для кожного $\varepsilon > 0$ існують множини $F, G \in \mathcal{E}$ такі, що $F$ - замкнена, $G$ - відкрита, $F\subset A\subset G$ і \begin{equation} \phi(G)-\varepsilon\le\phi(A)\le\phi(F)+\varepsilon.\label{eq:setfun_reg} \end{equation}
Означення
Нехай $\mu$ адитивна, регулярна, невід'ємна і скінченна на $\mathcal{E}$. Розглянемо зліченне покриття будь-якої множини $E\subset R^p$ відкритими елементарними множинами $A_n$: $$E\subset\cup_{n=1}^{\infty}A_n.$$ Визначимо $$\mu^*(E) = \inf\sum_{n=1}^{\infty}\mu(A_n),$$ $\inf$ береться по всіх можливих покриттях $E$ відкритими елементарними множинами. $\mu^*(E)$ називається зовнішня міра $E$, для відповідного $\mu$.
Очевидно, що $\mu^*(E)\ge 0$ для всіх $E$ і якщо $E_1 \subset E_2$, то $$\mu^*(E_1) < \mu^*(E_2).$$
Теорема
  1. \begin{equation} \forall A\in\mathcal{E},\mu^*(A)=\mu(A).\label{eq:mu_equal_elemset} \end{equation}
  2. Якщо $E = \cup_1^{\infty}E_n$, тоді \begin{equation} \mu^*(E)\le\sum_{n=1}^{\infty}\mu^*(E_n).\label{eq:single_le_cup} \end{equation}
Означення
Для будь-яких $A \subset R^p, B \subset R^p$, визначимо \begin{align} S(A,B) = (A-B)\cup(B-A),\label{eq:def_symdiff}\\ d(A,B) = \mu^*(S(A,B)).\label{eq:def_dist} \end{align}
Ми записуватимемо $A_n\to A$ якщо $\lim_{n\to\infty}d(A,A_n) = 0.$

Якщо існує послідовність $\{A_n\}$ елементарних множин таких, що $A_n \to A$, то ми кажемо, що $A$ скінченно $\mu-$вимірна і пишемо $A \in \mathfrak{M}_F(\mu)$.

Якщо $A$ - це об'єднання зліченної колекції скінченно $\mu-$вимірних множин, ми кажемо, що $A$ - це $\mu-$вимірна і пишемо $A \in \mathfrak{M}(\mu)$.

$S(A,B)$ - це так звана симетрична різниця $A$ і $B$. Далі буде видно, що $d(A,B)$ по суті є функцією відстані.

Перед наступною теоремою нам потрібно розробити декілька властивостей $S(A, B)$ і $d(A, B)$. \begin{equation} S(A,B) = S(B,A),\quad S(A,A) = 0.\label{eq:s_p_sym} \end{equation} \begin{equation} S(A,B) \subset S(A,C)\cup S(C,B).\label{eq:s_p_tri} \end{equation} \begin{equation} S(A_1,B_1)\cup S(A_2,B_2)\supset \begin{cases} S(A_1\cup A_2, B_1\cup B_2)\\ S(A_1\cap A_2, B_1\cap B_2)\\ S(A_1 - A_2, B_1 - B_2). \end{cases} \label{eq:s_p_cases} \end{equation} \begin{equation} d(A,B) = d(B,A),\quad d(A,A) = 0.\label{eq:d_p_sym} \end{equation} \begin{equation} d(A,B) \le d(A,C) + d(C,B).\label{eq:d_p_tri} \end{equation} \begin{equation} d(A_1,B_1)+d(A_2,B_2)\ge \begin{cases} d(A_1\cup A_2, B_1\cup B_2)\\ d(A_1\cap A_2, B_1\cap B_2)\\ d(A_1 - A_2, B_1 - B_2). \end{cases} \label{eq:d_p_cases} \end{equation} І якщо хоча б одни з $\mu^*(A), \mu^*(B)$ - скінченне, то \begin{equation} |\mu^*(A) - \mu^*(B)|\le d(A,B).\label{eq:d_p_mudiff} \end{equation}
Теорема
$\mathfrak{M}(\mu)$ - це $\sigma-$кільце і $\mu^*$ зліченно-адитивна на $\mathfrak{M}(\mu)$.
Ми доводитемо теорему в чотири кроки: $\mathfrak{M}_F(\mu)$ - кільце, $\mu^*$ - адитивна, $\mu^*$ - зліченно-адитивна, $\mathfrak{M}(\mu)$ - $\sigma-$кільце.
  1. Розглянемо дві множини $A, B \in \mathfrak{M}_F(\mu)$. Це означає, що існують послідовності $\{A_n\}, \{B_n\}$ елементарних множин, такі що $A_n\to A, B_n\to B$. Але з \eqref{eq:d_p_cases} випливає, що $$d(A_n\cup B_n, A\cup B) \le d(A_n, A) + d(B_n, B),$$ тобто $A_n\cup B_n \to A\cup B$. Аналогічно, $A_n - B_n \to A - B$.
  2. Згідно з \eqref{eq:prop_add_cupcap}, $$\mu(A_n) + \mu(B_n) = \mu(A_n\cup B_n) + \mu(A_n\cap B_n).$$ З \eqref{eq:d_p_mudiff} $\mu^*(A_n)\to\mu^*(A)$. Тоді з \eqref{eq:mu_equal_elemset} випливає, що $$\mu^*(A) + \mu^*(B) = \mu^*(A\cup B) + \mu^*(A\cap B).$$ Якщо $A \cap B = 0$, то $\mu^*(A\cap B) = 0$, значить $\mu^*$ адитивна на $\mathfrak{M}_F(\mu)$.
  3. Припустимо, що $A \in \mathfrak{M}(\mu)$ і ми маємо покриття $A = \cap_{n=1}^{\infty} A'_n$, де $A'_n \in \mathfrak{M}_F$. Нам потрібно представити $A$ як зліченне об'єднання неперетинних множин з $\mathfrak{M}_F$. Запишемо $A_1 = A'_1$, $A_n = (A'_1 \cup \cdots \cup A'_n) - (A'_1 \cup \cdots \cup A'_{n-1})$. Ці множини взаємо-неперетинні. І тепер ми можемо використати \eqref{eq:single_le_cup}: $$\mu^*(A) \le \sum_{n=1}^{\infty}A_n.$$ Також $A \supset A_1 \cup \cdots \cup A_n$, тому $\mu^*(A) \ge \mu^*(A_1) + \cdots + \mu^*(A_n)$. Отже, $$\mu^*(A) = \sum_{n=1}^{\infty} \mu^*(A_n).$$ Але цього не досить, бо $A \in \mathfrak{M}(\mu)$, а $A_n \in \mathfrak{M}_F(\mu)$. Нехай $\mu^*(A) < \infty$. Позначимо $B_n = \cup_{k=1}^{n}A_k$. Тоді $\lim_{n\to\infty} B_n = A$ і $B_n \in \mathfrak{M}_F(\mu)$, отже $A \in \mathfrak{M}_F$. Тепер видно, що $\mu^*$ зліченно-адитивна на $\mathfrak{M}(\mu)$.
  4. Тут нам потрібно довести, що зліченне об'єднання і різниця множин з $\mathfrak{M}(\mu)$ є множинами з $\mathfrak{M}(\mu)$.
    • Зліченне об'єднання зліченної колекції скінченно $\mu-$вимірних множин, знов є зліченною колекцією скінченно $\mu-$вимірних множин. (Див. як порахувати раціональні числа)
    • В пункті 3 ми довели, що якщо $A \in \mathfrak{M}(\mu)$ і $\mu^*(A) < \infty$, то $A \in \mathfrak{M}_F(\mu)$.
      Розглянемо $A, B \in \mathfrak{M}(\mu)$, $A = \cup_{n=1}^{\infty} A_n$ і $B = \cup_{n=1}^{\infty} B_n$, де $A_n, B_n \in \mathfrak{M}_F(\mu)$. Маємо, що $$A_n \supset A_n \cap B = \cup_{i=1}^{\infty} (A_n \cap B_i) \in \mathfrak{M}(\mu).$$ Тому $\mu^*(A_n \cap B) \le \mu^*(A_n) < +\infty$, отже $A_n \cup B \in \mathfrak{M}_F(\mu)$. А їх зліченне об'єднання в $\mathfrak{M}(\mu)$. $$A-B = \cup_{n=1}^{\infty} (A_n \cap B).$$
Тепер ми заміняємо $\mu^*(A)$ на $\mu(A)$, якщо $A \in \mathfrak{M}(\mu)$. Отже, ми розширили $\mu$ з $\mathcal{E}$ до зліченно-адитивної функції на $\sigma-$кільці $\mathfrak{M}(\mu)$. Таку розширену функцію називають мірою. Особливий випадок коли $\mu = m$ називають мірою Лебега на $R^p$.

неділя, 21 січня 2018 р.

Збірка широковживаних тестів значущості пов'язаних з нормальним розподілом

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

$z$-тест

  • Використання: Чи рівні середні значення популяції і гіпотетичне середнє.
  • Дані: $x_1, x_2, \dots, x_n$.
  • Припущення: Дані - це незалежні нормальні проби: $x_i \sim N(\mu, \sigma^2)$, де $\mu$ невідоме, а $\sigma$ - відоме.
  • $H_0$: для певного $\mu_0$, $\mu = \mu_0$.
  • $H_A$: $\mu \ne \mu_0, \mu > \mu_0, \mu < \mu_0$.
  • Тестова статистика: $z = \frac{\bar{x}-\mu_0}{\sigma/\sqrt{n}}$
  • Нульовий розподіл: $f(z|H_0)$ - це густина ймовіності $Z\sim N(0,1)$.
  • $p$-значення:
    Двостороння:$p = P(|Z| > z) = 2 * (1 - \mbox{pnorm}(|z|, 0, 1))$
    Одностороння-більше:$p = P(Z > z) = 1 - \mbox{pnorm}(z, 0, 1)$
    Одностороння-менше:$p = P(Z > z) = \mbox{pnorm}(z, 0, 1)$

Одновибірковий $t$-тест

  • Використання: Чи рівні середні значення популяції і гіпотетичне середнє.
  • Дані: $x_1, x_2, \dots, x_n$.
  • Припущення: Дані - це незалежні нормальні проби: $x_i \sim N(\mu, \sigma^2)$, де $\mu$ і $\sigma$ невідомі.
  • $H_0$: для певного $\mu_0$, $\mu = \mu_0$.
  • $H_A$: $\mu \ne \mu_0, \mu > \mu_0, \mu < \mu_0$.
  • Тестова статистика: $t = \frac{\bar{x}-\mu_0}{s/\sqrt{n}}$,
    де $s^2$ - це дисперсія вибірки: $s^2 = \frac {1}{n-1} \sum_{i=1}^n \left(x_i - \overline{x} \right)^ 2$
  • Нульовий розподіл: $f(t|H_0)$ - це густина ймовіності $T\sim t(n-1)$. (t-розподіл Стьюдента з $n-1$ ступенем свободи)
  • $p$-значення:
    Двостороння:$p = P(|Z| > z) = 2 * (1 - \mbox{pt}(|z|, n-1))$
    Одностороння-більше:$p = P(Z > z) = 1 - \mbox{pt}(z, n-1)$
    Одностороння-менше:$p = P(Z > z) = \mbox{pt}(z, n-1)$

Двовибірковий $t$-тест

Випадок рівних дисперсій

  • Використання: Чи різняться середні значення двох популяцій на гіпотетичну величину.
  • Дані: $x_1, x_2, \dots, x_n$ і $y_1, y_2, \cdots, y_m$.
  • Припущення: Дані - це незалежні нормальні проби: $x_i \sim N(\mu_x, \sigma^2)$, $y_i \sim N(\mu_y, \sigma^2)$, де $\mu_x$ і $\mu_y$ невідомі, можливо різні і $\sigma$ також невідома.
  • $H_0$: для певного $\mu_0$, $\mu_x - \mu_y = \mu_0$.
  • $H_A$: $\mu_x - \mu_y \ne \mu_0, \mu_x - \mu_y > \mu_0, \mu_x - \mu_y < \mu_0$.
  • Тестова статистика: $z = \frac{\bar{x}-\bar{y} - \mu_0}{s_{\bar{x}-\bar{y}}}$,
    де $s_x^2, s_y^2$ - це дисперсі] двох вибірок, а $s_{\bar{x}-\bar{y}}$ - оцінкова дисперсія: \begin{equation} s_{\bar{x}-\bar{y}} = s_P\left(\frac{1}{n} + \frac{1}{m}\right)\label{eq:sdiff} \end{equation} де $s_P$ - об'єднана (pooled) дисперсія (cереднє зважене): \begin{equation} s_P = \frac{(n-1)s_x^2 + (m-1)s_y^2}{n+m-2} \end{equation} де кількість ступенів свободи $df = n+m-2$.
    Множник $\frac{1}{n}+\frac{1}{m}$ в \eqref{eq:sdiff} випливає з того факту, що $\mbox{Var}(\bar{x}) = \frac{\sigma_x^2}{n}$ і $\mbox{Var}(\bar{x}-\bar{y}) = \frac{\sigma_x^2}{n} + \frac{\sigma_y^2}{m} = \sigma^2\left(\frac{1}{n}+\frac{1}{m}\right)$, бо $\sigma_x = \sigma_y = \sigma$.
  • Нульовий розподіл: $f(t|H_0)$ - це густина ймовіності $T\sim t(df)$.
  • $p$-значення:
    Двостороння:$p = P(|Z| > z) = 2 * (1 - \mbox{pt}(|z|, n-1))$
    Одностороння-більше:$p = P(Z > z) = 1 - \mbox{pt}(z, n-1)$
    Одностороння-менше:$p = P(Z > z) = \mbox{pt}(z, n-1)$

Випадок відмінних дисперсій

Цей випадок майже повністю повторює випадок з рівними дисперсіями, нам лише треба змінити наші припущення і формулу для об'єднаної дисперсії.
  • Використання: Чи різняться середні значення двох популяцій на гіпотетичну величину.
  • Дані: $x_1, x_2, \dots, x_n$ і $y_1, y_2, \cdots, y_m$.
  • Припущення: Дані - це незалежні нормальні проби: $x_i \sim N(\mu_x, \sigma^2)$, $y_i \sim N(\mu_y, \sigma^2)$, де $\mu_x$ і $\mu_y$ невідомі і різні і $\sigma$ також невідома.
  • $H_0$: для певного $\mu_0$, $\mu_x - \mu_y = \mu_0$.
  • $H_A$: $\mu_x - \mu_y \ne \mu_0, \mu_x - \mu_y > \mu_0, \mu_x - \mu_y < \mu_0$.
  • Тестова статистика: $z = \frac{\bar{x}-\bar{y} - \mu_0}{s_{\bar{x}-\bar{y}}}$, де \begin{equation} s_{\bar{x}-\bar{y}} = \frac{s_x^2}{n} + \frac{s_y^2}{m},\ df = \frac{(s_x^2/n+s_y^2/m)^2}{\frac{(s_x^2/n)^2}{n-1}+\frac{(s_y^2/m)^2}{m-1}} \end{equation}
  • Нульовий розподіл: $f(t|H_0)$ - це густина ймовіності $T\sim t(df)$.
  • $p$-значення:
    Двостороння:$p = P(|Z| > z) = 2 * (1 - \mbox{pt}(|z|, n-1))$
    Одностороння-більше:$p = P(Z > z) = 1 - \mbox{pt}(z, n-1)$
    Одностороння-менше:$p = P(Z > z) = \mbox{pt}(z, n-1)$

субота, 6 січня 2018 р.

Арифметика рухомої коми

Відносні помилки і ulp'и

Обчисленням з рухомою комою властиві помилки заокруглення, тому важливо мати змогу вимірювати їх. Розглянемо формат рухомої коми з основою ($\beta$) 10 і мантисою ($p$) 4. Якщо результат обчислень з рухомою комою $3.143\times10^{-1}$,а результат із використанням нескінченної точності $0.3141$, тоді помилка рівна 2 одиницям в останній позиції (units in last place - ulp). \begin{equation} 1\mbox{ulp} = \beta\times \beta^{-p}\label{eq:ulp} \end{equation}

Якщо результат обчислень з рухомою комою лежить найближче до справжнього результата, то помилка все ще може становити $0.5$ ulp. Якщо дійсне число наближене за допомогою числа з рухомою комою, то помилка може бути завбільшки з $0.0005$ ($p=4, \beta=10$), або $\beta/2\times \beta^{-p}$.

Щоб обчислити відносну помилку, що відповідає $0.5$ ulp, зауважимо, що коли дійсне число наближається за допомогою числа з рухомою точкою $$d.\underbrace{d\dots d}_{p \mbox{ раз}}\times \beta^e,$$ то помилка може бути завбільшки з $$0.\underbrace{0\dots 0}_{p \mbox{ раз}}\beta'\times\beta^e,$$ де $\beta' = \beta/2$. Ця помилка рівна $(\beta/2)\times\beta^{-p}\times\beta^e$. Числа у вигляді $d.d\dots d\times\beta^e$ набувають значень в діапазоні $[\beta^e,\beta\times\beta^e)$. Отже, діапазон відносної помилки простягається від $\frac{(\beta/2)\times\beta^{-p}\times\beta^e}{\beta\times\beta^e}$ до $\frac{(\beta/2)\times\beta^{-p}\times\beta^e}{\beta^e}$. Скорочуючи, отримуємо такий діапазон для відносної помилки, що відповідає половині ulp \begin{equation} (\frac{1}{2}\times\beta^{-p}, \frac{\beta}{2}\times\beta^{-p}].\label{eq:re_ulp_diap} \end{equation} Таким чином відносна помилка, що відповідає одному ulp може відрізнятись на множник $\beta$. Встановлюючи $\varepsilon=(\beta/2)\times\beta^{-p}$ у значення більшої межі діапазону \eqref{eq:re_ulp_diap}, ми можемо сказати, що коли дійсне число заокруглене до найближчого числа з рухомою точкою, то відносна помилка завжди обмежена $\varepsilon$, яке ми називаємо машинним епсілон.

З того, що $\varepsilon$ може переоцінити ефект заокруглювання до найближчого числа з рухомою комою у $\beta$ разів, оцінки помилок у формулах будуть тугішими на машинах із маленьким $\beta.$

Відносна помилка при відніманні

У цьому прикладі ми використаємо той самий формат рухомої коми з основою ($\beta$) 10 і мантисою ($p$) 4. Припустимо нам потрібно порахувати різницю між $1$ і $0,9999$, в заявленому форматі вони виглядатимуть так $1.000\times 10^0$ та $9.999\times 10^{-1}$. Але для виконання арифметичних дій нам потрібно, щоб порядок у чисел був однаковим, тому друге число втратить одну цифру \begin{equation*}\begin{array}{c} \phantom{-}1.000\times 10^0\\ \underline{-0.999\times 10^0}\\ \phantom{-}0.001\times 10^0\\ \end{array}\end{equation*} Тоді як відповідь у випадку з нескінченною точністю становить $0.0001$, отже відносна помилка $\eta$, яка дорівнює різниці між двома значеннями поділеній на справжнє значення, така $$ \frac{|0.0001-0.001|}{0.0001} = \frac{0.0009}{0.0001} = 9 $$ Отже, відносна помилка може сягати $9 = \beta - 1.$

Вартові числа

Для подолання такої великої відностної помилки використовують вартові цифри (guard digits).
Теорема
Якщо $x$ і $y$ - це числа з рухомою комою у форматі з параметрами $\beta$ і $p$, і віднімання виконується із використанням одного вартового числа, тобто з $p+1$ цифрами, тоді помилка заокруглення результату менша ніж $2\varepsilon$.

Відкидання

Тут ми можемо підсумувати, що без використання вартових чисел, відносна помилка утворена через віднімання двох близьких величин може бути дуже великою. Інакше кажучи, обчислення будь-якого виразу, що містить віднімання (або додавання величин з різним знаком) може у висліді дати відносну помилку настільки вилику, що всі цифри беззмістовні. Ми розрізняємо два типи відкидання: катастрофічне і доброякісне.

Катастрофічне відкидання трапляється коли операнди були заокруглені. Наприклад, воно відбуваєтсья у виразі $b^2-4ac$. Покладемо $b = 3.345; a = 1.219; c = 2.294$. Точне значення $b^2-4ac$ рівне $0.003481$. Але $b^2 = 11.189025$ і $4ac = 11.185544$, а з урахуванням заокруглення, $b^2 = 11.19$ і $4ac = 11.19$, і їх різниця це $0$. Тобто, ми отримали помилку в розмірі 35 ulp, див. \eqref{eq:ulp}.

Доброякісне відкидання відбувається коли віднімаюємо дві точно відомі величини. Про це нам каже теорема з попереднього розділу.

Іншим прикладом формули де відбувається катастрофічне відкидання є $x^2 - y^2$, але ми можемо уникнути цього ефекту замінивши її на $(x-y)(x+y)$. Обчислення $(x-y)(x+y)$ і $x^2-y^2$ потребує приблизно однакової кількості дій, тому очевидно якій формі варто віддавати перевагу. Однак, на загал, заміна катастрофічного на доброякісне відкидання не варта зусиль якщо витрати на це занадто великі, бо входові дані часто (хоча й не завжди) є наближеними. Але повне виключення видкидання має сенс навіть якщо дані не точні.

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

\begin{equation} A = \sqrt{s(s-a)(s-b)(s-c)}, \mbox{де} s=(a+b+c)/2\label{eq:tr_area_catas} \end{equation} Припустимо, що трикутник дуже тонкий; тобто, $a\approx b+c$. Тоді $s\approx a$, а множник $s-a$ віднімає два майже однакових числа, одне з яких містить помилку заокруглення. Але ми можемо переписати формулу \eqref{eq:tr_area_catas} так, що вона повертає акуратні результати навіть для тонких трикутників: \begin{equation} A = \frac{\sqrt{(a+(b+c))(c-(a-b))(c+(a-b))(a+(b-c))}}{4}, a\ge b\ge c \label{eq:tr_area_benign} \end{equation}

неділя, 3 грудня 2017 р.

Розподіл порядкових статистик

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

В цьому дописі я розгляну лише випадкові вибірки отримані з неперервних розподілів (тобто функція розподілу ймовірностей неперервна). Нехай $X_1, \dots, X_n$ - це випадкова вибірка розміру $n$ з неперервного розподілу з функцією розподілу ймовірностей $F(x)$. Впорядкувавши цю вибірку у зростному порядку отримуємо $Y_1, \dots, Y_n$.

Порядкова статистика $Y_i$ називається $i$-тою порядковою статистикою. Тут ми припускаємо, що $Y_1 < Y_2 < \dots < Y_n$, тобто збігів бути не може.

Функція розподілу ймовірностей

Якщо має місце подія $Y_i \le y$ тоді ми знаємо, що щонайменше $i$ з $X_j$ лежать ліворуч від $y$. Розглянемо подію, що $Y_j \le y$ як успіх, тоді має місце $n$ проб Бернуллі з імовірністю успіху $F(X \le y)$. Нас цікавить імовірність мати не менше ніж $i$ успіхів: \begin{equation} F_{Y_i}(y) = P(Y_i \le y) = \sum_{k=i}^n\binom{n}{k}F(y)^k(1-F(y))^{n-k} \label{eq:cdf} \end{equation}

Густина ймовірності

Для знаходження густини ймовірності нам знадобиться така формула: \begin{equation} F_{Y_i}(y) = F_{Y_{i-1}}(y) - \binom{n}{i-1}F(y)^{i-1}(1-F(y))^{n-i+1} \label{eq:cdf_diff} \end{equation} Я одразу наведу формулу густину ймовірності, а потім ми доведемо її індукційно. \begin{equation} f_{Y_i}(y) = \frac{n!}{(i-1)!(n-i)!}F(y)^{i-1}(1-F(y))^{n-i}f_X(y)\label{eq:pdf} \end{equation} Почнемо з бази індукції. Розглянемо $i = 1$. Отже, нас цікавить іморність того, що хоча б одна проба менша ніж $y$. Доповняльною подією буде те, що усі проби більші ніж $y$. Тому, $F_{Y_1}(y) = (1 - F(y))^n$. Для отримання густини ймовірності нам потрібно обчислити похідну: $$f_{Y_1}(y) = n(1 - F(y))^{n-1}f_X(y)$$ Припустимо, що ми довели \ref{eq:pdf} для $Y_{i-1}$: $$f_{Y_{i-1}}(y) = \frac{n!}{(i-2)!(n-i+1)!}F(y)^{i-2}(1-F(y))^{n-i+1}f_X(y)$$ Тепер розглянемо випадок для $Y_i$. Тут нам можна скористатись рівнянням \ref{eq:cdf_diff}: \begin{align*} f_{Y_i}(y) &= f_{Y_{i-1}}(y)\\ &\phantom{=}-\binom{n}{i-1}(i-1)F(y)^{i-2}f_X(y)(1-F(y))^{n-i+1}\\ &\phantom{=}-\binom{n}{i-1}F(y)^{i-1}(n-i+1)(1-F(y))^{n-i}f_X(y) \end{align*} Трошки алгебри і правий бік рівняння перетвориться в \ref{eq:pdf}.

Коментар

Розглянемо другу половину \ref{eq:pdf}: $$F(y)^{i-1}f_X(y)(1-F(y))^{n-i}$$ тут перший множник - це імовірність, що $i-1$ проб менші ніж $y$, другий - імовірність, що одна проба поблизу $y$ і третій - імовірність, що $n-i$ проб більші ніж $y$. Таким чином \ref{eq:pdf} - це такий мультиноміальний розподіл: \begin{equation} f_{Y_i}(y) = \frac{n!}{(i-1)!1!(n-i)!}F(y)^{i-1}f_X(y)(1-F(y))^{n-i}\label{eq:pdf_multi} \end{equation}

четвер, 2 листопада 2017 р.

Двоїстість у лінійному програмуванні

Згадаємо формулювання прямої і двоїстої програм у стандартній формі.
\begin{align} \text{максимізувати}\ &\vec{c}^T\vec{x}\\ \text{за умов}\ &A\vec{x}\le \vec{b}\\ &\vec{x}\ge 0 \end{align} \begin{align} \text{мінімізувати}\ &\vec{b}^T\vec{y}\\ \text{за умов}\ &A^T\vec{y} \ge \vec{c} \label{eq:dual_constr}\\ &\vec{y}\ge 0 \end{align}

Фізичне тлумачення

  • Застосувати поле $\vec{c}$ до кульки в середині обмеженого багатогранника $A\vec{x} \le \vec{b}$.
  • Зрештою, кулька зупиниться біля краю багатогранника (в точці $\vec{x}'$).
  • Грань $\vec{a}_i\vec{x}'\le b_i$ накладає певну силу $-y_i\vec{a}_i$ на кульку.
  • Кулька стабільна, тому $\vec{c}^T = \sum_i y'_i\vec{a}_i = \vec{y}'^TA$.
  • Двоїстість можна уявити як намагання мінімізувати "роботу" $\sum_i y_ib_i$, щоб перенести м'ячик у початок координат рухаючи багатогранник.
  • В оптимальній позиції тиснуть лише стіни суміжні з кулькою. (На рисунку лише $y_1, y_2$ не нульові.)
Звернемо увагу на те, що на кульку тиснуть лише грані до яких вона притуляється. Тобто для кожного $i$, або вона перебуває на лінії $\vec{a}_i^T \vec{x} = b_i$ або тиск дорівнює нулю: \begin{align} \vec{a}_i\vec{x} = b_i, y_i > 0&\ \text{торкається грані, тиск є} \label{eq:compslack1}\\ \vec{a}_i\vec{x} < b_i, y_i = 0&\ \text{не торкається грані, тиску немає}\label{eq:compslack2} \end{align} Тепер покажемо, що наш $\vec{x}'$ і є оптимальним розв'язком. Згідно з принципом слабкої двоїстості $\vec{c}^T\vec{x} \le \vec{b}^T\vec{y}$. Отже, якщо ми покажемо, що $\vec{c}^T\vec{x}' = \vec{b}^T\vec{y}'$, то це і є оптимальний розв'язок. Cтабільність кульки разом із третім законом Ньютона даюсть нам, що $A^T \vec{y}' = \vec{c}$. Завдяки рівнянням (\ref{eq:compslack1}, \ref{eq:compslack2}) ми можемо записати, що $\vec{a}_i\vec{x}'y'_i = b_iy'_i$, а отже $$\vec{c}^T\vec{x}' = \vec{x}'^T\vec{c} = \vec{x}'^TA^T\vec{y}' = (A\vec{x}')^T\vec{y}' = \sum_i\vec{a}_i\vec{x}'y'_i = \sum_i b_iy'_i = \vec{b}^T\vec{y}'$$ Отже, $\vec{x}'$ - оптимальний розв'язок, а $\vec{y}'$ - двоїстий до нього.

Економічне тлумачення

Розглянемо задачу оптимального виробництва:
  • $n$ продуктів, $m$ сирих матеріалів.
  • Продукт $j$ потребує $a_{ij}$ одиниць сирого матеріалу $i$.
  • Всього доступно $b_i$ одиниць матеріалу $i$.
  • Продукт $j$ дає $c_j$ прибутку на одиницю.
  • Майстерня хоче максимізувати прибуток за наявних матеріалів
\begin{matrix} \text{максимізувати}&\vec{c}^T\vec{x}&\\ \text{за умов}&\vec{a}_i^T\vec{x}\le b_i,&\ i = 1, \dots, m.\\ &x_j\ge 0,& j = 1,\dots,n. \end{matrix} Якщо задача розв'язна і розв'язок скінченний, то допустима область - це опуклий багатогранник. Завжди існує щонайменше одна вершина в якій цільова функція набуває оптимального значення. Для задачі оптимального виробництва це означає, що для $n$ продуктів і $m$ сирих матеріалів, існує оптимальний план з щонайбільше $m$ продуктами.

Тлумачення двоїстих змінних

Для кожної одиниці продукту $j, j \in \{1,n\}$, нам потріно $a_{ij}, i \in \{1, m\}$. Для того, щоб було вигідно випускати продукт, його ціна має бути не меншою ніж сумарна ціна всіх складових. Позначимо ціну продукту через $c_p(j)$ і ціну матеріалу через $c_m(i)$, тоді $$c_p(j) \ge \sum_i a_{ij} c_m(i)$$ Але ж з точністю до знаку це те саме, що і обмеження двоїстої програми (\ref{eq:dual_constr}). Поглянемо на це з іншого боку. За умови що $\vec{x}$ і $\vec{y}$ оптимальні розв'язки відповідно прямої і двоїстої задач згідно з принципом сильної двоїстості маємо: $$\vec{c}^T\vec{x}=\vec{b}^T\vec{y}$$
  • Ліворуч маємо максимальний виторг.
  • Праворуч - $\sum_i (\text{наявність матеріалу}\ i) \times (\text{виторг на одиницю матеріалу}\ i)$
Інакше кажучі: Значення $y_i$ в оптимумі - це двоїста ціна ресурсу $i$.

Якщо ж значення не оптимальні, тобто

$$\vec{c}^T\vec{x}<\vec{b}^T\vec{y}$$ Розв'язок не оптимальний, бо ресурси не використані повністю.