disorder.ru

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

Сбор данных

Как и в прошлый раз, начнем со сбора данных. Рассматриваемое хранилище обслуживает виртуальные машины в «облаке», предоставляя им образы дисков. Подавляющее большинство клиентов — веб-проекты. Это важное замечание, так как, как будет показано ниже, именно этот факт определяет профиль нагрузки данного хранилища. Не смотря на то, что конечные выводы в статье относятся только к хранилищам для веб-проектов или схожих систем, обслуживающих большие потоки случайных посетителей, описанная методика может быть применена для анализа работы хранилищ любого типа.

Хранилище представляет собой кластер из множества узлов, на которых хранятся данные. Доступ к узлам осуществляется через высокопроизводительную сеть; все пользовательские данные резервируются. Нагрузка между узлами распределена примерно поровну за исключением тех узлов, которые были введены в эксплуатацию недавно и еще не заполнились полностью. На нескольких произвольных узлах раз в секунду измерялось количество выполненных операций ввода-вывода (IOPS) и записывалось в файл. В прошлый раз мы замеряли отдельно операции чтения и записи, в этот раз учитывалась только их сумма. Это не имеет принципиального значения, но несколько упрощает обработку данных. Для лучшей репрезентативности замеры проводились как на полностью заполненных узлах, так и на полупустых. Выборка данных происходила в будний день во время пиковой нагрузки на хранилище (это около 16:00 по Москве) в течение часа, а точнее 3000 секунд. Столь небольшой промежуток времени связан с тем, что средняя загруженность хранилища быстро меняется в течение дня, и данные за больший интервал получились бы сильно «смазанными». Вот так выглядит один из полученных файлов:

$ head node1.dat
675
638
365
467
627
429
488
549
702
409

Введение в R

В предыдущей статье для анализа данных мы использовали самописные скрипты на Ruby. В этот раз воспользуемся специально предназначенным для этого инструментом — средой для статистических расчетов R. Столь странное название обусловлено тем, что это свободная реализация языка S, разработанного Bell Labs (создатели UNIX) для обработки данных. S, видимо, происходит от слова statistics и согласуется с названиями других языков, разработанных этой компанией — B и C. В сфере статистических расчетов R является бесплатной альтернативой таким коммерческим пакетам как MATLAB. Существуют версии для Unix, Windows и MacOS. Чтобы установить R на Ubuntu Linux, достаточно выполнить команду

sudo apt-get install r-mathlib

Изучение R можно начать с официального большого введения или с более простого учебника. Для понимания выкладок в этой статье достаточно знать несколько основных вещей.

Запускается R следующим образом:

$ R

R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R -- это свободное ПО, и оно поставляется безо всяких гарантий.
Вы вольны распространять его при соблюдении некоторых условий.
Введите 'license()' для получения более подробной информации.

R -- это проект, в котором сотрудничает множество разработчиков.
Введите 'contributors()' для получения дополнительной информации и
'citation()' для ознакомления с правилами упоминания R и его пакетов
в публикациях.

Введите 'demo()' для запуска демонстрационных программ, 'help()' -- для
получения справки, 'help.start()' -- для доступа к справке через браузер.
Введите 'q()', чтобы выйти из R.

>

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

> x <- 10
> x
[1] 10

Загрузим данные из файла

> x <- scan("node1.dat")
Read 3000 items

Сделаем что-нибудь с этими данными

> summary(x)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    215     834    1152    1331    1561    5471

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

> x <- c(3, 5, 7)
> x
[1] 3 5 7

Вектор последовательных значений можно задать, указав начальное и конечное значение через двоеточие.

> x <- 1:10
> x
 [1]  1  2  3  4  5  6  7  8  9 10

Для индексации элементов в векторе (и массиве) используются квадратные скобки.

> x[2]
[1] 2

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

> ?summary

В R можно создавать собственные функции с помощью функции function(). Определим например функцию возведения в квадрат.

> square <- function(x) { x^2 }
> square(5)
[1] 25

В интерактивном режиме можно вводить относительно небольшое количество команд. Большие порции кода удобно сохранять в файле и затем выполнять командой source().

> source("program.R")

Функции в R сгруппированы по библиотекам. Для загрузки необходимой библиотеки служит команда library().

> library(MASS)

Для R существует огромное количество сторонних пакетов на все случаи жизни. Все пакеты собраны в одном месте — CRAN. Установить нужный пакет можно прямо из интерпретатора командой install.packages().

> install.packages("actuar")

Символ «точка» в R не имеет специального значения и часто используется просто как разделитель в именах функций и переменных.

Предварительный анализ данных

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

> x <- scan("node1.dat")
Read 3000 items
> hist(x, freq = FALSE, breaks = "FD")

По оси абсцисс отложены IOPS, по оси ординат — плотность вероятности. Параметр freq = FALSE указывает, что нам нужна гистограмма плотностей, а не частот значений; параметр breaks = «FD» задает алгоритм разбивки данных на столбцы, удобный для просмотра.

Теперь построим график функции плотности вероятности. Непосредственно плотность распределения можно получить с помощью функции density(), а для ее визуализации воспользуемся функцией plot().

> plot(density(x))

И, наконец, построим график функции распределения вероятностей, используя plot() и ecdf().

> plot(ecdf(x))

По умолчанию все графические результаты выполнения операций отображаются на экране. Чтобы сохранить график в файле, например «hist.png», необходимо выполнить следующие команды:

> png("hist.png")
> hist(x, freq = FALSE, breaks = "FD")
> dev.off()

Так же легко можно вычислить и основные параметры распределения: математическое ожидание, дисперсию и стандартное отклонение.

> mean(x)
[1] 1331.184
> var(x)
[1] 591465.7
> sqrt(var(x))
[1] 769.0681

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

Подгонка распределения

Выводом математических законов, описывающих с той или иной степенью точности наблюдаемые случайные процессы, занимается математическая статистика; она в свою очередь опирается на теорию вероятностей. Типичной задачей в математической статистике является подгонка распределения (distribution fitting), т. е. выяснение, какое из известных в теории вероятностей распределений достаточно хорошо описывает наблюдаемую случайную величину и оценка параметров (parameter estimation) этого распределения на основании выборочных данных (sample data), т.е. данных, полученных в результате непосредственного наблюдения случайной величины. Методика решения этой задачи с помощью R хорошо описана в статье «Fitting distributions with R».

Прежде всего нужно выбрать предполагаемую функцию распределения. Это можно сделать либо на основании знаний о природе наблюдаемой случайной величины, либо анализируя выборку данных. Для начала выберем самый наивный и простой путь: посмотрим, форма кривой плотности какого из известных распределений наиболее похожа на кривую плотности наших выборочных данных. Это удобно делать, глядя на список популярных распределений из Википедии. Строго говоря, наше распределение должно быть дискретным, так как количество IOPS — целая величина. Однако мы вполне можем использовать и плотности непрерывных распределений, считая, что они определены только для целых значений аргумента.

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

Подобную форму имеют и другие распределения, но пока остановимся только на этом. Сама функция плотности имеет вид

, где

 — гамма-функция.

Запись означает, что речь идет о семействе функций , каждая из которых определяется двумя параметрами — и . В этом смысле гамма-распределение является двухпараметрическим семейством распределений. Параметр задает форму кривой, а — масштаб. Для оценки параметров распределения в R существует функция fitdistr() из библиотеки MASS, реализующая метод максимального правдоподобия (maximum-likelihood estimation, MLE). Это весьма популярный метод оценки, но он связан с некоторыми громоздкими вычислениями, поэтому наличие специальной функции в R сильно облегчает задачу. В качестве параметров в функцию передаются вектор выборочных данных и имя подгоняемого распределения.

> library(MASS)
> fitdistr(x, "gamma")
      shape           rate
  2.996031e+00   2.250652e-03
 (4.080286e-02) (2.374736e-05)

Результатом выполнения функции являются параметры подгоняемого распределения. Здесь shape = , а rate = , так как в R используется альтернативная параметризация гамма-распределения. В скобках указаны ошибки оценки параметров.

Чтобы оценить, насколько хорошо полученное распределение описывает изучаемую случайную величину, построим на одном графике гистограмму плотности выборочных данных и кривую плотности распределения (probability density function, PDF).

> hist(x, freq = FALSE, breaks = "FD", main = "Fitting gamma distribution, PDF")
> lines(dgamma(x = 0:6000, 2.996031e+00, 2.250652e-03), col = 2)

Функция lines() похожа на plot() и позволяет рисовать поверх уже существующего графика; ее параметр col задает цвет линии. Функция dgamma() возвращает набор значений плотности вероятности гамма-распределения. В качестве параметров ей передаются вектор аргументов (в нашем случае последовательность от 0 до 6000 IOPS), shape и rate.

Как видно, края распределения хорошо подходят под выборку, но середина и особенно пик — плохо. Можно также сравнить графики функций распределения (cumulative distribution function, CDF).

> plot(ecdf(x), main = "Fitting gamma distribution, CDF")
> lines(0:6000, pgamma(0:6000, 2.996031e+00, 2.250652e-03), col = 2)

Функция pgamma() аналогична dgamma(), но возвращает значения функции распределения.

Еще одним графическим методом оценки качества подгонки распределения является график квантилей (quantile). Квантиль — это такое число, что заданная случайная величина не превышает его лишь с указанной вероятностью. Можно рассматривать квантиль как функцию вероятности Q(p), обратную функции распределения вероятностей. Например, глядя на график распределения выборочных данных, приведенный выше, можно сказать, что квантиль 80% выборочных данных приблизительно равен 2000 IOPS. То есть с вероятностью 80% все выборочные значения не превышают 2000 IOPS. На графике квантилей по одной оси откладываются квантили выборочных значений, а по другой — подгоняемого распределения. Если теоретическое распределение хорошо описывает практические данные, то точки на графике должны лежать рядом с прямой y = x. Для построения графика квантилей служит функция qqplot(). В качестве параметров она принимает два набора данных. Чтобы получить набор данных для гамма-распределения, используется функция rgamma(), которая генерирует последовательность из n случайных чисел, подчиняющихся гамма-распределению с заданными параметрами. Дополнительно на график нанесена прямая y = x  при помощи функции abline().

> qqplot(x, rgamma(n = 3000, 2.996031e+00, 2.250652e-03),
+ main = "Fitting gamma distribution, QQ-plot")
> abline(0, 1)

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

Сложное распределение Пуассона

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

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

Известно, что величины, подобные количеству обращений к веб-сайту в единицу времени, подчиняются распределению Пуассона. Другим примером является число людей, подходящих к кассе в магазине в единицу времени. Распределение Пуассона является дискретным распределением, поэтому вместо функции плотности вероятности его характеристикой служит функция вероятности p(k), показывающая вероятность того, что случайная величина X примет значение k. График функции вероятности распределения Пуассона приведен ниже.

Сама функция задается формулой

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

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

Пусть  — случайная величина с распределением Пуассона, а  — одинаково распределенные случайные величины, не зависящие друг от друга и от .  Распределение случайной величины

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

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

Массовое обслуживание и страховые выплаты

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

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

Распределение отдельных выплат зависит от рассматриваемого портфеля страховых договоров. Для страхования здоровья может быть одно распределение, для страхования автомобилей — другое. Одним из распределений, хорошо описывающих разные формы рисков, является распределение Парето. Это распределение играет большую роль не только в страховании, но и в экономике, социологии, геофизике. Кстати, именно с этим распределением связан знаменитый принцип «20/80». График функции плотности вероятности распределения Парето выглядит следующим образом:

Сама функция имеет вид

Параметр задает масштаб и он же определяет «сдвиг» вправо от оси ординат. На рисунке выше . Параметр определяет форму кривой.

Чтобы понять, может ли распределение Парето описывать количество дисковых операций при обращении к веб-сайту, будем рассуждать следующим образом. Существует минимальное количество операций ввода-вывода, необходимых для отображения любой страницы веб-сайта. Это, скажем, запись в лог-файл веб-сервера и выборка каких-то статических или динамических данных. Это же количество будет иметь и максимальную плотность вероятности, так как наиболее популярные страницы веб-сайта стараются оптимизировать так, чтобы они отдавались за минимальное время. Большее количество операций ввода-вывода будет требоваться для менее частых обращений, например добавление комментария или поиск по сайту. Качественно, такие рассуждения соответствуют кривой плотности вероятности на графике. Примем за рабочую гипотезу, что количество операций ввода-вывода на одно обращение к сайту распределено по Парето и посмотрим, какие это даст результаты.

Для решения задач из области страхования, в том числе связанных с распределением Парето, существует пакет actuar (от англ. actuarial science — актуарная наука). Из него нам потребуются функции dpareto() для распределения Парето и aggregateDist() для построения сложного распределения. Чтобы построить сложное распределение Пуассона, в котором суммируемые величины подчиняются распределению Парето, необходимо выполнить следующие команды:

> library(actuar)
> fx <- c(dpareto(x = 1:6000, 10, 1000))
> Fs <- aggregateDist("recursive", model.freq = "poisson", lambda = 10,
+ model.sev = fx, maxit = 6000)
> plot(Fs)

Функция dpareto() аналогична dgamma(), первым аргументом идет вектор абсцисс, далее параметры распределения: . Параметры взяты «с потолка» исключительно для демонстрации. В функцию aggregateDist() передаются имя распределения частоты событий, в нашем случае это распределение Пуассона, параметр распределения и распределение суммируемых величин. В результате получается функция сложного распределения, график которой приведен выше. Чтобы получить плотность вероятности достаточно продифференцировать полученную функцию.

> plot(diff(Fs), main = "Poisson-Pareto PDF")

Форма кривой плотности вероятности вполне соответствует нашим ожиданиям. Теперь совместим гистограмму выборочных данных с кривой плотности вероятности. Надо сказать, что выбранные параметры обоих распределений взяты не совсем «с потолка». Они подобраны так, чтобы полученное распределение оказалось более или менее близко к реальным данным. Подбор делался простым перебором вручную нескольких значений.

> hist(x, freq = FALSE, breaks = "FD", main = "Fitting Poisson-Pareto distribution, PDF")
> lines(diff(Fs), col = 2)

Видно, что даже наспех подобранные параметры дают неплохой результат. Если подобрать их более тщательно, то подгонку распределения можно значительно улучшить. Кроме того, использованные параметры имеют вполне логичный физический смысл: показывает, что средняя интенсивность просмотра составляет 10 запросов к веб-сайту в секунду, говорит о том, что минимальное количество операций ввода-вывода (суммарно чтение и запись), требуемых для отображения страницы, равно 10. Это, в общем-то, вполне согласуется с практикой.

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

У актуариев часто для анализа также есть только распределение суммарных выплат. Поэтому существует множество теоретических работ, направленных на поиск приближений для конечных распределений без анализа отдельных выплат. В частности, в «Seal, H., 1977. Approximations to risk theories F(X, t) by the means of the gamma distribution» в качестве приближения используется гамма-распределение, которое мы уже пытались применить, а в «Dropkin, L.B., 1964. Size of loss distributions in Workmen’s Compensation insurance» и «Bickerstaff, D.R., 1972. Automobile collision deductibles and repair cost groups: The lognormal model» — логнормальное распределение, с которым мы еще не знакомы. Проверим, подойдут ли страховые теории для описания наших процессов.

Логнормальное распределение

Логнормальное распределение описывает случайную величину, логарифм которой подчиняется нормальному распределению; т. е. если — нормально распределенная случайная величина, то величина имеет логнормальное распределение. В отличии от нормального распределения, которое описывает сумму большого количества независимых случайных величин, логнормальное распределение описывает их произведение.

Ниже приведен график функции плотности вероятности логнормального распределения.

Сама функция имеет вид

Параметр задает масштаб, а — форму кривой плотности. В терминологи R они называются meanlog и sdlog соответственно по аналогии с нормальным распределением.

Прежде чем мы попробуем применить логнормальное распределения к нашим данным, познакомимся с еще одним методом оценки параметров распределения — методом моментов. Моменты случайной величины являются ее числовыми характеристиками. k-м начальным моментом случайной величины X называется математическое ожидание величины , k-м центральным моментом называется математическое ожидание величины , где E(X) — математическое ожидание X. Некоторые моменты имеют собственное название, например математическое ожидание E(X) — это первый начальный момент, а дисперсия D(X) — второй центральный момент. Суть метода моментов заключается в том, что если для данного распределения существуют выражения, связывающие параметры распределения с его моментами, то, приравнивая моменты выборочных данных к моментам подгоняемого распределения, получим систему уравнений, решив которую найдем искомые параметры. В отличии от метода максимального правдоподобия метод моментов очень прост и не требует громоздких вычислений.  Недостаток метода в том, что он не всегда работает. Например для того же гамма-распределения на некоторых выборочных данных получаются совсем неверные результаты.

Параметры логнормального распределения связаны с моментами следующими соотношениями:

Применим их в R.

> meanlog <- function(mean, var) { log(mean) - 0.5*log(1 + var/mean^2) }
> sdlog <- function(mean, var) { sqrt(log(1 + var/mean^2)) }
> meanlog(mean(x), var(x))
[1] 7.049817
> sdlog(mean(x), var(x))
[1] 0.5366686

Мы создали две новые функции для вычисления параметров логнормального распределения методом моментов. Сравним полученные параметры с тем, что выдает метод максимального правдоподобия.

> fitdistr(x, "lognormal")
     meanlog        sdlog
  7.056246398   0.518336033
 (0.009463478) (0.006691689)

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

> hist(x, freq = FALSE, breaks = "FD", main = "Fitting log-normal distribution, PDF")
> lines(dlnorm(x = 0:6000, 7.049817, 0.5366686), col = 2)

> plot(ecdf(x), main = "Fitting log-normal distribution, CDF")
> lines(0:6000, plnorm(0:6000, 7.049817, 0.5366686), col = 2)

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

Вернемся теперь к изначальной постановке задачи. Нам необходимо определить такой верхний порог производительности хранилища, чтобы все запросы клиентов в него укладывались с заданной вероятностью, другими словами вычислить квантиль. В прошлый раз в качестве вероятности мы взяли величину 99.7%, которая фигурирует в «правиле трех сигм»; поступим и сейчас также. Вычислить квантиль логнормального распределения можно с помощью функции qlnorm().

> qlnorm(0.997, 7.049817, 0.5366686)
[1] 5036.492

Сравним это значение с квантилем выборочных данных.

> quantile(x, 0.997)
   99.7%
5021.006

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

Окончательная обработка данных

Для автоматизации обработки данных с узлов хранилища была написана функция iostat.fitdistr(), получающая на входе имя файла с данными. Результатом ее работы являются графики логнормального и гамма-распределения, а также значения квантилей 99.7%. Параметры логнормального распределения оцениваются методом моментов, параметры гамма-распределения — методом максимального правдоподобия. Ниже приведен код функции и результат ее работы для нескольких узлов с разной заполненностью.

iostat.fitdistr <- function(file)
{
	p <- 0.997
	png(sprintf("%s.png", file))

	x <- scan(file)
	mean <- round(mean(x))
	var <- round(var(x))
	qs <- round(quantile(x, p))
	hist(x, freq = FALSE, breaks = "FD",
	    main = paste("Fitting distribution to", file), xlab = "IOPS",
	    xlim = c(0, 7000))

	meanlog <- log(mean) - 0.5*log(1 + var/mean^2)
	sdlog <- sqrt(log(1 + var/mean^2))
	ql <- round(qlnorm(p, meanlog, sdlog))
	lines(dlnorm(x = 0:7000, meanlog, sdlog), lwd = 3, col = 3)

	d <- fitdistr(x, "gamma")
	shape <- d[[1]][[1]]
	rate <- d[[1]][[2]]
	qg <- round(qgamma(p, shape, rate))
	lines(dgamma(x = 0:7000, shape, rate), lwd = 3, col = 2)

	legend("topright", c(
	    sprintf("Sample data\tQ(%.1f%%) = %d IOPS", 100*p, qs),
	    sprintf("Log-normal\tQ(%.1f%%) = %d IOPS", 100*p, ql),
	    sprintf("Gamma\t\tQ(%.1f%%) = %d IOPS", 100*p, qg)),
	    col = c(1, 3, 2), lwd = 3)
	dev.off()

	c(mean, var)
}

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

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

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

> sapply(c("node1.dat", "node2.dat", "node3.dat"), iostat.fitdistr)
Read 3000 items
Read 3000 items
Read 3000 items
     node1.dat node2.dat node3.dat
[1,]      1331      1563      1910
[2,]    591466    507101    618331

Функция sapply() последовательно применяет функцию iostat.fitdistr() к элементам вектора, передаваемого в качестве первого аргумента. Возвращаемые значения заносятся в таблицу, после чего их можно анализировать.

Возьмем побольше узлов и посмотрим на зависимость дисперсии от математического ожидания.

> a <- sapply(c("node1.dat", "node2.dat", "node3.dat", "node4.dat", "node5.dat", "node6.dat"),
+ iostat.fitdistr)
Read 3000 items
Read 3000 items
Read 3000 items
Read 3000 items
Read 3000 items
Read 3000 items
> plot(a[1,], a[2,], main = "D = f(E)", xlab = "Mean", ylab = "Variance")

Запись a[1,] означает выборку первой строки из таблицы. Конечно, по столь малому количеству точек и большому их разбросу нельзя с уверенностью говорить о какой-либо явной зависимости. Но можно предположить линейную зависимость, во всяком случае на наблюдаемом интервале. Так как наше предположение о линейной зависимости весьма шатко, смысла вычислять наилучшие коэффициенты прямой нет, достаточно подобрать их вручную буквально за несколько попыток.

> lines(a[1,], sapply(a[1,], function(x) { 500*(x - 400) }))

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

> mem <- c(160, 146, 183, 210, 106, 130)

Построим аналогичным образом зависимость математического ожидания от объема памяти.

> plot(mem, a[1,], main = "E = f(memory)", xlab = "GB", ylab = "Mean")

Опять же с некоторой натяжкой можно предположить линейную зависимость. Подберем параметры прямой.

> lines(mem, sapply(mem, function(x) { 23*(x - 90) }))

Вот и все. Осталось последние две эмпирические формулы оформить в виде функции и метод вычисления необходимой производительности хранилища данных готов. Код функции iostat.quantile(), вычисляющей необходимую производительность хранилища в зависимости от суммарного объема памяти виртуальных машин и вероятности выполнения запросов приведен ниже.

iostat.quantile <- function(mem, p)
{
	mean <- 23*(mem - 90)
	var <- 500*(mean - 400)
	meanlog <- log(mean) - 0.5*log(1 + var/mean^2)
	sdlog <- sqrt(log(1 + var/mean^2))
	round(qlnorm(p, meanlog, sdlog))
}

Итак, допустим, нам нужно хранилище для обслуживания виртуальных машин с суммарным объемом физической памяти в один терабайт, при этом все операции ввода-вывода должны выполняться с вероятностью не менее 99.9%.

> iostat.quantile(1024, 0.999)
[1] 33796

Получается, что производительность такого хранилища должна быть около 34000 IOPS. Этого хватит на 500 средних виртуальных машин, каждая по 2ГБ ОЗУ. Нужно заметить, что в проанализированном хранилище все данные резервируются, и измеренная в самом начале нагрузка содержит в себе затраты на репликацию данных; соответственно вычисленная мощность хранилища уже имеет как минимум двухкратный запас.

Заключение

В результате всех наших исследований нам удалось показать, что распределение количества дисковых операций хранилища данных, обслуживающего большое количество веб-проектов, хорошо аппроксимируется логнормальным распределением; а также удалось вывести формулу (точнее написать функцию на языке R), позволяющую рассчитать производительность такого хранилища, необходимую для выполнения заданного SLA для заданного числа клиентов. Однако главной целью статьи было продемонстрировать некоторые практические методы математической статистики, позволяющие выявлять и использовать закономерности случайных процессов. Кроме того было показано, как один и тот же математический аппарат позволяет описывать такие разные сферы человеческой деятельности, как страхование рисков и проектирование хранилищ данных.