Библиотека Scipy

Гауссова плотность распределения вероятностей и кумулятивная функция распределения

Здравствуйте и вновь добро пожаловать на занятия по теме «Инструментарий Numpy на языке Python».

В этой заключительной статье мы начнём рассмотрение библиотеки Scipy. Вновь-таки, сугубо для повторения: целью этого курса является не показать до конца все возможности каждой из библиотек, а отобрать те из них, которые чаще всего используются в моих курсах по машинному и глубокому обучению. Смысл в том, что необходимо помнить только то, что часто используется; если же что-то используется редко, вы всё равно это забудете, а попытка заучить абсолютно всё – лишь пустая трата времени, которое можно было бы потратить на изучение более полезных вещей. Не забывайте, что если (и когда) вам понадобится определённая функция, вы всегда можете заглянуть в документацию. Этот курс поможет вам познакомиться и освоиться со всеми этими библиотеками, так что когда в конце концов в будущем вам понадобится что-то в документации, это будет не страшно.

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

Унас уже есть всё необходимое, чтобы вычислить значение этой функции: мы знаем,как возводить в квадрат, как делить, как вычислять экспоненту и как извлекатьквадратный корень. Возникает вопрос: зачем тогда вообще это обсуждать? Конечноже, потому что Scipy предоставляетлучший способ вычислить плотность распределения вероятностей. На самом деле по возможностинужно всегда использовать Scipy, потому как этабиблиотека работает очень быстро.

Вначале импортируем модуль norm:

from scipy.stats import norm

Найдём,к примеру, плотность вероятности нуля для стандартного нормальногораспределения:

norm.pdf(0)

Разумеется,вы можете работать с гауссовым распределением, у которого среднее значениеотлично от нуля, а дисперсия не равна единице. В таком случае можно подставитьэти значения в функцию pdf:

norm.pdf(0, loc=5, scale=10)

Обратитевнимание, что параметр scale относится кстандартному отклонению, а не дисперсии, то есть стандартное отклонение, равное10, соответствует дисперсии, равной 100. Как и ожидалось, плотность вероятностипри этом стала значительно меньше.

Частобывает необходимо вычислить плотности распределения вероятностей сразу многихразных значений, содержащихся в массиве. Первое, что может прийти в голову, –использовать цикл for для вычислениякаждого из значений плотности распределения вероятностей по отдельности, но,конечно же, это не лучший вариант. Как и функции Numpy,рассматриваемая нами функция также позволяет проводить поэлементные вычисления.Так, имея случайный массив, мы можем одновременно вычислить все значенияплотности распределения вероятностей:

import numpy as np

r = np.random.randn(10)

norm.pdf(r)

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

Унас уже есть всё необходимое, чтобы вычислить значение этой функции: мы знаем,как возводить в квадрат, как делить, как вычислять экспоненту и как извлекатьквадратный корень. Возникает вопрос: зачем тогда вообще это обсуждать? Конечноже, потому что Scipy предоставляетлучший способ вычислить плотность распределения вероятностей. На самом деле по возможностинужно всегда использовать Scipy, потому как этабиблиотека работает очень быстро.

Вначале импортируем модуль norm:

from scipy.stats import norm

Найдём,к примеру, плотность вероятности нуля для стандартного нормальногораспределения:

norm.pdf(0)

Разумеется,вы можете работать с гауссовым распределением, у которого среднее значениеотлично от нуля, а дисперсия не равна единице. В таком случае можно подставитьэти значения в функцию pdf:

norm.pdf(0, loc=5, scale=10)

Обратитевнимание, что параметр scale относится кстандартному отклонению, а не дисперсии, то есть стандартное отклонение, равное10, соответствует дисперсии, равной 100. Как и ожидалось, плотность вероятностипри этом стала значительно меньше.

Частобывает необходимо вычислить плотности распределения вероятностей сразу многихразных значений, содержащихся в массиве. Первое, что может прийти в голову, –использовать цикл for для вычислениякаждого из значений плотности распределения вероятностей по отдельности, но,конечно же, это не лучший вариант. Как и функции Numpy,рассматриваемая нами функция также позволяет проводить поэлементные вычисления.Так, имея случайный массив, мы можем одновременно вычислить все значенияплотности распределения вероятностей:

import numpy as np

r = np.random.randn(10)

norm.pdf(r)

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

Каквы знаете, сложение – менее трудоёмкая операция, нежели умножение, поэтомутакая замена позволяет сэкономить время. Поскольку гауссово распределение имеетформу экспоненты, взятие её логарифма позволяет убрать экспоненциальную функцию,что также экономит время, ведь вычисление экспоненты также является трудоёмкойоперацией:

Поэтомув большинстве случаев целесообразнее работать с логарифмом распределениявероятностей, а не самой функцией распределения вероятностей. Очень удобно то,что и для этого есть соответствующая функция:

norm.logpdf(r)

Ипоследнее, о чём мы здесь поговорим, – это кумулятивная функция распределения.Напомню, что она равна интегралу функции плотности распределения вероятностейот минус бесконечности до x. Также напомню, что этот интеграл нельзя взятьнапрямую, то есть найти математическое выражение для него. Можно лишь вычислитьего численно, и, конечно же, Scipy предлагает намфункцию и для этого:

norm.cdf(r)

И,конечно же, есть функция для логарифма кумулятивной функции распределения:

norm.logcdf(r)

Выборка из одномерного гауссового распределения

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

r = np.random.randn(10000)

plt.hist(r, bins=100)

plt.show()

Вэтот раз мы использовали больше сегментов, поскольку у нас больше точек.

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

r = 10*np.random.randn(10000) + 5

Этодаёт нам выборку из гауссового распределения со стандартным отклонением 10 исредним значением 5. Обратите внимание, что данное выражение представляет собойлишь обратное преобразование из нестандартной нормальной выборки в стандартнуюнормальную. Вновь построим гистограмму, чтобы убедиться, что мы получилижелаемое:

plt.hist(r, bins=100)

plt.show()

Похоже, всё работает.

Выборка из сферического и эллиптического гауссового распределения с привязкой к осям координат

Вмашинном обучении часто приходится работать с многомерными данными.Предположим, нам нужна выборка из двухмерного гауссового распределения сосредним значением 0 и дисперсией 1. Это называется сферическим гауссовымраспределением, поскольку каждая размерность не коррелирует и не зависит отдругой, но все они имеют одну и ту же дисперсию. При отображении данных онираспределяются по кругу. Разумеется, в этом случае мы можем использовать ту жефункцию, что и ранее, лишь добавить больше размерностей:

r = np.random.randn(10000, 2)

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

plt.scatter(r[:,0], r[:,1])

plt.show()

Наэкране выглядеть это может не совсем похожим на круг, если оси координатпо-разному отмасштабированы, так что, возможно, придётся немного сузить окно.

Теперьпредположим, что нам нужно эллиптическое гауссово распределение. Это когдадисперсия своя для каждой из размерностей. Предположим, например, чтостандартное отклонение второй размерности равно 5, а среднее значение равно 2.Тогда можно просто отмасштабировать нужный столбец данных:

r[:,1] = 5*r[:,1] + 2

Ещёраз создадим диаграмму рассеяния, чтобы убедиться, что всё работает:

plt.scatter(r[:,0], r[:,1])

plt.show()

Заметьте,что выглядит всё это по-прежнему похожим на круг, но это потому, что Matplotlib масштабирует оси координат такимобразом, чтобы данные были сбалансированы по каждой оси. Если же мы попробуемотмасштабировать график так, чтобы всё выглядело как на самом деле, тополучится нечто вроде эллипса. Это можно сделать перед вызовом функции plt.show. Так и сделаем:

plt.scatter(r[:,0], r[:,1])

plt.axis(‘equal’)

plt.show()

Теперь видим данные с равными масштабами по осям координат. Тут у вас может возникнуть вопрос: а возможно ли создать эллипсы, не привязанные к осям координат? Это значит, что у нас будет полная ковариационная матрица с ненулевыми элементами, лежащими вне диагонали, и это значит, что размерности будут коррелировать друг с другом. Мы разберём этот вопрос на следующей лекции по Scipy.

Выборка в общем многомерном нормальном случае

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

Итак,определим ковариационную матрицу:

cov = np.array([[1, 0.8], [0.8, 3]])

Этозначит, что дисперсия для первой размерности равна 1, для второй размерностиравна 3, а ковариация между этими двумя размерностями равна 0,8. Импортируеммодуль для многомерного нормального распределения и установим среднее значение,равное 0,2:

from scipy.stats import multivariate_normal as mvn

Заглянувв документацию по Scipy, можно увидетьфункцию rvs для выборок с многомерным нормальнымраспределением. При этом в документации по старой версии Scipy есть только два параметра – среднее значение иковариация. Но, разумеется, мы хотим иметь возможность обрабатывать более чемодин пример за раз. В документации по новейшей версии Scipyесть также параметр size, благодарякоторому можно определить нужное количество примеров. К сожалению, в топерезультатов поиска Google до сих порнаходится документация именно по старой версии библиотеки, поэтому будьтевнимательны.

Создадимряд случайных величин:

r = mvn.rvs(mean=mu, cov=cov, size=1000)

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

plt.scatter(r[:,0], r[:.1])

plt.axis(‘equal’)

plt.show()

Каквидим, есть некоторая корреляция. Отметим, что есть возможность отобразить этипримеры с помощью Numpy:

r = np.random.multivariate_normal(mean=mu, cov=cov,size=1000)

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

plt.scatter(r[:,0], r[:.1])

plt.axis(‘equal’)

plt.show()

Всё правильно.

Другие интересные функции Scipy

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

Во-первых, если у вас нет соответствующей теоретической подготовки, то выполнение их в Scipy не имеет большого смысла. Во-вторых, они используются не слишком часто, лишь при необходимости, так что на самом деле нет особой необходимости испытывать их, пока не появится нужда в них. Цель этой лекции – рассказать об их существовании и где их найти, а также подчеркнуть, что это не слишком сложно и что вам не нужен ещё кто-то, чтобы показать, как пользоваться этими функциями. Если вы освоили остальную часть этого курса, то у вас не должно возникнуть проблем с самостоятельным их освоением, они должны работать, как предусмотрено.

Итак, первой мы рассмотрим функцию loadmat. В академических кругах MATLAB является очень популярным компьютерным языком для вычислений, с которым вы периодически будете сталкиваться, а если учитесь, то, наверное, уже сталкивались. В MATLAB есть собственный формат файлов – .mat; он не является текстовым, поэтому его нельзя просто прочитать строка за строкой. Данная функция Scipy позволяет корректно открывать и работать с этими файлами, если у вас появится необходимость работать с данными, хранящимися в формате MAT.

Вдополнение к изображениям существует другой распространённый источник «физических»данных – аудио. Аудио – это звук, а звук, в простейшем формате, хранится в WAV-файлах. Данные файлы содержат амплитуду сигнала вкаждый момент времени; обычно частота дискретизации составляет 44,1 килогерц,что значит, что в каждую секунду имеется 44100 целых чисел, представляющихзвук. Эти файлы можно считывать с помощью функции scipy.io.wavfile.read и записывать с помощью функции scipy.io.wavfile.write.

Оченьтесно с машинным обучением связана обработка сигналов. Одной из популярнейшихфункций для этого является свёртка. Если вы слышали о свёрточных нейронныхсетях, то как раз в них этот метод и работает. В Scipyесть несколько функций для свёртки. Функция scipy.signal.convolve используетсядля любых сигналов любого размера, а для изображений, являющихся двухмернымимассивами (если они чёрно-белые), можно использовать scipy.signal.convolve2d.

Кроме того, в модуле signal есть ряд интересных функций для создания фильтров. Фильтры – это, по сути, формы колебаний, тем или иным образом изменяющие входной сигнал, например ревербация, эффект эха и так далее. Но, как ни странно, преобразование Фурье, являющееся одной из наиболее распространённых функций обработки сигналов, содержится не в Scipy, а в Numpy, так что воспользоваться ею можно, даже не будучи знакомым с библиотекой Scipy. Преобразование Фурье преобразует сигнал из временной области в частотную. Другими словами, оно показывает частотные составляющие исходного сигнала.

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

x = np.linspace(0, 100, 1000)

y = np.sin(x) + np.sin(3*x) + np.sin(5*x)

plt.plot(y)

plt.show()

Вотона – периодическая волна с несколькими частотными составляющими. Теперьсделаем быстрое преобразование Фурье. Не забывайте, что преобразование даёт намсигнал с комплексными числами, поэтому прежде, чем отображать его, необходимонайти модуль:

Y = np.fft.fft(y)

plt.plot(np.abs(Y))

plt.show()

Иполучаем преобразование Фурье. Если увеличить изображение, то можно увидетьпики, представляющие исходные частотные компоненты. Так, мы видим пики в районе80, 48 и 16. Закроем этот график и рассчитаем, какие частоты они представляют:

2*np.pi*16/100

Этодаёт 1.

2*np.pi*48/100

Этодаёт 3.

2*np.pi*80/100

Иэто даёт нам 5.

Это и есть оригинальные частоты. Не переживайте особо, если не понимаете, что это значит, – я лишь хотел показать возможности инструментария Numpy.

Вот и подошел к концу наш курс по инструментарию Numpy. Спасибо, что изучили его и увидимся на следующих занятиях!

Понравилась статья? Поделить с друзьями:
Добавить комментарий

;-) :| :x :twisted: :smile: :shock: :sad: :roll: :razz: :oops: :o :mrgreen: :lol: :idea: :grin: :evil: :cry: :cool: :arrow: :???: :?: :!: