|
GNU R: Вопросы и ответы
[table 1 3 0]3^http://aspirantura.spb.ru/forum/pict...&pictureid=979 |
Внимание![BR]Автор темы совместно с Администрацией портала просит писать в эту тему только относящееся к GNU R. Благодарности, разговоры о погоде будут безжалостно удаляться. Все "чмоки", пожалуйста, во флейме. Вопросы по R прошу задавать в личных сообщениях, чтобы тема была удобна для восприятия. Надеюсь на понимание. [/table]Итак, поскольку я показал как пользоваться R для решения задачи подбора распределения, а добро не должно остаться безнаказанным, в личные сообщения посыпались вопросы. Чтобы два раза не вставать, ответы публикую сюда. Рассмотрено в других темах (пример кода) 1. Подгонка распределения. fitdistr из MASS 2. Подгонка распределения. fitdist из fitdistrplus. Лучше! 3. Нелинейный метод наименьших квадратов 4. Гамма-распределение. Вывод графика. 5. Диаграмма "Ящик с усами" для нескольких векторов 6. Упрощение себе жизни с помощью написания функции в R 7. Подгонка МНК Q: Где взять R? A: Лучше всего в специально предназначенном для этого месте, называется CRAN. Русское зеркало находится вот тут. Выбираете операционную систему, далее скачиваете и устанавливаете. Q:C чего начать? A: Почитать статьи вот тут Q: Как попробовать то, что показано в примерах? A: Очень просто. Запускаем R, копируем команду из примера в консоль интерпретатора. Команды начинаются с ">", сам символ ">" приглашения интерпретатора R не копируем. Например, скопируйте все строки из примера ниже (я оставил только команды из вышеприведенной заметки, а символ ">" стёр) и запустите в R (Не забудьте в конце нажать Enter, чтобы выполнилась последняя команда). Код:
# Тестовый пример. Скопировать все строки и вставить в окно интерпретатора R http://aspirantura.spb.ru/forum/pict...&pictureid=973 Q: Замечательно, получилось. А как теперь скопировать картинку? A: В MS Windows -- правая кнопка мыши на рисунке и выбрать нужный Вам вариант. Более продвинутые способы рассмотрим позже. Q: А как бы мне открыто новую картинку на экране, да не поверх предыдущей? A: Воспользуйтесь командой windows() Код:
>windows() Код:
>?windows Q: У меня не подключилась библиотека MASS. A: Понимаю, нужно её установить. В меню "Пакеты" выберите "Установить пакеты", дальше разберетесь сами. Q: Я ввел неверную команду или хочу повторить предыдущую, опять заново набивать надо? A: Нет. Достаточно нажать стрелку вверх. Это позволит "путешествовать" по истории команд, а если проскочили лишнее, то стрелка вниз. Далее можно редактировать строку. Потом нажмите Enter и она выполнится. Либо наберите команду history(), в отдельном окне покажется список последних введенных команд (25 штук). Если нужно больше, то поставьте соответствующую цифру в скобках, например history(100). Из открывшегося окна можно копировать нужные команды. |
Вывод в файл
Вложений: 1
Q; Как корректно вывести графики с русскими буквами в PDF (SVG)?
A: Путем долгих экспериментов я пришел к выводу, что для меня лучше всего использовать библиотеку cairoDevice, которая использует библиотеку GTK+ . Пример (использованы данные из темы про параметрические методы в филологии) Код:
> library(cairoDevice) Q: Можно ли редактировать результирующий график? A: Да. Достаточно вывести его в SVG и воспользоваться графическим редактором, который понимает это формат. Например, Inkscape. |
Получение данных
Q: Как удобнее редактировать скрипт в R.
A: Меню Файл->Новый скрипт. В новом окне пишите команды подряд. Любую команду можно запустить в окне интерпретатора, нажав Ctrl-R. Если нужно запустить весь скрипт (отдельные строки), то выделите нужные строки и нажмите Ctrl-R. Выделение всего скрипта -- Ctrl-A. Q: Как отредактировать данные. A: Например, воспользоваться командой edit или fix. Команда edit не меняет исходные данные, а создает копию. Поэтому её следует использовать так Код:
>LT1<-edit(LT) Код:
>fix(LT) Допустим, у вас имеются данные, которые нужно загрузить в R. Существует несколько способов это сделать. Самые полезные. 1) Загрузить из csv файла Код:
>Data2012<-read.csv("D://R//data2012.csv", header=TRUE, sep = ";", quote="\"", dec=".") 2) Получить из Clipboard Код:
> xx<-read.table("clipboard",header=T) Как это работает
3) Получить прямо из Интернет Для этого нужна библиотека XML. Вот таким нехитрым способом скачиваем 5 таблицу из статьи в Википедии. Код:
> library(XML) 4) Ввести вручную Код:
> mydata <- data.frame(age=numeric(0), gender=character(0), weight=numeric(0)) Подробнее см. тут |
Коллеги, приглашаю посетить мой блог "R: Анализ и визуализация данных": r-analytics.blogspot.com
Надеюсь, он окажется кому-то полезным в поисках ответов на вопросы, касающиеся R. |
Цитата:
|
Вопрос
Есть ли в R 1)Фурье-анализ 2)подбор нелинейного уравнения (из фиксированного перечня+пользовательск.) для данного набора 2-х или 3-х мерных данных , а потом отобразить графически исходные данные и аппроксимирующую кривую \поверхность + повращать (для трехмерной графики)? |
Цитата:
fft
fft {stats} R Documentation Fast Discrete Fourier Transform Description Performs the Fast Fourier Transform of an array. Usage fft(z, inverse = FALSE) mvfft(z, inverse = FALSE) Arguments z a real or complex array containing the values to be transformed. inverse if TRUE, the unnormalized inverse transform is computed (the inverse has a + in the exponent of e, but here, we do not divide by 1/length(x)). Value When z is a vector, the value computed and returned by fft is the unnormalized univariate Fourier transform of the sequence of values in z. When z contains an array, fft computes and returns the multivariate (spatial) transform. If inverse is TRUE, the (unnormalized) inverse Fourier transform is returned, i.e., if y <- fft(z), then z is fft(y, inverse = TRUE) / length(y). By contrast, mvfft takes a real or complex matrix as argument, and returns a similar shaped matrix, but with each column replaced by its discrete Fourier transform. This is useful for analyzing vector-valued series. The FFT is fastest when the length of the series being transformed is highly composite (i.e., has many factors). If this is not the case, the transform may take a long time to compute and will use a large amount of memory. Source Uses C translation of Fortran code in Singleton (1979). References Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole. Singleton, R. C. (1979) Mixed Radix Fast Fourier Transforms, in Programs for Digital Signal Processing, IEEE Digital Signal Processing Committee eds. IEEE Press. Код:
> x <- 1:4 Цитата:
Классический эффектный пример нормального распределения в двух измерениях [1]. Код:
mu1<-0 # setting the expected value of x1 |
Допиливаем пакет fitdistrplus
Вложений: 1
Рассмотренная ранее диаграмма Каллена и Фрея (Cullen and Frey graph) всем хороша, но она, во-первых, на языке эвентуального противника, во-вторых, нет возможности менять заголовок. Такое безобразие долго продолжаться не могло и, в процессе работы над статьёй, пришлось подправить функцию и русифицировать её вывод. Для этого идём на CRAN, скачиваем "Package source", распаковываем, в папке R находим нужную функцию и редактируем.
Результат во вложении. Пользуйтесь на здоровье. Если работаете в стандартной оболочке R, то достаточно распаковать архив, в меню "Файл" выбрать "Загрузить код R..." и дальше спокойно пользоваться оной функцией, например так, как показано ниже. Обратите внимание, параметр "main" меняет заголовок диаграммы. Этого нет в стандартном функционале descdistr. Код:
> descdist(LT,boot=100,main="Диграмма Каллена и Фрея\n(Редакция Hogfather'а для Портала Аспирантов)") Все хорошо на картинке, но буржуи из эксцесса тройку не вычитают. Вот и думаю: пилить дальше или и так неплохо. |
Обживаем RStudio
Прочитав про RStudio решил её обжить и делюсь решением первых проблем.
Проблемы с кодировкой Кодировка внутри файлов воспринимается весьма неплохо, если скрипт прочитался неверно, его можно прочитать в другой кодировке с помощью команды Файл->Reopen with Encoding... Но! Если вы на компьютере имеете логин на русском языке, и в папке Users (Пользователи) у Вас папка на русском языке, соответствующая логину, то в Win7 возникают проблемы. Главная, RStudio категорически не видит установленные пакеты, поскольку они устанавливаются в C:/Users/Логин/Documents/R/... Вместо вот этого "Логин", который у Вас вполне может быть "Милый суслик", программа пойдет в папку C:/Users/Логин/Documents/R/ (Логин = попытка прочитать в кодировке Windows-1251 UTF-8 слово "Логин") Поскольку ни считать её, ни создать она не может, у программы возникает печалька, и она не видит установленных пакетов, а также не может поставить новые. Вот такая ерунда. Выход есть! Нужно сделать ссылку типа "связь каталогов" с таким идиотским именем в папке Users Для этого (как вариант). 1. Скачиваем и устанавливаем Far. Хорошему грустному человеку он завсегда пригодится. 2. Заходим в каталог C:/Users/ 3. Нажимаем Alt-F6 и в меню выбираем тип ссылки "Связь каталогов" или "directory junction", зависит от языка в Far. 4. На противоположной панели появляется ссылка. 5. Переименовываем её в кракозябровое имя, которое мы подглядим в ошибке RStudio, выделим, скопируем, нажмем в Far F6 и вставим его в нужное место окна, дополнив строку пути нашими кракозябрами 6. Перезапускаем RStudio и радуемся жизни. |
Множественная линейная регрессия
Наши постоянные читатели спрашивают
Цитата:
Код:
> x <- c(75, 75, 75, 75, 100, 125, 125, 150, 150, 175, 175, 200, 200, 225, 250, 250, 250, 275, 300, 300, 300, 300) Код:
> library(lattice) Кого-как, а меня не впечатлило. Икс у нас явно не при делах, а Игрек на линейную модель никак не тянет. Напрашивается экспонента. Проверим эту нехитрую мысль. Введем фиктивную переменную lnZ равную натуральному логарифму от Z. В R это делается так. Код:
> myData$lnZ<-log(myData$z) Поясняю, мы от модель вида http://www.texify.com/img/%5Cnormals...%5Cepsilon.gif плавно перешли к http://www.texify.com/img/%5Cnormals...ilon%29%7D.gif. В данном случае эпсилон у нас изображает остатки модели, которые, в идеальном варианте, должны иметь нормальное распределение. Третий вариант мы можем получить логарифмированием всего и вся. В таком случае имеем http://www.texify.com/img/%5Cnormals...scr%7BE%7D.gif Код:
> splom(~log(myData[,c("x","y","z")]),cex.labels=1.2,varnames=c("ln(x)","ln(y)","ln(z)"),col="dark green") После такой работы уже можно и к моделированию переходить. Чисто чтоб поржать сделаем все три модели. Код:
> summary(lm1<-lm(z~x+y)) Код:
> summary(lm4<-lm(log(z)~y)) Код:
> summary(lm5<-lm(log(z)~-1+y)) Моделька имеет вид http://www.texify.com/img/%5Cnormals...scr%7BE%7D.gif Едем дальше. Впереди дисперсионный анализ. Код:
> anova(lm5) Код:
> par(mfrow=c(2,2)) http://aspirantura.spb.ru/forum/pict...pictureid=1003 Итак, у нас есть модель и остатки. Если все хорошо, они должны иметь нормальное распределение с матожиданием 0. Судя по картинке (график квантилей) мы должны получить неплохой результат в тесте не нормальность распределения остатков модели. Делаем тест Шапиро-Уилка. Смотрим и балдеем. Код:
> shapiro.test(lm5$residuals) Код:
> sd(lm5$residuals) Понятно, что в форумной заметке нельзя рассказать весь курс, но основные идеи я, надеюсь, рассказал. P.S. Вот тут спрашивают, а как нарисовать 3D-график этого безобразия. Отвечаю. Просто. Например, так. Код:
> library(scatterplot3d) |
Множественная регрессия. Краткий справочник.
Попробуем резюмировать.
Для подгонки модели МНК, используем код Код:
> fit <- lm(lnZ ~ x + y, data=myData) Код:
> coefficients(fit) # коэффициенты модели Код:
> fit1 <- lm(lnZ ~ y, data=myData) Автоматическая подгонка модели посредством исключения незначимых членов. Код:
> library(MASS) |
Краткий ликбез по моделям, или как читать то, что мы получаем в результате подгонки.
Хотя меня особо не спрашивают, но имея опыт работы со студентами постараюсь кратко и без лишнего занудства изложить, что же Пытливый Ум может увидеть в сборище букв на языке эвентуального противника, когда строит модель.
Рассмотрим пример, который был приведен выше. Код:
> summary(fit) F-statistic: 76.08 on 2 and 19 DF, p-value: 8.53e-10 В данном случае нулевая гипотеза выглядит так: H0: β1=β2=...βn=0, т.е. все коэффициенты кроме свободного члена равны 0. Это мы проводим проверку значимости модели. Если p-value у нас превышает 0.05, то дальнейшая работа с моделью просто бессмысленна, поскольку модель незначимая. Пример незначимой модели
Код:
> lm0<-lm(y~x,data=test) Второе, на что мы смотрим, это на остатки модели (Residuals). Код:
Residuals: Третье, на что мы смотрим -- коэффициенты. Код:
Coefficients: Дальше, Residual standard error: 1.924 on 18 degrees of freedom говорит нам о средней квадратической ошибке остатков (σ). Multiple R-squared: 0.889, Adjusted R-squared: 0.8773 говорит о коэффициентах детерминации и качестве модели. Чем они ближе к единице, тем лучше. В то же время, проверку адекватности модели следует произвести графически. Именно для этого все графики и строились.. Пример неадекватной модели
Код:
> test<-data.frame(x<-1:20) А теперь смотрите на график http://aspirantura.spb.ru/forum/pict...pictureid=1017 Что должно быть на графике: 1. Точки на графике Residuals vs Fitted должны быть случайно разбросаны, без явно выраженного тренда. Тут этого нет. 2. Точки на графике Normal Q-Q должны более-менее лежать на линии, показывая что остатки имеют нормальное распредение 3. Точки на двух остальных графиках должны быть в группах не слишком далеко от центра. Т.е. модель у нас неадекватна, не смотря на милый Эр Квадрат. Такие дела. Для сравнения посмотрите график на предыдущей странице Добавлено через 4 часа 3 минуты Использование метода Бокса-Кокса для преобразования регрессионной модели В предыдущем примере мы рассмотрели случай, когда модель выбрана неадекватной. Судя по параболе, которую мы наблюдаем, здесь необходимо какое-то степенное преобразование. Попробуем использовать преобразование Бокса-Кокса, которое входит в пакет MASS. Т.е. найдем такую λ, чтобы yλ сделал нам аж боже мой. Исходные данные у нас те же и та же модель. Код:
> test<-data.frame(x<-1:20) Код:
> library(MASS) http://aspirantura.spb.ru/forum/pict...pictureid=1018 К сожалению, конкретного результата, который нужно подставить у нас нигде нет. Но мы можем его получить. Функция which.max(bc$y) вернёт нам номер элемента, у которого Игрек наибольший, а bc$x[which.max(bc$y)]) -- собственно нужную лямбду. Код:
> (lambda<-bc$x[which.max(bc$y)]) Делаем преобразование и строим новую модель. Код:
> lm1<-lm(I(y^lambda)~x,data=test) Код:
> plot(lm1,which=1) Воот! Теперь совсем другое дело. Следует отметить особые случаи в преобразовании Бокса-Кокса. Если λ=1, то преобразование не требуется. Если λ=0, то следует использовать логарифмическое преобразование log(y). Пример. Код:
> test<-data.frame(x<-1:20) Скрипт для графика
Код:
> par(mfrow=c(1,2)) http://aspirantura.spb.ru/forum/pict...pictureid=1020 Добавлено через 18 часов 12 минут Кореллограммы В рассмотренной задаче (#10) рисовалась корреляционная матрица в виде соответствующих диаграмм рассеяния. Также можно получить корреляционную матрицу в виде коэффициентов корреляции в текстовом виде. Код:
> # Инициализация данных Код:
> symnum(cor(myData[,c("x","y","z")])) Код:
> symnum(cor(myData[,c("x","y","lnZ")])) Однако, есть и другой способ сделать это, да и к тому же красиво, нарисовав коррелограмму. Такой рисунок несказанно обогатит статью и поразит всех читателей. Для этого нужна установленная библиотека corrgram. Примеры вызовов функции приведены ниже. Существует 4 возможных отображения корреляции: panel.pts -- диаграмма разброса panel.pie -- круговая диаграмма, меняет цвет и заполненность в зависимость от коэффициента корреляции (отрицательная -- красная, положительная -- синяя). panel.ellipse -- рисует красоту, как показана на рисунке 3. В общем, границы и аппроксимацию. panel.shade -- закрашивает прямоугольник в зависимости от коэффициента корреляции. Их присваивают, соответственно lower.panel или upper.panel, в зависимости от того, где нужен рисунок: сверху или снизу. Также можно выводить название переменных и минимум/максимум значений, как это показано ниже. Код:
> Код:
> corrgram(myData[,c("x","y","lnZ")], lower.panel=panel.pts, Код:
> corrgram(log(myData[,c("x","y","z")]), order=TRUE, lower.panel=panel.shade, Понятно, что на трех переменных это не сильно показательно, но если их будет десяток, то картинка будет гораздо симпатичнее. |
Регресионные деревья
Чем хорош R, что в нём можно много всякого забавного сделать. например, захотелось поиграть в регрессионные деревья -- никто слова плохого не скажет.
Рассмотрим ещё раз нашу задачу. Код:
Результат, как видите следующий. Если y<912.5, то ln(z) равен 1.319, а если y>=912.5, то 3.499. Смысла для данной задачи никакого, разве что лишний раз убедились, что Икс у нас не при делах. Как это должно работать
В данном случае я немного модифицировал штатный пример из документации В качестве исходных данных берется таблица данных по кифозу у детей из 81 строк and 4 столбцов, описывающая детей у которых были корректирующие операции на позвоночнике Таблица содержит следующие колонки: Kyphosis - качественные данные (factor) с вариантами "absent" и "present", показывающая был ли кифоз после операции. Age- возраст в месяцах Number - число затронутых позвонков Start - номер верхнего позвонка, на котором происходила операция Код:
> # Первые строки таблицы |
Hogfather,
мне бывает нужен МНК на очень подробных сетках (1000\times 1000, например). Потянет это оболочка? (без ограничения производительности железа) |
Цитата:
Обратите внимание, я специально выводил время. 13 секунд на модель. Дольше таблицу инициализировал... Процесcор старый, Intel Core Duo T7250 2 ГГц. Озу 3 Гб Скрытый текст
Код:
> # Инициализируем матрицу 1000х1000 для x1..x1000, с матожидаением i |
Текущее время: 01:21. Часовой пояс GMT +3. |
|
Powered by vBulletin® Version 3.8.8
Copyright ©2000 - 2024, vBulletin Solutions, Inc. Перевод: zCarot
© 2001—2024, «Аспирантура. Портал аспирантов»