Язык программирования R/Математика
Основные сведения
правитьБазовые сведения по математическим операциям можно почерпнуть из встроенной справочной системы:
?Arithmetic ?Special
Линейная алгебра
правитьВекторы
правитьСкалярное произведение (inner product)
правитьПоследовательное произведение каждого члена вектора на соответствующий член другого вектора и сумма получившихся произведений.
В представленном примере вектор (3,3,3)
умножается на вектор (1,2,3)
, то есть выполняются следующие действия:
- 3*1=3 (запоминаем)
- 3*2=6 (запоминаем)
- 3*3=9 (запоминаем)
- (вспоминаем) 3+6+9=18
Таким образом векторы R совместимы с векторами обычной линейной алгебры.
> u <- rep(3,3) > v <- 1:3 > u%*%v # Скалярное произведение [,1] [1,] 18
Внешнее произведение (outer product)
правитьИз линейной алгебры известно:
Вот именно эту операцию R считает внешним произведением:
> u=1:3 > v=4:6 > u [1] 1 2 3 > v [1] 4 5 6 > u%o%v [,1] [,2] [,3] [1,] 4 5 6 [2,] 8 10 12 [3,] 12 15 18
Матрицы
правитьЕсли требуется создать матрицу, то нужно использовать функцию matrix()
. Для этого требуется:
- Передать вектор данных.
- Количество столбцов и/или строк.
- Задать способ обработки: по столбцам (по умолчанию) или строкам; при помощи опции byrow.
> matrix(data = NA, nrow = 5, ncol = 5, byrow = T) [,1] [,2] [,3] [,4] [,5] [1,] NA NA NA NA NA [2,] NA NA NA NA NA [3,] NA NA NA NA NA [4,] NA NA NA NA NA [5,] NA NA NA NA NA > X <- matrix(data = 1:15, nrow = 5, ncol = 5, byrow = T) > X [,1] [,2] [,3] [,4] [,5] [1,] 1 2 3 4 5 [2,] 6 7 8 9 10 [3,] 11 12 13 14 15 [4,] 1 2 3 4 5 [5,] 6 7 8 9 10
"Склеить" матрицу из нескольких векторов можно функциями cbind(v1,v2)
(по столбцам) или rbind(v1,v2)
(по строкам).
> v1 <- 1:5 > v2 <- 5:1 > cbind(v1,v2) v1 v2 [1,] 1 5 [2,] 2 4 [3,] 3 3 [4,] 4 2 [5,] 5 1 > rbind(v1,v2) [,1] [,2] [,3] [,4] [,5] v1 1 2 3 4 5 v2 5 4 3 2 1
Размерность матрицы может быть получена функцией dim()
. Узнать количество столбцов или строк можно, соответственно, функциями nrow()
или ncol()
.
> dim(X) [1] 5 5 > nrow(X) [1] 5 > ncol(X) [1] 5
Некоторые замечательные виды матриц
правитьДля работы с этими видами матриц требуется специальный набор функций из пакета matlab (в самом Matlab-е эти функции называются идентично).
> install.packages("matlab") > library("matlab")
Единичная матрица
правитьЕдиничная матрица в линейной алгебре выполняет функцию единицы из обычной алгебры и является замечательным элементом множества матриц из-за своего свойства: , где - любая матрица, а - единичная. Фактически, единичная матрица - это квадратная матрица подходящего, для совершения операции, размера, где по главной диагонали расположены единицы, а все остальные элементы равны нулю.
R позволяет сгенерировать единичную матрицу нужного размера при помощи функции eye()
, которой в качестве аргумента передаётся размер матрицы.
> eye(3) [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] 0 0 1
Матрица единиц
правитьКвадратную матрицу полностью состоящую из единиц можно получить при помощи функции ones(), которой в качестве параметра передаётся размер:
> ones(3) [,1] [,2] [,3] [1,] 1 1 1 [2,] 1 1 1 [3,] 1 1 1
Нулевая матрица
правитьТакже, замечательный элемент множества матриц, выполняющий в матричной алгебре функции нуля. Обладает свойством: , где - любая матрица, - нулевая матрица подходящего размера. Фактически, нулевая матрица полностью состоит из нулей.
Сгенерировать нулевую матрицу заданного размера можно при помощи функции zeros(), которой в качестве параметра передаётся размер:
> zeros(3) [,1] [,2] [,3] [1,] 0 0 0 [2,] 0 0 0 [3,] 0 0 0
Операции над матрицами
правитьПолучение диагонали
правитьДля получения диагонали матрицы нужно отдать матрицу в качестве аргумента функции diag()
, и она вернёт вектор составленный из элементов диагонали:
> (x <- matrix(0:15,ncol=4)) [,1] [,2] [,3] [,4] [1,] 0 4 8 12 [2,] 1 5 9 13 [3,] 2 6 10 14 [4,] 3 7 11 15 > diag(x) [1] 0 5 10 15
Вычленение верхнего угла или нижнего
правитьДля получения верхнего угла можно использовать функцию upper.tri()
, нижнего - lower.tri()
. Эти функции принимают матрицу в качестве аргумента и возвращают матрицу идентичного размера с элементами типа Boolean. На месте тех членов, которые должны быть включены, стоит TRUE, на месте прочих - FALSE. Дополнительный аргумент diag
(по умолчанию FALSE) определяет включать диагональ или нет. Более наглядно работа этих функций видна на примерах:
> (x <- matrix(0:15,ncol=4)) [,1] [,2] [,3] [,4] [1,] 0 4 8 12 [2,] 1 5 9 13 [3,] 2 6 10 14 [4,] 3 7 11 15 > lower.tri(X) [,1] [,2] [,3] [,4] [,5] [1,] FALSE FALSE FALSE FALSE FALSE [2,] TRUE FALSE FALSE FALSE FALSE [3,] TRUE TRUE FALSE FALSE FALSE [4,] TRUE TRUE TRUE FALSE FALSE [5,] TRUE TRUE TRUE TRUE FALSE > upper.tri(X, diag=T) [,1] [,2] [,3] [,4] [,5] [1,] TRUE TRUE TRUE TRUE TRUE [2,] FALSE TRUE TRUE TRUE TRUE [3,] FALSE FALSE TRUE TRUE TRUE [4,] FALSE FALSE FALSE TRUE TRUE [5,] FALSE FALSE FALSE FALSE TRUE
Зная логику работы функций несложно "отрезать" от имеющейся матрицы, например, нижний угол без диагонали:
> Y <- X > Y[lower.tri(X)] <- 0 > Y [,1] [,2] [,3] [,4] [,5] [1,] 1 2 3 4 5 [2,] 0 7 8 9 10 [3,] 0 0 13 14 15 [4,] 0 0 0 4 5 [5,] 0 0 0 0 10
Произведение матриц
правитьПример ниже демонстрирует:
> a=matrix(nrow=2,ncol=2,c(1,0,0,-1)) > b=matrix(nrow=2,ncol=2,c(1,2,3,4)) > a [,1] [,2] [1,] 1 0 [2,] 0 -1 > b [,1] [,2] [1,] 1 3 [2,] 2 4 > a%*%b [,1] [,2] [1,] 1 3 [2,] -2 -4 > b%*%a [,1] [,2] [1,] 1 -3 [2,] 2 -4
Тензорное произведение (Произведение Кронекера)
правитьТензорное произведение векторов [1]:
Тензорное произведение матриц:
В R эта операция реализуется оператором %x%
или функцией kron()
из библиотеки fUtilities:
> M <- matrix(rep(2,4),nrow = 2) > M [,1] [,2] [1,] 2 2 [2,] 2 2 > I <- eye(2) > I [,1] [,2] [1,] 1 0 [2,] 0 1 > I %x% M [,1] [,2] [,3] [,4] [1,] 2 2 0 0 [2,] 2 2 0 0 [3,] 0 0 2 2 [4,] 0 0 2 2 > library(fUtilities) > kron(I,M) [,1] [,2] [,3] [,4] [1,] 2 2 0 0 [2,] 2 2 0 0 [3,] 0 0 2 2 [4,] 0 0 2 2
Транспонирование
правитьДля транспонирования матрицы нужно применить функцию t()
:
> (M <- matrix(1:9, ncol=3)) [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > t(M) [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 [3,] 7 8 9
След матрицы
правитьДля получения следа матрицы применяется функция tr()
из пакета fUtilities:
> install.packages("fUtilities") > library(fUtilities) > (M <- matrix(1:9, ncol=3)) [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > tr(M) [1] 15
Ранг матрицы
правитьДля определения ранга матрицы используется функция qr()$rank:
> (M <- matrix(1:9, ncol=3)) [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > qr(M)$rank [1] 2
Определитель матрицы
правитьДля вычисления определителя матрицы применяется функция det()
:
> (M <- matrix(c(1,2,3,5), ncol=2)) [,1] [,2] [1,] 1 3 [2,] 2 5 > det(M) [1] -1
Получение обратной матрицы
правитьДля получения обратной матрицы можно воспользоваться функциями solve()
, inv()
(пакет fUtilities) или ginv()
(пакет MASS):
> M <- cbind(c(1,0,1),c(0,1,2),c(0,0,1)) > solve(M) [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] -1 -2 1 > inv(M) [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] -1 -2 1 > ginv(M) [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] -1 -2 1
По определению обратной матрицы: , где - любая матрица, а - единичная. Проверяем:
> solve(M)%*%M [,1] [,2] [,3] [1,] 1 0 0 [2,] 0 1 0 [3,] 0 0 1
Решение системы линейных уравнений
правитьИз линейной алгебры известно, что система:
эквивалентна:
где:
Следовательно:
Этот алгоритм в коде R будет выглядеть следующим образом:
> a=matrix(nrow=2,ncol=2,c(1,-.8,1,.2)) > a [,1] [,2] [1,] 1.0 1.0 [2,] -0.8 0.2 > > b=matrix(c(1.0+25.0/18,25.0/18.0)) > b [,1] [1,] 2.388889 [2,] 1.388889 > > x=solve(a,b) > x [,1] [1,] -0.9111111 [2,] 3.3000000 > > a%*%x # Проверяем ответ [,1] [1,] 2.388889 [2,] 1.388889
Собственные вектора, собственные числа и корневые подпространства
правитьТеория этих терминов не совсем тривиальна, а потому рекомендуется ознакомиться с со статьёй Википедии "Собственные векторы, значения и пространства".
Для вычисления собственных элементов (чисел и векторов) следует применять функцию eigen()
.
> (M <- matrix(1:9, ncol=3)) [,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9 > eigen(M) $values [1] 1.611684e+01 -1.116844e+00 -4.054214e-16 $vectors [,1] [,2] [,3] [1,] -0.4645473 -0.8829060 0.4082483 [2,] -0.5707955 -0.2395204 -0.8164966 [3,] -0.6770438 0.4038651 0.4082483
Прочее
править- Вычислить норму можно функцией
norm()
(пакет fUtilities). - Проверить, что матрица положительно определена можно функцией
isPositiveDefinite()
(пакет fUtilities). - Привести матрицу к положительно определённой форме можно функцией
makePositiveDefinite()
(пакет fUtilities). - Вычленить верхний или нижний треугольник матрицы можно, соответственно, функциями
Triang()
иtriang()
(пакет fUtilities).
Также полезными могут оказать другие функции из пакетов matrix, matlab, matrixcalc и matrixStats.
Математический анализ
правитьЛогарифмы и экспоненты
править- Вычислить степень числа можно двумя операторами (например, ):
10^3
10**3
- Взять логарифм можно функцией
log()
, где первый аргумент - это подлогорифменное значение, а второй, опциональный, равный по умолчанию - база. Таким образом, функцияlog()
без указания базы, по умолчанию, берёт натуральный логарифм. - Взять десятичный логарифм можно функцией
log10()
. - Вычислить можно функцией
exp()
.
Примеры:
> 10^3 # Степень [1] 1000 > 10**3 # Другой способ возведения в степень [1] 1000 > exp(1) # e в первой степени [1] 2.718282 > log(2.71) # Взятие натурального логарифма [1] 0.9969486 > log10(1000) # Взятие десятичного логарифма [1] 3 > log(1000,base = 10) # Взятие произвольного логарифма по основанию 10 [1] 3 > log(1024,base = 2) # Взятие произвольного логарифма по основанию 2 [1] 10
Решение многочленов
правитьЧто бы решить уравнение , где - данные числа, следует применить функцию polyroot()
:
> polyroot(c(n,...,b,a))
Например, вычислим корни уравнения :
> polyroot(c(-3,-5,2)) [1] -0.5+0i 3.0-0i
ответ следует читать как: .
Смотрите также пакеты polynom и multipol.
Производные
правитьСимвольные вычисления производных
правитьR может вычислять производные от выражений. Нужно сконвертировать имеющееся выражение при помощи функции expression()
, после чего можно будет дифференцировать его.
Вот несколько примеров:
> D(expression(x^n),"x") x^(n - 1) * n > D(expression(exp(a*x)),"x") exp(a * x) * a > D(expression(1/x),"x") -(1/x^2) > D(expression(x^3),"x") 3 * x^2 > D(expression(pnorm(x)),"x") dnorm(x) > D(expression(dnorm(x)),"x") -(x * dnorm(x))
Численные приближения
правитьСмотрите документацию по пакету numDeriv(англ.).
Интегралы
правитьИнтегрирование
правитьR может производить одноразмерное интегрирование. Например, мы можем проинтегрировать по плотности нормального распределения от до .
> integrate(dnorm,-Inf,Inf) 1 with absolute error < 9.4e-05 > integrate(dnorm,-1.96,1.96) 0.9500042 with absolute error < 1.0e-11 > integrate(dnorm,-1.64,1.64) 0.8989948 with absolute error < 6.8e-14 # Можно сохранить результат в объект > ci90 <- integrate(dnorm,-1.64,1.64) > ci90$value [1] 0.8989948 > integrate(dnorm,-1.64,1.64)$value [1] 0.8989948
Для взятия интегралов от многих переменных обратитесь к пакету cubature (так как пакет adapt был удалён из CRAN по лицензионным соображениям). Ранее механизм интегрирования применялся следующим образом (сейчас предполагается заменить функцию adapt()
из пакета adapt на функцию adaptIntegrate()
из пакета cubature, но лично я код не тестировал):
> library(adapt) > ?adapt > ir2pi <- 1/sqrt(2*pi) > fred <- function(z) { ir2pi^length(z) * exp(-0.5 * sum(z * z))} > > adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred) value relerr minpts lenwrk ifail 1.039222 0.0007911264 231 73 0 > adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred, eps = 1e-4) value relerr minpts lenwrk ifail 1.000237 1.653498e-05 655 143 0 > adapt(2, lo = c(-5,-5), up = c(5,5), functn = fred, eps = 1e-6) value relerr minpts lenwrk ifail 1.000039 3.22439e-07 1719 283 0
Вероятности
править- Сочетание вычисляется с помощью функции
choose()
:
> choose(3,2) [1] 3
- Объединения и пересечения вычисляются, соответственно, при помощи функций
union()
иintersect()
.
> union(1:10,5:7) [1] 1 2 3 4 5 6 7 8 9 10 > intersect(1:10,5:7) [1] 5 6 7
Арифметика
правитьВычисление факториала
правитьФункция factorial()
возвращает факториал натурального числа. Конечно, также, можно применить функцию prod()
к вектору последовательных натуральных чисел от единицы.
> factorial(3) [1] 6 > prod(1:3) [1] 6
Внимание: имейте ввиду, что по соглашению , и функция factorial()
возвращает 1 от нулевого аргумента, однако, функция prod()
об этом "не знает" и вернёт 0.
> factorial(0) [1] 1 > prod(0) [1] 0
Значения факториала могут быть очень большими, и достаточно легко превысить лимиты ОС:
> factorial(170) [1] 7.257416e+306 > factorial(171) [1] Inf Предупреждение In factorial(171) : значение в 'gammafn' -- за пределами
Целочисленное деление
правитьИмеется возможность производить целочисленное деление:
> # Остаток от деления > 5%%2 [1] 1 > # Деление нацело >5%/%2 [1] 2
Геометрия
править- Число Пи - встроенная константа
pi
. - Тригонометрические функции:
cos()
,sin()
,tan()
.
Символьные вычисления
правитьrSymPy (rsympy(англ.)) предоставляет пакет функций sympy для R (link(англ.)).
Если вы желаете производить множество символьных вычислений, то имеет смысл обратить внимание на пакеты Maxima[2], SAGE[3], Maple[4], Mathematica[5].
См. также
правитьКоманда ?Special
позволяет получить справку по специальным математическим функциям, например, таким как Бетта и Гамма функции.
Ссылки
править- ↑ материал Википедии о тензорных произведениях
- ↑ Maxima - проект с открытым исходным кодом http://maxima.sourceforge.net/(англ.)
- ↑ SAGE - проект с открытым исходным кодом включающий R и Maxima: http://www.sagemath.org/(англ.)
- ↑ Maple - коммерческий закрытый, но популярный продукт http://www.maplesoft.com/products/maple/(англ.)
- ↑ Mathematica - коммерческий закрытый, но, также, популярный продукт http://www.wolfram.com/products/mathematica/index.html(англ.)