diff --git a/docs/index.md b/docs/index.md index 652281d..06782f4 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,16 +1,3 @@ # Обо мне -Проверка работы *GitHub Actions*. - -## Commands - -* `mkdocs new [dir-name]` - Create a new project. -* `mkdocs serve` - Start the live-reloading docs server. -* `mkdocs build` - Build the documentation site. -* `mkdocs -h` - Print help message and exit. - -## Testing - -$$\Huge -\lim_{N \rightarrow \infty}\sum_{i=1}^N \frac{1}{i} = \infty -$$ +**TBD** diff --git a/docs/Математика/assets/baes/mode.jpg b/docs/Математика/assets/baes/mode.jpg new file mode 100644 index 0000000..cdab31b Binary files /dev/null and b/docs/Математика/assets/baes/mode.jpg differ diff --git a/docs/Математика/assets/baes/video.gif b/docs/Математика/assets/baes/video.gif new file mode 100644 index 0000000..176fa90 Binary files /dev/null and b/docs/Математика/assets/baes/video.gif differ diff --git a/docs/Математика/baes.md b/docs/Математика/baes.md new file mode 100644 index 0000000..8e6e0be --- /dev/null +++ b/docs/Математика/baes.md @@ -0,0 +1,221 @@ +# Немного про Байесовскую статистику +## Задача + +Представим, что у нас есть монета, честность которой нам неизвестна (если быть +точным, мы не знаем с каким шансом выпадает орёл, с каким — решка). Поэтому мы +бы хотели оценить вероятность $p$ — шанс выпадения орла [^1]. + +В [классической статистике](https://ru.wikipedia.org/wiki/%D0%A1%D1%82%D0%B0%D1%82%D0%B8%D1%81%D1%82%D0%B8%D0%BA%D0%B0) +эта задача бы решалась бы так: + +- Монета подбрасывается $n$ раз. +- Из них $m$ — количество выпавших орлов. +- Отношение $\frac{m}{n}$ будет оценкой $p$. + +В [Байесовской +статистике](https://ru.wikipedia.org/wiki/%D0%91%D0%B0%D0%B9%D0%B5%D1%81%D0%BE%D0%B2%D1%81%D0%BA%D0%B0%D1%8F_%D1%81%D1%82%D0%B0%D1%82%D0%B8%D1%81%D1%82%D0%B8%D0%BA%D0%B0) +подход иной: + +1. Обозначим монету как [бернуллевскую случайную + величину](https://ru.wikipedia.org/wiki/%D0%A0%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5_%D0%91%D0%B5%D1%80%D0%BD%D1%83%D0%BB%D0%BB%D0%B8) + $\xi$ с параметром $\theta$, у которой $1$ — это выпадение орла, $0$ — + решки. +2. Предполагается априорное распределение $\pi(\theta)$ (т.е. распределение, + которое мы предполагаем, исходя из того, что нам известно о параметре + $\theta$), как правило, это [равномерное + распределение](https://ru.wikipedia.org/wiki/%D0%9D%D0%B5%D0%BF%D1%80%D0%B5%D1%80%D1%8B%D0%B2%D0%BD%D0%BE%D0%B5_%D1%80%D0%B0%D0%B2%D0%BD%D0%BE%D0%BC%D0%B5%D1%80%D0%BD%D0%BE%D0%B5_%D1%80%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5) + $U \left(0,1\right)$. +3. Монета подбрасывается. +4. Распределение $\theta$ уточняется по формуле[^2]: + $$\Large + \pi(\theta | \xi) = \frac{p(\xi|\theta)p(\theta)} + {\int\limits_{\Theta}p(\xi|\theta)p(\theta)d\theta} + $$ +5. Повторить пункты 2-4, предполагая $\pi(\theta) = \pi(\theta|\xi)$ + + +Тогда оценкой $\theta$ будет: +$$\Large +\hat \theta = \arg \max_\theta \pi ( \theta | \xi ) +$$ + +Иначе говоря, +[мода](https://ru.wikipedia.org/wiki/%D0%9C%D0%BE%D0%B4%D0%B0_%28%D1%81%D1%82%D0%B0%D1%82%D0%B8%D1%81%D1%82%D0%B8%D0%BA%D0%B0%29) +апостериорного распределения. + +Мы, в качестве априорного распределения, взяли [Бета-распределения](https://ru.wikipedia.org/wiki/%D0%91%D0%B5%D1%82%D0%B0-%D1%80%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5) $Beta(2,2)$. В +таком случае, если на $i$-ом шаге мы имеем: +$$\Large +\theta \sim Beta \left( \alpha_i, \beta_i \right) +$$ + +то апостериорное распределение будет: + +$$\Large +Beta \left( \alpha_i, \beta_i \right) = +\begin{cases} +Beta \left( \alpha_i + 1, \beta_i \right), &\xi_i = 1 \\ +Beta \left( \alpha_i, \beta_i + 1 \right), &\xi_i = 0 \\ +\end{cases} +$$ + +Тогда не нужно интегрировать на каждом шаге, что существенно упрощает вычисления. + +## Решение + +Приведем решение на `python`. + +### Библиотеки + +Сначала нужно импортировать и настроить библиотеки: + +```python +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import seaborn as sns + +from tqdm import tqdm + +plt.rcParams.update({'font.size': 14}) +``` + +### Константы и функции + +Обозначим константы: + +- `p` — истинная вероятность выпадения орла. +- `NMODEL` — количество бросков монеты. + +```python +p = 0.523 +NMODEL = 7501 +``` + +Введем функции: +- `beta` — объект для Бета-распределения +- `coin` — функция одиночного броска монеты + +```python +from scipy.stats import beta +coin = lambda: np.random.binomial(n=1, p=p) +``` + +### Моделирование + +В качестве априорного распределения тут используется Бета-распределение +$Beta(2,2)$. Это сделано потому, что расчет моды накладывает ограничение, что оба +параметра должны быть строго больше 1. + +Мода расчитывается по формуле: + +$$\Large +\text{Mode} = \frac{\alpha - 1}{\alpha + \beta - 2} +$$ + +```python +alphas = np.zeros(NMODEL) +betas = np.zeros(NMODEL) +alphas[0], betas[0] = 2,2 +for i in tqdm(range(1, NMODEL)): + A = coin() + if A == 0: + betas[i] = 1 + else: + alphas[i] = 1 + +betas = np.cumsum(betas) +alphas = np.cumsum(alphas) + +Es = (alphas - 1) / (alphas + betas - 2) +``` + +### График + +#### Мода + +Построим график зависимости моды от количества бросков. + +```python +fig, ax = plt.subplots(1, 1, figsize=(10,8), dpi=80) +sns.lineplot(x=np.arange(NMODEL), + y=Es, + ax=ax, + label=r"arg$\max_p \; f(p|\mathbb{X})$", ) +sns.lineplot(x=np.arange(NMODEL), + y=p, + ax=ax, + label=r"Истинное $p = {:.3f}$".format(p)) +ax.grid() +ax.set(ylabel=r"$E(p|\mathbb{X})$", xlabel=r"$i$"); +``` + +![mode.jpg](assets/baes/mode.jpg) + +Гифка + +Сделаем гифку на тему. Сначала нужно импортировать библиотеку `imageio` и зададим +шаг `STEP` (так как, потом нужно будет сохранять каждый график в опертивной +памяти, что проблематично при большом `NMODEL`): + +```python +import imageio +STEP = NMODEL//150 +``` + +После, создадим список функций плотности вероятности: + +```python +pi = [] +for i in tqdm(range(0, NMODEL)): + pi.append(beta(alphas[i], betas[i])) +``` + +И по ним сделаем gif-изображение: + +```python +def create_frame(i): + fig, ax = plt.subplots(1, 1, figsize=(10,8), dpi=80) + x = np.linspace(0,1,500) + y = pi[i].pdf(x) + i_max = np.argmax(y) + x_max = x[i_max] + y_max = y[i_max] + ax.grid() + ax.plot(x,y, label=r"$f(p|\mathbb{X})$") + ax.axvline(x = p, color = 'red', label = f"Истинное $p={p}$") + plt.scatter([x_max], + [y_max], + color="green", + marker="o", + label=r"Мода $f(p|\mathbb{X})$") + plt.xlim([0,1]) + plt.xlabel('x', fontsize = 14) + plt.ylim([0,100]) + plt.ylabel(r'$f(p|\mathbb{X})$', fontsize = 14) + plt.title(r'Распределение p, $i={}$'.format(i+1), + fontsize=14) + ax.legend() + plt.savefig(f'./img/img_{i:04d}.png', + transparent = False, + facecolor = 'white' + ) + plt.close() + +for i in tqdm(range(0, NMODEL, STEP)): + create_frame(i) +frames = [] +for i in tqdm(range(0, NMODEL, STEP)): + image = imageio.v2.imread(f'./img/img_{i:04d}.png') + frames.append(image) +imageio.mimsave('./distributions.gif', + frames, + fps = 30) +``` + +Результат: + +![gif](assets/baes/video.gif) + +[^1]: Понятно, что вероятность выпадения решки равна $1 - p$ +[^2]: $p(\theta|\xi)$ -- это функция правдоподобия. diff --git a/mkdocs.yml b/mkdocs.yml index 0436b30..1a6ddac 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -1,7 +1,12 @@ site_name: Записки шизофреника site_url: https://rustbas.github.io/blog theme: - name: material + name: terminal + palette: dark + +plugins: + - search: + lang: ["ru", "en"] markdown_extensions: - pymdownx.arithmatex: @@ -9,6 +14,7 @@ markdown_extensions: - admonition - pymdownx.details - pymdownx.superfences + - footnotes extra_javascript: - javascripts/mathjax.js diff --git a/requirements.txt b/requirements.txt index 4c8f017..3fdd471 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1 +1,2 @@ mkdocs-material +mkdocs-terminal