Многомерный анализ

Нейрон

Нейроны, используемые в большинстве нейронных сетей, состоят из сумматора и нелинейного преобразователя (рис. 7.2). Сигналы, приходящие от нейронов предыдущего слоя, умножаются на весовые коэффициенты связей и суммируются в сумматоре. Полученное значение нелинейно преобразуется таким образом, чтобы результат снова оказался в интервале [-1, 1], например, с помощью функции активации

или

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

Рис. 7.2. Нейрон

Обучение нейронных сетей

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

Примеры применений нейронных сетей

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

– оптическое распознавание символов, включая распознавание подписи (в частности, разработана система идентификации подписи, учитывающая не только окончательный ее рисунок, но и скорость ручки на различных участках, что значительно затрудняет подделку чужой подписи);

– обработка изображений (например, система сканирует видеоизображения станций лондонского метро и определяет, насколько станция заполнена народом, причем работа системы не зависит от условий освещенности и движения поездов);

– прогнозирование финансовых временных рядов (компания LBS Capital Management объявила о значительных успехах в финансовых операциях, достигнутых за счет прогнозирования цен акций с помощью многослойных персептронов);

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

Нейросети используются фирмой Amoco для выделения характерных пиков в показаниях сейсмических датчиков. Надежность распознавания пиков - 95% по каждой сейсмо-линии. По сравнению с ручной обработкой скорость анализа данных увеличилась в 8 раз;

– медицинская диагностика (например, прогнозирование эпилептических припадков, определение размеров опухоли простаты).

– синтез речи (знаменитая экспериментальная система Nettalk, способная произносить фонемы из написанного текста);

– прогнозирование хаотических временных рядов (целый ряд исследований продемонстрировал хорошие способности нейронных сетей к прогнозированию хаотических временных данных);

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

Ford Motors Company внедрила у себя нейросистему для диагностики двигателей после неудачных попыток построить экспертную систему, т.к. хотя опытный механик и может диагностировать неисправности, он не в состоянии описать алгоритм такого распознавания. На вход нейро-системы подаются данные от 31 датчика. Нейросеть обучалась различным видам неисправностей по 868 примерам. "После полного цикла обучения качество диагностирования неисправностей сетью достигло уровня наших лучших экспертов, и значительно превосходило их в скорости";

– лингвистический анализ (пример: сеть с неконтролируемым обучением используется для идентификации ключевых фраз и слов в языках туземцев Южной Америки).

Набор медицинских, но достаточно разнообразных примеров приведен в статье Д.А.Россиева (Горбань и др., 1998):

– экспертная система для диагностики и лечения артериальной гипертонии;

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

– диагностика заболеваний периферических сосудов, например, определение форм артериита;

– диагностика клапанных шумов сердца и оценка систолической и диастолической фаз сердечного сокращения с постановкой предварительного диагноза;

– диагностика заболеваний печени по лабораторным данным исследования функций печени, а также дифференциальной диагностика заболеваний печени и желчного пузыря по УЗИ.

– ранняя диагностика с 80%-й точностью меланом кожи – одного из самых злокачественных заболеваний;

– прогнозирование осложнений инфаркта миокарда;

– система назначения оптимальной стратегии лечения больных облитерирующим тромбангиитом и прогнозирования его непосредственных исходов;

– система дифференциальной диагностики “острого живота”;

Группа НейроКомп из Красноярска (под руководством Александра Николаевича Горбаня) совместно с Красноярским межобластном офтальмологическом центром им. Макарова разработали систему ранней диагностики меланомы сосудистой оболочки глаза. Этот вид рака составляют почти 90% всех внутриглазных опухолей и легко диагностируется лишь на поздней стадии. Метод основан на косвенном измерении содержания меланина в ресницах. Полученные данные спектрофотометрии, а также общие характеристики обследуемого (пол, возраст и др.) подаются на входные синапсы 43-нейронного классификатора. Нейросеть решает, имеется ли у пациента опухоль, и если да, то определяет ее стадию, выдавая, кроме этого, процентную вероятность своей уверенности (Горбань и др., 1998);

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

Исследователь вправе выбрать любое расстояние (меру сходства или различия), которое считает нужным, исходя из содержательных соображений. Например, в зоогеографических исследованиях часто применяется индекс сходства Жаккара-Наумова между вариантами населения. Большой список индексов сходства и мер различия приведен в работе Ю.А.Песенко (1982). Уместно заметить, что мера сходства между признаками – коэффициент корреляции Браве-Пирсона – тоже неевклидова, если рассматривать их как объекты в двойственном пространстве. Однако методы работы с неевклидовыми расстояниями разработаны гораздо хуже.

Термокарты (heatmaps) и иерархическая кластеризация

Пусть имеется таблица “объект – признак”. Простейший способ получить визуальное представление о всей таблице сразу – это ее раскрасить (рис. 7.3, 7.4, пример условный). Раскраска осуществляется следующим образом. Каждому значению таблицы сопоставляется отдельная клетка. Клетка раскрашивается в зеленый цвет, если значение меньше среднего (по столбцу), и в красный, если значение больше. Причем, чем больше значение по абсолютной величине, тем цвет ярче. В черный (или серый, или белый) красятся клетки, значении в которых близки к среднему. Иногда вместо зеленого и красного используются синий и желтый цвета – для лиц с ограниченным цветовосприятием (дальтоников). Раскраска таблиц широко применяется в работах молекулярных генетиков, однако нет никаких причин не применять ее в других областях биологии, где требуется кластерный анализ.

Но раскрашенная таблица выглядит очень пестро, если ее не структурировать. Для этого используется кластерный анализ. Кластерный анализ – это разбиение исходного множества объектов на классы таким образом, чтобы близкие объекты попали в одни и те же классы, а далекие – в разные. Мера сходства или различия может быть измерена в количественной или даже ранговой шкале. Один из самых популярных способов структурирования – иерархическая классификация. “Иерархическая” означает, что каждый класс вложен в некоторый другой. Самый известный и часто используемый алгоритм иерархической классификации – алгоритм ближайшего соседа или единственной связи. Вначале каждый объект считается отдельным классом. На следующем шаге ищется пара самых близких объектов, которая объединяется в новый класс. Расстояния (или меры сходства) для нового класса со старыми пересчитываются по следующему правилу: расстоянием между классами считается расстояние между ближайшими объектами в этих классах (отсюда и название). Далее все повторяется до тех пор, пока не останется ровно один класс, содержащий все объекты. Если за расстояние между классами принять расстояние между самыми далекими объектами в этих классах, то получим метод дальнего соседа или полной связи. Можно также за расстояние между классами принять среднее расстояние между объектами этих классов, тогда получим метод UPGMA или средней связи. Несколько особняком стоит метод Уорда, в котором учитывается еще и разброс объектов внутри кластера.

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

Рис. 7.3. Термокарта и дендрограммы (условный пример)

Рис. 7.4. Термокарта и дендрограммы

после перестановки строк

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

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

В отличие от иерархической классификации, все классы равноправны и находятся на одном уровне. Еще одна особенность – классы не обязательно удовлетворяют так называемому условию “компактности”, т.е. не являются “хорошими” или “естественными” в том смысле, что ближайшими к некоторым объектам одного класса, могут быть объекты из другого класса. Поэтому некоторые авторы предпочитают называть его алгоритмом группировки, а не классификации. Есть критерии, позволяющие оценить удачность разбиения на классы (аналогичные методу Уорда). Если разбиение оказалось не очень удачным, необходимо изменить и весь процесс повторить с другим К. Рекомендуется для начала брать равным квадратному корню из числа объектов, однако это сильно зависит от исследуемого множества.
Многомерное шкалирование

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

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

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

частным случаем которой является евклидова метрика (при ). Критерием служит соответствие между расстояниями в этом пространстве и исходной матрицей сходства-различия между объектами.

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

Поэтому последние несколько десятков лет, в основном, используется неметрическое шкалирование в квазиметрическом варианте, восходящем к Крускалу (Kruskal, 1964а, 1964б), хотя оно требует очень много машинного времени и поэтому число объектов, которое можно обработать этим методом на персональных компьютерах с помощью профессиональных статистических пакетов не превышает сотни.

Пусть имеется конечное множество объектов и матрица различий или мер сходства между ними, а также произвольное представление объектов этого множества в виде точек в метрическом пространстве размерности с метрикой . Определим критерий различия между множеством и его представлением (“стресс” по Крускалу) в виде

где – некоторое монотонное преобразование.

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

Популярность этого метода объясняется исключительно тем, что ему не было альтернативы. Ситуация радикально изменилась после появления работы Й.Тагучи и Й.Ооно (Taguchi, Oono, 2005), в которой произошел возврат к первоначальной идее Р.Шепарда (Shepard, 1962) и неметрическому шкалированию, образно выражаясь, вернули права гражданства. Теперь речь идет об обработке тысяч и десятков тысяч объектов без потери качества метризации, что открывает огромные перспективы для исследователей во всех областях знаний.

Алгоритм Шепарда-Тагучи-Ооно работает следующим образом. Исходные оценки различия ранжируются. Выбирается размерность и метрика результирующего пространства. В этом пространстве случайным образом помещается совокупность точек, каждая из которых соответствует одному объекту. Между ними вычисляется матрица расстояний, которая также ранжируется. Каждой из пар объектов соответствует два ранга, в одной и другой ранжировке. Если ранжировки полностью соответствуют друг другу, то первый этап работы алгоритма закончен. Если нет, то имеется пара объектов, для которых ранги в двух ранжировках различны. Если ранг расстояния в результирующем пространстве больше ранга различия той же пары объектов в исходной матрице, то точки, представляющие объекты, чуть-чуть сдвигаются друг к другу, если меньше – раздвигаются. После прохождения всех пар объектов расстояния между точками результирующего пространства пересчитываются и ранжируются заново. Процесс продолжается до тех пор, пока сходство между ранжировками, например, ранговый коэффициент корреляции Спирмена, не перестанет расти. Если оно слишком мало, размерность пространства увеличивается на единицу и весь процесс повторяется. Скорость этого алгоритма оказалась, по меньшей мере, на порядок больше, чем алгоритма Крускала, что позволяет обрабатывать значительно большее число исходных данных.

Почему ранговые оценки сходства различий позволяют с такой большой точностью восстановить метрическую структуру данных? На этот вопрос лучше всего ответил сам автор неметрического шкалирования. “Парадоксальная возможность восстановления количественной структуры из качественных данных связана с тем обстоятельством, что число пар точек и, следовательно, число порядковых ограничений на их расстояния возрастает приблизительно как квадрат числа определяемых количественных координат точек. Такие методы называются «неметрическими», поскольку в этом случае используются только порядковые свойства входных данных. Однако выход может достигать большой метрической точности и всегда будет метричным в смысле соответствия аксиомам расстояния.” (Шепард, 1980).

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

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

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

Содержательная интерпретация полученных результатов, в силу наглядности представления, обычно не вызывает особенных затруднений. Рассмотрим, например, матрицу дивергенции (%) 402 пн участка гена цитохрома b мтДНК лесных и полевых мышей рода Apodemus, приведенную в статье Челоминой и др. (Генетика, 1998, т.34, №5, 650-661) .

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

Существуют и другие, более сложные, модели неметрического шкалирования, основанные на одновременном рассмотрении целого ряда матриц (Caroll, Chang, 1970; Caroll, 1976), на которых мы не будем останавливаться.

Бутстреп

Bootstrap (англ.) – ремешок на заднике ботинка, облегчающий его надевание. В английском языке существует идиома – lift oneself by one's own bootstrap – дословно, поднять самого себя за ремешок собственного ботинка. В переносном смысле – выбиться в люди благодаря собственным усилиям; самому пробить себе дорогу; быть всем обязанным самому себе. В статистике так называется процедура, предложенная Б. Эфроном (Efron, 1979, 1982; Диаконис, Эфрон, 1983). Предположим, что у нас есть данные и некоторая последовательность вычислительных действий, например, с использованием нейронных сетей или многомерного шкалирования или любых других эвристических алгоритмов. Мы хотим иметь представление о статистической устойчивости результатов расчета.

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

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

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

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

ЛЕКЦИЯ 8. Временные ряды
Устойчивость статистических связей

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

(1)

где – показатель численности в момент t,

– внешний по отношению к популяции фактор,

– остаток или “шум” без автокорреляции (Moran, 1953).

Кроме того, выбор класса схем и подгонка коэффициентов должны осуществляться на одной части статистического материала, а оценка соответствия – на другой (Колмогоров, 1933, 1986; Дрейпер, Смит, 1987). Для временных рядов обычно "зажимают" несколько последних элементов ряда, оценивают коэффициенты по оставшимся и проверяют степень расхождения на последней части ряда. По существу, мерой качества статистической модели является ее предсказательная сила. В качестве критериев адекватности избираются различные показатели: минимум дисперсии прогноза, близость спектральных и автоковариационных функций и т.д. (Кашьяп, Рао, 1983). Одним из возможных критериев является коэффициент корреляции.

Будем использовать следующую процедуру. Прогнозируемый ряд и предполагаемый предиктор (ряд, используемый для предсказания) разобьем на две части: обучающую и проверочную. Длину проверочной части будем брать не меньше двух десятков отсчетов, чтобы можно было воспользоваться аппроксимацией Р.Фишера для выборочного коэффициента корреляции (Большев, Смирнов, 1983). Подгонку параметров линейной регрессионной модели будем проводить на обучающей части рядов по минимуму суммы квадратов отклонений, что влечет максимизацию коэффициента корреляции между исходными и расчетными данными. Адекватность модели будем оценивать величиной коэффициента корреляции между проверочной частью исходных данных и соответствующими расчетными данными, вычисленными по подогнанным на обучающей части параметрам. В случае одной независимой переменной речь идет просто о вычислении коэффициентов корреляции между прогнозируемый рядом и предиктором на обучающей и проверочной частях отдельно. Предиктор проходит статистический тест и остается в списке для содержательного рассмотрения, если коэффициенты корреляции с прогнозируемым рядом по обеим частям отдельно достаточно велики и имеют одинаковый знак. Если каждый из коэффициентов корреляции превышает величину, соответствующую некоторому уровню значимости , то вероятность для обоих коэффициентов превысить его одновременно, имея одинаковые знаки, равна . Выбирая стандартный уровень значимости (0.05; 0.01; 0.001), получим, что (0.3162; 0.1414; 0.0447). Само собой разумеется, что значимой на уровне должна быть и корреляция между предиктором и прогнозируемым рядом на всем интервале наблюдений. Уровень значимости одного коэффициента корреляции при заданном числе наблюдений будем оценивать по уровню значимости нормально распределенной случайной величины (Большев, Смирнов, 1983).

В исследованиях динамик численности животных принято логарифмировать данные учетов или заготовок шкурок (Уильямсон, 1975). Это связано с тем, что при отсутствии внутренних и внешних ограничений рост популяции описывается уравнением:

где – коэффициент воспроизводства, – производная по времени.

Это же уравнение можно записать в виде
где (Бигон и др., 1989). Дискретный аналог этого уравнения имеет вид

,

что совпадает с уравнением (1), если положить

Полученные ряды лучше описываются нормальным распределением.

Фазовые портреты. Теорема Такенса. Метод главных компонент для временных рядов

Другим способом анализа временных рядов для выявления внутренне присущих им закономерностей является разложение Карунена-Лоэва (метод главных компонент, разложение на естественные ортогональные составляющие, сингулярный спектральный анализ) (Ефимов и др., 1988; Главные компоненты …, 1997; Бобрецов и др., 2000). Этот метод применим к любому временному ряду, не требует его стационарности, как, например, спектральный анализ, автоматически выявляет тренды, если они имеются, без каких-либо предположений об их природе и форме, и (последнее по счету, но не по важности) позволяет получать многомерные представления временного ряда – фазовые портреты – дающие возможность визуального изучения траектории ряда в многомерном пространстве его состояний.

Сущность метода заключается в следующем. Если временной ряд порождается некоторой динамической системой с конечным числом параметров, то совокупность его отрезков можно рассматривать как точки многомерного фазового пространства. Соединяя их последовательно, например, сплайнами, получим траекторию ряда в этом пространстве, которая, как следует из знаменитой теоремы Такенса (Takens, 1981), воспроизводит многомерный фазовый портрет динамической системы, если длина отрезков превышает удвоенное число параметров. Для редукции размерности динамических систем обычно используют преобразование Пуанкаре (Балеску, 1978). Однако можно применить метод главных компонент, заключающийся в поиске координатных осей, в проекции на которые дисперсия траектории ряда максимальна (Ефимов, Галактионов, 1983; Ефимов и др., 1988; Главные компоненты …, 1997). Максимизация автоковариации вместо дисперсии приводит к методу гладких компонент. Оба метода оказываются исключительно полезны при анализе внутренних закономерностей и прогнозе динамики численности и структуры популяций животных и влияющих на них факторов.

Кроме того, в последнее время получили некоторую популярность вэйвлет-методы (wavelet methods), которые близки по своим принципиальным подходам к фильтрации рядов в методе главных компонент. В этих методах сначала выбирается так называемая «материнская волна», например, «мексиканская шляпа», зависящая от параметров сдвига и сжатия, а потом эта волна применяется в качестве фильтра к исходному ряду при всех возможных значениях этих параметров. Получающаяся при этом поверхность над двумерной плоскостью анализируется визуально. Литературы на русском языке по вэйвлет-методам практически нет, но их описание и матобеспечение на английском языке доступно через Интернет (поиск по ключевому слову “wavelet”).

Обработка одного временного ряда методом главных

компонент

Пусть имеется последовательность наблюдений некоторого показателя в равноотстоящие моменты времени . Выберем в качестве многомерной характеристики процесса в момент времени вектор , именуемый предысторией процесса за время . Параметр называется лагом (запаздыванием). Сведем полученные векторы в таблицу, имеющую строк (объектов) и столбец (признак) (табл. 8.1), которая называется матрицей Ганкеля.

Таблица 8.1

Сдвиг временного ряда на отсчетов

Обработка полученной матрицы методом главных компонент приводит к появлению новой матрицы тех же размеров. Новые признаки (компоненты) являются линейными комбинациями старых

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

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

Кроме того, компоненты имеет смысл использовать в качестве предикторов, так как за каждой компонентой предположительно стоит порождающая ее самостоятельная и статистически независимая от других причина (Ефимов и др., 1988).

В исследования по динамике численности животных фазовые портреты с впервые применил Моран (Уильямсон, 1975). Более сложный случай рассмотрен в работе (Schaffer, 1984). С помощью компьютерной графики исследовалась траектория заготовок шкур канадской рыси в трехмерном пространстве, где выбиралось таким образом, чтобы выйти за пределы значимой корреляции между и . Многомерная траектория динамики заготовок шкурок водяной полевки и ее представление в виде фазового портрета с помощью метода главных компонент впервые рассмотрены нами в публикациях (Ефимов, Галактионов, 1982, 1983; Галактионов и др., 1987) и монографиях (Ефимов и др., 1988) и (Бобрецов и др., 2000).

Очевидно, что сфера применения метода выходит далеко за пределы динамики численности животных и он может быть применен к временным рядам любой природы. Однако история его появления достаточно запутана. Первыми публикациями, относящимися к этому методу, считаются статьи (Colebrook, 1978; Broomhead, King, 1986a, 1986b). Однако его идеи неоднократно и независимо появлялись (и появляются до сих пор) в различных областях знаний, связанных с обработкой временных рядов. Один из обзоров публикаций на эту тему можно найти в сборнике (Главные компоненты … , 1997). Имеются две монографии: элементарное введение в метод (Elsner, Tsonis, 1996) и содержащая его теоретическое обоснование книга (Golyandina et al, 2001).
Обработка нескольких временных рядов методом

главных компонент

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

Так как каждая компонента отражает существующую по какой-либо причине общность временных рядов, постоянную, временную или даже случайную, проявляющуюся в коррелированности рядов, и компоненты не коррелируют между собой, то очень часто является осмысленным предположение, что эти причины также независимы. Если совокупность временных рядов представляет собой регистрацию одного показателя, относящегося к различным точкам или районам территории, то метод главных компонент можно использовать для районирования этой территории по каждой компоненте и, следовательно, отдельно по каждой причине, порождающей общность временных рядов (Ефимов, Галактионов, 1983; Гусев, Ефимов, 1985; Ефимов и др., 1988). Если эта совокупность объединяет группу близких по смыслу показателей, например, динамики урожайности нескольких видов культур, то с помощью главных компонент можно выявить, во-первых, общие для всех или частные для некоторых подгрупп факторы, а, во-вторых, расположить виды в соответствии с чувствительностью и направлением действия этих факторов. Правда, необходимо отметить, что метод главных компонент не предоставляет автоматической интерпретации получаемых факторов и об их смысле приходится догадываться отдельно, что в некоторых случаях представляет непростую задачу и предъявляет довольно высокие требования к квалификации интерпретатора.

Если обрабатывается транспонированная матрица, то временными рядами являются собственные векторы, а вклады признаков отражены в компонентах.

Метод гладких компонент

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

Пусть ,

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

где

Очевидно, что точно так же

и, следовательно, где но матрица уже симметрична.

Максимизируем при условии

Обычными методами дифференциального исчисления (Кульбак, 1967) получим, что вектор удовлетворяет матричному уравнению

,

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

где – ортогональная матрица собственных векторов,

– диагональная матрица собственных значений.

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

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

Существует довольно глубокое и неожиданное сходство между методами главных и гладких компонент для анализа временных рядов и многомерным генетическим анализом. Одним из его приемов является многомерный анализ фенотипической ковариационной матрицы (Thorpe, Leamy, 1983; Falconer, 1989). Более двадцати лет назад была введена матрица генетических ковариаций (Lande, 1979). По смыслу – это коэффициенты корреляции (ковариации) между признаками родителей и их потомков. После ее введения оказалось возможным оценивать аддитивную наследуемость любой линейной комбинации признаков, в том числе, главных компонент матриц и (Atchley et al., 1981). Однако направленный поиск комбинированных признаков с максимальной аддитивной наследуемостью предложен и проведен только в недавнее время. Например, Ott&Rabinowitz (1999) для максимизации аддитивной наследуемости предложили разложение матрицы на собственные вектора. Klingenberg&Leamy (2001) с помощью такого разложения получили линейную комбинацию промеров нижней челюсти с аддитивной наследуемостью 0.73, не совпадающую ни с одной из главных осей матриц и . При этом наследуемость общего размера нижней челюсти на этом же материале равна 0.42, что по порядку совпадает с оценками наследуемостей других краниометрических признаков (Leamy, 1974; Atchley et al., 1981). Таким образом, разложение соответствующих ковариационных матриц может привести к комбинированным признакам с существенно более высокими коэффициентами наследуемости, чем у исходных признаков.

В случае временных рядов матрицы и играют роль матриц и , соответственно, так состояние системы на следующий год является «потомком» по отношению к ее текущему состоянию. Аналогом наследуемости является предсказуемость или прогнозируемость, аналогом ДНК – инерционность системы. В широком смысле и ДНК и инерцию можно рассматривать как формы памяти – нечто инвариантное, наследуемое следующим поколением. Продолжая аналогию дальше, можно ставить вопрос о поиске линейных комбинаций во временных рядах с максимальной прогнозируемостью через разложение матрицы и пытаться выяснить их содержательный смысл.
ЛЕКЦИЯ 9. Большие массивы. Достоверность
Ситуацию с большими массивами лучше всего рассмотреть на конкретном примере. К числу многообещающих направлений современной биологии относится анализ экспрессии генов с помощью микрочипов. Микрочипы – упорядоченные комплекты фрагментов ДНК или РНК (зондов), иммобилизованные на специальных носителях (пластинках из стекла, пластика или кремния, каплях геля), называемых платформами (вместе с сопутствующими технологиями, как правило, разных производителей). Технология микрочипов позволяет одновременно анализировать экспрессию десятков тысяч генов в нескольких десятках проб (образцов), обычно представляемую в виде матрицы уровней экспрессии, в которой столбцам отвечают образцы (объекты), а строкам – гены (признаки). Столбцы и строки этой матрицы называются профилями экспрессии, соответственно, образцов или генов. (Иногда возникает терминологическая путаница: термин “профиль экспрессии генов” в некоторых работах может означать “профиль экспрессии образца по всем генам”.) Технология используется для поиска генов-кандидатов, экспрессия которых различается на нескольких группах образцов, например, принадлежащих пациентам с разными формами туберкулеза, рака или страдающим различными нейродегенеративными заболеваниями (болезни Альцгеймера, Хантингтона и т.д.) (Borovecki et al., 2005; Ефимов, Катохин, 2009). После проверки эффективности найденных генов-кандидатов другими методами (ПЦР) на их основе могут быть сконструированы микрочипы для ранней диагностики этих заболеваний.

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

Как измерять величину различий? Для нескольких групп обычно используется критерий ANOVA (analysis of variance). Идея заключается в том, что считаются среднее всей выборки, средние по группам, отклонения от этих средних внутри каждой группы, отклонения самих средних от среднего всей выборки. Дальше вычисляются межвыборочная и объединенная внутривыборочная дисперсия и берется их отношение F=V_m/V_in. (Межвыборочная – это средний квадрат отклонения групповых средних от среднего всей выборки. Объединенная внутривыборочная – это средний квадрат всех внутривыборочных отклонений от своих групповых средних.) Если справедлива нулевая гипотеза, то отношение межгрупповой к объединенной внутривыборочной не должно быть слишком большим и можно вычислить вероятность случайного получения такого , которое получилось в нашей выборке. Если эта вероятность достаточно мала , то мы отклоняем нулевую гипотезу на первом, втором или третьем уровне значимости и считаем, что в нашем профиле (выборке) есть реальные различия между группами. Если групп всего две, то критерий ANOVA сводится к -критерию Стьюдента.

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

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

Но этот критерий слишком жесткий. Если, например, , то . А может быть и больше (микрочипы с более, чем сотней тысяч зондов, уже в повестке дня). Через критерий Бонферрони проходит слишком мало полезной информации. Особенно это видно на массивах микрочиповых данных. Поэтому в 1995 году Бенджамини и Хохберг предложили критерий FDR (false discovery rate). Упорядочим все профили по возрастанию , рассчитанных для случая одного профиля. Для -го по порядку профиля и вероятности ожидаемое число случайных профилей с вероятностью равной или меньше равно , а реальное число, очевидно, равно . То есть, среди первых по порядку профилей случайных должно быть не больше, чем . Если достаточно мало по сравнению с , например , то мы этим пренебрегаем и отклоняем нулевую гипотезу для всех профилей. Этот критерий работает более разумно и сейчас является основным в микрочиповых исследованиях. Главная работа: (Benjamini, Hochberg, 1995). Есть библиотека для R-пакета: http://strimmerlab.org/notes/fdr.html

<< предыдущая страница следующая страница >>