WWW.BOOK.LIB-I.RU
БЕСПЛАТНАЯ  ИНТЕРНЕТ  БИБЛИОТЕКА - Электронные ресурсы
 
s

Pages:     | 1 | 2 || 4 | 5 |

«В.К. Шитиков, С.Э. Мастицкий Классификация, регрессия и другие алгоритмы Data Mining с использованием R Тольятти, Лондон - 2017 Шитиков В.К., Мастицкий С.Э. (2017) Классификация, ...»

-- [ Страница 3 ] --

Все переменные разделились на две группы: значимые для классификации (x2, x3, x6, x7, x10, x11, x14, x15) со значением t = ±1.032 и незначимые (x1, x4, x5, x8, x9, x12, x13, x16) с t = 0.

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

К сожалению, использование функции train() из пакета caret как с параметрами method = "glm", family=binomial, так и method = "binda" привело к неудаче. Видимо при n = 16 еще не достигнут порог репрезентативности исходных данных, при котором эта функция начинает стабильно выдавать адекватные результаты.

library(caret) (binda.train - train(Xtrain, Class, method = "binda"))

Resampling results across tuning parameters:

lambda.freqs Accuracy Kappa Accuracy SD Kappa SD 0.0 0.3435714 -0.11575021 0.1753883 0.2124598 0.5 0.3435714 -0.10959636 0.1753883 0.2135358 1.0 0.3649048 -0.08130732 0.1929750 0.2404408 Accuracy was used to select optimal model using largest value.

The final value used for the model was lambda.freqs = 1.

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

Однако итоговая модель подобных неожиданностей не несет:

binda.train$finalModel pred.train - predict(binda.train, Xtrain) CTab - table(Факт=DFace$Class, Прогноз=pred.train) Acc = mean(DFace$Class == pred.train) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт 1 2 [1] "Точность=87.5%"

5.2. Бинарные деревья решений Как было отмечено в разделе 4.3, использование иерархических древовидных структур (decision trees) является одним из наиболее универсальных и эффективных методов классификации и регрессии (Hunt et al., 1966; Breiman at al., 1984; Loh and Shih, 1997). Классификационные модели деревьев рекурсивно делят набор данных на подмножества, которые являются все более и более гомогенными относительно определенных признаков. Создается решающее правило классификации иерархического типа и формируется ассоциативный логический ключ, дающий возможность выполнять распознавание объектов из новых выборок.

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

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

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

Мера информационного выигрыша IG (Information Gain) для разбиения A определяется как разность энтропии Шеннона H(S) для исходного набора данных и суммы энтропий всех фрагментов данных после разбиения Н(Sv) с соответствующим весом в зависимости от размера каждой части:





| Sv | | S | Н(Sv), IG (S, А) = H(S) vV где v – выбранное опорное значение переменной V, использованное для разбиения А. Если разделение на классы выполнено оптимально, то каждое подмножество Sv будет иметь небольшую энтропию, и, следовательно, значение IG будет максимальным.

На языке R этот критерий может быть рассчитан следующим образом:

Entropy - function( vls ) { res - vls/sum(vls)*log2(vls/sum(vls)) ; res[vls == 0] - 0

-sum(res) } InformationGain - function( tble ) { entropyBefore - Entropy(colSums(tble)) s - rowSums(tble) entropyAfter - sum(s/sum(s)* apply(tble, MARGIN = 1, FUN = Entropy )) informationGain - entropyBefore - entropyAfter return (informationGain) } Для реализации алгоритма ID3 используем компоненты пакета data.tree: объект класса tree (дерево) и метод AddChild, добавляющий к выстраиваемому дереву очередной дочерний узел (http://ipub.com/wp-content):

IsPure - function(data) { length(unique(data[,ncol(data)])) == 1 } TrainID3 - function(node, data) { node$obsCount - nrow(data) # если текущий набор данных принадлежит к одному классу, то if (IsPure(data)) { #создается лист дерева с экземплярами этого класса child - node$AddChild(unique(data[,ncol(data)])) node$feature - tail(names(data), 1) child$obsCount - nrow(data) child$feature - '' } else { # рассчитывается вектор информационных выигрышей IG ig - sapply(colnames(data)[-ncol(data)], function(x) InformationGain( table(data[,x], data[,ncol(data)]))) # выбирается значение признака с наибольшей величиной IG feature - names(which.max(ig)) node$feature - feature # создается подмножества данных на основе этого значения childObs - split(data[, names(data) != feature, drop = FALSE], data[,feature], drop = TRUE) # создаются дочерние узлы дерева c именем признака for(i in 1:length(childObs)) { child - node$AddChild(names(childObs)[i]) # осуществляется рекурсия алгоритма для дочернего узла TrainID3(child, childObs[[i]]) } } } Обратим внимание, что наша функция TrainID3 является рекурсивной, т.е. для выполнения итераций вызов функции осуществляется из тела ее самой. Функция завершает работу, если все ветви дерева завершаются листьями с объектами одного класса.

Для иллюстрации работы алгоритма ID3 обратимся к нашему примеру по классификации лиц избирателей (см. рис. 5.1 в разделе 5.1). Для получения осмысленной визуализации дерева перейдем от численных обозначений к предметным названиям признаков и их градаций.

В результате получаем следующее дерево решений:

DFace - read.delim(file="Faces.txt",header=TRUE, row.names=1) DFaceN -DFace[,-17]; Class - DFace$Class DFaceN[DFaceN==1]- "Да" ; DFaceN[DFaceN==0] - "Нет" Class[Class==1] - "Патриот" ; Class[Class==2] - "Демократ" DFaceN - cbind(DFaceN,Class) colnames(DFaceN) - c("голова_круглая","уши_оттопырен", "нос_круглый","глаза_круглые","лоб_морщины", "носогубн_складка","губы_толстые","волосы","усы","борода", "очки","родинка_щеке","бабочка","брови_подняты","серьга", "курит_трубка","Группа") library(data.tree) tree - Node$new("DFaceN") # Создаем «пустое» дерево TrainID3(tree, DFaceN) print(tree, "feature", "obsCount")

–  –  –

Predict - function(tree, features) { if (tree$children[[1]]$isLeaf) return (tree$children[[1]]$name) child - tree$children[[features[[tree$feature]]]] return ( Predict(child, features)) } pred - apply(DFaceN[,-17],1, function(x) Predict(tree,x)) table(Факт=DFaceN$Группа, Прогноз=pred) Прогноз Факт Демократ Патриот Демократ 8 0 Патриот 0 8 Предсказание электоральной группы по всей обучающей выборке, как и в случае логистической регрессии, также оказалось безошибочным. Однако нет уверенности, что построенные деревья будут столь же успешны при практической реализации на "свежих" данных.

Для оценки эффективности модели выполним перекрестную проверку:

Nerr - 0 for (i in 1:nrow(DFaceN)) { tree - Node$new("DFaceN") TrainID3(tree, DFaceN[-i,]) Nerr = Nerr + sum(DFaceN[i,17]!=Predict(tree,DFaceN[i,-17])) } (Nerr/nrow(DFaceN)) [1] 0.25 При выполнении скользящего контроля четыре объекта из 16 оказались неверно предсказанными с помощью бинарного дерева, построенного после исключения тестируемого примера. Считается (Дюк, Самойленко, 2001), что алгоритм ID3 в случае взаимозависимых признаков нередко направляет ход логического вывода по ложному пути и создает лишь иллюзию правильного рассуждения. Но для более обоснованных выводов следует рассматривать результаты тестирования на более серьезных примерах, чем наш.

5.3. Поиск логических закономерностей в данных Основное требование к математическому аппарату Data Minning заключается не только в эффективности классификаторов, но и в интерпретируемости результатов. Методы поиска логических закономерностей в бинарных матрицах наблюдений (как и в некоторых классических методах многомерного анализа) апеллируют к информации, заключенной не в отдельных признаках, а в различных сочетаниях их значений. Для признаков, измеренных в альтернативных шкалах, их значения xi = 0 или xi = 1 рассматриваются как элементарные события, что позволяет сформулировать решающие правила на простом и понятном человеку языке логических высказываний, таких как ЕСЛИ {(событие 1) и (событие 2) и... и (событие N)} ТО...

Так, классификация лиц в рассмотренном примере рис.

5.1 может быть произведена, например, с помощью логических высказываний 1 и 2:

1. ЕСЛИ {(голова овальная) и (есть носогубная складка) и (есть очки) и (есть трубка)} ИЛИ {(глаза круглые) и (лоб без морщин) и (есть борода) и (есть серьга)} ТО ("Патриот")

2. ЕСЛИ {(нос круглый) и (лысый) и (есть усы) и (брови подняты кверху)} ИЛИ {(оттопыренные уши) и (толстые губы) и (нет родинки на щеке) и (есть бабочка)} ТО ("Демократ").

Математическая запись этих правил выглядит следующим образом:

[(x1=0) ^ (x6=1) ^ (x11=1) ^ (x16=1)] [(x4=1) ^ (x5=0) ^ (x10=1) ^ (x15=1)] = 1 [(х3=1) ^ (х8=0) ^ (х9=1) ^ (х14=1)] [(х2=1) ^ (x7=1) ^ (х12=0) ^ (х13=1)] = 2.

Здесь значки – конъюнкция ("и"), ^ – дизъюнкция ("или"), = – импликация ("если, то").

Будем использовать логический метод распознавания, известный под названием алгоритм "Кора", который в свое время широко использовался в геологоразведочных работах и скрининге лекарственных препаратов (Бонгард, 1967; Вайнцвайг, 1973). В детерминистической версии алгоритма обучающая выборка x, предварительно разделенная на два класса (1 и 2), многократно просматривается и из всего множества высказываний выделяются так называемые непротиворечивые логические высказывания (x, ), покрывающие все множество примеров. Непротиворечивым высказыванием для каждого класса считается конъюнкция, которая встречается с вероятностью не менее только в одном классе и ни разу не встречается в другом.

Представленный ниже скрипт выполняет перебор всех возможных комбинаций столбцов хi, i = 1, …, m=16, по 2, 3 и maxKSize = 4 с использованием функции combn(). Для каждой комбинации столбцов перебираются все возможные варианты событий (xi = 0 или xi = 1). Для сокращения объема вычислений подмножества исходной матрицы трансформировались в векторы символьных бинарных переменных, а комбинации значений xi выражались наборами "битовых масок" (например, "000", "001", "010", "011", "100", "101", "110", "111").

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

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

# Преобразование числа в набор битов (5 - "0101") number2binchar - function(number, nBits) { paste(tail(rev(as.numeric(intToBits(number))), nBits), collapse="") } # Поиск конъюнкций по заданной битовой маске MaskCompare - function(Nclass, KSize, BitMask, vec_pos, vec_neg, ColCom) { nK - sapply(BitMask, function(x) { if (sum(x==vec_neg)0) return (0) if (minNum(countK= sum(x==vec_pos))) return (0) # Cохранение конъюнкции в трех объектах list Value.list[[length(Value.list)+1]] list(Nclass=Nclass, KSize=KSize, countK=countK, Bits=x) ColCom.list[[length(ColCom.list)+1]] - list(ColCom) RowList.list[[length(RowList.list)+1]] list(which(vec_pos %in% x)) return (countK) } ) } Зададим минимальную частоту встречаемости конъюнкций minNum = 4 (т.е. = 50%) и выполним формирование всех логических правил для рассматриваемого примера:

DFace - read.delim(file="Faces.txt",header=TRUE, row.names=1) maxKSize - 4 minNum - 4 # Списки для хранения результатов Value.list - list() # Nclass, KSize, BitMask, countK ColCom.list - list() # Наименования переменных ColCom RowList.list - list() # Номера индексов строк RowList # Перебор конъюнкций разной длины for (KSize in 2:maxKSize) { BitMask - sapply(0:(2^KSize-1), function(x) number2binchar(x,KSize)) cols - combn(colnames(DFace[,-17]),KSize) for (i in 1:ncol(cols)) { SubArr - DFace[,(names(DFace) %in% cols[,i])] vec1 - apply(SubArr[DFace$Class==1,],1, function (x) paste(x,collapse="")) vec2 - apply(SubArr[DFace$Class==2,],1, function (x) paste(x,collapse="")) MaskCompare(1, KSize, BitMask, vec1, vec2, cols[,i]) MaskCompare(2, KSize, BitMask, vec2, vec1, cols[,i]) } } # Создание результирующей таблицы DFval = do.call(rbind.data.frame, Value.list) nrow = length(Value.list) DFvar - as.data.frame(matrix(NA, ncol=maxKSize+1, nrow=nrow, dimnames = list(1:nrow, c( paste("L", 1:maxKSize, sep=""),"Объекты:")))) for (i in 1:nrow) { Varl - unlist(ColCom.list[[i]]) DFvar[i, 1:length( Varl)] - Varl Objl - unlist(RowList.list[[i]]) DFvar[i, maxKSize+1] - paste(Objl,collapse=" ") } DFvar[is.na(DFvar)] - " " DFout - cbind(DFval, DFvar) # Вывод результатов print ("Конъюнкции класса 1") ; DFout[DFout$Nclass==1,] print ("Конъюнкции класса 2") ; DFout[DFout$Nclass==2,]

–  –  –

Дальнейшая работа с выделенными логическими высказываниями сводится к следующему (реализация алгоритма в R находится в стадии доработки).

Из генерируемого списка необходимо исключить подчиненные (или дочерние) конъюнкции, включающие более короткие претенденты:

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

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

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

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

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

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

Алгоритм "Кора", как и другие логические методы распознавания образов, является достаточно трудоемким, поскольку при отборе конъюнкций необходим полный или частично направленный перебор. Здесь может быть полезным использование процедур эволюционного поиска: генетический алгоритм (genetic algorithm), случайный поиск с адаптацией (СПА – см. книгу Г.С. Лбова, 1981 на сайте math.nsc.ru) и др.

5.4. Алгоритмы выделения ассоциативных правил Ассоциативные правила представляют собой механизм нахождения логических закономерностей между связанными элементами (событиями или объектами). Пусть имеется A = {a1, a2, a3, …, an} – конечное множество уникальных элементов (list of items). Из этих компонентов может быть составлено множество наборов T (sets of items), т.е. T A.

Ассоциативные правила A C имеют следующий вид:

если условие то результат, где, в отличие от деревьев классификации, условие – не логическое выражение, а набор объектов из множества А, с которыми связаны (ассоциированы) объекты того же множества, включенные в результат данного правила. Например, ассоциативное правило если (смородина, тля) то (муравьи) означает, что если на кусте смородины встретилась тля, то ищи поблизости и муравьев.

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

Например, при прогнозировании погоды одно из ассоциативных правил может выглядеть так:

направление ветра = NNW завтра будет дождь = TRUE

Выделяют три вида правил:

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

тривиальные правила, содержащие действительную и легко объяснимую информацию, отражающую известные законы в исследуемой области, и поэтому не приносящие какой-либо пользы;

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

Поиск ассоциативных правил обычно выполняют в два этапа:

в пуле имеющихся признаков A находят наиболее часто встречающиеся комбинации элементов T;

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

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

Правило A C имеет поддержку s, если оно справедливо для s% взятых в анализ случаев:

support(A C) = P(A C).

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

из A следует C ):

confidence(A C) = P(A C)/P(A) = support(A C)/support(A).

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

В пакете arules для R используются и другие показатели – подъемная сила, или лифт (lift), которая показывает, насколько повышается вероятность нахождения C в анализируемом случае, если в нем уже имеется A:

lift(A C) = confidence(A C)/support(C), и усиление (leverage), которое отражает, насколько интересной может быть более высокая частота A и C в сочетании с более низким подъемом leverage(A C) = support(A C) - support(A) * support(C).

Первый алгоритм поиска ассоциативных правил был разработан в 1993 г. сотрудниками исследовательского центра IBM, что сразу возбудило интерес к этому направлению. Каждый год появлялось несколько новых алгоритмов (DHP, Partition, DIC и др.), из которых наиболее известным остался алгоритм "Apriori" (Agrawal, Srikant, 1994).

Пакет arules позволяет находить часто встречающиеся сочетания элементов в данных (frequent itemsets) и отбирать ассоциативные правила, обеспечивая интерфейс к модулям на языке C, которые реализуют алгоритмы "Apriori" и "Eclat". Так как обычно обрабатываются большие множества наборов и правил, то для уменьшения объёмов требуемой памяти пакет содержит развитый инструментарий преобразования разреженных входных матриц в компактные наборы транзакций (Hahsler et al., 2016; Огнева, 2012).

Для работы с алгоритмами выделения ассоциаций в arules реализованы специальные типы данных, относящиеся к объектам трех классов: входной массив транзакций (transactions) и на выходе – часто встречающиеся фрагменты данных (itemsets) и правила (rules).

Объекты класса transactions представляют собой специально организованные бинарные матрицы со строками-наборами и столбцамипризнаками, содержащие значения элемента 1, если соответствующий признак есть в транзакции, и 0, если он отсутствует. В зависимости от типа данных и способа их загрузки, эти объекты могут иметь разные способы организации и состав дополнительных слотов. В частности, подкласс itemMatrix является одновременно средством представления разреженных матриц с использованием функционала пакета Matrix. Другим способом формирования экземпляров класса transactions является загрузка данных из файла функцией read.transactions().

Для выделения ассоциативных правил вновь обратимся к нашему примеру по классификации лиц избирателей (см. рис. 5.1). Метод "basket" функции read.transactions() предполагает, что каждая строка в файле представляет собой одну транзакцию, в которой признаки (т.е. их метки) разделены символами sep (по умолчанию – запятая).

Переформируем исходную таблицу данных в файл необходимого формата:

DFace - read.delim(file="Faces.txt", header=TRUE, row.names=1) Class - DFace$Class; DFaceN -DFace[,-17] colnames(DFaceN) - c("голова_круглая","уши_оттопырен", "нос_круглый","глаза_круглые","лоб_морщины", "носогубн_складка","губы_толстые","волосы","усы","борода", "очки","родинка_щеке","бабочка", "брови_подняты","серьга", "курит_трубка") Class[Class==1] - "Патриот" Class[Class==2] - "Демократ" items_list - sapply(1:nrow(DFaceN),function(i) paste(c(Class[i], colnames(DFaceN[i,DFaceN[i,]==1])), collapse=",", sep="\n")) head(items_list) write(items_list, file = "face_basket.txt") Патриот,уши_оттопырен,лоб_морщины,носогубн_складка,усы,борода,очки,бабочка,брови_подняты,курит_трубка Патриот,голова_круглая,нос_круглый,глаза_круглые,губы_толстые,волосы,борода,очки,родинка_щеке,серьга Патриот,глаза_круглые,лоб_морщины,носогубн_складка,волосы,усы,очки,родинка_щеке,бабочка,курит_трубка Патриот,уши_оттопырен,нос_круглый,носогубн_складка,губы_толстые,борода,очки,брови_подняты,серьга,курит_трубка Патриот,голова_круглая,уши_оттопырен,глаза_круглые,носогубн_складка,волосы,борода,родинка_щеке,брови_подняты,серьга Патриот,нос_круглый,лоб_морщины,носогубн_складка,губы_толстые,усы,очки,бабочка,серьга,курит_трубка Загрузим данные из файла и создадим объект itemMatrix.

Информацию о сформированном массиве транзакций можно получить, выполнив команды inspect() и summary().

library("arules") trans = read.transactions("face_basket.txt", format = "basket", sep=",") inspect(trans) # Выводимые данные не показаны summary(trans) # Выводимые данные показаны частично transactions as itemMatrix in sparse format with 16 rows (elements/itemsets/transactions) and 18 columns (items) and a density of 0.5555556 Для объектов класса itemMatrix можно осуществить быструю визуализацию данных или построить график распределения частот встречаемости признаков (рис. 5.4а, б).

image(trans) itemFrequencyPlot(trans, support = 0.1, cex.names=0.8)

–  –  –

Рис. 5.4б. Частотное распределение встречаемости признаков в транзакциях Поиск ассоциативных правил является не вполне тривиальной задачей, т.к. с ростом числа элементов в А экспоненциально растет число их потенциальных комбинаций. Алгоритм "Apriori" является итерационным: при этом сначала выполняются действия для одноэлементных наборов, затем для 2-х, 3-х элементных и т.д. (т.е. он во многом напоминает алгоритм "Кора").

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

сколько раз набор встречается в имеющемся наборе данных. При последующем поиске k-элементных наборов генерация претендентов состоит из двух фаз – формирование кандидатов нового уровня на основе (k - 1)элементных наборов, которые были определены на предыдущей итерации алгоритма, и удаление избыточных кандидатов. После того как найдены все часто встречающиеся наборы элементов, выполняют процедуру непосредственного извлечения правил из построенного хеш-дерева (Зайцев, 2009).

Результатом анализа транзакций с помощью пакета arules являются объекты класса associations, включающие описания множества отношений между признаками (в виде часто встречающихся фрагментов, или правил), которые отбираются в соответствии с различными перечисленными выше мерами качества. Подкласс rules состоит из двух объектов itemMatrix, представляющих левую lhs (left-hand-side) и правую rhs (right-hand-side) сторону правила A C, т.е. A – lhs, C – rhs.

Формирование правил осуществляется функцией apriori() с указанием пороговых значений поддержки и достоверности.

Функция summary() обеспечивает частотный анализ правил по их длине и достигнутым мерам качества:

rules - apriori(trans, parameter = list(support = 0.01, confidence = 0.6)) summary(rules) # Результаты показаны частично set of 48306 rules rule length distribution (lhs + rhs):sizes

summary of quality measures:

support confidence lift Min. :0.06250 Min. :0.6000 Min. :0.960 1st Qu.:0.06250 1st Qu.:1.0000 1st Qu.:1.600 Median :0.06250 Median :1.0000 Median :1.600 Mean :0.07788 Mean :0.9781 Mean :1.750 3rd Qu.:0.06250 3rd Qu.:1.0000 3rd Qu.:1.600 Max. :0.62500 Max. :1.0000 Max. :2.667 Всего было отобрано 48306 правил, длина которых большей частью составляла от 5 до 7 элементов. Функция plot() из пакета arulesViz позволяет получать различные формы визуализации синтезированных правил, в том числе, анализ изменчивости их мер качества (рис.

5.5):

library("arulesViz") plot(rules, measure=c("support","lift"), shading="confidence")

Рис. 5.5. Поддержка, лифт и достоверность сгенерированных правил

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

5.6):

rulesPat - subset(rules, subset = rhs %in% "Патриот" & lift 1.8) inspect(head(rulesPat, n = 10, by = "support")) plot(head(sort(rulesPat, by="support"), 10), method="paracoord")

–  –  –

Рис. 5.7. Визуализация в форме графа 10 лучших правил для демократов Метод "graph" функции plot() показывает правила и составляющие их признаки в виде графа, размер узлов которого пропорционален уровню поддержки каждого представленного правила (рис. 5.7).

5.5. Анализ последовательностей знаков или событий Проблема статистического анализа последовательностей возникла, в первую очередь, с появлением быстрых методов секвенирования ДНК, в результате чего объем баз первичных нуклеотидных последовательностей растет экспоненциально. Параллельно интерес к этим методам возник при проведении маркетинговых исследований. Появилось желание выяснить, какова спонтанная траектория передвижения покупателя в торговом центре или в каком порядке пользователь нажимает кнопки и ссылки в Интернетмагазине, когда покупает товар (Горбач, Цейтлин, 2011).

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

–  –  –

Функция seqdef() обеспечит нам создание объекта stslist, включающего последовательности с необходимыми метками.

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

library(TraMineR) ; data(biofam) summary(biofam) # Результаты не приведены biofam.lab - c("Родит", "Отд", "Сем+Род", "Сем.", "Один+Род", "Один+Отд", "Сем+Дет", "Развод") biofam.seq - seqdef(biofam[, 10:25], labels = biofam.lab) par(mfrow = c(2, 2)); seqiplot(biofam.seq, withlegend = FALSE) seqdplot(biofam.seq, withlegend = FALSE) seqfplot(biofam.seq, withlegend = FALSE) seqlegend(biofam.seq) Рис. 5.8. Визуализация 10 первых (а), наиболее часто встречающихся последовательностей (б) и график частотного распределения градаций (в) В этом наборе данных содержится также ряд дополнительных показателей об анкетируемых респондентах: пол (sex), год рождения от 1909 до 1957 г. (birthyr), национальность (nat_1_02), родной язык (plingu02), религиозность (p02r01, p02r04), профессиональный статус отца и матери (cspfaj, cspmoj). С использованием семейства функций seq*plot() можно визуализировать данные с разбивкой по группам. Например, можно сделать группировку по полу и установить, что только 20% 25-летних швейцарских девушек продолжают жить с родителями, тогда как таких юношей – более 33% (рис.

5.9):

seqdplot(biofam.seq,group=biofam$sex, withlegend = FALSE)

Рис. 5.9. График частотного распределения градаций с группировкой по полу

Для набора исходных данных можно оценить два вида энтропии Шеннона. Вертикальная, или "кросс-секторальная", энтропия рассчитывается на основе распределения состояний в каждый момент времени для всех анализируемых последовательностей (рис.

5.10):

seqHtplot(biofam.seq)

Рис. 5.10. Изменение относительной энтропии по градациям возрастов

Очевидно, что в 15 лет все респонденты находятся в одном и том же состоянии – проживают в доме родителей и энтропия равна 0. Возраст от 24 до 27 лет связан с активной сменой семейного статуса.

Для каждой из последовательностей можно рассчитать набор разнообразных метрик, таких как длина, число подпоследовательностей (number of subsequences) или различных состояний (distinct states), длина пребывания в стабильном состоянии (durations) и число внутренних переходов (transitions). Важной обобщенной мерой нестабильности последовательности является ее горизонтальная энтропия (longitudinal, или within-sequence entropy). В нормированном выражении она принимает значение 0, если последовательность состоит из одинаковых состояний, и 1, если каждое состояние из общего словаря встречается с равной вероятностью (т.е. напимер, каждый из 8 кодов семейного состояния встречается точно два раза).

Рассмотрим, как изменяется средняя энтропия последовательностей за исторический период середины ХХ века (рис.

5.11):

Entropy - seqient(biofam.seq) ageg = cut(biofam$birthy, c(1909,1918,1928,1938,1948,1958), label = c("1909-18","1919-28","1929-38","1939-48","1949-58"), include.lowest = TRUE) boxplot(Entropy ~ ageg, data = biofam, xlab = "Диапазон годов рождения", ylab = "Энтропия последовательностей", col = "cyan") Рис. 5.11. Изменение относительной энтропии семейного состояния в течение ХХ века Дотошный аналитик наверняка сможет усмотреть на рис. 5.11, что молодые швейцарцы 70-80-х годов существенно активнее меняли свой семейный статус, чем в 30-40 годы.

Индекс сложности (complexity) последовательности рассчитывается как геометрическое среднее энтропии и числа внутренних переходов, отнесенных к их максимальным значениям. По еще более сложным формулам, учитывающим длины фрагментов с постоянным состоянием и разнообразие переходов, рассчитывается мера турбулентности (turbulence). Достоинства и недостатки отдельных индексов постоянно дискутируются (Gabadinho et al., 2011а), поэтому интересно проследить зависимость между ними и убедиться в том, что каждый индекс характеризует последовательность несколько поразному (рис.

5.12):

Turbulence - seqST(biofam.seq) plot(Turbulence, Entropy, xlab = "Турбулентность", ylab = "Энтропия") Рис. 5.12.

Зависимость энтропии от турбулентности Поскольку любой из этих индексов является численной характеристикой последовательности, они могут использоваться как компоненты обычных линейных моделей:

m.turb - lm(Turbulence ~ sex + birthyr, data = biofam) summary(m.turb)

Coefficients:

Estimate Std. Error t value Pr(|t|) (Intercept) -50.534580 7.355288 -6.871 8.52e-12 *** sexwoman 0.373527 0.079528 4.697 2.82e-06 *** birthyr 0.028381 0.003786 7.495 9.88e-14 ***

--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Нетрудно сделать вывод, что турбулентность семейной жизни статистически значимо выше у женщин и увеличивается по мере приближения к ХХI веку.

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

1. Вычисление расстояния по Хеммингу, или числа несовпадающих позиций (nomatching positions).

2. Выделение длины совпадающей части с начала последовательности Sp (LCP, Longest Common Prefix, Dp = 1 - Sp) или максимальной длины общей подпоследовательности (LCS, Longest Common Subsequence).

3. Оценка метрики Левенштайна (Горбач, Цейтлин, 2011), т.е. суммарной "стоимости" преобразования одной последовательности в другую.

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

В пакете TraMineR функция seqsubm() с использованием method = "TRATE" выполняет оценку матрицы стоимостей переходов от а к б, а другая функция - seqdist() - с использованием этих коэффициентов рассчитывает матрицу дистанций "оптимального" редактирования (OM, optimal matching distances) между каждой парой взятых в анализ последовательностей:

couts - seqsubm(biofam.seq, method = "TRATE") biofam.om - seqdist(biofam.seq, method = "OM", indel = 3, sm = couts) round(biofam.om[1:10, 1:10], 1) [] 2000 sequences with 8 distinct events/states [] 537 distinct sequences [] min/max sequence length: 16/16 [] computing distances using OM metric [] total time: 0.99 secs # --------------------Мы привели таблицу парных расстояний между первыми 10 последовательностями, и читатель может сопоставить ее с результатом их визуализации на рис. 5.8а. Отметим, что для расчета матрицы размерностью 2000 x 2000 и объемом 30 Мб на стандартном компьютере понадобилось меньше 1 сек.

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

Используем функцию agnes(), на вход которой подается рассчитанная выше матрица расстояний:

library(cluster) clusterward - agnes(biofam.om, diss = TRUE, method = "ward") plot(clusterward, which.plots = 2) Рис. 5.13.

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

# "Распилим" дерево на три части cluster3 - cutree(clusterward, k = 3) cluster3 - factor(cluster3, labels = c("Кластер 1", "Кластер 2", "Кластер 3")) table(cluster3) par(mfrow = c(2,2)) seqmtplot(biofam.seq, group = cluster3) cluster3 Кластер 1 Кластер 2 Кластер 3 Рис. 5.14. Распределение частот семейных состояний по выделенным кластерам Матрица расстояний между последовательностями может быть использована для непараметрического дисперсионного анализа по методу Андерсона (Anderson, 2001), подробному описанию сущности которого посвящен раздел 9.3. В общих чертах, если группировка задается в соответствии с уровнями изучаемых факторов, то алгоритм осуществляет разложение многомерной дисперсии, заключенной в матрице расстояний, на общую SST и внутригрупповую SSW суммы квадратов расстояний dij. Пакет включает несколько эффективно работающих функций TraMineR непараметрического ANOVA, которые с успехом могут использоваться в любом анализе произвольных матриц дистанций.

Функция dissassoc() оценивает зависимость между степенью изменчивости элементов матрицы парных расстояний и дискретной ковариатой переменной). Рассчитываются обычные (номинальной компоненты таблицы дисперсионного анализа (суммы квадратов и их средние по числу степеней свободы), дисперсионное отношение F и коэффициент детерминации R2. Функция выполняет также тест на однородность дисперсий в выделенных группах с использованием статистик Левене (Levene) и Бартлетта (Bartlett). Статистическая значимость тестов оценивается путем рандомизации (т.е. элементы матрицы дистанций случайным образом перемешиваются R раз относительно уровней фактора).

Оценим, как влияет пол респондента (sex) на взаимные расстояния между последовательностями:

dissvar(biofam.om) # Общая дисперсия парных расстояний da - dissassoc(biofam.om, group=biofam$sex, R = 1000) print(da) [1] 8.316977

–  –  –

6.1. Дискриминантный анализ Линейный дискриминантный анализ (Linear Discriminant Analysis, LDA) является разделом многомерного анализа, который позволяет оценивать различия между двумя и более группами объектов по нескольким переменным одновременно (Афифи, Эйзенс, 1982; Айвазян и др., 1989). Он реализует две тесно связанные между собой статистические процедуры:

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

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

В основе дискриминантного анализа лежит предположение о том, что описания объектов каждого k-го класса представляют собой реализации многомерной случайной величины, распределенной по нормальному закону Nm(µk; k) со средними µk и ковариационной матрицей 1 nk Ck = (xik µ k )T (xik µ k ) nk 1 i =1 (индекс m указывает на размерность признакового пространства).

Рассмотрим несколько упрощенную геометрическую интерпретацию алгоритма LDA для случая двух классов. Пусть дискриминантные переменные x – оси m-мерного евклидова пространства. Каждый объект (наблюдение) является точкой этого пространства с координатами, представляющими собой фиксируемые значения каждой переменной. Если оба класса отличаются друг от друга по наблюдаемым пртзнакам, их можно представить как скопления точек в разных областях рассматриваемого пространства, которые могут частично перекрываться. Для определения положения каждого класса можно вычислить его "центроид", который является воображаемой точкой, координатами которой являются средние значения переменных в данном классе.

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

Ее положение задается линейной дискриминантной функцией (linear discriminant, LD) с весовыми коэффициентами j, определяющими вклад z (x) = 1x1+ 2x2+…+ mxm.

каждой исходной переменной xj:

Если сделать предположение, что ковариационные матрицы объектов классов 1 и 2 равны, т.е. C = C1 = C2, то вектор коэффициентов {1,…, m} линейного дискриминанта z (x) может быть вычислен по формуле = С-1(µ1 - µ2), где С-1 – матрица, обратная к ковариационной, µк – вектор средних k-го класса. Полученная ось совпадает с уравнением прямой, проходящей через центроиды двух групп объектов классов, а обобщенное расстояние Махаланобиса, равное дистанции между ними в многомерном пространстве признаков, оценивается как D2 = T (µ1 - µ2).

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

Если между ними нет серьезных отличий, их объединяют в расчетную ковариационную матрицу C = [C1(n1 – 1) + C2(n2 – 1)]/ (n1 + n2 – 2).

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

Для проверки гипотезы о многомерном нормальном распределении данных используется многомерная версия критерия согласия Шапиро-Уилка, которая реализована в функции mshapiro.test() из пакета mvnormtest. На вход этой функции подается матрица, строки которой соответствуют переменным, а столбцы – наблюдениям.

В разделах 2.4-2.5 мы подробно рассматривали пример анализа зависимости между двумя различными способами производства листового стекла (флеш-стекло и по принципу вертикального вытягивания) и составом его химических ингредиентов. С использованием статистики Пиллая и критерия Хотеллинга была показана статистическая значимость такой связи.

Применим многомерный критерий Шапиро-Уилка к оценке характера распределения этой выборки:

DGlass - read.table(file="Glass.txt", sep=",", header=TRUE, row.names=1) DGlass$F - as.factor(ifelse(DGlass$Class==2, 2, 1)) library(mvnormtest) mshapiro.test(t(DGlass[DGlass$F == 1, 1:9])) mshapiro.test(t(DGlass[DGlass$F == 2, 1:9])) Shapiro-Wilk normality test data: F == 1 W = 0.1965, p-value 2.2e-16 data: F == 2 W = 0.1368, p-value 2.2e-16 Для обеих групп выявлены значимые отклонения от многомерного нормального распределения.

Для проверки гипотезы о гомогенности матриц ковариаций используется так называемый M-критерий Бокса, который реализован в функции boxM() из пакета biotools:

library(biotools) boxM(as.matrix(DGlass[, 1:9]), DGlass$F ) Box's M-test for Homogeneity of Covariance Matrices data: as.matrix(DGlass[, 1:9]) Chi-Sq (approx.) = 443.1846, df = 45, p-value 2.2e-16 Гетерогенность матриц ковариаций также статистически значима.

Дискриминантный анализ реализован в нескольких пакетах для R, но мы рассмотрим применение функции lda() из базового пакета MASS.

Поскольку важной характеристикой прогнозирующей эффективности модели является ее ошибка при перекрестной проверке, то в функции lda() пакета MASS заложена реализация скользящего контроля (leave-one-out CV).

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

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

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

# Функция вывода результатов классификации Out_CTab - function (model, group, type="lda") { # Таблица неточностей "Факт/Прогноз" по обучающей выборке classified-predict(model)$class t1 - table(group,classified) # Точность классификации и расстояние Махалонобиса Err_S - mean(group != classified) ; mahDist - NA if (type=="lda") { mahDist - dist(model$means %*% model$scaling) } # Таблица "Факт/Прогноз" и ошибка при скользящем контроле t2 - table(group, update(model, CV=T)$class-LDA.cv) Err_CV - mean(group != LDA.cv) Err_S.MahD -c(Err_S, mahDist) Err_CV.N -c(Err_CV, length(group)) cbind(t1, Err_S.MahD, t2, Err_CV.N) } # --- Выполнение расчетов library(MASS) lda.all - lda(F ~., data=DGlass[,-10]) Out_CTab(lda.all, DGlass$F)

1 2 Err_S.MahD 1 2 Err_CV.N 1 69 18 0.2515337 66 21 0.3128834 2 23 53 1.2414682 30 46 163

Отметим существенный рост ошибки распознавания до 31% при выполнении скользящего контроля.

Естественно задаться вопросом, какие из имеющихся 9 признаков являются информативными при разделении, а какие – сопутствующим балластом.

Шаговая процедура выбора переменных при классификации, реализованная функцией stepclass() из пакета klaR, основана на вычислении сразу четырех параметров качества моделей-претендентов:

а) индекса ошибок (correctness rate), б) точности (accuracy), основанной на евклидовых расстояниях между векторами "факта" и "прогноза",

в) способности к разделимости (ability to seperate), также основанной на расстояниях, и г) доверительных интервалах центроидов классов. Все эти параметры оцениваются в режиме многократной перекрестной проверки.

library ( klaR) stepclass(F ~., data=DGlass[,-10], method="lda") `stepwise classification', using 10-fold cross-validated correctness rate: 0.69228; in: "Al"; variables (1): Al correctness rate: 0.76618; in: "Ca"; variables (2): Al, Ca final model : F ~ Al + Ca lda.step - lda(F ~ Mg + Al, data=DGlass[,-10])

Prior probabilities of groups:

0.5337423 0.4662577

Group means:

Mg Al 1 3.550690 1.

171149 2 3.002105 1.408158

Coefficients of linear discriminants:

LD1 Mg -0.8294482 Al 2.6894108 В результате получили компактную дискриминантную функцию z (x) = 2.69Al - 0.83Mg, зависящую только от двух переменных.

Найдем ошибку предсказания как на обучающей выборке, так и при скользящем контроле:

partimat(F ~ Mg + Al, data=DGlass[,-10], main='', method="lda") Out_CTab(lda.step, DGlass$F) 1 2 Err_S.MahD 1 2 Err_CV.N 1 72 15 0.2331288 72 15 0.2392638 2 23 53 1.0924355 24 52 163 По обоим критериям ошибка предсказания сокращенной модели существенно ниже аналогичных показателей при использовании полного набора переменных.

Это отражает общие методические закономерности, важные при выборе алгоритмов классификации и анализе результатов их работы:

ошибка перекрестной проверки ECV всегда превышает внутреннюю ошибку модели ES на самой обучающей выборке и является объективной характеристикой качества распознавания на внешнем дополнении;

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

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

Рис. 6.1. Диаграмма дискриминации двух классов 1 и 2 (ошибочное распознавание выделено шрифтом красного цвета); показаны линейная дискриминантная функция z и жирными точками – положение центроидов.

–  –  –

The top 3 variables (out of 3): Al, K, Fe Для того чтобы уточнить, какую из трех построенных моделей следует предпочесть, выполним их тестирование на основе 10-кратной перекрестной проверки с 5 повторностями:

# Модель на основе всех 9 предикторов lda.full.pro - train(DGlass[, 1:9], DGlass$F, data= DGlass,method="lda", trControl = trainControl(method="repeatedcv",repeats=5, classProbs = TRUE), metric = "Accuracy") # Модель на основе 2 предикторов stepclass lda.step.pro - train(F1 ~ Mg + Al, data= DGlass,method="lda", trControl = trainControl(method="repeatedcv",repeats=5, classProbs = TRUE), metric = "Accuracy") # Модель на основе 3 предикторов rfe lda.rfe.pro - train(F1 ~ Al + K + Fe, data= DGlass,method="lda", trControl = trainControl(method="repeatedcv",repeats=5, classProbs = TRUE), metric = "Accuracy") plot(varImp(lda.full.pro)) Linear Discriminant Analysis 9 predictor - полная модель Accuracy Kappa 0.6974559 0.3875972 2 predictor - модель stepclass() Accuracy Kappa 0.7613382 0.5175581 3 predictor - модель rfe() Accuracy Kappa 0.7253725 0.4456987

–  –  –

Функция rfe() провела отбор переменных в полном соответствии с рейтингом их важности на рис. 6.2, но это решение оказалось неоптимальным. Модель lda.step, полученная с использованием функции stepclass(), оказалась существенно эффективней.

6.2. Метод опорных векторов Метод опорных векторов (support vector), называемый ранее алгоритмом "обобщенного портрета", был разработан советскими математиками В. Н. Вапником и А. Я. Червоненкисом (1974) и с тех пор приобрел широкую популярность.

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

Если имеется два класса наблюдений и предполагается линейная форма границы между классами, то возможны два случая. Первый из них связан с возможностью идеального разделения данных при помощи некоторой p гиперплоскости zk ( x) = i =1i xi + 0 (двумерный вариант представлен слева на рис. 6.3). Поскольку таких гиперплоскостей может быть множество, то оптимальной является разделяющая поверхность, которая максимально удалена от обучающих точек, т.е. имеющая максимальный зазор M (margin).

–  –  –

На рис. 6.3 справа показан другой случай, когда облако точек перекрывается и оба класса линейно неразделимы. Собственно опорными векторами называются наблюдения, лежащие непосредственно на границе разделяющей полосы, либо на неправильной для своего класса стороне относительно границ зазора (такие точки отмечены *j на рис. 6.3). Для граничных и всех остальных точек принимается *j = 0.

Оптимальную разделяющую гиперплоскость такого классификатора p zk ( x) = i =1i xi + 0 также находят из условия максимизации ширины зазора M, но при этом разрешается неверно классифицировать некоторую небольшую группу наблюдений, относящихся к опорным векторам. Для этого задаются дополнительным условием оптимизации *j С, где С – допустимое число j нарушений границы зазора и их выраженность, которое обычно выбирается с использованием перекрестной проверки. Математически поиск решения сводится к удобной задаче квадратичной оптимизации с линейными ограничениями, которая гарантировано сходится к одному глобальному минимуму (Vapnik, 1995; Джеймс и др., 2016).

Поскольку на расположение гиперплоскости оказывают влияние только те наблюдения, которые лежат на границах зазора или нарушают его, то решающее правило такого классификатора довольно устойчиво к выбросам большинства точек, расположенных вне "критической зоны" разделения. Это свойство отличает его от свойств других классификаторов. Например, разделяющая плоскость LDA на рис. 6.3 проводится перпендикулярно дискриминационной функции z и зависит от средних и ковариационной матрицы, вычисляемых по всему множеству имеющихся наблюдений.

Для реализации метода опорных векторов в среде R обычно используется функция svm() из пакета e1071. Рассмотрим построение с помощью этой функции линейного классификатора для прогнозирования способа производства стекла по его химическому составу (см.

разделы 2.4Для подбора параметров модели зададим перекрестную проверку с делением исходной выборки на 10 равных частей (cross=10):

DGlass - read.table(file="Glass.txt", sep=",", header=TRUE, row.names=1) DGlass$F - as.factor(ifelse(DGlass$Class==2,2,1)) svm.all - svm(formula = F ~., data=DGlass[,-10], cross=10, kernel = "linear") table(Факт=DGlass$F, Прогноз=predict(svm.all)) Acc = mean(predict(svm.all) == DGlass$F ) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт 1 2 [1] "Точность=74.23%" Для подтверждения того, что svm() действительно проводила построение модели с применением кросс-проверки, выполним дополнительно скользящий контроль с использованием самостоятельно составленной функции:

# Функция для вычисления ошибки перекрестной проверки library(e1071) CVsvm - function(x, y) { n - nrow(x) ; Err_S - 0 for(i in 1:n) { svm.temp - svm(x=x[-i,], y = y[-i], kernel = "linear") if (predict(svm.temp, newdata=x[i,]) != y[i]) Err_S - Err_S +1 } Err_S/n } Acc - 1-CVsvm(DGlass[,1:9],DGlass$F) paste("Точность=", round(100*Acc, 2), "%", sep="") [1] "Точность=74.17%" Ошибка предсказания изменилась лишь незначительно.

Определенным резервом повышения эффективности модели является выбор оптимального набора предикторов.

Селекция переменных может быть выполнена несколькими способами:

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

с использованием алгоритма рекурсивного извлечения переменных SVMRFE (Recursive Feature Extraction Algorithm; Guyon et al., 2002);

на основе функций penalizedSVM или хорошо знакомого нам пакета caret.

Для реализации функций stepclass() и svmlight() из пакета klaR необходимо скачать (http://svmlight.joachims.org) и поместить в рабочую директорию бинарные исполняемые файлы метода SVMlight (Joachims, 1999), написанные на языке C (эти файлы включены в приложения к этой книге):

stepclass(F ~., data = DGlass[, -10], method = "svmlight", pathsvm = "D:/R/SVMlight") Однако выполнение этой команды оказалось в наших условиях столь продолжительным, что мы не смогли продемонстрировать читателям эту возможность.

Рекурсивное извлечение переменных SVM-RFE происходит с использованием специального критерия, оценивающего релевантность каждого предиктора по отношению к качеству полученного классификатора.

Авторы метода (Guyon et al., 2002) написали и предоставили в открытый доступ скрипт небольшой функции (http://citeseer.ist.psu.edu) возвращающей ранжированный список svmrfeFeatureRanking(), переменных.

Мы также включили файл с командами этой функции в наш набор скриптов, сопутствующих книге:

source("SVM-RFE.r") featureRankedList = svmrfeFeatureRanking(DGlass[,1:9], DGlass$F) ErrSvm - sapply(1:9, function(nf) { svmModel = svm(DGlass[,featureRankedList[1:nf]], DGlass$F, kernel="linear") mean( predict(svmModel) != DGlass$F ) } ) data.frame(Index=featureRankedList, NameFeat = names(DGlass[,featureRankedList]), ErrSvm = ErrSvm) Index NameFeat ErrSvm 1 4 Al 0.2699387 2 6 K 0.2883436 3 9 Fe 0.2515337 4 8 Ba 0.2576687 5 3 Mg 0.2515337 6 5 Si 0.2576687 7 7 Ca 0.2331288 8 2 Na 0.2331288 9 1 RI 0.2576687 Отметим, что метод опорных векторов слабо зависит от эффектов коллинеарности предикторов и, по сравнению с LDA, значительно меньше переменных было исключено из модели как малоинформативные.

Построим гиперплоскость на основе 7 отобранных признаков:

DGlass.sel - DGlass[, c(featureRankedList[1:7], 11)] (svm.sel - svm(formula = F ~., data=DGlass.sel, cross=10, kernel = "linear", prob=TRUE)) table(Факт=DGlass.sel$F, Прогноз=predict(svm.sel)) Acc = mean(predict(svm.all) == DGlass$F ) paste("Точность=", round(100*Acc, 2), "%", sep="") SVM-Kernel: linear cost: 1 gamma: 0.1428571 Number of Support Vectors: 109 Прогноз Факт 1 2 [1] "Точность=76.69%" Действительно, число правильно классифицированных образцов оказалось на 4 единицы больше, чем при полном наборе признаков.

Рассмотрим процедуры селекции признаков и тестирования моделей с использованием пакета caret.

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

DGlass$F1 - as.factor(ifelse(DGlass$Class==2,"No","Flash")) svmProfile - rfe(DGlass[,1:9], DGlass$F1, sizes = 2:9, rfeControl = rfeControl(functions = caretFuncs, method = "cv"), method = "svmLinear") Recursive feature selection Variables Accuracy Kappa AccuracySD KappaSD Selected 2 0.5702 0.08002 0.03966 0.09278 3 0.6180 0.19280 0.07887 0.16994 4 0.7108 0.41206 0.11048 0.21729 5 0.7285 0.44707 0.12608 0.25065 * 6 0.6991 0.38878 0.09543 0.18624 7 0.6991 0.39014 0.10365 0.20313 8 0.7175 0.42935 0.11787 0.23245 9 0.6670 0.32381 0.11461 0.22449

The top 5 variables (out of 5):

Al, K, Fe, Si, Ba Информативный набор предикторов, полученный методом рекурсивного исключения, состоит из 5 признаков.

Выполним тестирование моделей на основе разного числа предикторов (9, 7 и 5) с использованием показателя AUC (площади под ROC-кривой – см.

раздел 2.4), для чего в trainControl() зададим параметр summaryFunction = twoClassSummary:

ctrl - trainControl(method="repeatedcv",repeats=5, summaryFunction=twoClassSummary, classProbs=TRUE) print(" Модель на основе всех 9 предикторов ") train(DGlass[,1:9], DGlass$F1, method = "svmLinear", metric="ROC", trControl=ctrl) print(" Модель на основе 7 предикторов ") train(DGlass[,c(4,6,9,8,3,5,7)], DGlass$F1, method = "svmLinear", metric="ROC", trControl=ctrl) print(" Модель на основе 5 предикторов ") train(DGlass[,c(4,6,9,8,5)], DGlass$F1, method = "svmLinear", metric="ROC", trControl=ctrl) Модель на основе всех 9 предикторов 163 samples 9 predictor 2 classes: 'Flash', 'No' ROC Sens Spec 0.7564236 0.7730556 0.5703571

–  –  –

Модель на основе 5 предикторов ROC Sens Spec 0.6986111 0.8169444 0.5396429 Последняя модель на основе RFS-метода оказалась наименее качественной по AUC, но имеет наибольшую чувствительностью Sens.

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

Машину опорных векторов SVM (Support Vector Machine) можно рассматривать как нелинейное обобщение линейного классификатора, представленного в разделе 6.2, основанное на расширении размерности исходного пространства предикторов с помощью специальных ядерных функций. Это позволяет строить модели с использованием разделяющих поверхностей самой различной формы.

Рис. 6.4. Линейная и радиальная границы классов

Собственно ядром является любая симметричная, положительно полуопределенная матрица K, которая составлена из скалярных произведений пар векторов xi и xj: K(xi, xj) = (xi), (xj), характеризующих меру их близости. Здесь – произвольная преобразующая функция, формирующая ядро.

В качестве таких функций чаще всего используют следующие:

линейное ядро: K(xi, xj) = xiт xj, что соответствует классификатору на опорных векторах в исходном пространстве (см. рис. 6.3);

полиномиальное ядро со степенью р: K(xi, xj) = (1 + xiт xj)р;

гауссово ядро с радиальной базовой функцией (RBF):

K(xi, xj) = exp{||xi - xj||2};

сигмоидное ядро: K(xi, xj) = tanh(xiт xj + 0).

Каждое ядро характеризуется параметрами (p,, 0 и т.д.), которые подлежат оптимизации.

Основная идея использования ядер заключается в том, что при отображении данных в пространство более высокой размерности исходное множество точек может стать линейно разделимым (Cristianini, Shawe-Taylor, 2000). Тогда естественно предположить, что оптимальная разделяющая гиперплоскость машины опорных векторов может быть найдена путем p подбора коэффициентов i выражения z k ( x) = i =1 i K ( x, xi ) + 0. Таким образом, как и в линейном случае, решается математически удобная задача квадратичной оптимизации с использованием множителей Лагранжа.

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

Построение классификатора на опорных векторах с использованием перечисленных выше ядер можно осуществить с помощью знакомой нам функции svm() из пакета e1071. Отметим, что по умолчанию в этой функции установлено значение параметра kernel = "radial", т.е. принято RBF-ядро.

Для подгонки модели распознавания типа стекла (разделы 6.1-6.2) нам необходимо предварительно оценить значения двух параметров: С (cost) – допустимый штраф за нарушение границы зазора и (gamma) – параметр радиальной функции. При увеличении С происходит уменьшение зазора и количества используемых опорных векторов: граница между классами становится более извилистой, а классификатор более точным на обучающей выборке и менее точным – на контрольной выборке.

В состав пакета входит функция tune.svm(), которая по информации о строящейся модели и с использованием перекрестной проверки осуществляет предварительный отбор необходимых параметров:

DGlass - read.table(file="Glass.txt", sep=",", header=TRUE,row.names=1) DGlass$F - as.factor(ifelse(DGlass$Class==2,2,1)) library(e1071) tune.svm(F ~., data=DGlass[,-10], gamma = 2^(-1:1), cost = 2^(2:4))

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation

- best parameters:

gamma cost 0.5 4

- best performance: 0.1886029 svm.rbf - svm(formula = F ~., data=DGlass[,-10], kernel = "radial",cost = 4, gamma = 0.5) table(Факт=DGlass$F, Прогноз=predict(svm.rbf)) Acc - mean(predict(svm.rbf) == DGlass$F ) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт 1 2 [1] "Точность=94.48%" Мы получили классификатор, значительно превышающий по точности рассмотренные выше модели линейного дискриминантного анализа и опорных векторов с линейным ядром. Полученная разделяющая поверхность имеет многомерный характер, но мы можем оценить ее внешний вид, создавая двумерные проекции любых двух исходных переменных из 9 (например, концентрации алюминия и магния):

svm2.rbf - svm(formula = F ~ Al+Mg, data=DGlass[,-10], kernel = "radial",cost = 4, gamma = 2.5) Acc - mean(predict(svm2.rbf) == DGlass$F ) paste("Точность=", round(100*Acc, 2), "%", sep="") [1] "Точность=85.89%" plot(svm2.rbf, DGlass[,c(3,4,11)], svSymbol = 16, dataSymbol = 2, color.palette = terrain.colors) legend("bottomleft", c("Опорные кл.1", "Опорные кл.2", "Данные кл.1", "Данные кл.2"), col=c(1,2,1,2), pch=c(16,16,2,2)) Рис. 6.5. Пример классификатора на опорных векторах с использованием в качестве ядра радиальной базисной функции Выполним аналогичные действия с использованием функции train() из пакета caret, где для построения модели опорных векторов используется функция svm() из пакета kernlab. Здесь также осуществляется подбор двух параметров – C и sigma (вероятно, то же самое, что и gamma).

Определим диапазоны их варьирования в таблице grid: C – от 2 до 5, а sigma – от 0.3 до 0.6:

ctrl - trainControl(method="repeatedcv",repeats=5, summaryFunction=twoClassSummary, classProbs=TRUE) grid - expand.grid(sigma = (3:6)/10, C = 2:5) (svmRad.tune - train(DGlass[,1:9], DGlass$F1, method = "svmRadial", metric="ROC", tuneGrid = grid, trControl=ctrl))

sigma C ROC Sens Spec 0.3 2 0.8808209 0.8211111 0.7957143 0.3 3 0.8832316 0.8300000 0.8110714 0.3 4 0.8780878 0.8372222 0.8053571 0.3 5 0.8761682 0.8300000 0.8078571 0.4 2 0.8864906 0.8300000 0.8242857 0.4 3 0.8766642 0.8297222 0.8189286 0.4 4 0.8692634 0.8322222 0.8164286 0.4 5 0.8648983 0.8344444 0.8192857 0.5 2 0.8797693 0.8300000 0.8221429 0.5 3 0.8730729 0.8247222 0.8189286 0.5 4 0.8693080 0.8180556 0.8192857 0.5 5 0.8623041 0.8158333 0.8085714 0.6 2 0.8826364 0.8230556 0.8192857 0.6 3 0.8721453 0.8063889 0.8192857 0.6 4 0.8663715 0.8066667 0.8167857 0.6 5 0.8638814 0.7950000 0.8142857

ROC was used to select optimal model using the largest value.

The final values used for model were sigma = 0.4 and C = 2.

pred - predict(svmRad.tune,DGlass[,1:9]) table(Факт=DGlass$F1, Прогноз=pred) Acc - mean(pred == DGlass$F1) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт Flash No Flash 83 4 No 9 67 [1] "Точность=92.02%" При рассмотрении приведенного примера (рис. 6.5) может создаться впечатление, что полученные нелинейные плоскости хорошо зарекомендовали себя при предсказании, но трудно интерпретируемы с точки зрения объяснения. Так и есть. Однако основная задача подобных моделей – научиться распознавать образы (pattern recognition), что нашло широкое применение в вычислительной биологии, генетике, спектроскопии, анализе данных, полученных с использованием микрочипов (microarray) и т.д.

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

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

library('kernlab') ; library('ggplot2') set.seed(2335246L) ; data('spirals') # Спектральная функция выделяет две различные спирали sc - specc(spirals, centers = 2) s - data.frame(x=spirals[,1],y=spirals[,2], class=as.factor(sc)) # Данные делятся на обучающую и проверочную выборки s$group - sample.int(100,size=dim(s)[[1]],replace=T) sTrain - subset(s,group10) sTest - subset(s,group=10) # Снова используем Гауссово или радиальное ядро mSVMG - ksvm(class~x+y,data=sTrain,kernel='rbfdot') sTest$predSVMG - predict(mSVMG,newdata=sTest,type='response') table(Факт=sTest$class, Прогноз=sTest$predSVMG) Прогноз Факт 1 2 ggplot() + geom_text(data=sTest, aes(x=x, y=y, label=predSVMG), size=12) + geom_text(data=s, aes(x=x, y=y, label=class, color=class), alpha=0.7) + coord_fixed() + theme_bw() + theme(legend.position='none') В результате использования радиальной ядерной функции была построена разделяющая поверхность между двумя вложенными спиралями таким образом, что ошибочное предсказание правильного класса наблюдалось лишь для 3 объектов проверочной выборки из 29.

–  –  –

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

Как упоминалось ранее, деревья решений CART могут эффективно применяться и для прогнозирования бинарного отклика. При построении дерева важно установить оптимальный уровень ветвления и выбрать значение относительного параметра стоимости сложности cp. Визуализацию дерева DGlass - read.table(file=»Glass.txt», sep=",", header=TRUE,row.names=1) DGlass$F1 - as.factor(ifelse(DGlass$Class==2,"No","Flash")) control-trainControl(method="repeatedcv",number=10,repeats=3) set.seed(7) fit.cart - train(DGlass[,1:9], DGlass$F1, method="rpart", trControl=control) library(rpart.plot) rpart.plot(fit.cart$finalModel) выполним с помощью пакета rpart.plot (рис.

6.7):

Resampling: Cross-Validated (10 fold, repeated 3 times) cp Accuracy Kappa 0.05263158 0.7822549 0.5586296 0.13157895 0.7577124 0.5022827 0.46052632 0.6075490 0.1805564 Accuracy was used to select optimal model using largest value.

The final value used for the model was cp = 0.05263158.

Рис. 6.7. Дерево rpart для прогнозирования типа стекла

В узлах и листьях дерева на рис. 6.7 указаны доля примеров в этой группе от объема всей обучающей выборки (%) и вероятность прогноза "не флеш-стекло".

Оценим адекватность модели по точности прогноза на обучающей выборке:

pred - predict( fit.cart,DGlass[,1:9]) table(Факт=DGlass$F1, Прогноз=pred) Acc - mean(pred == DGlass$F1) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт Flash No Flash 77 10 No 21 55 [1] "Точность=80.98%" Метод случайного леса (random forest) осуществляет бутстрепагрегирование некоторого множества деревьев решений.

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

set.seed(7) fit.rf - train(DGlass[,1:9], DGlass$F1, method="rf", trControl=control) fit.rf$finalModel pred - predict(fit.rf, DGlass[,1:9]) table(Факт=DGlass$F1, Прогноз=pred) Acc - mean(pred == DGlass$F1) paste("Точность=", round(100*Acc, 2), "%", sep="") Random Forest mtry Accuracy Kappa 2 0.8801879 0.7575337 5 0.8556618 0.7075469 9 0.8295588 0.6551760 Accuracy was used to select optimal model using largest value.

The final value used for the model was mtry = 2.

Call:

randomForest(x = x, y = y, mtry = param$mtry) Type of random forest: classification Number of trees: 500 No. of variables tried at each split: 2 OOB estimate of error rate: 12.27%

Confusion matrix:

Flash No class.error Flash 81 6 0.06896552 No 14 62 0.18421053 Прогноз Факт Flash No Flash 87 0 No 0 76 [1] "Точность=100%" Из протокола расчетов видно, что: а) на каждом шаге ветвления выбор идет только из двух случайных переменных, б) в итоге был составлен ансамбль из 500 деревьев, в) средняя ошибка по наблюдениям, не попавшим в "сумку" (OOB error), равна 12.27%, г) модель правильно предсказывает класс всех объектов обучающей выборки.

Приступим теперь к построению логистической регрессии LR, как это было сделано в разделе 5.1:

DGlass$F1 - ifelse(DGlass$Class==2,1,0) logit.all - glm(F1 ~.,data=DGlass[,-10], family=binomial) mp.all - predict(logit.all, type="response") Acc - mean(DGlass$F1 == ifelse(mp.all0.5,1,0)) paste("Точность=", round(100*Acc, 2), "%", sep="") [1] "Точность=73.62%" Получим теперь ошибку перекрестной проверки по методу скользящего контроля с использованием функции cv.glm() из пакета boot. Обратим внимание на необходимость функции cost, которая пересчитывает значение

delta из ошибки оценки вероятностей в ошибку предсказания класса:

library(boot) cost - function(r, pi = 0) mean(abs(r-pi) 0.5) Acc - 1 - cv.glm(DGlass[,-10], logit.all, cost)$delta[1] paste("Точность=", round(100*Acc, 2), "%", sep="") [1] "Точность=68.72%" Выполним пошаговую процедуру селекции независимых переменных и оценим ошибки предсказания:

# Регрессия на информативные переменные logit.step - step(logit.all) summary(logit.step) Start: AIC=185.56

Coefficients:

Estimate Std. Error z value Pr(|z|) (Intercept) -328.3768 352.0813 -0.933 0.350989 RI 472.1567 235.8428 2.002 0.045285 * Na -3.6626 1.1071 -3.308 0.000939 *** Mg -5.9264 1.1531 -5.140 2.75e-07 *** Si -3.7897 0.9963 -3.804 0.000143 *** K -3.1497 2.1582 -1.459 0.144455 Ca -4.9433 1.1783 -4.195 2.73e-05 *** Ba -5.6440 3.0682 -1.840 0.065838.

--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Null deviance: 225.22 on 162 degrees of freedom Residual deviance: 166.92 on 155 degrees of freedom AIC: 182.92 # Ошибка на обучающей выборке mp.step - predict(logit.step, type="response") Acc - mean(DGlass$F1== ifelse(mp.step0.5,1,0)) paste("Точность=", round(100*Acc, 2), "%", sep="") [1] "Точность=71.78%" # Ошибка при скользящем контроле Acc - 1 - cv.glm(DGlass[,-10], logit.step, cost)$delta[1] paste("Точность=", round(100*Acc, 2), "%", sep="") [1] "Точность=69.33%" Исключение нескольких переменных не привело к улучшению характеристик модели логистической регрессии. Если использовать метод RFE, то можно придти к тому же подмножеству предикторов, а тестирование моделей с применением функции train() приведет примерно к такому же соотношению эффективности обоих моделей (читатели могут проверить это самостоятельно).

6.5. Процедуры сравнения эффективности моделей классификации Итак, мы использовали шесть методов для построения моделей бинарной классификации. Осуществим их сравнение для данного примера с полным сохранением всех канонов тестирования моделей прогнозирования (Kuhn, Johnson, 2013) и в соответствии с методикой, представленной в сообщении http://machinelearningmastery.com:

1. Проводим разделение имеющейся выборки на обучающую и проверочную (использование фактора в качестве аргумента функции createDataPartition обеспечивает необходимый баланс его уровней при разделении). Для обучения выделим 70% (115 элементов) исходной выборки:

set.seed(7) train - unlist(createDataPartition(DGlass$F1, p=0.7)) # определение схемы тестирования control - trainControl(method="repeatedcv", number=10, repeats=3, classProbs = T)

2. Выполняем обучение шести моделей:

# LDA - линейный дискриминантный анализ set.seed(7) fit.lda - train(DGlass[train,3:4], DGlass$F1[train], method="lda", trControl=control) # SVML - метод опорных векторов с линейным ядром set.seed(7) fit.svL - train(DGlass[train,1:9], DGlass$F1[train], method="svmLinear", trControl=control) # SVMR - метод опорных векторов с радиальным ядром set.seed(7) fit.svR - train(DGlass[train,1:9], DGlass$F1[train], method="svmRadial", trControl=control, tuneGrid = expand.grid(sigma = 0.4, C = 2)) # CART - дерево классификации set.seed(7) fit.cart - train(DGlass[train,1:9], DGlass$F1[train], method="rpart", trControl=control) # RF - случайный лес set.seed(7) fit.rf - train(DGlass[train,1:9], DGlass$F1[train], method="rf", trControl=control) # GLM - Логистическая регрессия set.seed(7) fit.glm - train(DGlass[train,-c(1,4,10,11)],DGlass$F1[train], method="glm", family=binomial, trControl=control)

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

caret.models - list(LDA=fit.lda, SVML=fit.svL, SVMR=fit.svR, CART=fit.cart, RF=fit.rf, GLM=fit.glm) # ресэмплинг коллекции моделей results - resamples(caret.models) # обобщение различий между моделями summary(results, metric = "Accuracy") # Оценка доверительных интервалов и построение графика scales - list(x=list(relation="free"), y=list(relation="free")) dotplot(results, scales=scales)

–  –  –

Рис. 6.8. Сравнительный статистический анализ двух метрик эффективности шести моделей прогнозирования типа стекла

4. Наконец, можно вычислить статистическую значимость различий между моделями на основе распределения оценок критериев качества. Для этого воспользуемся функцией diff():

diffs - diff(results) # summarize p-values for pair-wise comparisons summary(diffs) p-value adjustment: bonferroni Accuracy

LDA SVML SVMR CART RF GLM

LDA 0.032440 -0.049573 -0.029584 -0.164472 0.035509 SVML 1.00000 -0.082012 -0.062024 -0.196911 0.003069 SVMR 0.08830 0.019988 -0.114899 0.085082 0.04908 CART 0.22316 0.42928 1.00000 -0.134887 0.065093 RF 0.199981

5.143e-09 2.525e-06 7.359e-05 1.335e-07 GLM 1.00000 1.00000 0.20047 0.03557 1.763e-07 Результаты выше главной диагонали таблицы статистической значимости содержат разности между центрами распределения. Эти значения показывают, как модели соотносятся друг с другом по абсолютному значению точности. Данные ниже диагонали таблицы представляют собой pзначения для нулевой гипотезы, которая утверждает, что параметры распределения одинаковы. Например, очевидно, что модель Random Forrest статистически значимо точнее всех остальных.

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

pred.fm - data.frame( LDA=predict(fit.lda,DGlass[-train,3:4]), SVML=predict(fit.svL, DGlass[-train,1:9]), SVMR=predict(fit.svR, DGlass[-train,1:9]), CART=predict(fit.cart, DGlass[-train,1:9]), RF=predict(fit.rf, DGlass[-train,1:9]), GLM=predict(fit.glm, DGlass[-train,-c(1,4,10,11)]) ) Резерв повышения надежности моделей прогнозирования многими исследователями видится в объединении отдельных моделей в "ансамбль" (ensemble) и построении системы коллективного распознавания.

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

Добавим в таблицу прогнозов столбец с результатами голосования (методу Random Forrest, как наиболее эффективному предсказателю, отдано два голоса):

CombyPred - apply(pred.fm, 1, function (voice) { voice2 - c(voice, voice[5]) # У RF двойной голос ifelse(sum(voice2=="Flash") 3,"Flash", "No") } ) pred.fm - cbind(pred.fm, COMB=CombyPred) head(pred.fm)

–  –  –

(classProbs = T), построим ROC-кривые для трех используемых методов распознавания (LDA, SVMR и RF) по полной таблице исходных данных, включая обучающую и тестовую последовательности. Будем анализировать вероятности отнесения к классу Flash, т.е. 1-й столбец таблицы вероятностей, возвращаемой функцией predict() при значении параметра type = "prob":

# Извлекаем вероятности классов по всей выборке pred.lda.roc - predict(fit.lda,DGlass[,3:4],type = "prob") pred.svmR.roc - predict(fit.svR,DGlass[,1:9],type = "prob") pred.rf.roc - predict(fit.rf,DGlass[,1:9],type = "prob") library(pROC) # Строим три ROC-кривые m1.roc - roc(DGlass$F1, pred.lda.roc[,1]) m2.roc - roc(DGlass$F1, pred.svmR.roc[,1]) m3.roc - roc(DGlass$F1, pred.rf.roc[,1]) plot(m1.roc, grid.col=c("green", "red"), grid=c(0.1, 0.2), print.auc=TRUE,print.thres=TRUE) plot(m2.roc, add = T, col="green", print.auc=T, print.auc.y=0.45,print.thres=TRUE) plot(m3.roc, add = T, col="blue", print.auc=T, print.auc.y=0.40,print.thres=TRUE) legend("bottomright", c("LDA","SVM","RF","COMBY"),lwd=2, col=c("black","green","blue","red"))

Рис. 6.9. ROC-кривые для четырех моделей классификации

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

Например, классификатор SVM на опорных векторах имеет определенное преимущество при предсказании объектов в критической зоне разделения (р = 0.4 0.6) и при пороговом значении p = 0.53 он достигает максимальной точности (чувствительность SE = 93.1% и специфичность SP = 88.2%). В других диапазонах p эффективность модели SVM несколько уступает остальным классификаторам, что имеет свои объяснения, если вспомнить механизм распознавания на основе опорных векторов.

Максимум значения AUC на рис. 6.8 принадлежит модели случайного леса, но использование функций ci.auc() и roc.test() даст нам доверительные интервалы для AUC и позволит сделать вывод о значимости различий между классификаторами по этому показателю.

# Доверительные интервалы для параметров ROC-анализа:

ci.auc(m2.roc); ci.auc(m3.roc) roc.test(m2.roc, m3.roc) 95% CI: 0.8949-0.9814 (DeLong) 95% CI: 0.9619-0.9988 (DeLong) DeLong's test for two correlated ROC curves data: m2.roc and m3.roc Z = 2.7057, p-value = 0.006815 Результаты свидетельствуют, что по величине AUC метод случайного леса (m3.roc) на полной выборке значимо превосходит метод опорных векторов (m2.roc). Аналогичные функции ci.se() и ci.sp() дадут возможность сделать такой же вывод относительно чувствительности и специфичности.

Вернемся к нашим попыткам осуществить коллективный прогноз.

Кроме голосования моделей, вторым простым способом комбинирования результатов бинарной классификации является расчет взвешенной средней апостериорной вероятности отнесения объекта к тому или иному классу pc = w j p j, где pj – частные вероятности, полученные g различными j классификаторами, wj – веса, нормирующие "компетентность" отдельных моделей. В качестве весов могут использоваться оценки экспертов или любые другие статистики. Будем использовать, например, значения AUC по трем выбранным методам – см. рис.

6.3, предварительно возведенные в квадрат для увеличения "контраста" и нормированные на их сумму:

weight - c(as.numeric(m1.roc$auc)^2,as.numeric(m2.roc$auc)^2, as.numeric(m3.roc$auc)^2) weight - weight/sum(weight) Pred=data.frame(lda=pred.lda.roc[,1], svm=pred.svmR.roc[,1], rf=pred.rf.roc[,1]) CombyPred - apply(Pred,1, function (x) sum(x*weight)) Построим ROC-кривую для комбинированного прогноза по этому методу и добавим ее на рис. 6.9 – COMBY. С помощью функции coord() можно вывести полную информацию об оптимальной пороговой точке (threshold) полученного классификатора.

m.Comby -roc(DGlass$F1, CombyPred) plot(m.Comby, add = T, col="red", print.auc=T, print.auc.y=0.35) coords(m.Comby, "best", ret = c("threshold", "specificity", "sensitivity", "accuracy")) threshold specificity sensitivity accuracy 0.4965734 0.9655172 0.9210526 0.9447853 Коллективные прогнозы, полученные как усреднением вероятностей, так и голосованием, имели более низкие показатели точности предсказания по сравнению с лучшими индивидуальными моделями. Видимо, такова объективная закономерность процедуры, когда решения в коллективе принимаются механическим усреднением. Далее мы снова вернемся к обсуждению этой проблемы в главе 8 при рассмотрении функций пакета ForecastCombinations.

7. МОДЕЛИ КЛАССИФИКАЦИИ ДЛЯ НЕСКОЛЬКИХ КЛАССОВ

7.1. Ирисы Фишера и метод k-ближайших соседей Ирисы Фишера – самый популярный в статистической литературе набор данных, часто используемый для иллюстрации работы различных алгоритмов классификации. При всем желании мы не смогли без него обойтись, поскольку в современных реальных приложениях редко встречаются такие компактные наборы данных, позволяющие построить хороший классификатор при минимуме исходных признаков.

Выборка состоит из 150 экземпляров ирисов трех видов (рис. 7.1), для которых измерялись четыре характеристики: длина и ширина чашелистика (Sepal.Length и Sepal.Width), длина и ширина лепестка (Petal.Length и Petal.Width). Поставим своей задачей на основе этого набора данных, включенного в базовый дистрибутив R, построить различные модели многоклассовой классификации, прогнозирующие каждый из трех видов растения по данным проведенных измерений.

–  –  –

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

data(iris) ; library(ggplot2) qplot(Sepal.Length, Sepal.Width, data = iris) + facet_grid(facets = ~ Species)+ geom_smooth(color="red", se = FALSE) qplot(Petal.Length, Petal.Width,data = iris) + facet_grid(facets = ~ Species)+ geom_smooth(color="red", se = FALSE)

Рис. 7.2. Категоризованная диаграмма выборки ирисов

Выполним теперь визуализацию анализируемых данных в пространстве первых двух главных компонент:

library(vegan) mod.pca - rda(iris[,-5], scale = TRUE) scores - as.data.frame(scores(mod.pca, display ="sites", scaling = 3)) scores$Species - iris$Species # Включаем в названия осей доли объясненных дисперсий axX - paste("PC1 (", as.integer(100*mod.pca$CA$eig[1]/sum(mod.pca$CA$eig)),"%)") axY - paste("PC2 (", as.integer(100*mod.pca$CA$eig[2]/sum(mod.pca$CA$eig)),"%)") # Составляем таблицу для "каркаса" точек l - lapply(unique(scores$Species), function(c) { f - subset(scores,Species==c); f[chull(f),]}) hull - do.call(rbind, l) # Выводим ординационную диаграмму ggplot() + geom_polygon(data=hull, aes(x=PC1,y=PC2, fill=Species), alpha=0.4, linetype=0) + geom_point(data=scores, aes(x=PC1,y=PC2,shape=Species, colour=Species), size=3) + scale_colour_manual(values = c('purple', 'green', 'blue'))+ xlab(axX) + ylab(axY) + coord_equal() + theme_bw()

–  –  –

Как видно на рис. 7.3, один из классов (Iris setosa) линейно отделим от двух остальных.

В разделах 3.4 и 4.2 мы использовали простейший, но достаточно эффективный алгоритм k ближайших соседей (k nearest neighbors), или kNNрегрессию. В основе kNN-классификатора лежит гипотеза компактности (Загоруйко, 1999) которая предполагает, что тестируемый объект d будет иметь такую же метку класса, как и обучающие объекты в локальной области его ближайшего окружения. В варианте 1NN анализируемый объект относят к определенному классу в зависимости от информации о его ближайшем соседе. В варианте kNN каждый объект относится к преобладающему классу ближайших соседей, где k – параметр алгоритма.

Решающие правила в методе kNN определяются границами смежных сегментов диаграммы Вороного (Voronoi resselation), разделяющей плоскость на |n| выпуклых многоугольников, каждый из которых содержит один и только один объект обучающей выборки (рис. 7.4). В p-мерных пространствах границы решений состоят уже из сегментов (p - 1)-мерных полуплоскостей, образованных выпуклыми многогранниками Вороного.

Рис. 7.4. Иллюстрация работы алгоритма k ближайших соседей Алгоритм предсказания строится по принципу "большинства голосов", т.е. в качестве результата объявляется метка класса-победителя. На рис. 7.4 тестируемый объект "звёздочка" попадает в ячейку объекта класса "треугольников" и при k = 1 будет отнесён к этому классу. Однако при k = 3 по голосам двух ближайших соседей из трех экзаменуемая точка будет отнесена к классу "кружочков". Вероятностный вариант метода kNN использует для ранжирования предполагаемых классов сумму "голосов" соседей с учетом их весов, в частности, евклидовой меры расстояния между тестируемым объектом и каждым из соседей.

Вариант 1NN всегда обеспечивает 100% правильное распознавание примеров обучающей выборки (самый ближайший сосед – он сам), однако часто ошибается на неизвестных ему данных. При увеличении k 1 до некоторых пределов качество распознавания на контрольной выборке будет возрастать. Оптимальное в смысле точности предсказаний значение k может быть найдено с использованием перекрестной проверки. Для этого по фиксированному значению k строится модель k-ближайших соседей и оценивается CV-ошибка классификации. Эти действия повторяются для различных k и значение, соответствующее наименьшей ошибке распознавания, принимается как оптимальное.

В среде R классификация методом ближайших соседей может быть выполнена функциями knn() из базового пакета class и kknn() из одноименного пакета с расширенными возможностями настройки меры близости и формы используемой ядерной функции.

Существует несколько возможностей оценить оптимальное значение k с помощью кросс-проверки:

1. Воспользоваться функцией train.kknn() из пакета kknn (по умолчанию задается скользящий контроль):

library(kknn) # -------- Метод к-ближайших соседей kNN train.kknn(Species ~., iris, kmax = 50, kernel="rectangular") Type of response variable: nominal Minimal misclassification: 0.02666667 Best kernel: rectangular Best k: 16

2. Составить собственный скрипт для выполнения перекрестной проверки, например (рис. 7.5):

max_K=20 ; gen.err.kknn - numeric(max_K) mycv.err.kknn - numeric(max_K) ; n - nrow(iris) # Рассматриваем число возможных соседей от 1 до 20 for (k.val in 1:max_K) { pred.train - kknn(Species ~., train=iris,test=iris,k=k.val,kernel="rectangular") gen.err.kknn[k.val] - mean(pred.train$fit != iris$Species) for (i in 1:n) { pred.mycv - kknn(Species~., train=iris[-i,],test=iris[i,], k=k.val,kernel="rectangular") mycv.err.kknn[k.val] - mycv.err.kknn[k.val] + (pred.mycv$fit != iris$Species[i]) } } ; mycv.err.kknn - mycv.err.kknn/n plot(1:20,gen.err.kknn,type="l",xlab='k', ylim=c(0,0.07), ylab='Ошибка классификации', col="limegreen", lwd=2) points(1:max_K,mycv.err.kknn,type="l",col="red", lwd=2) legend("bottomright", c("При обучении", "Скользящий контроль"), lwd=2, col=c("limegreen", "red")) Рис. 7.5. Нахождение оптимального числа k ближайших соседей при классификации видов ириса

3. Использовать функцию train() из пакета caret, с которой мы уже неоднократно сталкивались ранее (см. главу 3):

library(caret) ; library(class) ; set.seed(123) contrl - trainControl(method="repeatedcv", repeats = 3) train(Species ~., data = iris, method = "kknn", trControl = contrl, preProcess = c("center", "scale"), tuneLength = 20) Resampling: Cross-Validation (10 fold, repeated 3 times) k Accuracy Kappa Accuracy SD Kappa SD 7 0.96 0.94 0.0414 0.0621 9 0.96 0.94 0.0414 0.0621 11 0.962 0.943 0.0453 0.0679 13 0.967 0.95 0.042 0.063 15 0.964 0.947 0.0454 0.0681 17 0.958 0.937 0.0479 0.0718 19 0.958 0.937 0.0446 0.0669 21 0.951 0.927 0.0461 0.0691 Accuracy was used to select the optimal model using the largest value.The final value used for the model was k = 13.

Здесь была осуществлена трехкратная перекрестная проверка с разбивкой на 10 блоков, т.е. случайным образом исходная выборка делилась на 10 частей, после чего поочередно 9 частей использовались для оценки k, а по одной части оценивалась точность классификации, и такая процедура повторялась три раза.

Выполним теперь стандартную процедуру тестирования модели при k = 13.

Исходную выборку разбиваем на две части: обучающую (train) и проверочную (test) в соотношении 3:1:

set.seed(123) ; (samp.size - floor(nrow(iris) *.75)) train.ind - sample(seq_len(nrow(iris)), size = samp.size) train - iris[train.ind,] ; test - iris[-train.ind, ] [1] 112 Функция knn() из пакета class сразу на выходе дает предсказываемые значения для тестовой последовательности, которые можно использовать для построения матрицы неточностей:

knn.iris - knn(train = train[,-5], test = test[,-5], cl = train[,"Species"], k = 13, prob = T) table(Факт=test$Species,Прогноз=knn.iris) Acc = mean(knn.iris == test$Species) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт setosa versicolor virginica setosa 11 0 0 versicolor 0 13 0 virginica 0 1 13 [1] "Точность=97.37%" Метод дал одну ошибку при предсказании 38 "свежих" примеров, что является вполне приличным результатом.

–  –  –

Остается только позаботиться, чтобы значения логарифмируемых вероятностей были не слишком близки к 0 (чтобы этого не случилось, можно применить лапласовское сглаживание – http://stats.stackexchange.com).

В среде R расчеты выполняются обычно с использованием функций NaiveBayes() из пакета klaR или naiveBayes() из пакета e1071.

Рассмотрим их использование на примере классификации цветков ириса:

library(klaR) naive_iris - NaiveBayes(iris$Species ~., data = iris) naive_iris$tables$Petal.Width [,1] [,2] setosa 0.246 0.1053856 versicolor 1.326 0.1977527 virginica 2.026 0.2746501 Для каждой метрической независимой переменной были выведены средние значения Petal.Width (первый столбец) и их стандартные отклонения (второй столбец) для каждого выделенного класса. Можно установить, что ширина лепестка у виргинского ириса существенно выше, чем у остальных двух видов.

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

7.6):

plot(naive_iris, lwd=2)

–  –  –

Выполним прогноз вида ирисов для объектов из обучающей выборки:

pred - predict(naive_iris, iris[,-5])$class table(Факт=iris$Species, Прогноз=pred) Acc = mean(pred == iris$Species) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт setosa versicolor virginica setosa 50 0 0 versicolor 0 47 3 virginica 0 3 47 [1] "Точность=96.00%" Осуществим 10-кратную перекрестную проверку построенной модели, чтобы установить качество ее предсказаний на контрольных примерах.

Одновременно уточним значения некоторых гиперпараметров процедуры вычисления вероятностей по частотам: надо ли использовать ядерное сглаживание (usekernel=TRUE) или можно предположить нормальное распределение, а также следует ли проводить коррекцию вероятностей по Лапласу fL. Для этого снова воспользуемся функцией train() из пакета

caret:

library(caret) # Определим условия перекрестной проверки:

train_control - trainControl(method='cv',number=10) Test - train(Species~., data = iris, trControl=train_control, method="nb") print(Test) Acc = mean(predict(Test$finalModel,iris[,-5])$class == iris$Species) paste("Точность=", round(100*Acc, 2), "%", sep="") Resampling: Cross-Validation (10 fold) usekernel Accuracy Kappa Accuracy SD Kappa SD FALSE 0.96 0.94 0.0562 0.0843 TRUE 0.96 0.94 0.0466 0.0699 Tuning parameter 'fL' was held constant at a value of 0 Accuracy was used to select optimal model using largest value.

The final values used for the model were fL = 0 and usekernel = FALSE.

[1] "Точность=96.00%" Более качественной модели в результате перекрестной проверки построить не удалось, но точность классификатора на обучающей и контрольной выборках осталась стабильно высокой.

7.3. Классификация в линейном дискриминантном пространстве Исходные предпосылки линейного дискриминантного анализа, изложенные нами в разделе 6.1, естественным образом распространяются на случай произвольного числа классов k. Если вновь обратиться к аналогам евклидовой геометрии, то центроиды двух классов определяют положение прямой линии, трех – плоскости, четырех – четырехмерной гиперплоскости и т.д. Принцип сводится к тому, что тем самым определяется "натянутое на центроиды" дискриминантное пространство, размерность которого на единицу меньше, чем число классов.

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

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

–  –  –

где k – номер линейного дискриминанта, xj – значение j-й переменной, x j – среднее значение j-й переменной (по всем группам), kj – весовой коэффициент j-й переменной в k-ом линейном дискриминанте; суммирование ведется по m переменным. Первый дискриминант подбирается таким образом, чтобы максимизировать основные различия между группами, второй и последующие также максимизируют оставшуюся межгрупповую изменчивость, но подбираются ортогонально первому и всем предыдущим.

Отсюда следует, что любая каноническая дискриминантная функция LDk имеет нулевую корреляцию с LD1, …, LDg-1. Если число групп равно g, то число канонических дискриминантных функций будет на единицу меньше числа групп.

Один из методов поиска оптимальных дискриминантных функций заключается в максимизации отношения межгрупповой вариации к внутригрупповой, т.е. B / W max, где B – межгрупповая, а W – внутригрупповая матрица рассеяния наблюдаемых переменных относительно средних. Тогда нахождение канонических коэффициентов дискриминантных функций сводится к нахождению собственных значений k и векторов uk: если спроецировать g групп m–мерных выборок на (g – 1)мерное пространство, порожденное собственными векторами uk, то отношение = B / W будет максимальным, т.е. разброс между группами будет наиболее контрастным при заданной внутригрупповой вариации.

Обратимся вновь к распознаванию видов цветков ириса, выборка которых носит имя сэра Р.

Фишера, поскольку именно на ней он и продемонстрировал возможности дискриминантного анализа:

require(MASS) LDA.iris - lda(formula = Species ~.,data = iris) LDA.iris$scaling # Коэффициенты линейных дискриминантов LD1 LD2 Sepal.Length 0.8293776 0.02410215 Sepal.Width 1.5344731 2.16452123 Petal.Length -2.2012117 -0.93192121 Petal.Width -2.8104603 2.83918785 LDA.iris$svd (prop = LDA.iris$svd^2/sum(LDA.iris$svd^2)) [1] 48.642644 4.579983 [1] 0.991212605 0.008787395 Последними действиями мы вывели значения, которые являются отношением межгрупповой дисперсии к внутригрупповой дисперсии, и вычислили долю межгрупповой дисперсии, объясненную каждым линейным дискриминантом. На рис. 7.7 представлена ординационная диаграмма наблюдений в осях LD1-LD2, подобная графику главных компонент на рис.

7.3. Можно обратить внимание на то, что трансгрессия точек между видами практически исчезла, а главная ось дискриминации объясняет 99% межгруппового разброса:

prop - percent(prop) pred -predict(LDA.iris, newdata = iris) scores - data.frame(Species = iris$Species,pred$x) # Выводим ординационную диаграмму require(ggplot2) ggplot() + geom_point(data=scores, aes(x=LD1,y=LD2,shape=Species,colour=Species),size=3) + scale_colour_manual( values = c('purple', 'green', 'blue'))+ labs(x = paste("LD1 (", prop[1], ")", sep=""), y = paste("LD2 (", prop[2], ")", sep="")) + theme_bw()

–  –  –

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

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

lda(scale(iris[,1:4]), gr = iris$Species)$scaling LD1 LD2 Sepal.Length 0.6867795 0.01995817 Sepal.Width 0.6688251 0.94344183 Petal.Length -3.8857950 -1.64511887 Petal.Width -2.1422387 2.16413593 В результате применения линейного дискриминантного анализа к примерам из обучающей выборки можно сформировать решающее правило, позволяющее отнести новый объект к классу, имеющему наибольшую плотность вероятности ("уровень правдоподобия") в точке х m-мерного пространства признаков. Сама процедура классификации основана на вычислении апостериорных вероятностей Pr(Gk | x) принадлежности произвольного вектора наблюдений х к группе Gk, k = 1, 2,…, g, с учетом анализа функций плотности распределения fk(х).

Величина Pr(Gk | x) связана с дистанцией между точкой x и k-м центроидом, и в случае выполнения предпосылок дискриминантного анализа может быть оценена с использованием значений классифицирующих функций:

dk(x) = µkT С-1x - 0.5 µkT С-1µk + lnk, где С-1 – матрица, обратная к ковариационной, µк и lnk – вектор средних и априорная вероятность наблюдения объектов k-го класса соответственно.

Объект x относится к тому классу k из g, для которого значение dk, а следовательно и апостериорная вероятность, максимальны.

Если выделить из обучающей выборки 10 объектов для проверки качества предсказаний, то получим следующие (безошибочные) результаты:

train - sample(1:150, 140) LDA.iris3 - lda(Species ~., iris, subset = train) plda = predict(LDA.iris3, newdata = iris[-train, ]) data.frame(Species=iris[-train,5], plda$class, plda$posterior) Species plda.class setosa versicolor virginica 32 setosa setosa 1.000000e+00 2.075304e-20 1.232723e-39 38 setosa setosa 1.000000e+00 1.274233e-24 1.270493e-45 45 setosa setosa 1.000000e+00 1.185866e-18 1.103632e-36 56 versicolor versicolor 7.594352e-24 9.982990e-01 1.700985e-03 81 versicolor versicolor 1.751434e-18 9.999970e-01 2.952626e-06 94 versicolor versicolor 1.020832e-14 9.999999e-01 9.062344e-08 97 versicolor versicolor 4.091110e-20 9.998948e-01 1.052484e-04 116 virginica virginica 1.386213e-41 3.491571e-05 9.999651e-01 119 virginica virginica 5.658468e-63 5.997144e-10 1.000000e+00 146 virginica virginica 8.347791e-41 9.898162e-05 9.999010e-01 Подчеркнуты значения максимальных апостериорных вероятностей, на основании которых предсказывался вид цветка plda.class.

При анализе всей обучающей выборки цветков ириса мы имеем три ошибочных результата классификации из 150.

Acc = mean(pred$class == iris$Species) paste("Точность=", round(100*Acc, 2), "%", sep="") [1] "Точность=98.00%"

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

LDA.irisCV - lda(Species ~., data = iris, CV = TRUE) table(Факт=iris$Species,Прогноз=LDA.irisCV$class) Acc = mean(LDA.irisCV$class==iris$Species) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт setosa versicolor virginica setosa 50 0 0 versicolor 0 48 2 virginica 0 1 49 [1] "Точность=98.00%"

7.4. Нелинейные классификаторы в R Квадратичный дискриминантный анализ (QDA) является нелинейным обобщением метода LDA, который не использует предположения об однородности ковариационной матрицы. В качестве решающего правила применяется квадратичная функция dk(x) = - 0.5 (x - µk)TCk-1(x - µk) – 0.5 ln|Ck| + lnk, где |Ck| – детерминант ковариационной матрицы k-го класса, Ck-1 – ее обратная матрица; k – априорная вероятность наблюдения объектов k-го класса. Экзаменуемый объект, как и в линейном случае, относится к классу с максимальным значением dk.

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

Однако он сохраняет большинство недостатков LDA: использует предположение о нормальности распределения и не работает, когда матрицы ковариаций вырождены (например, при большом числе переменных). Другим недостатком QDA является то, что уравнение разделяющей гиперповерхности выражено в неявном виде и не может быть использовано для "объяснения" результатов.

Для проведения квадратичного дискриминантного анализа можно воспользоваться функцией qda() из пакета MASS:

require(MASS) QDA.iris = qda(Species~ Petal.Length + Petal.Width, data = iris) pred = predict(QDA.iris)$class table(Факт=iris$Species,Прогноз=pred) Acc = mean(pred ==iris$Species) paste("Точность=", round(100*Acc, 2), "%", sep="")

–  –  –

Рис. 7.8.

Разделяющие границы, полученные в результате применения квадратичного дискриминантого анализа для классификации цветков ириса (ошибочное распознавание выделено шрифтом красного цвета); жирными точками показаны центроиды классов C использованием функции train() из пакета caret сравним точность классификации по полной модели и модели с двумя предикторами на основании 10-кратной перекрестной проверки:

library(caret) set.seed(123) # Используем только размеры лепестка train(Species ~ Petal.Length + Petal.Width, data = iris, method = "qda", trControl = trainControl(method = "cv")) Accuracy Kappa Accuracy SD Kappa SD 0.973 0.96 0.0344 0.0516 # Используем весь набор признаков train(Species ~., data = iris, method = "qda", trControl = trainControl(method = "cv")) Accuracy Kappa Accuracy SD Kappa SD 0.973 0.96 0.0344 0.0516 Очевидно, что ботаникам вполне можно было обойтись без утомительной работы по фиксации размеров чашелистника.

Регуляризованный дискриминантный анализ (RDA – Regularized Discriminant Analysis), предложенный Дж. Фридманом, является своеобразным промежуточным звеном между LDA и QDA. Пусть – положительный параметр, подобранный таким образом, чтобы управлять степенью кривизны разделяющей поверхности.

При его использовании регуляризованная ковариационная матрица Ck () оказывается как бы более "притянутой" к общей ковариационной матрице C0, чем индивидуальная Ck для k-го класса:

Ck() = (1 - ) Ck + C0.

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

Более точная регуляризация получается при искусственном увеличении числа доступных объектов обучения за счет добавления к данным шума, интенсивность которого связана со вторым параметром – :

Ck (, ) = (1 - ) Ck() + trace[Ck()] I / m, где trace[Ck()] – диагональные элементы матрицы Ck(), или вектор дисперсий исходных m переменных, описывающих объекты k-го класса; I – единичная матрица.

При экстремальных значениях и ковариационная структура сводится к четырем особым случаям:

(=0, =0): QDA – индивидуальная ковариация в каждой группе;

(=0, =1): LDA – общая матрица ковариации;

(=1, =0): исходные переменные условно независимы (подобно тому, что предполагает наивный байесовский классификатор), но их внутригрупповые дисперсии однородны;

(=1, =1): дисперсии равны для всех групп и классификация объектов осуществляется по евклидовым расстояниям до центроидов.

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

library(klaR); set.seed(123) # 1 - этап грубой оптимизации train(Species ~., data = iris, method = "rda", trControl = trainControl(method = "cv")) gamma lambda Accuracy Kappa Accuracy SD Kappa SD 0 0 0.973 0.96 0.0344 0.0516 0 0.5 0.98 0.97 0.0322 0.0483 0 1 0.98 0.97 0.0322 0.0483 0.5 0 0.967 0.95 0.0471 0.0707 0.5 0.5 0.967 0.95 0.0471 0.0707 0.5 1 0.96 0.94 0.0466 0.0699 # 2 - этап с диапазоном lambda = 0.1:0.5 и gamma = 0.02:0.1 RDAGrid - expand.grid(.lambda = (1:5)/10,.gamma = (1:5)/50) RDA.iris - train(Species ~., data = iris, method = "rda", tuneGrid =RDAGrid, trControl = trainControl(method = "cv")) gamma lambda Accuracy Kappa Accuracy SD Kappa SD 0.02 0.3 0.98 0.97 0.0322 0.0483 0.02 0.4 0.98 0.97 0.0322 0.0483 0.02 0.5 0.98 0.97 0.0322 0.0483 0.04 0.1 0.98 0.97 0.0322 0.0483 0.04 0.2 0.98 0.97 0.0322 0.0483 Accuracy was used to select optimal model using the largest value.

The final values used for model were gamma = 0.02 and lambda = 0.5.

# Оценка точности итоговой модели:

pred = predict(RDA.iris) Acc = mean(pred ==iris$Species) paste("Точность=", round(100*Acc, 2), "%", sep="") [1] "Точность=98.00%" Обратим внимание, что при выполнении скрипта нужно отключить загруженный пакет vegan, где также есть функция rda().

Машина опорных векторов SVM, описанная в разделах 6.2-6.3, является бинарным классификатором и естественным образом на случай с несколькими классами не переносится. Однако имеется возможность использовать этот метод, применяя стратегии "каждый против каждого" (one against one) или "один против всех" (one against all).

Применим метод опорных векторов к данным по ирисам для нахождения разделяющих поверхностей, которые по умолчанию в функции svm() формируются с использованием радиального ядра RBF:

library(e1071) SVM.iris - svm(Species~., data=iris) pred.DV - predict(SVM.iris, iris[,1:4], decision.values = TRUE) ind - sort(sample(150,9)) data.frame(attr(pred.DV, "decision.values")[ind,], iris$Species[ind])

–  –  –

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

plot(cmdscale(dist(iris[,-5])), col = as.integer(iris[,5])+1, pch = c("o","+")[1:150 %in% SVM.iris$index + 1], font=2, xlab="Шкала 1",ylab="Шкала 2" ) legend (0,1.2, c("setosa","versicolor","virginica"),pch = "o", col =2:4) Здесь функция dist() строит матрицу евклидовых расстояний размерностью 150150 для всех пар наблюдений по всей совокупности измеренных показателей. Функция cmdscale() проецирует эту матрицу на плоскость таким образом, чтобы минимизировать искажения взаимных дистанций между объектами при сокращении исходного многомерного пространства до двумерной ординации (рис. 7.9). Опорные векторы, представленные на диаграмме крестиками, задают границы разделяющих RBF-поверхностей, аккуратно окаймляющих области каждого класса.

Рис. 7.9. Ординационная диаграмма наблюдений из набора данных по ирисам, построенная методом многомерного шкалирования (крестиками соответствующего цвета показаны опорные векторы)

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

Рассмотрим сначала модель логита (multinomial logit) для отклика с несколькими классами. Пусть J обозначает число классов для случайной независимой величины Y (требование независимости всех J альтернатив тут играет важное значение). Многомерное распределение вероятности наблюдения каждого из ее классов может быть представлено вектором {1,..., J}, j j = 1. Тогда для каждого из n независимых наблюдений оценка j будет соответствовать вероятности того, насколько предпочтительно присвоение произвольному тестируемому объекту метки класса j.

Отметим, что для оценки всех вероятностей достаточно построить J моделей, поскольку недостающее значение легко получить из условий нормировки (сумма j всегда равна 1).

Поэтому один из классов можно принять за базовый (примем для определенности, что это класс с номером 1), и тогда совокупность логит-моделей (baseline-category logit) будет иметь вид:

log(j / 1) = j + j x, j = 2,..., J.

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

Для анализа категориальных данных в R существуют мощные специализированные пакеты VGAM и mlogit. Мы же воспользуемся более простой в освоении функцией multinom() из пакета nnet, которая осуществляет оценку коэффициентов системы логит-моделей с использованием алгоритмов построения искусственных нейронных сетей (см.

раздел 7.6).

Для данных по ирисам Фишера получаем:

library(nnet) MN.iris - multinom(Species~., data=iris) summary(MN.iris)

Coefficients:

(Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width versicolor 18.69037 -5.458424 -8.707401 14.24477 -3.097684 virginica -23.83628 -7.923634 -15.370769 23.65978 15.135301

Std. Errors:

(Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width versicolor 34.97116 89.89215 157.0415 60.19170 45.48852 virginica 35.76649 89.91153 157.1196 60.46753 45.93406 Residual Deviance: 11.89973 AIC: 31.8997 Анализируя коэффициенты, можно заметить, например, что ширина чашелистика Sepal.Width уменьшается при переходе от setosa к virginica, причем логарифм отношения шансов log(vir / set) уменьшается на величину 23 = -15.3 при увеличении ширины на 1 см.

Метод fitted(), применяемый к объекту multinom, возвращает матрицу апостериорных вероятностей принадлежности каждого наблюдения к соответствующим классам.

Предсказание класса осуществляется с учетом максимума такой вероятности:

Probs - fitted(MN.iris) setosa versicolor virginica 1 1.000000e+00 1.526406e-09 2.716417e-36 2 9.999996e-01 3.536476e-07 2.883729e-32 3 1.000000e+00 4.443506e-08 6.103424e-34 pred=apply(Probs,1,function(x) colnames(Probs)[which(x==max(x))]) table(Факт=iris$Species,Прогноз=pred) Acc = mean(pred ==iris$Species) paste("Точность=", round(100*Acc, 2), "%", sep="")

–  –  –

# p-значения на основе t-статистики pt(z, df = nrow(iris) - 5, lower=FALSE) (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width versicolor 0.2969240 0.5241679 0.5220705 0.4066286 0.5270993 virginica 0.7469061 0.5350513 0.5388982 0.3480821 0.3711264 И мы с недоумением видим, что, несмотря на приличную точность предсказаний модели, все отношения log(k / 1) согласно вычисленным рзначениям статистически не отличаются от 0, т.е. все шансы равны. Ответ видится в резковатой реплике проф. Б. Рипли (B. Ripley) при дискуссии в Интернете: «Попробуйте почитать наши книги. Там объясняется, почему тесты Вальда выполнять нельзя, и что асимптотическая теория может здесь дико вводить в заблуждение». Читатели, желающие убедиться в правоте сказанного могут пересчитать р-значения коэффициентов с использованием теста множителей Лагранжа или бутстрепа, в частности, функцией cluster.im.mlogit() из пакета clusterSEs.

7.6. Классификаторы на основе искусственных нейронных сетей Нейросетевые модели, родившиеся в процессе развития концепции искусственного интеллекта, имеют две вполне прозрачные аналогии – биологическая нейронная система мозга и компьютерная сеть. Их основная парадигма состоит в том, что решение в сети формируется множеством простых нейроно-подобных элементов, образующих граф с взвешенными синаптическими (информационными) связями, которые совместно и целенаправленно работают на получение общего результата. Таким образом, искусственная нейронная сеть (ИНС) является (Дук, Самойленко, 2001):

параллельной системой, поскольку в любой момент времени в активном состоянии могут находиться несколько процессов;

распределенной системой, поскольку каждый из процессов может независимо обрабатывать локальные данные;

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

Главным строительным блоком ИНС является формальный нейрон, структура которого имеет вид, представленный на рис.

7.10:

Рис. 7.10. Структура искусственного нейрона (слева) и вид некоторых функций активации F (справа); обозначения по тексту Основная функция искусственного нейрона – cформировать выходной сигнал y в зависимости от сигналов x1...xn, поступающих на его входы. Эти значения могут усиливаться или "тормозиться" в зависимости от знака весов синапсов w1...wn. Входные сигналы обрабатываются адаптивным n сумматором S = wi xi T, где T – порог нейрона, а затем выходной сигнал i =1 сумматора поступает в нелинейный преобразователь F с некоторой функцией активации, после чего результат подается на выход (в точку ветвления).

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

Например, преобразование может осуществляться функцией, определенной как единичный скачок (а на рис. 7.10 справа), линейный порог (гистерезис, б), гиперболический тангенс (в) или сигмоид (г). Наиболее распространена нелинейная функция с насыщением – так называемая логистическая функция, или сигмоид, y =, которая имеет много ценных свойств (монотонность 1 + e cS и дифференцируемость, устойчивость к выбросам, простое выражение для ее производной и т.д.). Выходное значение сигмоидального нейрона лежит в диапазоне [0, 1]. При уменьшении коэффициента с сигмоид становится более пологим, вырождаясь в пределе при с = 0 в горизонтальную линию на уровне

0.5. При увеличении с сигмоид приближается по внешнему виду к функции единичного скачка с порогом T в точке x = 0.

По способу соединения нейронов выделяют сети с разной архитектурой: персептроны, сети адаптивного резонанса, рециркуляционные, рекуррентные, встречного распространения, ИНС с обратными связями Хэмминга и Хопфилда, с двунаправленной ассоциативной памятью, с радиально-базисной функцией активации, самоорганизующиеся ИНС Кохонена и т.п.

Мы будем рассматривать возможности обучения ИНС на примере многослойного персептрона с прямым распространением сигнала (и обратным распространением ошибки), структура которого представлена на рис 7.11:

Рис. 7.11. Структура многослойного персептрона

–  –  –

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



Pages:     | 1 | 2 || 4 | 5 |

Похожие работы:

«ПАРАЗИТОЛОГИЯ, 55, 2007 УДК 576.895.775 О С О Б Е Н Н О С Т И С Т Р О Е Н И Я ГРУДНЫХ И Б Р Ю Ш Н Ы Х К Т Е Н И Д И Е В БЛОХ (SIPHONAPTERA) © С.Г.Медведев Методами световой и растровой электронной микроскопии изучено строение псевдосет, зубчиков и зубцов гребней (ктенидиев) у представителей 80 % родов и подродов мировой фауны о...»

«ДЕПАРТАМЕНТ СМОЛЕНСКОЙ ОБЛАСТИ ПО ЗДРАВООХРАНЕНИЮ 214008, г. Смоленск, пл. Ленина, д. 1 Тел.: (4812) 38-67-58, e-mail: zdrav@admin.smolensk.ru окпо: 00097063, огрн: 1026701426069 инн: 6730009960, кпп: 673001001 от № на № от Выписка из приказа Департамента Смоленской области по здравоохранению от 10.09.2013 № 1221 В соответствии со ст. 18 Федерально...»

«Российская библиотека Холокоста МЫ НЕ МОЖЕМ МОЛЧАТЬ ШКОЛЬНИКИ И СТУДЕНТЫ О ХОЛОКОСТЕ Выпуск 11 Составители М.В. Гилёва, Д.В. Прокудин Москва Центр и Фонд "Холокост" Издательство "МИК" УДК 63.3(0)62 ББК 94(100)"1939/1945" М 94 "РОССИЙСКАЯ БИБЛИОТЕКА ХОЛОКОСТА" И. А. Альтман (отв. редактор), А. Е. Гербер, Ю. А. Домбр...»

«Брось сигарету! (Петрова Зинаида Валерьевна врач психиатр-нарколог ГУЗ "РНД") Пристрастие к сигаретам все травите себя табаком Не чаще становится причиной конфликтов в общественных местах, на производстве, в семье. Сегодня курение сталоГУЗ "РНД") не толь...»

«Сеть социальных контактов: мобилизация социального окружения детей и семей в кризисной ситуации Редакторы: Наталия Власова, Бритта Хольмберг, Наталия Снурникова МОСКВА, КИРИЛЛ 30 августа. 16/00. еще день и в школу. Мама снова пьяна. Бабушка с утр...»

«Клод Дюран Франция НЕКОТОРЫЕ ЗАМЕТКИ ОБ ИСПОЛЬЗОВАНИИ ПОСЛОВИЦ В ЛИТЕРАТУРЕ, В ОСОБЕННОСТИ В РУССКОЙ ЛИТЕРАТУРЕ, В ЧАСТНОСТИ — В "КРАСНОМ КОЛЕСЕ" АЛЕКСАНДРА СОЛЖЕНИЦЫНА1 Жизнь и творчество Александра Солженицына: на пути к "...»

«О ПРЕДПРИЯТИИ Северодонецкое научно-производственное объединение Импульс – разработчик, производитель и поставщик высоконадежных информационных и управляющих систем для атомной и тепловой энергетики, нефтегазового комплекса, железнодорожного транспорта. Фирма основана в 1956 году как базово...»

«М. Хайдеггер. Что значит мыслить? Мы попадаем в то, что называется мышлением, когда мыслим сами. Чтобы нам это удалось, мы должны быть готовы учиться мыслить. Как только мы принимаемся за это учение, мы сразу по...»

«А К А Д Е М И Я Н А У К СССР т •Р У Д JM ИНСТИТУТА ГЕОЛОГИЧЕСКИХ НАУК ВЫП. 123. ГЕОЛОГИЧЕСКАЯ СЕРИЯ (№ 44). 1950 В. В. Т И Х О М И Р О В МАЛЫЙ КАВКАЗ В ВЕРХНЕМЕЛОВОЕ ВРЕМЯ (ОСНОВНЫЕ ТИПЫ ОТЛОЖЕНИЙ И УСЛОВИЯ ИХ ОБРАЗОВАНИЯ) ИЗДАТЕЛЬСТВО АКАДЕМИИ НАУК СССР...»

«КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ 2009 Т. 1 № 3 С. 321–336 АНАЛИЗ И МОДЕЛИРОВАНИЕ СЛОЖНЫХ ЖИВЫХ СИСТЕМ Ангармонические колебательные резонансы в малых водных ассоциатах А. В. Карговский1,a Московский госу...»

«Дифф. геометрия и векторные расслоения М. Вербицкий Векторные расслоения, лекция 2: связности Миша Вербицкий 16 сентября, 2013 матфак ВШЭ и НМУ Дифф. геометрия и векторные расслоения М. Вербицкий Векторные расслоения (повторение) ОПРЕДЕЛЕНИЕ: Векторное р...»

«Титульный лист рабочей учебной программа Форма Ф СО ПГУ 7.18.3/30 Министерство образования и науки Республики Казахстан Павлодарский государственный университет им. С. Торайгырова Кафедра "Исполнительское искусство" РАБОЧАЯ УЧЕБНАЯ ПРОГРАММА дисциплины "Основы теории музыки и гармонии" для студентов специальности 5В01060...»

«Электрически обогреваемые печи для работы в среде защитного газа и вакуума Авторы: Р. Вайтц, д-р. П. Вюббен, Б. Гайс, В. Мюллер В этой статье описаны виды печей для термообработки материалов в среде защитного газа и вакуума. В начале дается краткое описание об основе использования защитного газа и...»

«Измерение шума приборами серии АССИСТЕНТ. Краткое описание прибора. Рекомендации по использованию. Методики выполнения измерений. Москва 2009 Описаны эффективные приемы работы при измерении шума приборами серии АССИСТЕНТ. Измерение и частотный анализ...»

«УДК 821.161.1-312.4 ББК 84(2Рос=Рус)6-44 А94 Афанасьев, Александр. А94 Воздушные головорезы / Александр Афанасьев. — Москва : Эксмо, 2016. — 320 с. — (Спецназ. Группа Антитеррор). ISBN 978-5-699-89513-7 С таким коварным замыслом мир еще не сталкивался: террористы запланировали атаки на крупные российские города, в которы...»

«59 Теоретическая и прикладная лингвистика, 2017, 3 (1), 59–66 УДК 811’111 UDC 811’111 Матвеева Елена Владимировна Амурский государственный университет г. Благовещенск, Российская Федерация Elena V. Matveeva Amur State University Blagoveshchensk, Russian Federation e-mail: matveeva_555@mail....»

«Author: Огородников Вадим Зиновьевич Сахалин: Люди. Жители. Временщики.                                                Люди. Жители. Временщики. Начало всех начал и начальник и отец командир в армиях всех государств и народов, командир части. Командиром час...»

«DONIDO ОБОРУДОВАНИЕ И ТЕХНОЛОГИИ ДЛЯ ПЕРЕРАБОТКИ МОЛОКА ПРОЕКТИРОВАНИЕ ПРОИЗВОДСТВА Поиск оптимального решения Диалог с Заказчиком Анализ исходных данных и пожеланий Заказчика • ассортимент производимой продукции;• количество производимой продукции;• энергетические и человеческие ресурсы;• ожидаемая себестоимость готового продукта; • ожидаем...»

«Оборудование слива-налива (ОСН) Стояк Верхнего Налива ОСН-СВН-100-Н-А-4-НН с Ду=100 для Автоналива сырой Нефти с зоной обслуживания 4 м (корневой фланец вНиз, патрубок №2 Низ) Сертификат соответствия ГОСТ Р № РОСС RU/ME60.В00629. Разрешение на применение Госгортехнадзора России № РРС 02-5223...»

«Московский государственный университет им. М.В. Ломоносова Географический факультет Научно-исследовательская лаборатория эрозии почв и русловых процессов им. Н.И. Маккавеева МАККАВЕЕВСКИЕ ЧТЕНИЯ – 2005 Научный редактор – профессор Р.С. Чалов Москва – 2006 УДК 6.31.4: 55.3 М...»

«Антон ФАРБ ГЛИФ ЧАСТЬ ПЕРВАЯ. КАК СКУЧНО ЖИТЬ В ПРОВИНЦИИ. С приходом весны с дорог сошел снег, а вместе со снегом сошел и асфальт. Трасса от Киева была еще куда ни шло, но после Коростышева трясти начало так, что даже кино стало смотреть невозможно. Ника закрыла ноутбук, вытащила наушн...»

«Карта установки охранного комплекса систем Scher-Khan Magicar 11 и Falcon CI-20 на автомобиль Mazda 6 (2012) (ver.2) Оглавление Устанавливаемые компоненты: Общая информация Подключение питания системы Подключение к CAN шине автомобиля Управление центральным замком автомобиля Подключение...»










 
2017 www.book.lib-i.ru - «Бесплатная электронная библиотека - электронные ресурсы»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.