[R] нубские вопросы

2

>x[x<0]=0

1 http://rss.acs.unt.edu/Rdoc/library/TTR/html/MovingAverages.html

О как - знал же что просто...

Спасибо, добрый человек!
Вот только вопрос попутно - почему такая простая функция загнана в какую-то трейдерскую библиотеку?

Более менее:

stats::decompose        Classical Seasonal Decomposition by Moving
                        Averages

По первому:

> ??approximation

ЗЫ ниасилил и запилил аппроксимацию по методу наименьших квадратов

3 Из help(lm):


na.action: a function which indicates what should happen when the data
contain ‘NA’s. The default is set by the ‘na.action’ setting
of ‘options’, and is ‘na.fail’ if that is unset. The
‘factory-fresh’ default is ‘na.omit’. Another possible value
is ‘NULL’, no action. Value ‘na.exclude’ can be useful.

Да нет - NA здесь не причём - они нужны только для отсчёта Y при построении прямой линейной регрессии. Прямая нужна в виде вектора значений равного по длине вектору по которому она строится.

Вот это я как раз видел - сам находил, но так и не понял как использовать.

В принципе поставил TTR - всё работает на ура.

То есть просто сами считали «руками»?

Ну да - примерно это я и имел ввиду.

Спасибо, что помогаете. Необычный язык. Да и язык ли? С большим бы удовольствием юзал Tcl и в нём функции R, а не наоборот...

Я нихрена не понял.

Необычный язык.

Считай, что это такой необычный лисп, с инфиксным синтаксисом и неимоверным количеством сахара.

mean_my <- function(x) mean(window(1:20, start =x, end = x+3))

sapply(1:10, mean_my)

считается что настоящий .... ну и так далее :)

считается что хватает window

шутка :)

на самом деле все называется ?filter как более общий случай

Про window не знал, про filter читал, но до сих пор не понял как это работает.

Не могу понять не только прицип но даже частный случай.

Допустим имеем:

x

[1] 1 2 3 4 5 6 7 8 9 10

filter(x, rep(1,3))

Time Series:
Start = 1
End = 10
Frequency = 1
[1] NA 6 9 12 15 18 21 24 27 NA

Согласну функции из мана (y[i] = f[1]*x[i+o] + … + f[p]*x[i+o-(p-1)]) третий элемент должен быть:

1*x[3]+1*x[3-(2-1)]+1*x[3-(3-1)] == 1*3+1*2+1*1 == 6

А он у нас 9. А откуда 2 элемент?

> filter(1:10,c(1,1,1), method="convolution",sides=1)
Time Series:
Start = 1 
End = 10 
Frequency = 1 
 [1] NA NA  6  9 12 15 18 21 24 27
> filter(1:10,c(1,1,1), method="convolution",sides=2)
Time Series:
Start = 1 
End = 10 
Frequency = 1 
 [1] NA  6  9 12 15 18 21 24 27 NA

в первом случае первые два NA следствие неполноты окна в котором свертка делается, третье значение сумма первых трех значений первоначального ряда 6=1+2+3; следующее значение 9=2+3+4

далее меняем окно на двустороннее...

The convolution filter is

y = f[1]*x[i+o] + ... + f

*x[i+o-(p-1)] where ‘o’ is the offset: see ‘sides’ for how it is determined.

sides: for convolution filters only. If ‘sides=1’ the filter coefficients are for past values only; if ‘sides=2’ they are centred around lag 0. In this case the length of the filter should be odd, but if it is even, more of the filter is forward in time than backward.

То есть при sides=1 offset устанавливается равным 0. При sides=2 offset устанавливается так, что бы подгруппа сворачивалась в центр. Правильно?

Другой offset задать видимо нельзя?

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

filter(x,c(0.25,0.25,0.25,0.25), method=«convolution»,sides=1)

а треугольное сглаживание вот так:

(filter(1:10,c(1,1,1), method=«convolution»,sides=2))/3

К чему тогда все манипуляции с mean и window в прошлом посте?

я рад что вы освоили filter :)

К чему тогда все манипуляции с mean и window в прошлом посте?

годный вброс :)

поступив так как написано вы можете получить ответ на

Другой offset задать видимо нельзя?

простым путем --- изменяя смещение window в функции помещаемой в sapply

Треугольное сглаживание применяется к трем последовательным ячейкам таблицы, причем значение в средней ячейке заменяется на новое сглаженное значение. Оператор сглаживания относится к типу 25-50-25, то есть сглаженное значение рассчитывается как сумма 25% первого значения, 50% второго значения и 25% третьего значения. Затем фильтр смещается на одну ячейку вниз и т.д., пока не достигается конец выделенного диапазона ячеек. Первая и последняя ячейки выделенного диапазона являются особыми ячейками. Сглаженное значение для первой ячейки вычисляется как сумма 2/3 первого значения и 1/3 второго значения. Сглаженное значение для последней ячейки вычисляется как сумма 2/3 последнего значения и 1/3 предпоследнего значения. Это сглаживание не преобразует пустые ячейки. Если одна из трех ячеек пустая, то к этим трем ячейкам применяется оператор 2/3-1/3 вместо оператора 25-50-25. Можно использовать несколько итераций треугольного сглаживания. На каждой итерации сглаживание применяется к обновленным значениям (т.е. к значениям после предыдущей итерации).

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

но судя по описанию

> filter(rep(1,10),c(0.25,0.5,0.25), method="convolution",sides=2)
Time Series:
Start = 1 
End = 10 
Frequency = 1 
 [1] NA  1  1  1  1  1  1  1  1 NA

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

> я рад что вы освоили filter :)

А я как рад - не передать ;)

Про window теперь понял - спасибо. Но уже сделал так:

movingAverageBwd <- rev(filter(rev(temp),              rep(1,fMA[1]),              method=«convolution»,              sides=1)           )/fMA[1]

Потом может переделаю...

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

Тут такое дело - у каждого своё представление о треугольном сглаживании, есть даже какое-то такое:

filter(rep(1,10),c(0.5,0,0.5), method=«convolution»,sides=2)

Повторенное трижды.

Ну это на столько тривиально, что даже у меня нет вопросов ;)

а почему усредняется «будущее» временного ряда? новые веяния в «стране вечно скользящих средних»? :)

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

P.S. Это не реальный проект - пытаюсь разобраться в R. Возможно в феврале придётся заняться проектом где много как статистики, так и инженерных расчётов. И ещё целый ряд интересных требований ТЗ. Такчто учу параллельно R и Tcl...

Очень рекомендую посмотреть на Sweave (редактировать естественно в emacs с ess и модами для LaTeX). Разрабатывать эталонное решение сложной задачи очень продуктивно.

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

PS а весь ряд реверсить не обязательно можно например задать вот так в фильтре c(rep(0,6),rep(1,7)/7) при sides=2

Спасибо.
Задача как раз несложная, просто ресурсоёмкая. Куча решений СЛАУ - там кода - дюжена строк, но много раз - важна скорость.
Возможно даже и не будем использовать R есть пакет frame3dd устраивающий нас в части инженерных расчётов, а может и будем. Но нужно быть готовым ко всему - задачка нетривиальная...
А так как сейчас есть чуть-чуть свободного времени - тренируюсь.
Да и вообще интересно осилить годный калькулятор - пользуются же виндоюзеры Фотошопом для просмотра фотографий. А чем мы хуже? Я уже исправил пункт вызова калькулятора в меню на 'gnome-terminal -x R' :)

Но посмотрю Sweave обязательно.

PS
Логично - крутилась какая-то такая мысль, но думать её не стал - так и перепишу.
Ещё раз - огромное спасибо за помощь. С моим уровнем английского (да и образования вообще) трудно въезжать в это самому...

> Я уже исправил пункт вызова калькулятора в меню на 'gnome-terminal -x R' :)

В таком случае очень должен понравится J http://www.jsoftware.com/

У него хороший учебник в комплекте. Кроме того он наследник APL, а именно APL среда предшественник S(R). Техника вычисления на нем частично применима в R.

Достоинство: J работает на смартофонах основанных на покетпц, превращая их в неслабые вычислители. График с 32 тыс. точек не представляет проблемы.

Ну и очень немногословен.

Вот ещё один вопрос - при экспорте таблицы так:

write.table(basis,
            file = pathes["pathToBasis"],
            quote = TRUE,
            append = FALSE,
            sep=";", dec=",")

Десятичный разделитель устанавливается как "," только для колонок которые содержат только числа. Колонки содержащие NA экспортируются с разделителем "."

Это баг, я чего-то недопонимаю, или в этом есть какой-то тайный смысл?

Есть у меня покетпц от хп - посмотрю :)

?write.table

The ‘dec’ argument only applies to columns that are not subject to conversion to character because they have a class or are part of a matrix-like column (or matrix), in particular to columns protected by ‘I()’. Use ‘options(«OutDec»)’ to control such conversions.

?options

‘OutDec’: character string containing a single-byte character. The character to be used as the decimal point in output conversions, that is in printing, plotting and ‘as.character’ but not deparsing.

ага - не баг, и то хорошо...

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

Например, есть таблица, где в колонке value:
1, 2, 4, 1, 4, 3, 2, 1, 5, 4, 2, 4, 1, 1, 3, 4, 6, 1
А N==9.
Тогда разбивка нужна такая:
1, 2, 4, 1
4, 3, 2
1, 5, 4
2, 4, 1, 1
3, 4
6, 1

Так как в этом случае сумма разности сумм каждого ряда и N минимальна.
Значит функция типа endpoints, должна вернуть что то типа:
0, 4, 7, 10, 14, 16, 18

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

В манах ничего не нашёл (или пропустил непоняв) - искал по ??interval и ??devide

> data<-c(1, 2, 4, 1, 4, 3, 2, 1, 5, 4, 2, 4, 1, 1, 3, 4, 6, 1)
> cumsum(data)
 [1]  1  3  7  8 12 15 17 18 23 27 29 33 34 35 38 42 48 49
> cumsum(data)-9
 [1] -8 -6 -2 -1  3  6  8  9 14 18 20 24 25 26 29 33 39 40
> cumsum(data)-9-9
 [1] -17 -15 -11 -10  -6  -3  -1   0   5   9  11  15  16  17  20  24  30  31
> cumsum(data)-9-9-9
 [1] -26 -24 -20 -19 -15 -12 -10  -9  -4   0   2   6   7   8  11  15  21  22
> cumsum(data)-9-9-9-9
 [1] -35 -33 -29 -28 -24 -21 -19 -18 -13  -9  -7  -3  -2  -1   2   6  12  13
> cumsum(data)-9-9-9-9-9
 [1] -44 -42 -38 -37 -33 -30 -28 -27 -22 -18 -16 -12 -11 -10  -7  -3   3   4
> cumsum(data)-9-9-9-9-9-9
 [1] -53 -51 -47 -46 -42 -39 -37 -36 -31 -27 -25 -21 -20 -19 -16 -12  -6  -5
> 

как то так я думаю.

Идею понял.
Кстати, не знал про кумулятивные команды. Вообще люто не хватает какого-нибудь справочника команд, где была бы команда и краткое описание в 2-3 слова...

Кстати такой подход может дать только интервалы где сумма меньше или равна N. Ну да ладно я подумаю и разовью идею.
Спасибо.

по моему неплохой справочник строит команда help.start()

и то что он показывает в base надо прочитать, я считаю, этот раздел формировали очень компетентные люди

в каждой статье есть раздел на что еще похожее можно опереться

например очень неплохо поиск делать через RSiteSearch(«что ищем»), только галочки добавить, доступ к спискам рассылки дорогого стоит.

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

Да, я видел, но это не очень удобно - много букв и к тому же не русских. У меня плохо с ними. Но и по русски такое читать иногда затруднительно.

Мне нравится когда просто тупо перечень команд. Ведь о многих можно вообще по их написанию догадаться. Ну вот тот же cumsum, например.

в base оглавлении несколько экранов всего, всё кратко и по английски... перевода мануалов нет, есть перевод введения в R, но там естественно не все команды.

я думаю в процессе чтения выучите, там всё на уровне cumsum :)

... что бы легче читалось надо почитать что либо по анализу данных

Да, такую тему испортили, советами читать base. После такого чтения, нубских вопросов почти и не остаётся... :(

Ну ничего - я тот ещё нуб! :)

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

Т.е. Есть вектор V например длинной 25. И фактор, например такой: 1,1,1,1,1,2,2,2,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,1,1. Нужно применить функцию соотвенно к его участкам 1:5,6:8,9:15,16:23 и 24:25.

Не могу сообразить как...

Лучше конечно вообще получить вектор 1,5,6,8,9,15,16,23,24,25. Тогда и задачу решить можно с помощью period.apply и в дальнейшем понадобится, но моожно и другой способ...

получить список с разбиением by(1:25,f,сbind)

получить первый элемент в виде вектора by(1:25,f,cbind)[[1]][,1]

или tapply(1:25,f,cbind)

надо прочитать Введение в R

на русском http://m7876.wiki.zoho.com/Introduction-to-R.html&pid=78912000000002015

Спасибо - я знаю и про by и про tapply, но они дают разбиение по уровням. Т.е.в итоге я получу вектор разбитый на два подвектора - для уровня 1 и для уровня 2:

> x <- c(1:25)
> f <- factor(c(rep(1,7),))
> f <- factor(c(rep(1,7),rep(2,4),rep(1,5),rep(2,9)))
> f
 [1] 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 2 2 2 2 2 2 2 2 2
Levels: 1 2
> tapply(x,f,cbind)
$`1`
      [,1]
 [1,]    1
 [2,]    2
 [3,]    3
 [4,]    4
 [5,]    5
 [6,]    6
 [7,]    7
 [8,]   12
 [9,]   13
[10,]   14
[11,]   15
[12,]   16

$`2`
      [,1]
 [1,]    8
 [2,]    9
 [3,]   10
 [4,]   11
 [5,]   17
 [6,]   18
 [7,]   19
 [8,]   20
 [9,]   21
[10,]   22
[11,]   23
[12,]   24
[13,]   25

То есть групировка по уровням. А надо как-то так:

> f <- factor(c(rep(1,7),))
> f <- factor(c(rep(1,7),rep(2,4),rep(1,5),rep(2,9)))
> f
 [1] 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 2 2 2 2 2 2 2 2 2
Levels: 1 2
> ????(x,f,cbind)
$`1`
      [,1]
 [1,]    1
 [2,]    2
 [3,]    3
 [4,]    4
 [5,]    5
 [6,]    6
 [7,]    7

$`2`
      [,1]
 [1,]    8
 [2,]    9
 [3,]   10
 [4,]   11

$`3`
      [,1]
 [1,]   12
 [2,]   13
 [3,]   14
 [4,]   15
 [5,]   16

$`4`
      [,1]
[1,]   17
[2,]   18
[3,]   19
[4,]   20
[5,]   21
[6,]   22
[7,]   23
[8,]   24
[9,]   25

> split(1:25,f)[[1]]-c(split(1:25,f)[[1]][-1],NA)
 [1] -1 -1 -1 -1 -4 -1 -1 -1 -1 -1 -1 -9 -1 NA
> split(1:25,f)[[2]]-c(split(1:25,f)[[2]][-1],NA)
 [1] -1 -1 -8 -1 -1 -1 -1 -1 -1 -1 NA

> split(1:25,f)[[1]]-c(split(1:25,f)[[1]][-1],NA) +1
 [1]  0  0  0  0 -3  0  0  0  0  0  0 -8  0 NA
> cumsum(split(1:25,f)[[1]]-c(split(1:25,f)[[1]][-1],NA) +1)
 [1]   0   0   0   0  -3  -3  -3  -3  -3  -3  -3 -11 -11  NA

вполне фактор появляется, им можно получить индекс в исходном векторе

еще нашел ?rle

Compute the lengths and values of runs of equal values in a vector - or the reverse operation.

> rle(f)
Run Length Encoding
  lengths: int [1:5] 5 3 7 8 2
  values : num [1:5] 1 2 1 2 1

Интересно я один такой тупой, или это язык такой не для средних умов? Раньше не замечал за собой такого ни в асме, ни в basic, ни в tcl...

Спасибо!

Кстати вектор не обязательно упорядоченный брать - можно из фактора расчитать сразу:

# плучаем факторы для каждого диапазоны
		# здесь вектор вида:
# [1]  1  1  1  1  1  1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1  1  1  1
# [26] 1  1  1  1  1  1  1  1  1  1  1 -1 -1 -1 -1  1  1  1  1  1  1  1  1 -1 -1
		# превращается в вектор вида:
# [1]  0  0  0  0  0  0  0  2  0  0  0  0  0  0  0  0  0 -2  0  0  0  0  0  0  0
# [26] 0  0  0  0  0  0  0  0  0  0  2  0  0  0 -2  0  0  0  0  0  0  0  2  0  0
		uniqPrdFactor <- dyrtiDirFactor -
				 c(dyrtiDirFactor[-1],last(dyrtiDirFactor));
		# далее делаем его положительным чтобы не обнулять cumsum
# [1]  0  0  0  0  0  0  0  2  0  0  0  0  0  0  0  0  0  2  0  0  0  0  0  0  0
# [26] 0  0  0  0  0  0  0  0  0  0  2  0  0  0  2  0  0  0  0  0  0  0  2  0  0
		uniqPrdFactor <- abs(uniqPrdFactor);
		# создаём уникальные факторы с помощью cumsum
		uniqPrdFactor <- cumsum(uniqPrdFactor);
		# на выходе имеем
# [1]  0  0  0  0  0  0  0  2  2  2  2  2  2  2  2  2  2  4  4  4  4  4  4  4  4
# [26] 4  4  4  4  4  4  4  4  4  4  6  6  6  6  8  8  8  8  8  8  8  8 10 10 10
		# остаётся сдвинуть его вправо и преобразовать в фактор
	uniqPrdFactor <- factor(c(0,uniqPrdFactor[length(uniqPrdFactor)*-1]))

Ага - посмотрю сейчас. Пока сделал как в сообщении выше - работает!

Другая парадигма. «Векторный» язык, наследие APL,J,K.

> надо прочитать Введение в R ... на русском http://m7876.wiki.zoho.com/Introduction-to-R.html&pid=78912000000002015

А я то думал, у меня уже есть все материалы по R на великом и могучем. А тут такая книжечка. Огромный респект переводчикам.

Материал очень полно и главное доходчиво изложен, с объяснением принципов, что очень важно.

> Огромный респект переводчикам.

Спасибо! :)

> Спасибо! :)

Я так и думал ;)

Есть фактор и вектор, надо получить из вектора экстремумы по факторам, но не значения, а их индексы в векторах. Как?