DataWondering Mar 28 at 15:17

Байесовская собака: анализ пёсьего компаса

Medium

23 min

5.5K

Open Data Science corporate blogPython*Data visualization*

Tutorial

✏️ Technotext 7

tl;dr

Ориентируются ли собаки по компасу, когда делают свои грязные дела? Оказывается — да! Если вам интересно, как можно это подтвердить в домашних условиях, используя компас, Байесовскую статистику и собаку (собака не включена), то добро пожаловать под кат.

Я очень люблю заниматься нестандартными исследованиями данных и периодически пишу о них статьи и заметки (вроде этой). Если вам такое тоже нравится — приглашаю подписаться на мой телеграм-канал: Data Wondering.

Это моя собака, Аури. Ему пять лет и по документам он породистый Кавалер-кинг-чарльз-спаниель, хотя по повадкам до кавалера бывает далеко.

Как и многие владельцы собак, во время наших с ним прогулок я начал замечать, что у Аури есть очень специфический туалетный ритуал. Каждый раз, когда он находит идеальное место, он начинает кружится вокруг своей оси, пока не выстроится поудобнее, прямо как стрелка компаса.

Поначалу, это просто выглядело забавно. В конце концов, кто знает, что там происходит в пёсьей голове. Однако спустя какое-то время я вспомнил, что в далеком 2013-м году мне попалась на глаза исследовательская статья под названием "Собаки чувствительны к небольшим вариациям магнитного поля Земли" (оригинал: Dogs are sensitive to small variations of the Earth’s magnetic field). Исследование было проведено на достаточно большой выборке собак, а главный вывод подтверждал, что "собаки предпочитают испражняться, ориентируя своё тело вдоль оси Север-Юг".

"Исследование как раз в моём стиле, надо повторить!", подумал я. И, к счастью, у меня как раз под рукой оказался идеальный подопытный. Я решил воспроизвести эту работу и подтвердить (либо опровергнуть!) гипотезу на Аури, моём ничего не подозревающем N=1 испытуемом.

С этого момента начался мой увлекательный и очень длительный процесс сбора данных, по результатам которого я заполучил более 150 "сеансов ориентирования", если вы понимаете, о чем я.

Сбор данных

Для моего исследования, мне нужно было считывать показания компаса каждый раз, когда Аури ходил в туалет. Благодаря последним достижениям современных технологий, у нас не только есть калькулятор на iPad, но и довольно точный компас на телефоне. Его-то я и стал использовать.

Методология была очень проста. Каждый раз, когда мой пёс усаживался поудобнее, я открывал приложение с компасом, выстраивал его вдоль собачьей тушки и делал скриншот. В оригинальной статье, авторы очень лаконично назвали этот процесс "определением направления компаса по грудному отделу позвоночника (между лопатками) в сторону головы". Очень научный способ сказать "стрелка компаса и голова собаки должны смотреть в одну сторону".

Прекрасная иллюстрация из оригинальной статьи.

В итоге, эту увлекательную разметку данных я повторил, в общей сложности, 150 с лишним раз на протяжении нескольких месяцев. Особенно увлекательно было ловить взгляды прохожих, которые явно не понимали, почему я тщательно фотографирую такой интимный момент собачьей жизни. Стоило ли оно того? Сейчас всё узнаем.

Анализ

Я кратко расскажу про извлечение и предобработку данных, после чего мы сразу перейдем к работе с круговыми распределениями и проверке гипотез.

Весь код и данные можно найти по ссылке на мой гитхаб: Data Wondering.

Как обработать скриншоты приложения?

После сбор данных у меня на руках оказалась серия скриншотов приложения компаса:

Будучи в меру ленивым, я не особо хотел вручную перебивать все данные со скриншотов. Вместо этого, я решил быстренько обработать все картинки при помощи простого скрипта.

Источник вдохновения при любой автоматизации

Задача была несложной: единственная информация, которая мне была нужна — это большие числа внизу экрана. К счастью, есть огромное количество предобученных моделей, которые легко справляются с базовым оптическим распознаванием символов (OCR). Я остановился на easyocr, библиотеке, которая, может, и не самая шустрая, но зато бесплатная и с ней легко работать.

Я вкратце приведу пример использования easyocr и еще одной библиотеки для работы с изображениями, opencv для извлечения показаний компаса из одного скриншота.

Для начала, давайте подгрузим изображение и посмотрим на него:

import cv2
import os

image_dir = '../data/raw/'
img = cv2.imread(image_dir + 'IMG_6828.png')

plt.imshow(img)
plt.show()

Далее я конвертирую изображение в черно-белое, чтобы убрать все ненужные цвета и сделать скриншот менее шумным:

gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
plt.imshow(gray, cmap='gray')
plt.show()

После этого я обрезаю все лишнее, чтобы сфокусироваться только на нужной мне области изображения:

Наконец, используя easyocr я извлекаю число с картинки. Помимо самого числа, модель возвращает еще и оценку вероятности (confidence score), но для моей задачи это было не нужно.

import easyocr

reader = easyocr.Reader(['en'])
result = reader.readtext(gray[1850:2100, 200:580])
for bbox, text, prob in result:
    print(f"Detected text: {text} with confidence {prob}")

>> Detected text: 340 with confidence 0.999995182215476

Вот и всё! Я написал простой цикл, чтобы пройтись по всем скриншотам, и сохранил результаты в CSV файл.

Здесь можно найти полный ноутбук с препроцессингом: Data Preprocessing.

Крутим колесо: круговые распределения

Обычно, я не работаю с круговыми распределениями, так что пришлось немного почитать и поизучать. В отличие от привычных нам распределений, круговые обладают одной интересной особенностью: "концы" распределения соединены друг с другом.

Например, если мы возьмем распределение часов в дне, то мы увидим, что расстояние между 23:00 и 00:00 точно такое же, как между 00:00 и 01:00. Или в случае с компасом, расстояние между 359° и 0° будет таким же, как и между 0° и 1°.

Даже вычисление выборочного среднего будет отличаться. Стандартное арифметическое среднее между 360° и 0° даст нам 180°, несмотря на то, что и 360°, и 0° указывают в одном направлении.

Аналогичным образом, я получал практически противоположные результаты на своих данных, оценивая среднее арифметически и используя специальный метод. Я перевел градусы компаса в радианы, используя вспомогательную функцию из еще одной хорошей библиотеки: pingouin, и посчитал круговое среднее через circ_mean:

from pingouin import circ_mean

arithmetic_mean = data['radians'].mean()
circular_mean = circ_mean(data['radians'])

print(f"Arithmetic mean: {arithmetic_mean:.3f}; Circular mean: {circular_mean:.3f}")

>> Arithmetic mean: 0.082; Circular mean: 2.989

Дальше мне захотелось визуализировать распределение моих данных из компаса. Для этого я использовал распределение фон Мизеса, которое моделирует круговые распределения.

Распределение фон Мизеса — это круговой аналог нормального распределения. Оно определяется двумя параметрами: средней локацией $\mu$ и концентрацией $\kappa$ . Концентрация определяет разброс данных вокруг среднего и аналогична обратной дисперсии. Когда параметр $\kappa$ равен нулю, распределение становится равномерным, при увеличении концентрации, распределение сосредотачивается вокруг среднего.

Импортируем необходимые библиотеки и напишем несколько вспомогательных функций:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

from scipy.stats import vonmises
from pingouin import convert_angles
from typing import Tuple, List


def vonmises_kde(series: np.ndarray, kappa: float, n_bins: int = 100) -> Tuple[np.ndarray, np.ndarray]:
    """
    Estimate a von Mises kernel density estimate (KDE) over circular data using scipy.
    
    Parameters:
    series: np.ndarray 
        The input data in radians, expected to be a 1D array.
    kappa: float
        The concentration parameter for the von Mises distribution.
    n_bins: int
        The number of bins for the KDE estimate (default is 100).
    
    Returns:
    bins: np.ndarray
        The bin edges (x-values) used for the KDE.
    kde: np.ndarray
        The estimated density values (y-values) for each bin.
    """
    bins = np.linspace(-np.pi, np.pi, n_bins)
    kde = np.zeros(n_bins)
    
    for angle in series:
        kde += vonmises.pdf(bins, kappa, loc=angle)
    
    kde = kde / len(series)
    return bins, kde


def plot_circular_distribution(
    data: pd.DataFrame,
    plot_type: str = 'kde',
    bins: int = 30,
    figsize: tuple = (4, 4),
    **kwargs
) -> None:
    """
    Plot a compass rose with either KDE or histogram for circular data.
    
    Parameters:
    -----------
    data: pd.DataFrame
        DataFrame containing 'degrees'and 'radians' columns with circular data
    plot_type: str
        Type of plot to create: 'kde' or 'histogram'
    bins: int
        Number of bins for histogram or smoothing parameter for KDE
    figsize: tuple
        Figure size as (width, height)
    **kwargs: dict
        Additional styling arguments for histogram (color, edgecolor, etc.)
    """
    plt.figure(figsize=figsize)
    ax = plt.subplot(111, projection='polar')

    ax.set_theta_zero_location('N')
    ax.set_theta_direction(-1)

    # add cardinal directions
    directions = ['N', 'E', 'S', 'W']
    angles = [0, np.pi / 2, np.pi, 3 * np.pi / 2]
    for direction, angle in zip(directions, angles):
        ax.text(
            angle, 0.45, direction,
            horizontalalignment='center',
            verticalalignment='center',
            fontsize=12,
            weight='bold'
        )

    if plot_type.lower() == 'kde':
        x, kde = vonmises_kde(data['radians'].values, bins)
        ax.plot(x, kde, color=kwargs.get('color', 'red'), lw=2)
    
    elif plot_type.lower() == 'histogram':
        hist_kwargs = {
            'color': 'teal',
            'edgecolor': 'black',
            'alpha': 0.7
        }
        hist_kwargs.update(kwargs) 
        
        angles_rad = np.deg2rad(data['degrees'].values)
        counts, bin_edges = np.histogram(
            angles_rad, 
            bins=bins, 
            range=(0, 2*np.pi), 
            density=True
        )
        widths = np.diff(bin_edges)
        ax.bar(
            bin_edges[:-1],
            counts,
            width=widths,
            align='edge',
            **hist_kwargs
        )
    
    else:
        raise ValueError("plot_type must be either 'kde' or 'histogram'")
    
    ax.xaxis.grid(True, linestyle='--', alpha=0.5)
    ax.yaxis.grid(True, linestyle='--', alpha=0.5)
    ax.set_yticklabels([]) 
    plt.show()

Теперь подгрузим данные и построим графики:

data = pd.read_csv('../data/processed/compass_degrees.csv', index_col=0)
data['radians'] = convert_angles(data['degrees'], low=0, high=360)

plot_circular_distribution(data, plot_type='histogram', figsize=(6, 6))
plot_circular_distribution(data, plot_type='kde', figsize=(5, 5))

Уже из гистограммы легко заметить, что у Аури есть определенные предпочтения при выборе направления. Есть очень четкий пик, указывающий на Север, и явный сдвиг в сторону Юга. Неплохо!

На круговом распределении плотности (KDE) мы получаем более сглаженное представление данных. И всё ещё легко заметить, что это распределение очень далеко от равномерного круга.

Настало время подтвердить это статистически!

Статистически значимые собачьи дела

Так же как для круговых данных нужны особые распределения, для них требуются и специальные статистические тесты.

Я буду использовать несколько тестов из библиотеки pingouin, которую я уже упоминал. Первый из них — это Rayleigh test, который проверяет равномерность круговых данных. Нулевая гипотеза утверждает, что данные распределены равномерно по кругу, альтернативная — что это не так.

from pingouin import circ_rayleigh

z, pval = circ_rayleigh(data['radians'])
print(f"Z-statistics: {z:.3f}; p-value: {pval:.6f}")

>> Z-statistics: 3.893; p-value: 0.020128

Хорошие новости! Значение p-value оказалось меньше 0.05, и мы можем отвергнуть нулевую гипотезу. Стратегическое позиционирование Аури во время туалета не случайно!

Единственный минус — в предпосылках теста говорится, что распределение должно быть унимодальным (то есть в распределении должен быть только один пик), а данные должны быть выборкой из распределения фон Мизеса. Ну что ж, попробуем что-то другое. У Аури явно есть несколько любимых направлений.

Следующий на очереди — V-test. Он проверяет, является ли распределение неравномерным с заданным средним направлением. Из документации к тесту:

V-тест обладает большей мощностью, чем тест Rayleigh, и предпочтителен, если есть предположение об определенном среднем направлении.

Отлично! Попробуем.

Из распределений выше очевидно, что Аури особенно тяготеет к южному направлению. Я задам среднее направление как $\pi$ радиан (что соответствует югу) и прогоню тест.

from pingouin import circ_vtest

v, pval = circ_vtest(data['radians'], dir=np.pi)
print(f"V-statistics: {v:.3f}; p-value: {pval:.6f}")

>> V-statistics: 24.127; p-value: 0.002904

Вот мы и подобрались к чему-то интересному! Значение p-value практически равно нулю, так что мы снова вправе отвергнуть нулевую гипотезу. Аури — статистически значимый любитель справлять нужду мордой на юг!

Байесовская собака: математическая часть

А теперь — кое-что совершенно другое. Попробуем байесовский подход.

Для начала я решил посмотреть, как оценка среднего направления меняется с увеличением выборки. Идея проста: я начну с равномерного кругового априорного распределения и буду обновлять его с каждым новым наблюдением.

Придётся ввести немного математики, так что давайте разберёмся с основами. Если вы не фанат уравнений, можете сразу перейти к следующему разделу с классными визуализациями!

1. Распределение фон Мизеса
Функция плотности вероятности распределения фон Мизеса $p(\theta | \mu, \kappa)$ выглядит следующим образом:

$p(\theta | \mu, \kappa) = \frac{\exp(\kappa \cos(\theta - \mu))}{2 \pi I_0(\kappa)}$

где:

$\mu$ — среднее направление
$\kappa$ — параметр концентрации (аналогичен обратному от дисперсии в нормальном распределении)
$I_0(\kappa)$ — модифицированная функция Бесселя первого рода, отвечающая за нормализацию распределения

2. Априорное распределение и функция правдоподобия

Предположим, что у нас есть:

Априорное распределение $p(\theta | \mu_{\text{prior}}, \kappa_{\text{prior}})$ с параметрами $(\mu_{\text{prior}}, \kappa_{\text{prior}})$
Функция правдоподобия $p(\theta | \theta_n, \kappa_{\text{likelihood}})$ для нового наблюдения $\theta_n$ с параметрами $(\theta_n, \kappa_{\text{likelihood}})$

Мы хотим обновить наше априорное распределение с помощью теоремы Байеса:

$p(\theta | \text{data}) \propto p(\theta | \mu_{\text{prior}}, \kappa_{\text{prior}}) \cdot p(\text{data} | \theta)$

где $p(\text{data} | \theta) = p(\theta | \theta_n, \kappa_{\text{likelihood}})$

3. Перемножение априорного распределения и функции правдоподобия для распределения фон Мизеса

Произведение двух распределений фон Мизеса с параметрами $(\mu_1, \kappa_1)$ и $(\mu_2, \kappa_2)$ также даёт распределение фон Мизеса, но с обновлёнными параметрами. Давайте разберёмся, как именно это работает:

Имея:

$p(\theta | \mu_{\text{prior}}, \kappa_{\text{prior}}) = \frac{\exp(\kappa_{\text{prior}} \cos(\theta - \mu_{\text{prior}}))}{2 \pi I_0(\kappa_{\text{prior}})}$

$p(\theta | \theta_n, \kappa_{\text{likelihood}}) = \frac{\exp(\kappa_{\text{likelihood}} \cos(\theta - \theta_n))}{2 \pi I_0(\kappa_{\text{likelihood}})}$

апостериорное распределение будет пропорционально произведению:

$p(\theta | \text{data}) \propto \exp(\kappa_{\text{prior}} \cos(\theta - \mu_{\text{prior}})) \cdot \exp(\kappa_{\text{likelihood}} \cos(\theta - \theta_n))$

Используя тригонометрическое преобразование для суммы косинусов:

$\begin {align*} \cos(\theta - \mu_{\text{prior}}) + \cos(\theta - \theta_n) = &\cos(\theta) ( \cos(\mu_{\text{prior}}) + \cos(\theta_n)) \\ & + \sin(\theta) (\sin(\mu_{\text{prior}}) + \sin(\theta_n)) \end {align*}$

Итого получаем:

$\begin{align*} p(\theta | \text{data}) \propto & \exp\Big(\left(\kappa_{\text{prior}} \cos(\mu_{\text{prior}}) + \kappa_{\text{likelihood}} \cos(\theta_n)\right) \cos(\theta) \\ & + \left(\kappa_{\text{prior}} \sin(\mu_{\text{prior}}) + \kappa_{\text{likelihood}} \sin(\theta_n)\right) \sin(\theta)\Big) \end{align*}$

4. Преобразование в полярные координаты для апостериорного распределения

Финальный рывок! Оказывается, выражение выше — это распределение фон Мизеса, замаскированное в другой форме. Используя полярные координаты, мы можем переписать его, чтобы получить оценки обновленного среднего направления и параметра концентрации.

Пусть:

$C = \kappa_{\text{prior}} \cos(\mu_{\text{prior}}) + \kappa_{\text{likelihood}} \cos(\theta_n)$ $S = \kappa_{\text{prior}} \sin(\mu_{\text{prior}}) + \kappa_{\text{likelihood}} \sin(\theta_n)$

Теперь выражение для апостериорного распределения упрощается до:

$p(\theta | \text{data}) \propto \exp(C \cos(\theta) + S \sin(\theta))$

Давайте остановимся и внимательнее посмотрим на это упрощённое выражение.

Обратите внимание, что $C \cos(\theta) + S \sin(\theta)$ — это скалярное произведение двух векторов (C, S) и $(\cos(\theta), \sin(\theta))$ , которое можно записать как:

$C \cos(\theta) + S \sin(\theta) = |(C, S)| \cdot |(\cos(\theta), \sin(\theta))| \cdot \cos(\phi),$

где $\phi$ — угол между векторами (C, S) и $(\cos(\theta), \sin(\theta))$

Длины этих векторов:

$|(C, S)| = \sqrt{C^2 + S^2}$ $|(\cos(\theta), \sin(\theta))| = 1$

Угол между вектором $(\cos(\theta), \sin(\theta))$ и положительной осью— это просто $\theta$ , а угол между (C, S) и положительной осьюпо определению:

$\phi = \arctan2(S, C)$

Следовательно, угол между этими двумя векторами:

$\cos(\phi) = \cos(\theta - \arctan2(S, C))$

Теперь подставим это обратно в выражение для апостериорного распределения:

$p(\theta | \text{data}) \propto \exp\left(\sqrt{C^2 + S^2} \cos(\theta - \arctan2(S, C))\right)$

или,

$p(\theta | \text{data}) \propto \exp\left(\kappa_{\text{post}} \cos(\theta - \mu_{\text{post}})\right)$

где:

$\kappa_{\text{post}}$ — апостериорный параметр концентрации:

$\kappa_{\text{post}} = \sqrt{C^2 + S^2}$

$\mu_{\text{post}}$ — апостериорное среднее направление:

$\mu_{\text{post}} = \arctan2(S, C)$

Ура, мы сделали это! Апостериорное распределение также является распределением фон Мизеса с обновлёнными параметрами $(\mu_{\text{post}}, \kappa_{\text{post}})$ . Теперь мы можем обновлять априорное распределение с каждым новым наблюдением и следить, как изменяется среднее направление.

Байесовская собака: весёлая часть

Добро пожаловать всем тем, кто скипнул математическую часть и моё почтение тем, кто пробрался сквозь неё. Нам осталось закодить итерации байесовского обновления параметров и визуализировать результаты.

Для начала, определим вспомогательную функцию для визуализации апостериорного распределения. Она нам понадобиться для создания красивых анимаций.

import imageio
from io import BytesIO

def get_posterior_distribution_image_array(
    mu_grid: np.ndarray, 
    posterior_pdf: np.ndarray, 
    current_samples: List[float], 
    idx: int, 
    fig_size: Tuple[int, int], 
    dpi: int, 
    r_max_posterior: float
) -> np.ndarray:
    """
    Creates the posterior distribution and observed samples histogram on a polar plot, 
    converts it to an image array, and returns it for GIF processing.

    Parameters:
    -----------

    mu_grid (np.ndarray): 
        Grid of mean direction values for plotting the posterior PDF.
    posterior_pdf (np.ndarray): 
        Posterior probability density function values for the given `mu_grid`.
    current_samples (List[float]): 
        List of observed angle samples in radians.
    idx (int): 
        The current step index, used for labeling the plot.
    fig_size (Tuple[int, int]): 
        Size of the plot figure (width, height).
    dpi (int): 
        Dots per inch (resolution) for the plot.
    r_max_posterior (float): 
        Maximum radius for the posterior PDF plot, used to set plot limits.

    Returns:
        np.ndarray: Image array of the plot in RGB format, suitable for GIF processing.
    """
    fig = plt.figure(figsize=fig_size, dpi=dpi)
    ax = plt.subplot(1, 1, 1, projection='polar')
    ax.set_theta_zero_location('N')  
    ax.set_theta_direction(-1)  
    ax.plot(mu_grid, posterior_pdf, color='red', linewidth=2, label='Posterior PDF')

    # observed samples histogram
    n_bins = 48
    hist_bins = np.linspace(-np.pi, np.pi, n_bins + 1)
    hist_counts, _ = np.histogram(current_samples, bins=hist_bins)

    # normalize the histogram counts
    if np.max(hist_counts) > 0:
        hist_counts_normalized = hist_counts / np.max(hist_counts)
    else:
        hist_counts_normalized = hist_counts

    bin_centers = (hist_bins[:-1] + hist_bins[1:]) / 2
    bin_width = hist_bins[1] - hist_bins[0]

    # set the maximum radius to accommodate both the posterior pdf and histogram bars
    r_histogram_height = r_max_posterior * 0.9 
    r_max = r_max_posterior + r_histogram_height
    ax.set_ylim(0, r_max)

    # plot the histogram bars outside the circle
    for i in range(len(hist_counts_normalized)):
        theta = bin_centers[i]
        width = bin_width
        hist_height = hist_counts_normalized[i] * r_histogram_height
        if hist_counts_normalized[i] > 0:
            ax.bar(
                theta, hist_height, width=width, bottom=r_max_posterior, 
                color='teal', edgecolor='black', alpha=0.5
            )

    ax.text(
        0.5, 1.1, f'Posterior Distribution (Step {idx + 1})', 
        transform=ax.transAxes, ha='center', va='bottom', fontsize=18
    )
    ax.set_yticklabels([])
    ax.grid(linestyle='--')
    ax.yaxis.set_visible(False)
    ax.spines['polar'].set_visible(False)
    plt.subplots_adjust(top=0.85, bottom=0.05, left=0.05, right=0.95)

    # saving to buffer for gif processing
    buf = BytesIO()
    plt.savefig(buf, format='png', bbox_inches=None, pad_inches=0)
    buf.seek(0)
    img_array = plt.imread(buf)
    img_array = (img_array * 255).astype(np.uint8)
    plt.close(fig)
    return img_array

Теперь мы готовы написать цикл для обновления параметров. Не забываем, что сначала нужно задать априорное распределение. Я начну с равномерного кругового распределения, что эквивалентно распределению фон Мизеса с параметром концентрации, равным 0.

Для $\kappa_{\text{likelihood}}$ я задам фиксированное умеренное значение концентрации . Это сделает шаг обновления апостериорного распределения более наглядным.

# initial prior parameters
mu_prior = 0.0  # initial mean direction (any value, since kappa_prior = 0)
kappa_prior = 0.0  # uniform prior over the circle

# fixed concentration parameter for the likelihood
kappa_likelihood = 2.0

posterior_mus = []
posterior_kappas = []

mu_grid = np.linspace(-np.pi, np.pi, 200)

# vizualisation parameters
fig_size = (10, 10)
dpi = 100

current_samples = []
frames = []

for idx, theta_n in enumerate(data['radians']):

    # compute posterior parameters
    C = kappa_prior * np.cos(mu_prior) + kappa_likelihood * np.cos(theta_n)
    S = kappa_prior * np.sin(mu_prior) + kappa_likelihood * np.sin(theta_n)
    kappa_post = np.sqrt(C**2 + S**2)
    mu_post = np.arctan2(S, C)

    # posterior distribution
    posterior_pdf = np.exp(kappa_post * np.cos(mu_grid - mu_post)) / (2 * np.pi * i0(kappa_post))

    # store posterior parameters and observed samples
    posterior_mus.append(mu_post)
    posterior_kappas.append(kappa_post)
    current_samples.append(theta_n)

    # plot posterior distribution
    r_max_posterior = max(posterior_pdf) * 1.1
    img_array = get_posterior_distribution_image_array(
        mu_grid, 
        posterior_pdf, 
        current_samples, 
        idx, 
        fig_size, 
        dpi, 
        r_max_posterior
        )
    frames.append(img_array)

    # updating priors for next iteration
    mu_prior = mu_post
    kappa_prior = kappa_post

# Create GIF
fps = 10
frames.extend([img_array]*fps*3) # repeat last frame a few times to make a "pause" at the end of the GIF
imageio.mimsave('../images/posterior_updates.gif', frames, fps=fps)

Вот и всё! Этот код сгенерирует гифку, которая показывает последовательное обновление апостериорного распределения с каждым новым наблюдением. Вот она во всём своём великолепии:

С каждым новым наблюдением, апостериорное распределение становится всё более и более сконцентрированным вокруг среднего направления. Еще бы заменить красную линию силуэтом Аури, и визуализация была бы идеальной.

Мы можем также посмотреть на историю обновлений апостериорных параметров среднего и концентрации:

# Convert posterior_mus to degrees
posterior_mus_deg = np.rad2deg(posterior_mus) % 360
n_samples = data.shape[0]
true_mu = data['degrees'].mean()
# Plot evolution of posterior mean direction
fig, ax1 = plt.subplots(figsize=(12, 6))

color = 'tab:blue'
ax1.set_xlabel('Number of Observations')
ax1.set_ylabel('Posterior Mean Direction (Degrees)', color=color)
ax1.plot(range(1, n_samples + 1), posterior_mus_deg, marker='o', color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax1.axhline(true_mu, color='red', linestyle='--', label='Sample Distribution Mean Direction')
ax1.legend(loc='upper left')
ax1.grid(True)

ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
color = 'tab:orange'
ax2.set_ylabel('Posterior Concentration Parameter (kappa)', color=color)  # we already handled the x-label with ax1
ax2.plot(range(1, n_samples + 1), posterior_kappas, marker='o', color=color)
ax2.tick_params(axis='y', labelcolor=color)

fig.tight_layout()  # otherwise the right y-label is slightly clipped
sns.despine()
plt.title('Evolution of Posterior Mean Direction and Concentration Over Time')
plt.show()

Апостериорное среднее направление и концентрация

График показывает, как меняются апостериорные параметры распределения с каждым новым наблюдением. Среднее направление постепенно сходится к выборочному среднему, а параметр концентрации увеличивается по мере того, как оценка становится более уверенной.

Коэффициент Байеса: PyMC для собачьего компаса

Последняя вещь, которую я хочу попробовать — это применить коэффициент Байеса (он же Байесовский фактор) для тестирования гипотез. Идея простая — нужно посчитать отношение маргинальных правдоподобий для двух конкурирующих гипотез или моделей, а затем интерпретировать полученное значение.

В общем случае Байесовский фактор определяется как:

$BF_{ij} = \frac{p(D|M_i)}{p(D|M_j)}=\frac{p(M_i|D)}{p(M_j|D)} \cdot \frac{p(M_j)}{p(M_i)}$

где:

и — это маргинальные функци правдоподобия данных при условии гипотез и
и — это апостериорные вероятностные распределения моделей при условии наблюдаемых данных
и — это априорные вероятностные распределения моделей

Результат вычислений — это число, которое говорит нам, насколько более вероятна одна гипотеза по сравнению с другой. Есть множество разных подходов к интерпретации коэффициента Байеса, один из часто встречающихся — это шкала Джеффриса, предложенная Гарольдом Джеффрисом.

Коэффициент Байеса	Интерпретация
< 1	Отрицательный эффект
От 1 до 3	Анекдотичный эффект
От 3 до 10	Достаточный эффект
От 30 до 100	Очень сильный эффект
> 100	Решающий эффект

Вы можете спросить, а что же за модели мы используем? Всё просто! Это вероятностные распределения с различными параметрами. Я буду использовать PyMC, чтобы задать эти модели и насэмплировать апостериорные распределения из них.

Для начала, давайте заново определим нашу нулевую гипотезу. Я вновь буду предполагать круговое равномерное распределение фон Мизеса с параметром концентрации $\kappa=0$ , но на этот раз нам нужно оценить функцию правдоподобия данных при этой гипотезе. Чтобы упростить вычисления, здесь и далее мы будем использовать логарифмом правдоподобия.

# Calculate log likelihood for H0
log_likelihood_h0 = vonmises.logpdf(data['radians'], kappa=0, loc=0).sum()

Далее, давайте введем альтернативную модель-гипотезу. Начнем с простого сценария: унимодальный Юг, где я предполагаю, что Аури в основном смотрит на юг, а значит, распределение сконцентрировано вокруг 180° на компасе или $\pi$ в радианах.

import pymc as pm
import arviz as az
import arviz.data.inference_data as InferenceData
from scipy.stats import halfnorm, gaussian_kde

with pm.Model() as model_uni:
    # Prior for kappa 
    kappa = pm.HalfNormal('kappa', sigma=10)
    # Likelihood
    likelihood_h1 = pm.VonMises('angles', mu=np.pi, kappa=kappa, observed=data['radians'])
    # Sample from posterior 
    trace_uni = pm.sample(
        10000, tune=3000, chains=4, 
        return_inferencedata=True, 
        idata_kwargs={'log_likelihood': True})

В результате получаем простую модель, которую мы можем визуализировать в виде графа:

# Model graph
pm.model_to_graphviz(model_uni)

А вот и апостериорное распределение параметра концентрации:

az.plot_posterior(trace_uni, var_names=['kappa'])
plt.show()

Всё, что осталось сделать — это посчитать логарифм правдоподобия для альтернативной модели и оценить коэффициент Байеса.

# Posterior samples for kappa
kappa_samples = trace_uni.posterior.kappa.values.flatten()
# Log likelihood for each sample
log_likes = []
for k in kappa_samples:
    # Von Mises log likelihood
    log_like = vonmises.logpdf(data['radians'], k, loc=np.pi).sum()
    log_likes.append(log_like)
# Log-mean-exp trick for numerical stability
log_likelihood_h1 = np.max(log_likes) +\
       np.log(np.mean(np.exp(log_likes - np.max(log_likes))))
BF = np.exp(log_likelihood_h1 - log_likelihood_h0)
print(f"Bayes Factor: {BF:.4f}")
print(f"Probability kappa > 0.5: {np.mean(kappa_samples > 0.5):.4f}")

>> Bayes Factor: 32.4645
>> Probability kappa > 0.5: 0.0649

Так как мы делим правдоподобие альтернативной модели на правдоподобие нулевой гипотезы, то коэффициент Байеса говорит нам о том, насколько более вероятно пронаблюдать данные при условии, что верна альтернативная гипотеза. В нашем случае, мы получили значение коэффициента , что соответствует очень сильному эффекту. А значит, данные не распределены равномерно по кругу и есть основания полагать, что южное направление предпочтительнее.

Однако, мы дополнительно посчитали вероятность того, что параметр концентрации больше . Это простой способ проверить, значимо ли распределение отличается от равномерного. У нашей модели унимодального Юга эта вероятность составляет всего , что означает, что распределение всё ещё очень сильно размазано по кругу и не особо сконцентрировано.

Попробуем другую модель: бимодальную смесь Север-Юг.

Бимодальная смесь Северг-Юг

На этот раз я предположу, что собака предпочитает смотреть либо на северг, либо на юг, когда занимается своими важными делами. В рамках модели это означает концентрацию вокруг 0° и 180° градусов компаса.

Чтобы построить такую модель, мне потребуется смесь двух распределений фон Мизеса с разными фиксированными значениями средних направлений и общим параметром концентрации.

Давайте зададим несколько вспомогательных функций:

# Type aliases
ArrayLike = Union[np.ndarray, pd.Series]
ResultDict = Dict[str, Union[float, InferenceData.InferenceData]]

def compute_mixture_vonmises_logpdf(
    series: ArrayLike,
    kappa: float,
    weights: npt.NDArray[np.float64],
    mus: List[float]
) -> float:
    """
    Compute log PDF for a mixture of von Mises distributions
    
    Parameters:
    -----------
    series: ArrayLike 
        Array of observed angles in radians
    kappa: float
        Concentration parameter
    weights: npt.NDArray[np.float64],
        Array of mixture weights
    mus: List[float] 
        Array of means for each component
    
    Returns:
    --------
    float: Sum of log probabilities for all data points
    """
    mixture_pdf = np.zeros_like(series)
    
    for w, mu in zip(weights, mus):
        mixture_pdf += w * vonmises.pdf(series, kappa, loc=mu)
    
    return np.log(np.maximum(mixture_pdf, 1e-300)).sum()

def compute_log_likelihoods(
    trace: az.InferenceData, 
    series: ArrayLike,
    mus: List[float]
    ) -> np.ndarray:
    """
    Compute log likelihoods for each sample in the trace

    Parameters:
    -----------
    trace: az.InferenceData
        The trace from the PyMC3 model sampling.

    series: ArrayLike
        Array of observed angles in radians
    
    """

    kappa_samples = trace.posterior.kappa.values.flatten()
    weights_samples = trace.posterior.weights.values.reshape(-1, 2)
    # Calculate log likelihood for each posterior sample
    log_likes = []
    for k, w in zip(kappa_samples, weights_samples):
        log_like = compute_mixture_vonmises_logpdf(
            series, 
            kappa=k, 
            weights=w, 
            mus=mus
        )
        log_likes.append(log_like)
    
    # Calculate marginal likelihood using log-sum-exp trick
    log_likelihood_h1 = np.max(log_likes) + np.log(np.mean(np.exp(log_likes - np.max(log_likes))))
    return log_likelihood_h1

def posterior_report(
    log_likelihood_h0: float, 
    log_likelihood_h1: float, 
    kappa_samples: ArrayLike,
    kappa_threshold: float = 0.5
    ) -> str:

    """
    Generate a report with Bayes Factor and probability kappa > threshold

    Parameters:
    -----------
    log_likelihood_h0: float
        Log likelihood for the null hypothesis
    log_likelihood_h1: float
        Log likelihood for the alternative hypothesis
    kappa_samples: ArrayLike
        Flattened posterior samples of the concentration parameter
    kappa_threshold: float
        Threshold for computing the probability that kappa > threshold

    Returns:
    --------
    summary: str
        A formatted string containing the summary statistics.
    """
    BF = np.exp(log_likelihood_h1 - log_likelihood_h0)

    summary = (
        f"Bayes Factor: {BF:.4f}\n"
        f"Probability kappa > {kappa_threshold}: {np.mean(kappa_samples > kappa_threshold):.4f}"
    )

    return summary

А теперь обратно к модели:

mu1 = 0            # 0 degrees
mu2 = np.pi        # 180 degrees

with pm.Model() as model_mixture_bimodal_NS:
    # Priors for concentration parameters
    kappa = pm.HalfNormal('kappa', sigma=10) 
    # Priors for component weights
    weights = pm.Dirichlet('weights', a=np.ones(2))
    
    # Define the von Mises components
    vm1 = pm.VonMises.dist(mu=mu1, kappa=kappa)
    vm2 = pm.VonMises.dist(mu=mu2, kappa=kappa)
    
    # Mixture distribution
    likelihood = pm.Mixture(
        'angles',
        w=weights,
        comp_dists=[vm1, vm2],
        observed=data['radians']
    )
    
    # Sample from the posterior
    trace_mixture_bimodal_NS = pm.sample(
        10000, tune=3000, chains=4, return_inferencedata=True, idata_kwargs={'log_likelihood': True})
    
    # Get kappa samples
    kappa_samples = trace_mixture_bimodal_NS.posterior.kappa.values.flatten()

И снова мы можем визуализировать граф модели и апостериорное распределение параметра концентрации $\kappa$ :

# Model graph
pm.model_to_graphviz(model_mixture_bimodal_NS)

# Posterior Analysis
az.plot_posterior(trace_mixture_bimodal_NS, var_names=['kappa'])
plt.show()

И, наконец, посчитаем коэффициент Байеса и вероятность того, что параметр концентрации больше 0.5:

log_likelihood_h1 = compute_log_likelihoods(trace_mixture_bimodal_NS, 
                                            data['radians'], [mu1, mu2])
print(posterior_report(log_likelihood_h0, log_likelihood_h1, kappa_samples))

>> Bayes Factor: 214.2333
>> Probability kappa > 0.5: 0.9110

Критическая удача! Обе метрики указывают на то, что эта модель гораздо лучше подходит под наблюдаемые данные. Коэффициент Байеса говорит о решающем эффекте и больше 90% всех апостериорных значений параметра концентрации больше 0.5 со средним в 0.99, как мы можем убедиться на распределении выше.

Попробуем еще пару моделей перед тем, как перейти к результатам.

Бимодальная смесь Запад-Юг

В этой модели мы вновь предполагаем, что собака предпочитает смотреть на две стороны света, но в этот раз на запад и на юг или 270° и 180° — тоже частые направления на самом первом графике с распределением по компасу.

mu1 = np.pi          # 180 degrees
mu2 = 3 * np.pi / 2  # 270 degrees

with pm.Model() as model_mixture_bimodal_WS:
    # Priors for concentration parameters
    kappa = pm.HalfNormal('kappa', sigma=10)
    
    # Priors for component weights
    weights = pm.Dirichlet('weights', a=np.ones(2))
    
    # Define the four von Mises components
    vm1 = pm.VonMises.dist(mu=mu1, kappa=kappa)
    vm2 = pm.VonMises.dist(mu=mu2, kappa=kappa)
    
    # Mixture distribution
    likelihood = pm.Mixture(
        'angles',
        w=weights,
        comp_dists=[vm1, vm2],
        observed=data['radians']
    )
    
    # Sample from the posterior
    trace_mixture_bimodal_WS = pm.sample(
        10000, tune=3000, chains=4, return_inferencedata=True, 
      idata_kwargs={'log_likelihood': True})
    
    # Get kappa samples
    kappa_samples = trace_mixture_bimodal_WS.posterior.kappa.values.flatten()

# Posterior Analysis
az.plot_posterior(trace_mixture_bimodal_WS, var_names=['kappa'])
plt.show()

log_likelihood_h1 = compute_log_likelihoods(trace_mixture_bimodal_WS, 
                                            data['radians'], [mu1, mu2])
print(posterior_report(log_likelihood_h0, log_likelihood_h1, kappa_samples))

>> Bayes Factor: 20.2361
>> Probability kappa > 0.5: 0.1329

Ну не очень. Точно не на уровне предыдущей модели. Поехали дальше!

Смесь "на все четыре стороны"

Последняя попытка. Что если собака ориентируется строго по сторонам света? Попробуем распределение с четырьмя средними направлениями: 0°, 90°, 180° и 270°.

mu1 = 0            # 0 degrees
mu2 = np.pi / 2    # 90 degrees
mu3 = np.pi        # 180 degrees
mu4 = 3 * np.pi / 2  # 270 degrees

with pm.Model() as model_mixture_quad:
    # Priors for concentration parameters
    kappa = pm.HalfNormal('kappa', sigma=10)
    
    # Priors for component weights
    weights = pm.Dirichlet('weights', a=np.ones(4))
    
    # Define the four von Mises components
    vm1 = pm.VonMises.dist(mu=mu1, kappa=kappa)
    vm2 = pm.VonMises.dist(mu=mu2, kappa=kappa)
    vm3 = pm.VonMises.dist(mu=mu3, kappa=kappa)
    vm4 = pm.VonMises.dist(mu=mu4, kappa=kappa)
    
    # Mixture distribution
    likelihood = pm.Mixture(
        'angles',
        w=weights,
        comp_dists=[vm1, vm2, vm3, vm4],
        observed=data['radians']
    )
    
    # Sample from the posterior
    trace_mixture_quad = pm.sample(
        10000, tune=3000, chains=4, return_inferencedata=True, idata_kwargs={'log_likelihood': True}
    )
    # Get kappa samples
    kappa_samples = trace_mixture_quad.posterior.kappa.values.flatten()
# Posterior Analysis
az.plot_posterior(trace_mixture_quad, var_names=['kappa'])
plt.show()
log_likelihood_h1 = compute_log_likelihoods(trace_mixture_quad, data['radians'], [mu1, mu2, mu3, mu4])
print(posterior_report(log_likelihood_h0, log_likelihood_h1, kappa_samples))

>> Bayes Factor: 0.0000
>> Probability kappa > 0.5: 0.9644

Что ж, снова не очень. Хотя вероятность того, что параметр концентрации больше 0.5, коэффициент Байеса стал равен 0.0. А это означает, что модель объясняет наблюдаемые данные намного хуже даже нулевой гипотезы с равномерным распределением.

Отличная особенность коэффициента Байеса — встроенная регуляризация, которая наказывает слишком сложные модели.

Сравнение моделей

Давайте посмотрим на результаты моделирования под другим углом — используя информационные критерии. Мы возьмем два наиболее распространенных критерия, которые используют в Байесовском моделировании: Widely Applicable Information Criterion (WAIC) и Leave-One-Out Cross-Validation (LOO).

# Compute WAIC for each model
wail_uni = az.waic(trace_uni)
waic_quad = az.waic(trace_mixture_quad)
waic_bimodal_NS = az.waic(trace_mixture_bimodal_NS)
waic_bimodal_WS = az.waic(trace_mixture_bimodal_WS)

model_dict = {
    'Quadrimodal Model': trace_mixture_quad,
    'Bimodal Model (NS)': trace_mixture_bimodal_NS,
    'Bimodal Model (WS)': trace_mixture_bimodal_WS,
    'Unimodal Model': trace_uni 
}
# Compare models using WAIC
waic_comparison = az.compare(model_dict, ic='waic')
waic_comparison

# Compare models using LOO
loo_comparison = az.compare(model_dict, ic='loo')
loo_comparison

# Visualize the comparison
az.plot_compare(waic_comparison)
plt.show()

И у нас есть победитель! Бимодальная модель Север-Юг лучше всего объясняет наблюдаемые данные согласно обоим информационным критериям.

Выводы

Получилось неплохо! То, что началось, как простое наблюдение за странными туалетными привычками моей собаки, превратилось в полноценное байесовское исследование.

В этой статье я показал, как работать с круговыми данными, оценивать среднее направление и параметр концентрации, а также использовать байесовское моделирование для получения апостериорных оценок параметров распределения. Мы также посмотрели на применение коэффициентов Байеса для тестирование гипотез и на сравнение моделей при помощи информационных критериев.

И результаты тоже не разочаровали! Аури действительно имеет свои предпочтения и каким-то образом часто умудряется выстроиться вдоль оси Север-Юг. Если я когда-нибудь заблужусь в лесу со своим пёселем, то вместо мха на дереве смогу узнать примерное направление по-другому. Нужно только собрать достаточно большой размер выборки.

Надеюсь, вам тоже понравилось это путешествие. Если у вас есть вопросы или пожелания, смело пишите мне!

Еще больше интересных исследований и заметок про анализ данных, машинное обучение, HealthTech и много другое можно найти в моём телеграм-канале Data Wondering:

Data Wondering

t.me

Tags:

Hubs:

If this publication inspired you and you want to support the author, do not hesitate to click on the button

Байесовская собака: анализ пёсьего компаса

tl;dr

Сбор данных

Анализ

Как обработать скриншоты приложения?

Крутим колесо: круговые распределения

Статистически значимые собачьи дела

Байесовская собака: математическая часть

3. Перемножение априорного распределения и функции правдоподобия для распределения фон Мизеса

Байесовская собака: весёлая часть

Коэффициент Байеса: PyMC для собачьего компаса

Бимодальная смесь Северг-Юг

Бимодальная смесь Запад-Юг

Смесь "на все четыре стороны"

Сравнение моделей

Выводы

Articles

Information

zproxy.org