Метод наименьших квадратов и распределение Гаусса

Применение МНК для подбора параметров гауссиан выборки

Два колокола

Черным — частоты данных в выборке + два «колокола», использованных для ее генерации (пунктиром)
Красным — оба «колокола», подобранные для выборки при помощи нелинейного программирования

Функция нормального распределения (Гаусса) f(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(x-\mu)^2}{2\sigma^2}) называется гауссианой и представляет из себя ни что иное как «кривую колокола», верхушка которого будет располагаться в точке математического ожидания (медианы), а отклонение (разброс) определяет его ширину. Ее же мы увидим, если построим гистограмму по нормально-распределенной выборке: близкие к медиане величины встречаются чаще, чем величины с бОльшим отклонением. Когда колокол в гистограмме только один, получить его параметры очень легко: математическое ожидание и стандартное отклонение всей выборки и являются параметрами колокола.

Однако, предположим, что у выборки есть несколько источников случайных величин, каждый со своей конфигурацией колокола. Под катом рассмотрим применение метода наименьших квадратов для определения параметров всех колоколов при помощи языка программирования Python и библиотек NumPy, SciPy, matplotlib.

Генерация входных данных

Для генерации данных зададим параметры для каждого из колоколов- математическое ожидание (mean), стандартное отклонение (std) и еще один параметр- prob— доля колокола в выборке (интеграл по колоколу = его площадь = вероятность того что величина сгенерируется именно в этом колоколе). Сумма таких «весов»= интегралу итоговой функции плотности вероятности = 1.

Итак, исходный вектор параметров (real_a), генерация данных (data) на основе него и вычисление гистограммы данных (hist_xs, hist_vals):

import numpy

N = 40000  # число элементов в будущей выборке

real_a = [1.0, 2.0, 0.5,  # mean, std, prob
          12.0, 4.0, 0.5]
bell_cnt = len(real_a) / 3

data = numpy.concatenate(
        [numpy.random.normal(real_a[3 * i],
                             real_a[3 * i + 1],
                             (N * real_a[3 * i + 2], 1))
         for i in xrange(bell_cnt)
         ])
hist_vals, hist_bin_ranges = numpy.histogram(data, 
                                       bins=1000, 
                                       normed=True)
hist_xs = (hist_bin_ranges[:-1] 
                           + hist_bin_ranges[1:]) * 0.5

Выборка довольно большая  (40000 элементов), как и гистограмма (1000 интервалов одинаковой ширины).

Функция плотности вероятности нормального распределения

В дальшейшем нам понадобится код функции (если кто-нибудь предложит в комментариях стандартную реализацию — я поправлю)

def norm_prob_func(x, mean, std, prob):
    return  prob \
            * numpy.exp(numpy.square(x - mean) \
                        * (-0.5 / (std ** 2))) \
            / (numpy.sqrt(2 * numpy.pi) * std)

Функция принимает x — скалярную величину либо вектор (numpy array)  середин интервалов (речь об интервалах гистограммы) и возвращает соответственно либо скалярную плотность, либо вектор плотностей. Звучит запутанно, но пример все разъяснит:

>>> norm_prob_func(0.2, 0.0, 1.0, 1.0)
0.391042693975
>>> norm_prob_func(0.6, 0.0, 1.0, 1.0)
0.333224602892
>>> norm_prob_func(0.8, 0.0, 1.0, 1.0)
0.289691552761
>>> norm_prob_func(numpy.array([0.2, 0.6, 0.8]), 0.0, 1.0, 1.0)
[ 0.39104269  0.3332246   0.28969155]

NumPy — очень хорошая библиотека,  мастхэв для анализа данных, де-факто — не знаю, как еще расхвалить. 🙂 Она содержит оптимизированные алгоритмы для операций над огромными массивами, что нам как раз и нужно.

Нам пригодится еще одна небольшая функция: sum_bells, принимает вектор параметров a, описывающий все колокола, по три параметра на каждый и вектор переменных x (середины интервалов, для которых хотим узнать плотность), а возвращает суммарную плотность вероятности для каждого интервала.

def sum_bells(a, x):
    return sum(norm_prob_func(x, *a[3 * i:3 * i + 3])
                           for i in xrange(bell_cnt))
>>> sum_bells(real_a, numpy.array([-2.0, 1.0, 6.0, 12.0, 20.0]))
[ 0.03248848  0.10087227  0.02057177  0.04986781  0.00674887]

График сгенерированных данных

from matplotlib import pyplot as plt

def draw_bell(mean, std, prob, style='b--'):
    X = numpy.linspace(mean - 3 * std, mean + 3 * std, 100)
    Y = norm_prob_func(X, mean, std, prob)
    plt.plot(X, Y, style, alpha=0.3)

plt.plot(hist_xs, hist_vals, 'k-', alpha=0.3)
for i in xrange(bell_cnt):
    draw_bell(*real_a[3 * i:3 * i + 3], style='k--')
plt.plot(hist_xs, sum_bells(real_a, hist_xs), 'k--', 
                            linewidth=2, alpha=0.5)
plt.show()
Сгенерированные данные

Ломаный серый график — это 1000-интервальная гистограмма по выборке.
Бледно-серый пунктир — отдельные колокола
Темный пунктир — итоговая функция плотности распределения

 Задача

Собственно, сама задача теперь- используя только сгенерированную выборку (data) и ее гистограмму (hist_xs, hist_vals) вычислить параметры исходных колоколов. С одним колоколом все просто: мат. ожидание = numpy.mean(data), стандартное отклонение = numpy.std(data). А вот для двух и более колоколов придется решать задачу оптимизации.

Нелинейное программирование — метод наименьших квадратов

Один из подходов оптимизации — подбор параметров в направлении минимизации суммы квадратов ошибок. Причем это «направление» выбирается в соответствии со значением производной в текущей точке.

Я использовал пакет SciPy — расширение numpy с разными стандартными алгоритмами, а именно scipy.optimize.leastsq (least squares).

Для того, чтобы уменьшать квадраты ошибки, мы должны предоставить начальный вектор параметров (init_a) и саму функцию ошибки (err_func). Следует отметить, что во многих случаях успех оптимизации зависит именно от начального вектора параметров — чем ближе он к «правильным» значениям, тем лучше. Я сталкивался с тем, что иногда процесс оптимизации вообще не может нормально завершиться из-за того что начальный вектор слишком далек от искомого. Я все еще в поиске «правильного» способа задать этот начальный вектор. Вам привожу как есть.

mn = numpy.mean(data)
dev = numpy.std(data)
init_a = [((mn - dev + i * 2) if i % 3 == 0 else (1.0 / bell_cnt)) 
          for i in xrange(bell_cnt * 3)]

def err_func(a, x, y):
    for i in xrange(bell_cnt):
        if a[i * 3 + 2] < 0:  # отрицательные веса исключаем
            return 1.0e10  # Очень Большая Ошибка
    return sum_bells(a, x) - y  # разность оценки и 
                                # реальных данных (из гистограммы)

a — это уже привычные нам параметры колоколов, x — это значения величин, на которых подбираются параметры. y— реальные данные из гистограммы, соответствующие x. (вектора x и y имеют одинаковую длину)

Ну и, собственно, сам подбор:

from scipy.optimize import leastsq

result, success = leastsq(err_func, 
                          init_a, 
                          args=(hist_xs, hist_vals), 
                          maxfev=100000)  # лимит итераций

result— это и есть вектор параметров, которые оптимизатор подобрал, которые, по идее, сводят ошибку к минимуму. То есть, мы получаем мат.ожидания, отклонения и веса отдельных колоколов. Теперь мы можем отобразить их на графике, поверх данных и посмотреть, насколько хорошо эти параметры описывают данные.

Результаты

for i in xrange(bell_cnt):
    draw_bell(*result[3 * i:3 * i + 3], style='r--')
plt.plot(hist_xs, sum_bells(result, hist_xs), 'r--', 
linewidth=2, alpha=0.5)

Два колокола

>>> result
[ 0.98332766, 2.00368246, 0.49786907,
 11.98476499, 4.0236398, 0.50173233]

Обратите внимание, как точно были подобраны параметры: даже не верится, что это было вычислено лишь на основе выборки случайных данных. Мне остается лишь предложить вам не верить мне на слово и проверить все самим 🙂

Напоследок, еще несколько «показательных» картинок (кликабельно):

Три колокола

Три колокола

Один большой и четыре маленьких колокола поверх

«Лапа»: один большой колокол с большой дисперсией и четыре маленьких поверх

Пять разных колоколов

Последнее изображение, кстати, примечательно тем, что на нем видно как ошибки на соседних колоколах нивелируют друг друга: один колокол выше, зато соседний ниже — сумма примерно та же.

Резюме

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

Однако есть и разочаровавшие меня моменты: процесс оптимизации в scipy.optimization.leastsq  стохастический: в зависимости от условий, далеко не всегда с первого раза выдает достойный результат. Кроме того, leastsq внутри использует «жадный» алгоритм Левенберга — Марквардт (метод градиента), который имеет свойство проваливаться в локальный минимум. Что это значит? Да то, что если на «пути» по «поверхности» суммарной ошибки от начального вектора параметров к искомому встретится «ямка» (такие параметры, отступление от которых в любую сторону повлечет за собой увеличение ошибки) то он алгоритм в ней и останется, приняв ее за конечное решение.  (Хотя можно предоставить свою функцию производной, в которой учитывать что-то вроде инерции- тогда, думаю, результаты можно улучшить, заставив его «перескакивать» через «ямки»)

В целом, работой оптимизатора я доволен, однако на production я бы не рискнул его использовать — слишком часто все зависит от магии и рандома. Думаю, должны быть более надежные реализации (пусть даже медленные — leastsq, кстати, стоит упомянуть, довольно быстрый)

Если есть вопросы — задавайте в комментариях, я с радостью постараюсь ответить.

Спасибо, что читаете!

Весь код целиком:

from scipy.optimize import leastsq
import itertools
from matplotlib import pyplot as plt
import numpy

N = 40000
real_a = [15.0, 20.0, 0.8,
          1.0, 1.0, 0.05,
          10.0, 1.0, 0.05,
          20.0, 1.0, 0.05,
          30.0, 1.0, 0.05]

bell_cnt = len(real_a) / 3

data = numpy.concatenate(
        [numpy.random.normal(real_a[3 * i],
                             real_a[3 * i + 1],
                             (N * real_a[3 * i + 2], 1))
         for i in xrange(bell_cnt)])
hist_vals, hist_x_ranges = numpy.histogram(data, bins=1000, normed=True)
hist_xs = (hist_x_ranges[:-1] + hist_x_ranges[1:]) * 0.5

def norm_prob_func(x, mean, std, prob):
    return  prob * \
            numpy.exp(numpy.square(x - mean) * (-0.5 / (std ** 2))) / \
            (numpy.sqrt(2 * numpy.pi) * std)

def sum_bells(a, x):
    return sum(norm_prob_func(x, *a[3 * i:3 * i + 3]) for i in xrange(bell_cnt))

def draw_bell(mean, std, prob, style='b--'):
    X = numpy.linspace(mean - 3 * std, mean + 3 * std, 100)
    Y = norm_prob_func(X, mean, std, prob)
    plt.plot(X, Y, style, alpha=0.3)

plt.plot(hist_xs, hist_vals, 'k-', alpha=0.3)
plt.plot(hist_xs, sum_bells(real_a, hist_xs), 'k--', linewidth=2, alpha=0.5)
[draw_bell(*real_a[3 * i:3 * i + 3], style='k--') for i in xrange(bell_cnt)]

mn = numpy.mean(data)
dev = numpy.std(data)

init_a = list(itertools.chain(*[[mn - dev + i * 8, 2.0, 0.1] for i in xrange(bell_cnt)]))
print init_a

def err(a, x, y):
    for i in xrange(bell_cnt):
        if a[i * 3 + 2] < 0:
            return 1e6
    return sum_bells(a, x) - y

result, success = leastsq(err, init_a, args=(hist_xs, hist_vals))

print success

for i in xrange(bell_cnt):
    draw_bell(*result[3 * i:3 * i + 3], style='r--')
plt.plot(hist_xs, sum_bells(result, hist_xs), 'r--',
         linewidth=2, alpha=0.5)
print result.reshape((len(result) / 3, 3))

plt.show()

 

 

Запись опубликована в рубрике Алгоритмы с метками , , , , . Добавьте в закладки постоянную ссылку.

2 комментария: Метод наименьших квадратов и распределение Гаусса

  1. Максим говорит:

    Что ты куришь? поделись!

  2. Alexandr говорит:

    Какие практические задачи могли бы быть решены таким путём?

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *