Язык программирования 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)

Из линейной алгебры известно:

$X=\left(\begin{array}{ccc}1& 2& 3\end{array}\right);Y=\left(\begin{array}{c}4\\ 5\\ 6\end{array}\right);X\cdot Y=\left(\begin{array}{ccc}1\cdot 4& 1\cdot 5& 1\cdot 6\\ 2\cdot 4& 2\cdot 5& 2\cdot 6\\ 3\cdot 4& 3\cdot 5& 3\cdot 6\end{array}\right)=\left(\begin{array}{ccc}4& 5& 6\\ 8& 10& 12\\ 12& 15& 18\end{array}\right)$

Вот именно эту операцию 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")

Единичная матрица

Единичная матрица в линейной алгебре выполняет функцию единицы из обычной алгебры и является замечательным элементом множества матриц из-за своего свойства: $A\cdot E=A$, где $A$ - любая матрица, а $E$ - единичная. Фактически, единичная матрица - это квадратная матрица подходящего, для совершения операции, размера, где по главной диагонали расположены единицы, а все остальные элементы равны нулю.

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

Нулевая матрица

Также, замечательный элемент множества матриц, выполняющий в матричной алгебре функции нуля. Обладает свойством: $A\cdot O=O$, где $A$ - любая матрица, $O$ - нулевая матрица подходящего размера. Фактически, нулевая матрица полностью состоит из нулей.

Сгенерировать нулевую матрицу заданного размера можно при помощи функции 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=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right);B=\left(\begin{array}{cc}1& 3\\ 2& 4\end{array}\right);A\cdot B=\left(\begin{array}{cc}1& 3\\ -2& -4\end{array}\right);B\cdot A=\left(\begin{array}{cc}1& -3\\ 2& -4\end{array}\right)$

> 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

Тензорное произведение (Произведение Кронекера)

Тензорное произведение векторов :

$𝐚\otimes 𝐛\to \left[\begin{array}{c}{a}_1\\ a_2\\ a_3\\ a_4\end{array}\right]\left[\begin{array}{ccc}{b}_1& b_2& b_3\end{array}\right]=\left[\begin{array}{ccc}{a}_1{b}_1& a_1{b}_2& a_1{b}_3\\ a_2{b}_1& a_2{b}_2& a_2{b}_3\\ a_3{b}_1& a_3{b}_2& a_3{b}_3\\ a_4{b}_1& a_4{b}_2& a_4{b}_3\end{array}\right]$

Тензорное произведение матриц:

$\mathrm{A}=\left[\begin{array}{ccc}{a}_{11}& \cdots & a_{1n}\\ ⋮& \ddots & ⋮\\ a_{m1}& \cdots & a_{mn}\end{array}\right]$

$\mathrm{B}=\left[\begin{array}{ccc}{b}_{11}& \cdots & b_{1q}\\ ⋮& \ddots & ⋮\\ b_{p1}& \cdots & b_{pq}\end{array}\right]$

$\mathrm{A}\otimes \mathrm{B}=\left[\begin{array}{ccc}{a}_{11}B& \cdots & a_{1n}B\\ ⋮& \ddots & ⋮\\ a_{m1}B& \cdots & a_{mn}B\end{array}\right]=$ $=\left[\begin{array}{cccccccccc}{a}_{11}{b}_{11}& a_{11}{b}_{12}& \cdots & a_{11}{b}_{1q}& \cdots & \cdots & a_{1n}{b}_{11}& a_{1n}{b}_{12}& \cdots & a_{1n}{b}_{1q}\\ a_{11}{b}_{21}& a_{11}{b}_{22}& \cdots & a_{11}{b}_{2q}& \cdots & \cdots & a_{1n}{b}_{21}& a_{1n}{b}_{22}& \cdots & a_{1n}{b}_{2q}\\ ⋮& ⋮& \ddots & ⋮& & & ⋮& ⋮& \ddots & ⋮\\ a_{11}{b}_{p1}& a_{11}{b}_{p2}& \cdots & a_{11}{b}_{pq}& \cdots & \cdots & a_{1n}{b}_{p1}& a_{1n}{b}_{p2}& \cdots & a_{1n}{b}_{pq}\\ ⋮& ⋮& & ⋮& \ddots & & ⋮& ⋮& & ⋮\\ ⋮& ⋮& & ⋮& & \ddots & ⋮& ⋮& & ⋮\\ a_{m1}{b}_{11}& a_{m1}{b}_{12}& \cdots & a_{m1}{b}_{1q}& \cdots & \cdots & a_{mn}{b}_{11}& a_{mn}{b}_{12}& \cdots & a_{mn}{b}_{1q}\\ a_{m1}{b}_{21}& a_{m1}{b}_{22}& \cdots & a_{m1}{b}_{2q}& \cdots & \cdots & a_{mn}{b}_{21}& a_{mn}{b}_{22}& \cdots & a_{mn}{b}_{2q}\\ ⋮& ⋮& \ddots & ⋮& & & ⋮& ⋮& \ddots & ⋮\\ a_{m1}{b}_{p1}& a_{m1}{b}_{p2}& \cdots & a_{m1}{b}_{pq}& \cdots & \cdots & a_{mn}{b}_{p1}& a_{mn}{b}_{p2}& \cdots & a_{mn}{b}_{pq}\end{array}\right]$

В 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

По определению обратной матрицы: $A\cdot A^{-1}=E$, где $A$ - любая матрица, а $E$ - единичная. Проверяем:

> solve(M)%*%M
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

Решение системы линейных уравнений

Из линейной алгебры известно, что система:

$\{\begin{array}{c}{a}_{11}{x}_1+a_{12}{x}_2+\dots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\dots +a_{2n}{x}_n=b_2\\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \\ a_{m1}{x}_1+a_{m2}{x}_2+\dots +a_{mn}{x}_n=b_m\end{array}$

эквивалентна:

$A\cdot X=B,$

где:

$A=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right);X=\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right);B=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_m\end{array}\right).$

Следовательно:

$X=A^{-1}\cdot B$

Этот алгоритм в коде 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
  • 10**3
• Взять логарифм можно функцией log(), где первый аргумент - это подлогорифменное значение, а второй, опциональный, равный по умолчанию $e$ - база. Таким образом, функция log() без указания базы, по умолчанию, берёт натуральный логарифм.
• Взять десятичный логарифм можно функцией log10().
• Вычислить $e^x$ можно функцией 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

Решение многочленов

Что бы решить уравнение $a{x}^k+b{x}^{k-1}+\cdots +n=0$, где $a,b,\dots ,n$ - данные числа, следует применить функцию polyroot():

> polyroot(c(n,...,b,a))

Например, вычислим корни уравнения $2x^2-5x-3=0$:

> polyroot(c(-3,-5,2))
 [1] -0.5+0i  3.0-0i

ответ следует читать как: $x_1=-0.5\wedge x_2=3$.

Смотрите также пакеты 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Шаблон:Ref-en.

Интегралы

Интегрирование

R может производить одноразмерное интегрирование. Например, мы можем проинтегрировать по плотности нормального распределения от $-\mathrm{\infty }$ до $+\mathrm{\infty }$.

> 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 

Вероятности

• Сочетание $C_n^k$ вычисляется с помощью функции 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

Внимание: имейте ввиду, что по соглашению $0!=1$, и функция 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Шаблон:Ref-en) предоставляет пакет функций sympy для R (linkШаблон:Ref-en).

Если вы желаете производить множество символьных вычислений, то имеет смысл обратить внимание на пакеты Maxima, SAGE, Maple, Mathematica.

См. также

Команда ?Special позволяет получить справку по специальным математическим функциям, например, таким как Бетта и Гамма функции.