R Math
Basic
Sequence
1:5 - 1 # 0 1 2 3 4 5:1 # 5 4 3 2 1 pi:10 # 3.14159265359 4.14159265359 5.14159265359 6.14159265359 # 7.14159265359 8.14159265359 9.14159265359 seq(10) seq_len(10) # 1 2 3 4 5 6 7 8 9 10 seq.int(pi:10) # 1 2 3 4 5 6 7 seq(1, 20, by=3) # 1 4 7 10 13 16 19 seq(from = 1, by = .2, length = 10) # 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 my_seq = seq(5, 10,length=30) #== my_seq = seq(from = 5, to = 10,length = 30) # [1] 5.000000 5.172414 5.344828 5.517241 5.689655 5.862069 6.034483 6.206897 # [9] 6.379310 6.551724 6.724138 6.896552 7.068966 7.241379 7.413793 7.586207 # [17] 7.758621 7.931034 8.103448 8.275862 8.448276 8.620690 8.793103 8.965517 # [25] 9.137931 9.310345 9.482759 9.655172 9.827586 10.000000 seq(along.with = my_seq) ##== 1:30 seq_along(my_seq) ##== 1:30 |
rep(4:5, 4) ##== rep(4:5,times=4) # 4 5 4 5 4 5 4 5 rep(1:5, each = 4) # 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 rep(1:5, len = 12) # 1 2 3 4 5 1 2 3 4 5 1 2 rep(letters[1:4], 1:4) # [1] "a" "b" "b" "c" "c" "c" "d" "d" "d" "d" rep.int(1:4, 2) ## at least twice as fast when x has names # 1 2 3 4 1 2 3 4 rep_len(seq(1, 100, by=3), 9) # 1 4 7 10 13 16 19 22 25 |
sequence(3:7) # [1] 1 2 3 1 2 3 4 1 2 3 4 5 1 2 3 4 5 6 1 2 3 4 5 6 7 ## gl: Generate factor levels: levels, times, total gl(3, 10) # 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 # Levels: 1 2 3 gl(3, 10, 60, labels = c("A", "B","C")) # A A A A A A A A A A B B B B B B B B B B C C C C C C C C C C A A A A A A A A A # A B B B B B B B B B B C C C C C C C C C C # Levels: A B C interaction(gl(2, 5), gl(5, 2)) # 1.1 1.1 1.2 1.2 1.3 2.3 2.4 2.4 2.5 2.5 # Levels: 1.1 2.1 1.2 2.2 1.3 2.3 1.4 2.4 1.5 2.5 expand.grid(c(1:3),c(2:4)) # Var1 Var2 # 1 1 2 # 2 2 2 # 3 3 2 # 4 1 3 # 5 2 3 # 6 3 3 # 7 1 4 # 8 2 4 # 9 3 4 |
Simple Calculation
Avaiable triangle functions:
sin cos tan cot sec csc asin acos atan sinh cosh tanh
.5 + -2:4 # -1.5 -0.5 0.5 1.5 2.5 3.5 4.5 2^(1/2) # 1.414214 sqrt(2)+abs(-sin(pi)) # 1.41421356237 sqrt(-2+0i) # 0+1.414214i factorial(0) # 1 factorial(10) # 3628800 choose(10,4) # 210 choose(10,1:10) # 10 45 120 210 252 210 120 45 10 1 combn(5,3) # 1 1 1 1 1 1 2 2 2 3 # 2 2 2 3 3 4 3 3 4 4 # 3 4 5 4 5 5 4 5 5 5 |
x <- -1:12 2 * x + x ^ 2 #== 2 * x + x**2 # -1 0 3 8 15 24 35 48 63 80 99 120 143 168 c(1,2) * x # -1 0 1 4 3 8 5 12 7 16 9 20 11 24 ## Quotient x %% 5 # 4 0 1 2 3 4 0 1 2 3 4 0 1 2 ## Modulus x %/% 5 # -1 0 0 0 0 0 1 1 1 1 1 2 2 2 min(x);max(x);range(x);median(x);length(x) # -1 # 12 # -1 12 # 5.5 # 14 which(x < 2) # 1 2 3 which.max(x);which.min(x);range(x) # 14 # 1 # -1 12 mean(x) weighted.mean(x, rep(1,length(x))) # 5.5 weighted.mean(x, 1:length(x)) # 7.666667 median(x) # 5.5 sum(x);prod(x[-2]) # 77 # -479001600 cumsum(x);cumprod(x);cummax(x);cummin(x) # [1] -1 -1 0 2 5 9 14 20 27 35 44 54 65 77 # [1] -1 0 0 0 0 0 0 0 0 0 0 0 0 0 # [1] -1 0 1 2 3 4 5 6 7 8 9 10 11 12 # [1] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 sd(x)^2 ## sample standard deviation # sum((x-mean(x))^2)/(length(x)-1) # 17.5 var(x) ## sample variance # 17.5 mad(x) ## Median absolute deviation # 1.4826*mean(abs(x-median(x))) # 5.1891 sign(x) # -1 0 1 1 1 1 1 1 1 1 1 1 1 1 |
Set
x <- 1:10 y <- c(3:7,11:13) z <- c(4:5,1:2,12) intersect(x, y) # 3 4 5 6 7 union(x, y) # 1 2 3 4 5 6 7 8 9 10 11 12 13 setdiff(x, y) # 1 2 8 9 10 setequal(x, c(1:2, y[-c(6:8)], 8:10)) # TRUE unique(c(x,y,z)) # 1 2 3 4 5 6 7 8 9 10 11 12 13 c(x,y,z)[duplicated(c(x,y,z))] # 3 4 5 6 7 4 5 1 2 12 y %in% x is.element(y, x) # TRUE TRUE TRUE TRUE TRUE FALSE FALSE FALSE |
Quantile
fivenum(x) ## (x[floor(length(x)/4)] + x[ceiling(length(x)/4)+1])/2 ## (x[floor(length(x)/4*3)] + x[ceiling(length(x)/4*3)+1])/2 # -1.0 2.0 5.5 9.0 12.0 quantile(x) ## quantile(x,type = 7) # 0% 25% 50% 75% 100% # -1.00 2.25 5.50 8.75 12.00 3/4 + 1;x[4]*0.75+x[5]*0.25;13/4*3 + 1;x[10]*0.25+x[11]*0.75 # [1] 1.75 # [1] 2.25 # [1] 10.75 # [1] 8.75 quantile(x, type=6) # 0% 25% 50% 75% 100% # -1.00 1.75 5.50 9.25 12.00 15/4;x[3]*0.25+x[4]*0.75;15/4*3;x[11]*0.75+x[12]*0.25 # [1] 3.75 # [1] 1.75 # [1] 11.25 # [1] 9.25 IQR(x) # 6.5 IQR(x,type = 6) # [1] 7.5 |
Logarithm
log(8, base = 2) # 3 log10(10)s # 1 log2(8) # 3 log1p(2) # log1p(x) := log(1+x) # 1.098612 log(exp(3), base = exp(1)) log(exp(3)) # 3 expm1(1) # expm1(x) := exp(x) - 1 # 1.718282 exp(log(10)) # 10 |
R caution
log(8, base = 0) # 0!!! (-8)^(1/3) # !!!NaN sqrt(-2) # NaN |
Floor Ceiling Trunc
nums <-c(pi,-pi,exp(1),1.5,-1.5,0.5) floor(nums) # 3 -4 2 1 -2 0 ceiling(nums) # 4 -3 3 2 -1 1 trunc(nums) # 3 -3 2 1 -1 0 # nagtive numbers will be truncated, its round toward zero round(nums,1) # 3.1 -3.1 2.7 1.5 -1.5 0.5 round(nums) # 3 -3 3 2 -2 0 ## really weird, but from round(0.15, 1) could be either 0.1 or 0.2,coz 0.15 is not represented exactly round(-5.5+1:10) # -4 -4 -2 -2 0 0 2 2 4 4 # thus with .5 values rounded to even number comp <- c(1.4+2118.8i,2118,pi) round(comp,-2) # 0+2100i 2100+0i 0+0i signif(comp, digits = 6) # 1.4+2118.8i 2118.0+0.0i 3.14159+0.0i 29*50/50 # 29 29/50*50 # 29 floor(29/50 * 50) # 28 floor(29 * 50/50) # 29 .5-.8+.3 # -5.551115e-17 op <- options(digits = 22) 0.58*50 # 28.999999999999996 options(op) |
Formatting
x <- c(3.1,2.71828,100) format(x) # " 3.10000" " 2.71828" "100.00000" format(x, trim = T) # "3.10000" "2.71828" "100.00000" format(x, nsmall = 3,digits = 2) ## digits(minimum) --> nsmall # " 3.100" " 2.718" "100.000" format(c(x,"Hi"), justify = "left") # "3.1 " "2.71828" "100 " "Hi " format(x, scientific = T) # "3.10000e+00" "2.71828e+00" "1.00000e+02" format(ISOdate(2000, 1:12, 1), "%B") # "一月" "二月" "三月" "四月" "五月" "六月" "七月" "八月" "九月" "十月" "十一月" "十二月" sprintf("%.5s", x) # "3.1" "2.718" "100" sprintf("%1.0f", x) # "3" "3" "100" sprintf("%05.1f", x) # "003.1" "002.7" "100.0" sprintf("%-10f", x) ## left justified # "3.100000 " "2.718280 " "100.000000" sprintf("%g", 1e6 * x) ## exponential # "3.1e+06" "2.71828e+06" "1e+08" |
Interval
cut(0:10,breaks = 4) # [1] (-0.01,2.5] (-0.01,2.5] (-0.01,2.5] (2.5,5] (2.5,5] (2.5,5] # [7] (5,7.5] (5,7.5] (7.5,10] (7.5,10] (7.5,10] # Levels: (-0.01,2.5] (2.5,5] (5,7.5] (7.5,10] cbind(0:10, findInterval(0:10, c(2, 5, 7.5))) # [,1] [,2] # [1,] 0 0 # [2,] 1 0 # [3,] 2 1 # [4,] 3 1 # [5,] 4 1 # [6,] 5 2 # [7,] 6 2 # [8,] 7 2 # [9,] 8 3 # [10,] 9 3 # [11,] 10 3 cbind(0:10, findInterval(0:10, c(2, 5, 7.5),left.open = T)) # [,1] [,2] # [1,] 0 0 # [2,] 1 0 # [3,] 2 0 # [4,] 3 1 # [5,] 4 1 # [6,] 5 1 # [7,] 6 2 # [8,] 7 2 # [9,] 8 3 # [10,] 9 3 # [11,] 10 3 |
Formula
## This returns a string: "y ~ x1 + x2" # "y ~ x1 + x2" ## This returns a formula: as.formula("y ~ x1 + x2") # y ~ x1 + x2 class(y ~ x1 + x2) # "formula" f <- as.formula(paste("Y", paste(letters[1:5], collapse=" + "), sep=" ~ ")) # Y ~ a + b + c + d + e ## Take a look at f str(f) # Class 'formula' language Y ~ a + b + c + d + e # ..- attr(*, ".Environment")=<environment: R_GlobalEnv> environment(f) # Get each part f[[1]] # `~` f[[2]] # Y f[[3]] # a + b + c + d + e ## View the whole thing as a list as.list(f) # [[1]] # `~` # # [[2]] # Y # # [[3]] # a + b + c + d + e as.character(c(f[[1]],f[[2]],f[[3]])) # [1] "`~`" "Y" "a + b + c + d + e" as.character(f[[3]]) # "+" "a + b + c + d" "e" ## You can use deparse() to get a more natural looking string deparse(f[[3]]) # "a + b + c + d + e" deparse(f) # "Y ~ a + b + c + d + e" |
Piecewise function
\[\begin{equation} f(x)=\left\{ \begin{aligned} \cos(x), x<=0\\ \sin(x), 0<x<=2\\ \frac x5, x>2 \end{aligned} \right. \end{equation}\]x <- -5:5 f <- function(x){ if(x<=0) y=cos(x) else if(x>0&x<=2) y=sin(x) else y=x/5 return(y) } sapply(x, f) # 0.2836622 -0.6536436 -0.9899925 -0.4161468 0.5403023 1.0000000 0.8414710 0.9092974 0.6000000 0.8000000 1.0000000 |
y <- numeric(length(x)) y[x<=0] <- cos(x[x<=0]) y[x>0 & x<=2] <- sin(x[x>0 & x<=2]) y[x>2] <- x[x>2]/5 y # 0.2836622 -0.6536436 -0.9899925 -0.4161468 0.5403023 1.0000000 0.8414710 0.9092974 0.6000000 0.8000000 1.0000000 |
Complex
- complex(length.out = 0, real = numeric(), imaginary = numeric(),modulus = 1, argument = 0)
- as.complex(x, …)
- is.complex(x)
- Re(z)
- Im(z)
- Mod(z)
- Arg(z)
- Conj(z)
1i*(9:1) # 0+9i 0+8i 0+7i 0+6i 0+5i 0+4i 0+3i 0+2i 0+1i Re(1 + 2i);Im(1 + 2i) # 1 # 2 Mod(1 + 2i) # 2.236068 Arg(1 + 2i) # 1.107149 Conj(1 + 2i) # 1-2i |
z <- (rep(1:4, len = 9) + 1i*(9:1))/10 # 0.1+0.9i 0.2+0.8i 0.3+0.7i 0.4+0.6i 0.1+0.5i 0.2+0.4i 0.3+0.3i 0.4+0.2i 0.1+0.1i z.shift <- complex(modulus = Mod(z), argument = Arg(z) + pi) # -0.1-0.9i -0.2-0.8i -0.3-0.7i -0.4-0.6i -0.1-0.5i -0.2-0.4i -0.3-0.3i -0.4-0.2i -0.1-0.1i plot(z, xlim = c(-1,1), ylim = c(-1,1), col = "red", main = expression(paste("Rotation by "," ", pi == 180^o))) abline(h = 0, v = 0, col = "blue", lty = 3) points(z.shift, col = "orange") points(complex(re=Re(z)*2, im=tan(Re(z)*2)),col="gray20",pch = 16) points(complex(re=-Re(z)*2, im=tan(-Re(z)*2)),col="gray20",pch = 16) |
Large Integers
library(gmp) |
- add.bigz(e1, e2)
- sub.bigz(e1, e2 = NULL)
- mul.bigz(e1, e2)
- div.bigz(e1, e2)
- divq.bigz(e1,e2) ## == e1 %/% e2
- mod.bigz(e1, e2) ## == e1 %% e2
2^1024 # 1.60693804426e+60 as.bigz(2)^200 pow.bigz(2,200) # Big Integer ('bigz') : # [1] 1606938044258990275541962092341162602522202993782792835301376 nchar(as.character(as.bigz(1024^1024)^1024)[1]) # 2466038 |
Examples 1
Simple Calculator
add <- function(x, y) { return(x + y) } subtract <- function(x, y) { return(x - y) } multiply <- function(x, y) { return(x * y) } divide <- function(x, y) { return(x / y) } # take input from the user print("Select operation.") print("1.Add") print("2.Subtract") print("3.Multiply") print("4.Divide") choice = as.integer(readline(prompt="Enter choice[1/2/3/4]: ")) num1 = as.integer(readline(prompt="Enter first number: ")) num2 = as.integer(readline(prompt="Enter second number: ")) operator <- switch(choice,"+","-","*","/") result <- switch(choice, add(num1, num2), subtract(num1, num2), multiply(num1, num2), divide(num1, num2)) print(paste(num1, operator, num2, "=", result)) |
Triangle
triangle <- function(n, unit="degrees"){ rad <- n/180*pi mat <- array(0, dim = c(1,13)) colnames(mat) <- c("triangle", "sin","cos","tan", "csc","sec","cot", "asin","acos","atan", "sinh","cosh","tanh") mat[1] <- n mat[2] <- sin(rad) mat[3] <- cos(rad) mat[4] <- tan(rad) mat[5] <- 1/sin(rad) mat[6] <- 1/cos(rad) mat[7] <- 1/tan(rad) mat[8] <- asin(rad) mat[9] <- acos(rad) mat[10] <- atan(rad) mat[11] <- sinh(rad) mat[12] <- cosh(rad) mat[13] <- tanh(rad) return(mat) } ## NOTES: # (sqrt(6)-sqrt(2))/4 0.258 # (sqrt(6)+sqrt(2))/4 0.965 # 2-sqrt(3) 0.267 # sqrt(6)+sqrt(2) 3.86 # sqrt(6)-sqrt(2) 1.03 # 2+sqrt(3) 3.73 triangle(15) # triangle sin cos tan csc sec cot asin acos # 15 0.258819 0.9659258 0.2679492 3.863703 1.035276 3.732051 0.2648861 1.30591 # atan sinh cosh tanh # 0.2560528 0.2648002 1.034466 0.2559778 triangle(30) # triangle sin cos tan csc sec cot asin acos atan sinh # 30 0.5 0.8660254 0.5773503 2 1.154701 1.732051 0.5510696 1.019727 0.4823479 0.5478535 # cosh tanh # 1.140238 0.4804728 |
isPrime
# num = as.integer(readline(prompt="Enter a number: ")) num = 1573 flag = 0 if(num > 1) { flag = 1 for(i in 2:num ** 0.5) { if ((num %% i) == 0) { flag = 0 break } } } if(num == 2) flag = 1 if(flag == 1) { print(paste(num,"is a prime number")) } else { print(paste(num,"is not a prime number:",num,"/",i,"=",num/i)) } # "1573 is not a prime number: 1573 / 11 = 143" |
Factorial
num = 10 factorial(num) # 3628800 factorial <- function(num){ factorial = 1 if(num < 0) { print("Sorry, factorial does not exist for negative numbers") } else if(num == 0) { print("The factorial of 0 is 1") } else { for(i in 1:num) { factorial = factorial * i } return(paste("The factorial of", num ,"is",factorial)) } } factorial(10) # "The factorial of 10 is 3628800" ## Find Factorial of a number using recursion recursive_factorial <- function(n) { if(n <= 1) { return(1) } else { return(n * recursive_factorial(n-1)) } } recursive_factorial(num) # 3628800 rbind(system.time(recursive_factorial(3000)),system.time(factorial(3000))) # user.self sys.self elapsed user.child sys.child # [1,] 0.01 0 0.02 NA NA # [2,] 0.00 0 0.00 NA NA |
"%!%" <- function(x, y) {factorial(x)/factorial(y)} 10%!%8 # 90 |
Armstrong number(narcissistic number)
is.Armstrong <- function(num){ sum = 0 temp = num while(temp > 0) { digit = temp %% 10 sum = sum + (digit ** nchar(num)) temp = temp %/% 10 } if(num == sum) { return(TRUE) }else{ return(FALSE) } } c(100:100000)[sapply(100:100000,is.Armstrong)] # [1] 153 370 371 407 1634 8208 9474 54748 92727 93084 |
Examples 2
Bisect
cubic <- function (x, a=0, b=0, c=0, d=0) a*x^3 + b*x^2 + c*x + d a <- 1 ; b <- -2 ; c <- 3 ; d <- -4 result <- uniroot(cubic, c(-5,5), a = a, b = b,c = c,d = d, tol = 0.000001) result # $`root` # [1] 1.650629 # # $f.root # [1] 5.33745e-09 # # $iter # [1] 8 # # $init.it # [1] NA # # $estim.prec # [1] 5e-07 |
bisect <- function(f, a, b, eps = 1e-5, ...){ if (!is.function(f)){ list(fail="f is not a function!") } else if(f(a, ...) * f(b, ...) > 0) list(fail="finding root is fail!") else{ repeat{ if (abs(b - a) < eps) break x <- (a + b)/2 if (f(a, ...) * f(x, ...) < 0) b <- x else a <- x } list(root = (a+b)/2, fun = f(x, ...)) } } bisect(cubic, -5, 5, eps = 1e-6, 1, -2 , 3, -4) # $`root` # [1] 1.650629 # # $fun # [1] 2.048755e-06 |
Newton
Jacobi matrix:
\[\begin{equation} \textbf{J}(x) = \begin{bmatrix} \dfrac{\partial f_1}{\partial x_1} & \dfrac{\partial f_1}{\partial x_2} & \cdots &\dfrac{\partial f_1}{\partial x_n} \\ \dfrac{\partial f_2}{\partial x_1} & \dfrac{\partial f_2}{\partial x_2} &\cdots& \dfrac{\partial f_2}{\partial x_n} \\ \vdots& \vdots& \ddots & \vdots \\ \dfrac{\partial f_m}{\partial x_1} &\dfrac{\partial f_m}{\partial x_2} & \cdots &\dfrac{\partial f_m}{\partial x_n} \\ \end{bmatrix} \end{equation}\]If we want to solve this equation (start = c(0,1), ep <= 1e-5):
\[\begin{equation} \left\{ \begin{aligned} \ x_1^2 + x_2^2 - 5 = 0\\ \ (x_1 + 1)x_2 - (3x_2 + 1) = 0\\ \end{aligned} \right. \end{equation}\]Newtons<-function (fun, x, ep = 1e-5, it_max = 100){ index <- 0; k <- 1 while (k <= it_max){ x1 <- x; obj <- fun(x); x <- x - solve(obj$J, obj$f); norm <- sqrt((x - x1) %*% (x - x1)) if (norm < ep){ success <- 1; break } k <- k + 1 } obj <- fun(x); list(root = x, val = obj$f, iteration = k, success = success) } |
funs <- function(x){ f <- c(x[1]^2 + x[2]^2 - 5, (x[1] + 1) * x[2] - (3 * x[1] + 1)) J <- matrix( c(2 * x[1], 2 * x[2], x[2] - 3, x[1] + 1), nrow = 2, byrow = T) list(f = f, J = J) } |
Newtons(funs, c(0,1)) |
## $root ## [1] 1 2 ## ## $val ## [1] 1.598721e-14 6.217249e-15 ## ## $iteration ## [1] 6 ## ## $success ## [1] 1 |
Differences
diff(x, lag = 1, differences(一阶差分) = 1, …) lag(x, k = 1, …)
diff(1:10) # 1 1 1 1 1 1 1 1 1 diff(1:10, lag=2,differences=1) # 2 2 2 2 2 2 2 2 diff(1:10, 2, 2) # 0 0 0 0 0 0 lag(1:10,5) # [1] 1 2 3 4 5 6 7 8 9 10 # attr(,"tsp") # [1] -4 5 1 lagdf <- function(x, k) { c(rep(NA, k), x)[1 : length(x)] } |
Integral
Newton-Leibniz Formula:
\[ \int_{a}^{b} f(x)dx = F(b) - F(a) \]
梯形公式(n=1):
\[ T= (b - a)[\frac {1}{2} f(a) + \frac {1}{2} f(b)] \]
辛普森公式(n=2, Simpson’s rule):
\[ S=(b-a)[\frac {1}{6} f(a) + \frac {4}{6} f(\frac{a+b}{2}) + \frac {1}{6} f(b) ] \] \(Divide [a,b] to\, n fractions, step h = \frac{b-a}{n}, then\, x_i = a + ih,i = 1,2,\cdots,x_n:\)
复化梯形公式: \[ \int_{a}^{b} f(x)dx \approx \sum_{i=1}^{n}\int_{x_i}^{x_{i+1}}dx=T_n= \frac{h}{2}[f(a) + 2\sum_{i=1}^{n-1}f(a+ih) + f(b)] \] 复化辛普森公式: \[ S_n = \frac{h}{6}\sum_{i=0}^{n-1}[f(x_i) + 4f(x_{i + \frac{1}{2}}) + f(x_{i+1}))] \]
area <- function(f, a, b, eps = 1.0e-06, lim = 10) { fun1 <- function(f, a, b, fa, fb, a0, eps, lim, fun) { d <- (a + b) / 2 h <- (b - a) / 4 fd <- f(d) a1 <- h * (fa + fd) a2 <- h * (fb + fd) if (abs(a0 - a1 - a2) < eps || lim == 0) return(a1 + a2) else { return( fun(f, a, d, fa, fd, a1, eps, lim - 1, fun) + fun(f, d, b, fd, fb, a2, eps, lim - 1, fun) ) } } fa <- f(a) fb <- f(b) a0 <- ((fa + fb) * (b - a)) / 2 fun1(f, a, b, fa, fb, a0, eps, lim, fun1) } f <- function(x) 1 / x^2 quad <- area(f, 1, 5) quad |
All articles in this blog are licensed under CC BY-NC-SA 4.0 unless stating additionally.
Related Articles
Comment
