Портал аспирантов
 

Вернуться   Портал аспирантов > Компьютер для аспирантов > Software (программное обеспечение)

Ответ
 
Опции темы
Старый 11.11.2012, 15:55   #1
Hogfather
Platinum Member
 
Регистрация: 22.07.2010
Сообщений: 3,023
По умолчанию GNU R: Вопросы и ответы

/
Внимание!
Автор темы совместно с Администрацией портала просит писать в эту тему только относящееся к GNU R. Благодарности, разговоры о погоде будут безжалостно удаляться. Все "чмоки", пожалуйста, во флейме. Вопросы по R прошу задавать в личных сообщениях, чтобы тема была удобна для восприятия. Надеюсь на понимание.


Итак, поскольку я показал как пользоваться 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
LT<-c(rep(1,9),rep(2,267),rep(3,2843),rep(4,5450),rep(5,6564),rep(6,7044),rep(7,7518),rep(8,7071),rep(9,5620),rep(10,4016),rep(11,2545),rep(12,1494),rep(13,854),rep(14,416),rep(15,214),rep(16,122),rep(17,53),rep(18,16),rep(19,7),rep(20,2),21,22)
summary(LT)
length(LT)
# строим график
old.par <- par(mfrow=c(2,2))
hist(LT,main="Распределение букв",ylab="Число наблюдений",xlab="Число букв")
hist(LT,freq=F,ylab="Вероятность",xlab="Число букв",main="Распределение букв")
plot(ecdf(LT),verticals=T,main="График функции распределения")
boxplot(LT,main="Диаграмма Ящик-с-Усами",xlab="Число букв",horizontal=T)
par(old.par)
Вот что получится в результате:


Q: Замечательно, получилось. А как теперь скопировать картинку?

A: В MS Windows -- правая кнопка мыши на рисунке и выбрать нужный Вам вариант. Более продвинутые способы рассмотрим позже.

Q: А как бы мне открыто новую картинку на экране, да не поверх предыдущей?

A: Воспользуйтесь командой windows()
Код:
>windows()
Можно задать размеры окна в дюймах и еще кучу всего. Подробнее в справке.
Код:
>?windows

Q: У меня не подключилась библиотека MASS.

A: Понимаю, нужно её установить. В меню "Пакеты" выберите "Установить пакеты", дальше разберетесь сами.

Q: Я ввел неверную команду или хочу повторить предыдущую, опять заново набивать надо?

A: Нет. Достаточно нажать стрелку вверх. Это позволит "путешествовать" по истории команд, а если проскочили лишнее, то стрелка вниз. Далее можно редактировать строку. Потом нажмите Enter и она выполнится. Либо наберите команду history(), в отдельном окне покажется список последних введенных команд (25 штук). Если нужно больше, то поставьте соответствующую цифру в скобках, например history(100). Из открывшегося окна можно копировать нужные команды.

Последний раз редактировалось Hogfather; 06.12.2012 в 08:32.
---------
"So Long, and Thanks for all the Fish"
Hogfather вне форума   Ответить с цитированием
Реклама
Старый 12.11.2012, 15:27   #2
Hogfather
Platinum Member
 
Регистрация: 22.07.2010
Сообщений: 3,023
По умолчанию Вывод в файл

Q; Как корректно вывести графики с русскими буквами в PDF (SVG)?

A: Путем долгих экспериментов я пришел к выводу, что для меня лучше всего использовать библиотеку cairoDevice, которая использует библиотеку GTK+ .

Пример (использованы данные из темы про параметрические методы в филологии)

Код:
> library(cairoDevice)
> #####  Формируем PDF
> Cairo_pdf("test.pdf", width = 7, height = 7, pointsize = 10)
> boxplot(LT,main="Диаграмма ящик-с усами")
> dev.off()
windows 
      2 
> ##### Формируем SVG
> Cairo_svg("test.svg", width = 7, height = 7, pointsize = 10)
> boxplot(LT,main="Диаграмма ящик-с усами")
> dev.off()
windows 
      2 

> ##### Формируем PNG
> Cairo_png("test.png", width = 7, height = 7, pointsize = 10)
> boxplot(LT,main="Диаграмма ящик-с усами")
> dev.off()
windows 
      2
Результаты во вложении. Недавно вышла новая библиотека Cairo 2012 года, проверил, тоже работает.

Q: Можно ли редактировать результирующий график?
A: Да. Достаточно вывести его в SVG и воспользоваться графическим редактором, который понимает это формат. Например, Inkscape.
Вложения
Тип файла: zip test.zip (18.6 Кб, 4 просмотров)

Последний раз редактировалось Hogfather; 12.11.2012 в 18:13.
---------
"So Long, and Thanks for all the Fish"
Hogfather вне форума   Ответить с цитированием
Старый 12.11.2012, 20:22   #3
Hogfather
Platinum Member
 
Регистрация: 22.07.2010
Сообщений: 3,023
По умолчанию Получение данных

Q: Как удобнее редактировать скрипт в R.
A: Меню Файл->Новый скрипт. В новом окне пишите команды подряд.
Любую команду можно запустить в окне интерпретатора, нажав Ctrl-R.
Если нужно запустить весь скрипт (отдельные строки), то выделите нужные строки и нажмите Ctrl-R. Выделение всего скрипта -- Ctrl-A.

Q: Как отредактировать данные.
A: Например, воспользоваться командой edit или fix. Команда edit не меняет исходные данные, а создает копию. Поэтому её следует использовать так
Код:
>LT1<-edit(LT)
Если нужно внести правки, без сохранения исходных данных, то тогда всё проще.
Код:
>fix(LT)
Q: Как загрузить данные в R.
Допустим, у вас имеются данные, которые нужно загрузить в R.
Существует несколько способов это сделать. Самые полезные.
1) Загрузить из csv файла
Код:
>Data2012<-read.csv("D://R//data2012.csv", header=TRUE, sep = ";", quote="\"", dec=".")
По команде читается файл по адресу D:/R/data2012.csv (Обратите внимание на двойные слеш-символы), первая строка файла содержит заголовки (английский, без пробелов), разделитель записей символ ";", символ кавычки ["], десятичная точка -- ".".

2) Получить из Clipboard
Код:
> xx<-read.table("clipboard",header=T)
Как это работает




3) Получить прямо из Интернет
Для этого нужна библиотека XML. Вот таким нехитрым способом скачиваем 5 таблицу из статьи в Википедии.
Код:
> library(XML)
> url<-"http://en.wikipedia.org/wiki/World_population"
> (tbls<-readHTMLTable(url,which=5))
   Rank Country / Territory    Population              Date
1     1     В*China[note 2] 1,354,210,000 November 12, 2012
2     2             В*India 1,210,193,422        March 2011
3     3     В*United States   314,750,000 November 12, 2012
4     4         В*Indonesia   238,400,000          May 2010
5     5            В*Brazil   197,268,000 November 12, 2012
6     6          В*Pakistan   181,201,000 November 12, 2012
7     7           В*Nigeria   170,123,740         July 2012
8     8        В*Bangladesh   161,083,804         July 2012
9     9            В*Russia   141,927,297   January 1, 2010
10   10             В*Japan   127,610,000       May 1, 2012
   В*% of world\npopulation Source
1                     19.2%   [73]
2                       17%   [74]
3                     4.46%   [75]
4                     3.32%   [76]
5                      2.8%   [77]
6                     2.57%   [78]
7                     2.41%   [79]
8                     2.28%   [80]
9                    2.013%   [81]
10                    1.81%   [82]
>

4) Ввести вручную

Код:
> mydata <- data.frame(age=numeric(0), gender=character(0), weight=numeric(0))
> mydata <- edit(mydata)
5) Получить из базы данных
Подробнее см. тут

Последний раз редактировалось Hogfather; 13.12.2012 в 17:05.
---------
"So Long, and Thanks for all the Fish"
Hogfather вне форума   Ответить с цитированием
Старый 18.11.2012, 16:42   #4
SergeyMastitsky
Newbie
 
Регистрация: 18.11.2012
Сообщений: 1
По умолчанию

Коллеги, приглашаю посетить мой блог "R: Анализ и визуализация данных": r-analytics.blogspot.com
Надеюсь, он окажется кому-то полезным в поисках ответов на вопросы, касающиеся R.
SergeyMastitsky вне форума   Ответить с цитированием
Старый 18.11.2012, 18:25   #5
Hogfather
Platinum Member
 
Регистрация: 22.07.2010
Сообщений: 3,023
По умолчанию

Цитата:
Сообщение от SergeyMastitsky Посмотреть сообщение
приглашаю посетить мой блог "R: Анализ и визуализация данных": r-analytics.blogspot.com
Подтверждаю. Блог полезный, кошерный, рекомендуется к посещению и вдумчивому прочтению.
---------
"So Long, and Thanks for all the Fish"
Hogfather вне форума   Ответить с цитированием
Старый 18.11.2012, 19:59   #6
will
Platinum Member
 
Аватар для will
 
Регистрация: 17.09.2011
Сообщений: 2,770
По умолчанию

Вопрос
Есть ли в R
1)Фурье-анализ
2)подбор нелинейного уравнения (из фиксированного перечня+пользовательск.) для данного набора 2-х или 3-х мерных данных , а потом отобразить графически исходные данные и аппроксимирующую кривую \поверхность + повращать (для трехмерной графики)?
will вне форума   Ответить с цитированием
Старый 18.11.2012, 20:28   #7
Hogfather
Platinum Member
 
Регистрация: 22.07.2010
Сообщений: 3,023
По умолчанию

Цитата:
Сообщение от will Посмотреть сообщение
Есть ли в R
1)Фурье-анализ
Конечно, входит. Для быстрого преобразования Фурье даже пакетов подключать дополнительно не надо (пакет stats)
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
> fft(x)
[1] 10+0i -2+2i -2+0i -2-2i
> fft(fft(x), inverse = TRUE)/length(x)
[1] 1+0i 2+0i 3+0i 4+0i
Цитата:
Сообщение от will Посмотреть сообщение
2)подбор нелинейного уравнения (из фиксированного перечня+пользовательск.) для данного набора 2-х или 3-х мерных данных , а потом отобразить графически исходные данные и аппроксимирующую кривую \поверхность + повращать (для трехмерной графики)?
Естественно есть, причем разными методами. Я уже рассматривал нелинейный МНК. Трехмерные графики тоже есть. Вот тут есть разбор пакетов.

Классический эффектный пример нормального распределения в двух измерениях [1].

Код:
mu1<-0 # setting the expected value of x1
mu2<-0 # setting the expected value of x2
s11<-10 # setting the variance of x1
s12<-15 # setting the covariance between x1 and x2
s22<-10 # setting the variance of x2
rho<-0.5 # setting the correlation coefficient between x1 and x2
x1<-seq(-10,10,length=41) # generating the vector series x1
x2<-x1 # copying x1 to x2
#
f<-function(x1,x2)
{
term1<-1/(2*pi*sqrt(s11*s22*(1-rho^2)))
term2<--1/(2*(1-rho^2))
term3<-(x1-mu1)^2/s11
term4<-(x2-mu2)^2/s22
term5<--2*rho*((x1-mu1)*(x2-mu2))/(sqrt(s11)*sqrt(s22))
term1*exp(term2*(term3+term4-term5))
} # setting up the function of the multivariate normal density
#
z<-outer(x1,x2,f) # calculating the density values
#
persp(x1, x2, z,
main="Two dimensional Normal Distribution",
sub=expression(italic(f)~(bold(x))==frac(1,2~pi~sqrt(sigma[11]~
sigma[22]~(1-rho^2)))~phantom(0)^bold(.)~exp~bgroup("{",
list(-frac(1,2(1-rho^2)),
bgroup("[", frac((x[1]~-~mu[1])^2, sigma[11])~-~2~rho~frac(x[1]~-~mu[1],
sqrt(sigma[11]))~ frac(x[2]~-~mu[2],sqrt(sigma[22]))~+~
frac((x[2]~-~mu[2])^2, sigma[22]),"]")),"}")),
col="lightgreen",
theta=30, phi=20,
r=50,
d=0.1,
expand=0.5,
ltheta=90, lphi=180,
shade=0.75,
ticktype="detailed",
nticks=5) # produces the 3-D plot
#
mtext(expression(list(mu[1]==0,mu[2]==0,sigma[11]==10,sigma[22]==10,sigma[12
]==15,rho==0.5)), side=3) # adding a text line to the graph
---------
"So Long, and Thanks for all the Fish"
Hogfather вне форума   Ответить с цитированием
Старый 21.11.2012, 08:43   #8
Hogfather
Platinum Member
 
Регистрация: 22.07.2010
Сообщений: 3,023
По умолчанию Допиливаем пакет fitdistrplus

Рассмотренная ранее диаграмма Каллена и Фрея (Cullen and Frey graph) всем хороша, но она, во-первых, на языке эвентуального противника, во-вторых, нет возможности менять заголовок. Такое безобразие долго продолжаться не могло и, в процессе работы над статьёй, пришлось подправить функцию и русифицировать её вывод. Для этого идём на CRAN, скачиваем "Package source", распаковываем, в папке R находим нужную функцию и редактируем.

Результат во вложении. Пользуйтесь на здоровье.

Если работаете в стандартной оболочке R, то достаточно распаковать архив, в меню "Файл" выбрать "Загрузить код R..." и дальше спокойно пользоваться оной функцией, например так, как показано ниже. Обратите внимание, параметр "main" меняет заголовок диаграммы. Этого нет в стандартном функционале descdistr.

Код:
> descdist(LT,boot=100,main="Диграмма Каллена и Фрея\n(Редакция Hogfather'а для Портала Аспирантов)")
summary statistics
------
min:  1   max:  22 
median:  7 
mean:  7.199839 
estimated sd:  2.628829 
estimated skewness:  0.519882 
estimated kurtosis:  3.143716 
>


Все хорошо на картинке, но буржуи из эксцесса тройку не вычитают. Вот и думаю: пилить дальше или и так неплохо.
Вложения
Тип файла: zip descdist.zip (2.9 Кб, 3 просмотров)
---------
"So Long, and Thanks for all the Fish"
Hogfather вне форума   Ответить с цитированием
Старый 21.11.2012, 22:25   #9
Hogfather
Platinum Member
 
Регистрация: 22.07.2010
Сообщений: 3,023
По умолчанию Обживаем 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 и радуемся жизни.
---------
"So Long, and Thanks for all the Fish"
Hogfather вне форума   Ответить с цитированием
Старый 01.12.2012, 21:20   #10
Hogfather
Platinum Member
 
Регистрация: 22.07.2010
Сообщений: 3,023
По умолчанию Множественная линейная регрессия

Наши постоянные читатели спрашивают
Цитата:
Будет замечательно, если Вы рассмотрите данную задачу как пример.

Задача: получить зависимость z(x,y)
Исходные данные:

x <- c(75, 75, 75, 75, 100, 125, 125, 150, 150, 175, 175, 200, 200, 225, 250, 250, 250, 275, 300, 300, 300, 300)
y <- c(375, 625, 1000, 1500, 500, 375, 1125, 750, 1250, 825, 1500, 375, 1000, 1125, 750, 1250, 1500, 1375, 375, 825, 1500, 1250)
z <- c(1.452, 3.705, 18.857, 50.434, 1.956, 2.272, 40.306, 3.368, 47.686, 4.464, 55.63, 2.611, 9.077, 15.304, 7.559, 27.257, 60.174, 35.965, 3.355, 22.491, 38.036, 54.237)
Ну, что ж. Начнем. Загоняем данные в R и делаем Data Frame (таблицу)
Код:
> x <- c(75, 75, 75, 75, 100, 125, 125, 150, 150, 175, 175, 200, 200, 225, 250, 250, 250, 275, 300, 300, 300, 300)
> y <- c(375, 625, 1000, 1500, 500, 375, 1125, 750, 1250, 825, 1500, 375, 1000, 1125, 750, 1250, 1500, 1375, 375, 825, 1500, 1250)
> z <- c(1.452, 3.705, 18.857, 50.434, 1.956, 2.272, 40.306, 3.368, 47.686, 4.464, 55.63, 2.611, 9.077, 15.304, 7.559, 27.257, 60.174, 35.965, 3.355, 22.491, 38.036, 54.237)
> myData<-data.frame(x,y,z)
Вторым нашим шагом будет полюбоваться на наши данные, чтобы увидеть зависимость, используя корреляционную матрицу.
Код:
> library(lattice)
> splom(~myData[,c("x","y","z")],cex.labels=1.2,varnames=c("x","y","z"))


Кого-как, а меня не впечатлило. Икс у нас явно не при делах, а Игрек на линейную модель никак не тянет. Напрашивается экспонента. Проверим эту нехитрую мысль.
Введем фиктивную переменную lnZ равную натуральному логарифму от Z. В R это делается так.
Код:
> myData$lnZ<-log(myData$z)
> splom(~myData[,c("x","y","lnZ")],cex.labels=1.2,varnames=c("x","y","ln(z)"),col="red")

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

Третий вариант мы можем получить логарифмированием всего и вся. В таком случае имеем

Код:
> 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))

Call:
lm(formula = z ~ x + y)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.717  -8.306   1.616   6.760  17.840 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -21.233747   6.918102  -3.069  0.00631 ** 
x             0.001910   0.028261   0.068  0.94683    
y             0.045646   0.005698   8.010 1.64e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 10.24 on 19 degrees of freedom
Multiple R-squared: 0.7846,	Adjusted R-squared: 0.762 
F-statistic: 34.61 on 2 and 19 DF,  p-value: 4.624e-07 

> summary(lm2<-lm(log(z)~x+y))

Call:
lm(formula = log(z) ~ x + y)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.61700 -0.30351 -0.03218  0.19850  0.81873 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.5939071  0.3067107  -1.936   0.0678 .  
x            0.0016541  0.0012529   1.320   0.2025    
y            0.0029021  0.0002526  11.488  5.4e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4538 on 19 degrees of freedom
Multiple R-squared: 0.889,	Adjusted R-squared: 0.8773 
F-statistic: 76.08 on 2 and 19 DF,  p-value: 8.53e-10 

> summary(lm3<-lm(log(z)~log(x)+log(y)))

Call:
lm(formula = log(z) ~ log(x) + log(y))

Residuals:
    Min      1Q  Median      3Q     Max 
-0.9323 -0.2865  0.1118  0.3488  0.6532 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -14.5313     1.6847  -8.626 5.38e-08 ***
log(x)        0.2550     0.2275   1.121    0.276    
log(y)        2.3263     0.2280  10.203 3.81e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.5058 on 19 degrees of freedom
Multiple R-squared: 0.8621,	Adjusted R-squared: 0.8476 
F-statistic: 59.41 on 2 and 19 DF,  p-value: 6.682e-09
Смотрим колонку Pr(>|t|) и количество звездочек сбоку. Начинаем модифицировать модель "a posteriori". Собственно говоря, Икс нам совершенно не нужен во всех его проявлениях. Непонятно зачем он вообще тут оказался. Возьмем вторую модель, она посимпатичнее и выкинем лишнее.
Код:
> summary(lm4<-lm(log(z)~y))

Call:
lm(formula = log(z) ~ y)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.66190 -0.35391 -0.08276  0.28610  1.01267 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.3660607  0.2582089  -1.418    0.172    
y            0.0029897  0.0002483  12.043 1.28e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4622 on 20 degrees of freedom
Multiple R-squared: 0.8788,	Adjusted R-squared: 0.8727 
F-statistic:   145 on 1 and 20 DF,  p-value: 1.277e-10
Смотрим дальше. Свободный член (не радуйтесь девушки, это не то, что вы подумали) "Intercept" нам тоже не нужен. Да и пёс с ним.
Код:
> summary(lm5<-lm(log(z)~-1+y))

Call:
lm(formula = log(z) ~ -1 + y)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.78397 -0.35742 -0.05766  0.18365  0.91500 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
y 2.664e-03  9.699e-05   27.47   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4732 on 21 degrees of freedom
Multiple R-squared: 0.9729,	Adjusted R-squared: 0.9716 
F-statistic: 754.6 on 1 and 21 DF,  p-value: < 2.2e-16
Обратите внимание какой у нас славный получился Эр квадрат. Дай бог всякому.
Моделька имеет вид
Едем дальше. Впереди дисперсионный анализ.
Код:
> anova(lm5)
Analysis of Variance Table

Response: log(z)
          Df  Sum Sq Mean Sq F value    Pr(>F)    
y          1 168.946 168.946  754.61 < 2.2e-16 ***
Residuals 21   4.702   0.224
И графическое отображение остатков модели.
Код:
> par(mfrow=c(2,2))
> plot(lm5)
> par(mfrow=c(1,1))
Ну, лично мне нравится.


Итак, у нас есть модель и остатки. Если все хорошо, они должны иметь нормальное распределение с матожиданием 0. Судя по картинке (график квантилей) мы должны получить неплохой результат в тесте не нормальность распределения остатков модели.
Делаем тест Шапиро-Уилка. Смотрим и балдеем.
Код:
> shapiro.test(lm5$residuals)

	Shapiro-Wilk normality test

data:  lm5$residuals 
W = 0.9628, p-value = 0.5488
Раз все так здорово, находим стандартное отклонение.
Код:
> sd(lm5$residuals)
[1] 0.4700086
> summary(lm5$residuals)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.78400 -0.35740 -0.05766 -0.05331  0.18360  0.91500
Sapienti sat.

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

P.S. Вот тут спрашивают, а как нарисовать 3D-график этого безобразия. Отвечаю. Просто. Например, так.
Код:
> library(scatterplot3d)
> scatterplot3d(x , y, z,highlight.3d=T,main="Вот такой график")

Последний раз редактировалось Hogfather; 05.12.2012 в 20:20.
---------
"So Long, and Thanks for all the Fish"
Hogfather вне форума   Ответить с цитированием
Ответ

Опции темы

Ваши права в разделе
Вы не можете создавать новые темы
Вы не можете отвечать в темах
Вы не можете прикреплять вложения
Вы не можете редактировать свои сообщения

BB коды Вкл.
Смайлы Вкл.
[IMG] код Вкл.
HTML код Выкл.



Текущее время: 05:55. Часовой пояс GMT +3.


Powered by vBulletin® Version 3.8.8
Copyright ©2000 - 2017, vBulletin Solutions, Inc. Перевод: zCarot
© 2001—2017, «Аспирантура. Портал аспирантов»
Рейтинг@Mail.ru