- Надо всё-таки готовить заранее текст и прогнать его хотя бы раз перед выступлением. А лучше два или три. Ну ладно, второй жертвой будет мама)
- Чётче сформулировать цель работы - наши преподы полагают будто мы собираемся заниматься какими-то серьезными научными исследованиями в бакланской-то работе. Ну и да... Цель - обработать полученные данные.
- Слушателям было не понятно, с какими данными я работаю и какие они => пожалуй, всё-таки стоит включить в слайды скрин эксель-таблички или как-то расписать их. И сказать им, что эти данные собраны из интернета и не мной.
- Надо связать гравитационную модель миграций и линейную регрессию. Яшунский сказал, что если прологарифмировать гравитац. модель, то получим сумму логарифмов регрессоров модели. (то есть независимых параметров, от которых зависит количество миграций)
- Сделать summary по различным тестам: названия, ограничения, сложность, количество обнаруженных выбросов, процент существенных выбросов. Ну и в заключении сказать, какой из тестов (или какая последовательность) была оптимальной.
- таки выделить список outlier-ов для каждого города с количеством мигрантов (пара строчек кода), потом обобщить и эти полученные данные.
- Показать, что если учитывать выбросы (несущественные), то наша модель хорошо аппроксимируется прямой (собственно, это нам и надо).
- А не запилить ли мне блок-схемы алгоритмов? Так и мне прогать будет проще, да и в диплом вставить полезно будет.
- P.S. На каждом графике по городу представлена зависимость реальной и предсказанной миграции из города i в города a, b, ..., g, h, j, ....z. Никогда бы не подумала, что слушатели могут этого не понять из контекста. Забыла, что мысли читать умеют далеко не все! XD
Мда-а-а-а... выглядит списочек внушительно. Ну что ж, начнём, пожалуй.
О данных и о цели
Привет,
сводная таблица статистики_перемещений_контакта х данные_из_википедии
(широта, долгота, коды ОКТМО, регионы) и, наверное, какого-нибудь
Росстата или откуда там берутся данные по населению? Выглядит она так:
Да уж, выглядит она действительно как жирный судак - более чем на 17К строк! Табличка сведена по городам FROM. Их всего 26:
[1] Большой Камень Братск Верхотурье Вилючинск Губкинский
[6] Димитровград Дубна Дудинка Заречный Игарка
[11] Качканар Лесной Магадан Мегион Муравленко
[16] Невьянск Новоуральск Норильск Ноябрьск Обнинск
[21] Озерск Саров Снежинск Сосновый Бор Трехгорный
[26] Усть-Илимск
А вот готоров TO уже дофига - аж 671! (Вылезло форматирование? Ну да и хрен с ним!)
Абаза Абакан Абинск Агидель Агрыз Азнакаево Азов Аксай Алагир Алапаевск Алатырь Алейск Александров Алексеевка Алексин Алзамай
Альметьевск Анадырь Анапа Ангарск Анжеро-Судженск Анива Апатиты АпшеронскАрамиль Аргун Ардатов Арзамас Армавир Арсеньев Арск Артем Артемовский Архангельск Асбест Астрахань Ахтубинск Ачинск Аша Бавлы Байкальск Баксан Балабаново Балаково Балахна Балашиха Балашов Балтийск Барабинск Барнаул Батайск Белая Калитва Белая Холуница Белгород Белебей Белинский Белово Белогорск Белокуриха Белорецк Белореченск Белоусово Белоярский Бердск Березники Березовский Бийск Бикин Билибино Биробиджан Бирск Благовещенск Богданович Богородск Богучар Бодайбо Бокситогорск Бологое Бор Борзя Борисоглебск Боровичи Боровск Бородино Братск Брянск Бугульма Бугуруслан Бузулук Буй Валуйки Великие Луки Великий Новгород Верхний Тагил Верхний Уфалей Верхняя Пышма Верхняя Салда Верхотурье Видное Вилючинск Вихоревка Вичуга Владивосток Владикавказ Владимир Волгоград Волгодонск Волжск Волжский Вологда Волоколамск Волосово Вольск Воркута Воронеж Воскресенск Воткинск Всеволожск Выборг Выкса Вышний Волочек Вяземский Вязники Вязьма
Вятские Поляны Гаджиево Гатчина Гвардейск Гдов Геленджик Георгиевск Глазов Голицыно Горно-Алтайск Горнозаводск Городец Горячий Ключ Грозный Грязи Губаха Губкин Губкинский Гудермес Гуково Гулькевичи Гусев Гусь-Хрустальный
Давлеканово Дальнереченск Данилов Дегтярск Дедовск Десногорск Дзержинск
Дзержинский Дивногорск Димитровград Дмитров Добрянка Долгопрудный Домодедово Донской Дубна Дудинка Дюртюли Егорьевск Ейск Екатеринбург
Елабуга Елец Енисейск Ермолино Ессентуки Ефремов Железноводск Железногорск Железногорск-Илимский Железнодорожный Жигулевск Жирновск Жуков Жуковка Жуковский Заводоуковск Заволжск Заволжье Заинск Закаменск Заозерск Заполярный ЗаречныйЗаринск Звенигород Зверево Зеленогорск Зеленоград Зеленоградск Зеленодольск
Зерноград Зея Зима Златоуст Знаменск Ивангород Иваново Ивантеевка Ивдель Игарка Ижевск Избербаш ] Иланский Инсар Ирбит Иркутск Исилькуль Искитим Истра Ишим Ишимбай Йошкар-Ола Казань Калининград Калуга Каменск-Уральский Каменск-Шахтинский Камень-на-Оби Камышин Камышлов Канаш Кандалакша Канск Карабаново Карталы Касимов Касли Каспийск Катав-Ивановск Качканар Кашира Кедровый Кемерово Кимовск Кимры Кингисепп Кинель Кинешма Киржач Кириши Киров Кировград Кирово-Чепецк Кировск
Кисловодск Климовск Клин Ковров Когалым Кодинск Коломна Колпашево Кольчугино Комсомольск-на-Амуре Конаково Кондрово Копейск Кореновск Коркино Королёв Коряжма Костомукша Кострома Котельники Котельнич Котлас КотовоКрасавино Красновишерск Красногорск Краснодар Краснознаменск Краснокаменск Краснокамск Краснотурьинск Красноуральск Красноуфимск Красноярск Красный Кут Красный Сулин Кронштадт Кропоткин Крымск Кстово Кубинка Кудымкар Кузнецк Куйбышев Кумертау Кунгур Курган Курск Куртамыш Курчатов Кушва Кызыл Кыштым Кяхта Лабинск Лабытнанги Лаишево Лангепас Лениногорск Ленинск-Кузнецкий
Ленск Лермонтов Лесной Лесозаводск Лесосибирск Ливны Липецк Лиски Лихославль Лобня Лодейное Поле Ломоносов Луга Лысково Лысьва Лыткарино Люберцы Лянтор Магадан Магнитогорск Майкоп Малая Вишера Малоархангельск Малоярославец Мариинский Посад Маркс Махачкала Мегион Медвежьегорск Медногорск Медынь Междуреченск Меленки Мелеуз Менделеевск Миасс Микунь МиллеровоМинеральные Воды Минусинск Мирный Михайловка Мичуринск Можайск Можга Моздок Мончегорск Морозовск Москва Муравленко Мураши Мурманск Муром Мценск Мыски Мытищи Набережные Челны Надым Назарово Назрань Нальчик Наро-Фоминск Находка Невель Невинномысск Невьянск Неман Нерюнгри Нефтекамск Нефтекумск Нефтеюганск Нижневартовск
Нижнекамск Нижнеудинск Нижние Серги Нижний Новгород Нижний Тагил Нижняя Салда
Нижняя Тура Николаевск Новая Ляля Новоалтайск Нововоронеж Новокубанск Новокузнецк Новокуйбышевск Новомосковск Новопавловск Новороссийск Новосибирск Новотроицк Новоуральск Новочебоксарск Новочеркасск Новошахтинск Новый Оскол Новый Уренгой Ногинск Норильск Ноябрьск Нурлат Нытва Нягань Обнинск Одинцово Озерск Озёры Октябрьск Октябрьский Окуловка Оленегорск Омск Омутнинск Орёл ОренбургОрехово-Зуево Орск Оса Осинники Осташков Остров Острогожск Отрадное Отрадный Павлово Павловск Павловский Посад Партизанск Певек Пенза ПервоуральскПеревоз Пересвет Переславль-Залесский Пермь Петергоф Петрозаводск Петропавловск-Камчатски
Петушки Печора Пионерский Подольск Подпорожье Покров Полевской Полярные Зори Полярный Почеп Правдинск Приморско-Ахтарск Приозерск Прокопьевск Протвино Прохладный Псков Пугачев Пурпе Пушкин Пушкино Пущино Пыть-ЯхПятигорск Радужный Раменское Ревда Реж Реутов Ржев Рославль Россошь Ростов Ростов-на-Дону Рошаль Рубцовск Руза Рузаевка Рыбинск Рыбное Рыльск Рязань Салават Салехард Самара Санкт-Петербург Саранск Сарапул Саратов Саров Сасово Сатка Сафоново Саяногорск Саянск Светлогорск Северобайкальск Северодвинск Североморск Североуральск Северск Сельцо Семенов Семикаракорск Семилуки Серафимович Сергач Сергиев Посад Сердобск
Серов Серпухов Сибай Сим Сковородино Скопин Славгород Славянск-на-КубаниСлободской Слюдянка Смоленск Снежинск Снежногорск Советск Советский Сокол Соликамск Солнечногорск Соль-Илецк Сорочинск Сосенский Сосновоборск Сосновый Бор Сочи Спасск-Дальний Среднеуральск Ставрополь Старая Купавна Старая Русса Старица Старый Оскол Стерлитамак Стрежевой Строитель Струнино Ступино Суворов Суздаль Сургут Сурск Сусуман Сухиничи Сухой Лог Сызрань Сыктывкар ТаганрогТайга Тайшет Талица Тамбов Тарко-Сале Таруса Татарск Тверь Темрюк Тимашевск Тихвин Тихорецк Тобольск Тольятти Томск ТоржокТрехгорный Троицк Туапсе Туймазы Тула Тулун Тутаев Тюмень Углич Удачный Удомля Ужур Узловая Улан-Удэ Ульяновск Усинск Усолье-Сибирское Уссурийск Усть-Илимск Усть-Катав Усть-Кут Усть-Лабинск Уфа Ухта Учалы Фокино Фролово Фрязино Фурманов Хабаровск Ханты-Мансийск
Харабали Хилок Химки Хотьково Цивильск Цимлянск Чайковский ЧапаевскЧебаркуль Чебоксары Челябинск Черемхово Череповец Черкесск Черноголовка
Черногорск Чернушка Черняховск Чехов Чистополь Чита Чкаловск ЧудовоЧусовой Чухлома Шадринск Шали Шарыпово Шахты Шебекино Шелехов Шумерля Шушенское Шуя Щекино Щёлково Щербинка Щучье Электрогорск Электросталь Электроугли Элиста Энгельс Юбилейный Югорск Южно-Сахалинск Южноуральск Юрга Юрюзань Юхнов Якутск Ялуторовск Янаул ЯранскЯрославль Ярцево Ясный Яхрома Большой Камень
С городами разобрались. Теперь об оставшихся колонках таблички.
MIGR - миграционный поток
TO_POP - население
FROM|TO_LAT и FROM|TO_LON - широта и долгота городов откуда и куда
FROM|TO_OKTMO - коды ОКТМО. В принципе, я на них, как и на долготу с шириной, наверное, забью
SAME_REG - находятся ли города в одном регионе
ADM1 - является ли "куда" столицей региона
ADM2 - является ли "куда" градом стольным (т.е. Москва и СПБ)
/* Таки да, запрос unique(f$TO[f$ADM1==1]) выдает Москву и Питер */
Моя задача сводится к тому, чтобы обработать эти данные - построив для каждого из горо
дов FROM графики зависимости предсказанного потока от реального и найти среди них
города-выбросы, которые не ложатся на линию линейной регрессии.
Существует несколько способов находить выбросы, о них я расскажу в дальнейшем.
А пока вдохновимся песенкой о юном натурастите, который отчаянно хотел знать, как
устроена лягушка, и препарируем мой код.
Алгоритм весьма прост:
- загружаем данные
- делаем проекцию по количеству миграций
среди них берем выборку по уникальным город
- упражняемся с полученными данными
Шапочка весьма проста:
f <- read.csv("C:/Users/MSI/Desktop/Diploma/raw_data_v71.csv",header = TRUE, sep=";", quote="\"", encoding = "UTF-8"); #download data
names(f) <- c("FROM","TO", "MIGR", "TO_POP", "FROM_LAT", "FROM_LON", "TO_LAT","TO_LON" , "DISTANCE", "FR_OKTMO", "TO_OKTMO", "SAME_REG", "ADM2", "ADM1", "M2PRICE"); #Это специальная строчка во избежание вот такой, #например, херни в имени колонки - F.UFFU.X.FROM . Откуда это возникает я хз, #проще переименовать колонки)
MigrMore5 <- f[f$MIGR >5, ];
citiesMore5 <- unique(MigrMore5$FROM) ;
citiesMore5 <- as.vector(citiesMore5); #why not? (просто так компилятор меньше #говнится;) )
varlen <- rep(0, length(citiesMore5)); #summary output length of all variences #Я не помню, нахрена мне эта строчка....
.........
Пока я не перешла, собственно, к самой отвественной части программы - обсчете каждого города, скажу пару слов про алгоритм построения регрессионной модели и тех данных, которые будут обсчитываться алгоритмами в последствии.
Линейная модель строится как линейно зависимая функция от нескольких независимых переменных. Конкретно сейчас и конкретно в этом коде, я буду строить линейную модель зависимости потока миграций от населения, расстояния и цены за жильё. Про линейную регрессию можно прочитать здесь ( http://www.machinelearning.ru/wiki/index.php?title=Линейная_регрессия_(пример) ) , здесь ( https://ru.wikipedia.org/wiki/Линейная_регрессия ) или здесь (https://en.wikipedia.org/wiki/Linear_regression).
В языке R для этого есть функция lm(predictor ~ regressor1 + regressor2 + ...). Возвращает она объект типа list, который включает в себя коэффициенты, остатки, аппроксимированные значения и многое другое. Посмотреть хелп по этой функции можно в Rstudio набрав запрос здесть (два окна поиска)
или набрав в командной строке любой из этих трех вариантов:
>help(lm)
>??lm
>?lm
Я брала только коэффиценты, потому что с остатками всё получалось слишком хорошо. Брала свободный член, коэффиценты для каждого из параметров, вектор-столбец каждого из параметров и у меня получались предсказанные значения.
Кстати, реализовывать это можно, конечно, и через for(...){....}, а можно просто одной строчкой - благо R это почти как матлаб, только удобнее ;) В коде приведу оба варианта, просто один закоменчу. Вернемся к скрипту... ..............varlen <- rep(0, length(citiesMore5)); #summary output length of all variences
for (i in 1:length(citiesMore5))
{
#удаляем лишние строки из фрейма, поэтому конструкция ", ]", символизирующая #отсуствие условий на столбцы - обязательна:
City <- MigrMore5[MigrMore5$FROM == citiesMore5[i], ];
#логарифмирование поможет сделать данные более однородными.
lnMIGR <- log(City$MIGR);
lnDISTANCE <- log(City$DISTANCE);
lnPOP <- log(City$TO_POP);
lnPRICE <- log(City$M2PRICE);
#Строим линейную модель
LinearModel <- lm(lnMIGR ~ lnDISTANCE + lnPOP + lnPRICE );
coef<- LinearModel[[1]]; #присваиваем все коэффициенты Л.Р. в отдельную переменную
print(coef);
residuals <- LinearModel[[2]]; #присваиваем все остатки Л.Р.
dimention<-dim(City); #Находим размерность
deep <-dimention[1]; #количество строк
varlen[i]<-deep;
if( any(is.na(coef[1:4]))) #если встречается хотя бы одно NA в коэффициентах - шлём
{
print(c(i, "GTFO! Grub a better city!"))
next;
}
countRes <- rep (0, deep); #'cause my compiler say me "object countRes not found...
prediction <- rep (0, deep);
lnMIGR <- as.vector(lnMIGR); #why not? #Да-да, считаем предиктор и остатки
#--------next---may---be---counted---in----for-cycle-----------------------------#
prediction <-coef[1] + coef[2]*lnDISTANCE + coef[3]*lnPOP + coef[4]*lnPRICE;
countRes <- prediction - lnMIGR;
#-----------------------------------------------------------------------------------#
abscountRes<- abs(countRes); #cause I can but I don't use them
print(c(i, "passed"));
Отлично! Теперь у нас есть данные, которые мы уже можем обрабатывать!
Это вектор-столбцы prediction, countRes, lnMIGR и столбцы дата-фрейма City$To, City$FROM. Запихнём всё это в новый фрейм!
TO <- City$TO; FROM <-City$FROM;
CityData <- data.frame(FROM,TO, lnMIGR, prediction, countRes);
А вот теперь пошла магия различных методов.... Чтобы методы сработали, у нас сумма остатков должна равняться нулю. Ну, поскольку такого нифига не бывает в реальной жизни, то она должна быть близка к нулю. Я тестила. Там можно ооооочень много нулей написать ;) Поэтому условие на возможность прогона подготовленных данных статметодами следующее
countres <- CityData$countRes; #шоб было, ей богу!
if (sum(countres)>0.000001)
{
print(c(i, "BAD"));
next;
}
Метод №0. Всё, что меньше 25% квантиля и больше 75% квантиля - выброс. Метод дает достаточно много выбросов и жрет всё, как голодный мужчина после работы ;)#quantile() - высчитывает квантили. Выдает проименованный вектор 1:5. #Разыменовывать - моя прихоть
quantile <- unname(quantile(countres));
quantile25 <- quantile[2];
quantile75 <- quantile[4];
outliers3 <- countres[countres < quantile25]
outliers4 <- countres[countres > quantile75]
alloutliers2 <- c(outliers3,outliers4 ) #concatinate 2 numeric vectors
print("countres are"); print(sort(countres));
#находим индексы всех выбросов в остатках
index <- which( countres %in% alloutliers2);
print(c(index, "index"));
Migr <- lnMIGR[index]; #находим миграции, соотвествующие этим индексам
PreOut <- prediction[index]; #находим предсказания соотвествующие индексам
Надо бы это всё попроще и получще переписать и распихать по разным файлам....
IQR - метод. В принципе, всё тоже самое, только высчитывается доверительный интервал. Выбросов выявляет существенно меньше.
quantile <- unname(quantile(countres));
quantile25 <- quantile[2];
quantile75 <- quantile[4];
cond1 <- quantile25 - 1.5*(quantile75-quantile25);
print(cond1);
cond2 <- quantile75 + 1.5*(quantile75-quantile25);
print(cond2);
outliers <- countres[countres < cond1]
outliers2 <- countres[countres > cond2]
alloutliers <- c(outliers,outliers2 )
print("countres are"); print(sort(countres));
pos <- which( countres %in% alloutliers);
MigrOutlier <- lnMIGR[pos];
predictedOutlier <- prediction[pos];
Для методов 0 и 1 запись в файл:
#-----print---results----into---png-file---------------------------------------# adress<- paste0 ("C:\\Users\\MSI\\Desktop\\Diploma\\", i, ".png");
print(adress)
png(adress);
plot(lnMIGR, prediction, xlab = "Real migration", ylab = "Predicted migration");
title(as.character(citiesMore5[i]));
abline(0,1);
points(Migr, PreOut, type = "p", pch = 19, col = 'green');
dev.off()
#--------------------------end of printing in file------
Тест диксона. Работает только на малых выборках. Как работет, описано здесь
( http://sebastianraschka.com/Articles/2014_dixon_test.html ) и здесь (https://en.wikipedia.org/wiki/Dixon%27s_Q_test). Таблицы для больших значений - легко гуглятся, но она у меня и в коде забита.
Кстати... Возможно, у меня в коде ошибка есть. Там же надо брать выборку и лишь одно из значений - кандидат в выбросы. А у меня ж их дофига после нулевого метода. Стоит ли мне исключать из выборки остальные кандидаты в выбросы, кроме исследуемого?
if (deep <30)
{
#----avalible--countres----alloutliers2---index----Migr----PreOut--------#
isoutlier = 0; #flag
outlier =0;
for( k in 1:length(index))
{
candidate <- alloutliers2[index[k]]; #may be it's outlier?
res<-sort(countres); #sort residuals
l = k+1;
while(res[l]==res[k]) #if the outlier is not unique)
{
l=l+1;
}
gap = res[l]-res[k]; #-candidate
range = max(res)-min(res);
Q = gap/range;
Qt <- c(0.970, 0.970, 0.970, 0.829, 0.710, 0.625, 0.568, 0.526, 0.493, 0.466, 0.444, 0.426, 0.410, 0.396, 0.384, 0.374, 0.365, 0.356, 0.349, 0.342, 0.337 , 0.331 , 0.326 , 0.321 , 0.317, 0.312, 0.308, 0.305, 0.301, 0.298 );
if (Q > Qt[length(res)])
{
outlier <- c(outlier, candidate);
isoutlier =1;
print("is outlier");
print (candidate);
}
}
outlier <- outlier[outlier != 0] #remove 0-elemets
indexD <- which( countres %in% outlier);
Migr <- lnMIGR[indexD];
PreOut <- prediction[indexD];
Отдельная тема про автоматическую запись в файл. Можно записывать в разные форматы: .pdf, .png, .tif, .jpg, .bmp. Делается это примерно одинаково:
Открывается пусть, дается команда на запись файла, выполняются графические команды, запись закрывается.
adress<- paste0 ("C:\\Users\\MSI\\Desktop\\Dixon's test\\", i, ".png");
print(adress)
png(adress);
plot(lnMIGR, prediction, xlab = "Real migration", ylab = "Predicted migration");
title(paste(as.character(citiesMore5[i])), "Тест Диксона. N=", length(res));
abline(0,1); #строится график вида y=x
points(Migr, PreOut, type = "p", pch = 5, col = 'red');
dev.off()Надо бы ещё приточить легенду.
Но мне леееееень. Я марал. Ну, то есть олень. И вообще пошла я спать.
А ещё надо приточить вывод городов . Скорее всего стартанет это примерно так:
print(CityData$TO[indexD]).
Ой, таки я была права. Сработало. Компилятор выдал:
print(CityData$TO[index])
[1] Абакан Анапа Белгород Волгоград Иркутск
[6] Канск Кемерово Киров Краснодар Красноярск
[11] Новокузнецк Новороссийск Новосибирск Пермь Санкт-Петербург
[16] Сочи Сургут Тюмень Улан-Удэ Уфа
[21] Хабаровск Чита
671 Levels: Абаза Абакан Абинск Агидель Агрыз Азнакаево Азов Аксай Алагир ... Яхрома
> CityData$TO
[1] Абакан Анапа Ангарск Барнаул Белгород
[6] Братск Владивосток Волгоград Воронеж Екатеринбург
[11] Иваново Иркутск Казань Калининград Канск
[16] Кемерово Киров Краснодар Красноярск Липецк
[21] Магнитогорск Москва Нижний Новгород Новокузнецк Новороссийск
[26] Новосибирск Омск Пермь Ростов-на-Дону Самара
[31] Санкт-Петербург Сочи Ставрополь Сургут Сыктывкар
[36] Томск Тюмень Улан-Удэ Уфа Хабаровск
[41] Челябинск Чита Шелехов Ярославль
Держите котика, господа читатели:
