Язык программирования R/Математика

Основные сведения

править

Базовые сведения по математическим операциям можно почерпнуть из встроенной справочной системы:

?Arithmetic
?Special

Линейная алгебра

править

Векторы

править

Скалярное произведение (inner product)

править

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

В представленном примере вектор (3,3,3) умножается на вектор (1,2,3), то есть выполняются следующие действия:

  1. 3*1=3 (запоминаем)
  2. 3*2=6 (запоминаем)
  3. 3*3=9 (запоминаем)
  4. (вспоминаем) 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(). Для этого требуется:

  1. Передать вектор данных.
  2. Количество столбцов и/или строк.
  3. Задать способ обработки: по столбцам (по умолчанию) или строкам; при помощи опции 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 позволяет получить справку по специальным математическим функциям, например, таким как Бетта и Гамма функции.

Ссылки

править
  1. материал Википедии о тензорных произведениях
  2. Maxima - проект с открытым исходным кодом http://maxima.sourceforge.net/(англ.)
  3. SAGE - проект с открытым исходным кодом включающий R и Maxima: http://www.sagemath.org/(англ.)
  4. Maple - коммерческий закрытый, но популярный продукт http://www.maplesoft.com/products/maple/(англ.)
  5. Mathematica - коммерческий закрытый, но, также, популярный продукт http://www.wolfram.com/products/mathematica/index.html(англ.)